# A preconditioned Hessian-like proximal gradient algorithm for sparse ML-DOA estimation

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

► To cite this version:

Thomas Aussaguès, Anne Ferréol, Alice Delmer, Pascal Larzabal. A preconditioned Hessian-like proximal gradient algorithm for sparse ML-DOA estimation. 10th IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP 2025), Dec 2025, Punta Cana, Dominican Republic. hal-05241987

HAL Id: hal-05241987

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

Submitted on 5 Sep 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 Attribution 4.0 International License
A PRECONDITIONED HESSIAN-LIKE
PROXIMAL GRADIENT ALGORITHM FOR
SPARSE ML-DOA ESTIMATION

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**—Maximum Likelihood (ML) Direction-of-Arrival (DoA) estimation under the Vectorized Covariance Matrix Model (VCMM) offers improved performance but remains computationally intensive. To alleviate this issue, sparse formulations have been adopted and shown to be asymptotically equivalent to the ML using a proper regularization parameter choice. Yet, the resulting sparse criterion remains intricate and is typically minimized using slow first-order methods such as the Proximal Gradient Algorithm (PGA). To improve convergence speed, an Hessian-like preconditioned PGA is introduced. As the associated proximal operator lacks a closed-form expression, an IRL1 scheme, known to converge to a critical point of the $l_0$-regularized criterion is employed, enabling its numerical evaluation via the Chambolle-Pock (CP) algorithm. This approach yields faster convergence and enhanced performance in low-SNR scenarios by exploiting the local geometry encoded in the Hessian.

**Index Terms**—Sparse optimization, Second-order, Direction-of-Arrival estimation, Maximum Likelihood, Preconditioning

I. INTRODUCTION

The advent of sparse estimation techniques over last decades marked a significant shift for Direction-of-Arrival (DoA) estimation [1], [2] offering superior performance over classical methods such as MUSIC [3] in adverse conditions [4], [5], [6], [7]. Sparse estimators are used in combination with the Vectorized Covariance Matrix Model (VCMM) [8] providing access to a Virtual Array (VA) of $N^2$ elements thereby enhancing angular resolution and increasing the number of identifiable sources [9]. Recent works [10] established equivalence between sparse and Maximum Likelihood (ML) estimators [11], for appropriate regularization parameter choice, thus enabling tractable ML implementation through sparse estimators.

Then, DoA estimation amounts to minimizing a regularized least-squares criterion, where sparsity is enforced through the $l_0$-norm, continuous relaxations such as CEL0 [12] (Continuous Exact $l_0$), or the convex $l_1$-norm. The optimization is typically handled by the Proximal Gradient Algorithm (PGA) also referred as the Forward-Backward Splitting algorithm [13] which reduces to Iterative Hard Thresholding (IHT) [14] for $l_0$ regularization, and to Iterative Soft-Thresholding (ISTA) [15], [16] in the convex $l_1$ case.

Despite its widespread use, PGA remains a first-order method which experiences slow convergence, particularly under ill-conditioning [17]. Inertial [18], [19], [20] and variable stepsize variants [21] have been proposed yielding moder-

ate improvements. To accelerate convergence, preconditioning techniques have been introduced [22], [23], [24]. Such strategies have been proven effective in numerous applications including tomography [25], hyperspectral imaging [26], radioastronomy [27], and low-rank matrix recovery [28]. When the proximal operator has no closed-form, as in the preconditioned case, it can be computed through inner iterations of the Chambolle-Pock (CP) algorithm [29], assuming a convex regularizer such as the $l_1$-norm.

This work introduces an Hessian-Like preconditioned PGA for ML-DoA estimation under CEL0 regularization. The preconditioning matrix is chosen as an approximation of the Hessian as it is not invertible. Since the CEL0 functional, whose interest for DoA estimation has been demonstrated [30], is non-convex [12], numerical computation of the proximal operator after preconditioning through the Chambolle-Pock (CP) algorithm, which requires convexity, is not feasible. To address this limitation, a convex weighted $l_1$ majorizer is found using the Iteratively Reweighted $l_1$ (IRL1) scheme [31], whose convergence to a critical point of the CEL0 criterion is proven [31]. Then, a CP-based evaluation of the proximal step is now possible. The resulting estimator achieves faster convergence and improved performance under low Signal-to-Noise Ratio (SNR) conditions by exploiting the criterion local geometry information encoded in the Hessian matrix.

II. SIGNAL MODELLING

Considering $M$ independent narrowband plane waves impinging on an array of $N$ antennas from directions $\theta = \{\theta_1, \dots, \theta_M\}$, the array output is:

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

where $\mathbf{A}(\theta) = [\mathbf{a}(\theta_1), \dots, \mathbf{a}(\theta_M)]$ is the steering matrix, $s_{\theta}(t) = [s_1(t), \dots, s_M(t)]^T$ the complex envelopes of the emitted signals which are considered unknown but deterministic and $\mathbf{n}(t)$ a complex white circular Gaussian noise vector, independent of $s_{\theta}(t)$, with covariance matrix $\mathbb{E}[\mathbf{n}(t)\mathbf{n}^H(t)] = \sigma^2\mathbf{I}_N$. The covariance matrix of (1) is then:

$$ \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}[s_{\theta}(t)s_{\theta}^H(t)]$ is the sources covariance matrix. Assuming temporally uncorrelated emitters, $\mathbf{R}_s$ is diagonal
such that $\mathbf{R}_s = \text{diag}(\gamma_\theta)$ with $\gamma_\theta$ the sources power vector. Then, vectorizing (2) yields the VCMM observation [8]:

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

with $\text{vec}(\cdot)$ the column-wise vectorisation operator and $\mathbf{B}(\theta)$ the VA steering matrix formed by the vectors $\mathbf{b}(\theta_m) = \mathbf{a}^*(\theta_m) \otimes \mathbf{a}(\theta_m)$. In practical applications, the true covariance matrix $\mathbf{R}_x$ is not accessible and is therefore replaced by its ML estimate $\hat{\mathbf{R}}_x$, computed from $K$ independent and identically distributed array snapshots $\mathbf{x}(t_k)$, with $1 \le k \le K$. Assuming temporally white noise, that is, $\mathbb{E}[\mathbf{n}^H(t_i)\mathbf{n}(t_j)] = 0$ for $i \ne j$, the estimated covariance matrix reads:

$$ \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 (4) $$

with $\Delta\mathbf{R}_x$ a random matrix following a complex Wishart distribution [32] yielding an estimation error $\delta$ in (3):

$$ \mathbf{r} = \text{vec}(\hat{\mathbf{R}}_x - \sigma^2\mathbf{I}_N) = \mathbf{B}(\theta)\gamma_\theta + \delta \quad (5) $$

The estimation error $\delta$ is asymptotically Gaussian with a covariance matrix $\mathbf{\Gamma} = \mathbb{E}[\delta\delta^H] = \frac{1}{K}(\mathbf{R}_x^T \otimes \mathbf{R}_x) \ne \mathbf{I}_{N^2}$ thereby revealing a non-white estimation error [33], [34]. Thus, a pre-whitening transform is performed on (5) to handle the correlated Gaussian estimation error $\delta$ [10]. This yields a novel observation $\mathbf{y}$ with a white Gaussian estimation error:

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

where $\hat{\mathbf{W}} = \sqrt{K}(\hat{\mathbf{R}}_x^T \otimes \hat{\mathbf{R}}_x)^{-1/2}$ is the estimated whitening matrix, $\mathbf{B}_w(\theta) = \hat{\mathbf{W}}\mathbf{B}(\theta)$ denotes the VA steering matrix after whitening and $\delta_w = \hat{\mathbf{W}}\delta$ a spatially white noise satisfying $\mathbb{E}[\delta_w\delta_w^H] = \mathbf{I}_{N^2}$.

### III. SPARSE ESTIMATION

#### A. Sparse problem formulation

To circumvent the ML implementation pitfalls, a sparse estimator based on (6) is employed, exploiting equivalence under white Gaussian estimation error between sparse and ML estimators [10]. Using a grid of $G$ directions $\varphi = \{\varphi_1, \dots, \varphi_G\}$ and assuming that the sources directions $\theta$ lie within $\varphi$, eq. (6) can be rewritten using the sparse framework:

$$ \mathbf{y} = \mathbf{B}_w\gamma_0 + \delta_w \quad (7) $$

where $\mathbf{B}_w = \hat{\mathbf{W}}\mathbf{B}$ with $\mathbf{B} = [\mathbf{b}(\varphi_1), \dots, \mathbf{b}(\varphi_G)]$ is an overcomplete dictionary of size $N^2 \times G$, $G \gg N^2$ and $\gamma_0$ is an $M$-sparse vector having $M$ non-zero components at directions satisfying $\varphi_g = \theta_m$.

Then, DoAs can be estimated from the sparse model (7) using the grid directions $\varphi_g$ associated to non-null entries of the sparse vector $\gamma_0$. The sparse vector estimation is then addressed by leveraging the sparsity prior of $\gamma_0$ through the addition of an $\ell_0$ regularizer [35]:

$$ \min_{\gamma} \{ \mathcal{J}_{\ell_0}(\lambda, \gamma) = \frac{1}{2} \|\mathbf{y} - \mathbf{B}_w \gamma\|_2^2 + \lambda \|\gamma\|_0 \} \quad (8) $$

where $\lambda > 0$ refers to the regularization parameter, balancing the solution sparsity towards data fidelity. To ensure equivalence between sparse and ML estimators, $\lambda$ is selected following the approach of [10] ie. $\lambda = \frac{1}{2}\mathcal{F}^{-1}(0.10)$ with $\mathcal{F}^{-1}$ the inverse $\chi_2$ cumulative distribution function with $N^2 - M$ degrees of freedom. Finally, the $\ell_0$-norm is replaced by the CEL0 penalty $\phi_{\text{CEL0}}$, [12]. The resulting objective function, denoted by $\mathcal{J}_{\text{CEL0}}$, exhibits fewer local minima while preserving the same global minimizer as $\mathcal{J}_{\ell_0}$, thereby facilitating its optimization [12], [36]. Therefore, the regularization parameter $\lambda$ previously derived remains valid in this relaxed formulation.

#### B. The Proximal Gradient Algorithm (PGA)

The minimization of (8) is performed using the PGA algorithm [13], namely:

$$ \gamma^{(i+1)} = \text{prox}_{\beta\phi_{\text{CEL0}}(\cdot, \lambda)}(\gamma^{(i)} - \beta\mathbf{B}_w^H(\mathbf{B}_w\gamma^{(i)} - \mathbf{y})) \quad (9) $$

with $\text{prox}_{\beta\phi_{\text{CEL0}}(\cdot, \lambda)}$ the CEL0 proximal operator [12], [37] and $\beta \in ]0, \frac{1}{L}$ [the stepsize where $\mathcal{L} = \lambda_{\text{max}}(\mathcal{H}_w)$ stands for the largest eigenvalue of the least squares term, in (8), Hessian matrix:

$$ \mathcal{H}_w = \mathbf{B}_w^H \mathbf{B}_w \quad (10) $$

### IV. A PRECONDITIONED HESSIAN-LIKE PROXIMAL GRADIENT ALGORITHM FOR SPARSE DOA ESTIMATION

#### A. Preconditioned PGA

Classical PGA (9) often experiences slow convergence, even using inertial extensions [18], [19]. To further enhance convergence rates, preconditioning strategies have been developed [38], [39], [22], [23], [24] involving a matrix $\mathbf{U}$ refered as a preconditioner. Formally, the preconditioned PGA reads [22]:

$$ \gamma^{(i+1)} = \text{prox}_{\beta\mathbf{U}\phi_{\text{CEL0}}(\cdot, \lambda)}(\gamma^{(i)} - \beta\mathbf{U}^{-1}\mathbf{B}_w^H(\mathbf{B}_w\gamma^{(i)} - \mathbf{y})) \quad (11) $$

with $\mathbf{U}$ the preconditioner and $\beta\mathbf{U} \in ]0, \frac{1}{L}$ [, where $\mathcal{L}_\mathbf{U} = \lambda_{\text{max}}(\mathbf{U}^{-1}\mathcal{H}_w)$ is the Lipschitz constant of the preconditioned gradient operator, a stepsize value ensuring convergence towards a critical point of $\mathcal{J}_{\text{CEL0}}$ [23], [25]. $\text{prox}_{\beta\phi_{\text{CEL0}}(\cdot, \lambda)}$ is the proximal operator of the CEL0 penalty relatively to the norm induced by $\mathbf{U}$ ie. $\|\mathbf{x}\|_\mathbf{U}^2 = \mathbf{x}^H\mathbf{U}\mathbf{x}$ [40]:

$$ \text{prox}_{\beta\mathbf{U}\phi_{\text{CEL0}}(\cdot, \lambda)}(\gamma) = \arg\min_{\mathbf{u}} \frac{1}{2\beta\mathbf{U}} \|\gamma - \mathbf{u}\|_{\mathbf{U}}^2 + \phi_{\text{CEL0}}(\mathbf{u}, \lambda) \quad (12) $$

#### B. Preconditioner design criteria

This section outlines the criteria for designing an efficient preconditioner. The preconditioner $\mathbf{U}$ should be designed such that: (i) its eigenvalues are clustered, yielding a small condition number and allowing for large gradient steps, and (ii) both $\mathbf{U}$ and its inverse $\mathbf{U}^{-1}$ should be easily computed [25]. Finally, to construct an efficient preconditioner, the preconditioned PGA (11) is interpreted as a Majorization-Minimization (MM) algorithm, with preconditioners chosen to ensure tight majorizing functions (iii). The following specifies this final design criterion. Let $f(\gamma) := \frac{1}{2}\|\mathbf{y} - \mathbf{B}_w\gamma\|^2$ the least squares term
in (8) and $g(\gamma) := \phi_{\text{CEL0}}(\gamma, \lambda)$ the penalty term. A majorizer of $f$ is then given by:

$$Q(\gamma, \gamma^{(i)}) = f(\gamma^{(i)}) + \langle \nabla f(\gamma^{(i)}) | \gamma - \gamma^{(i)} \rangle + \frac{1}{2\beta_U} \|\gamma - \gamma^{(i)}\|_{\text{U}}^2 \quad (13)$$

where $\langle x|y \rangle = x^H y$. Then, one PGA iteration (11) can be casted as the minimization of the following majorizer of $J_{\text{CEL0}}$ [19], [23]:

$$\gamma^{(i+1)} = \arg \min_{\gamma} \{ f(\gamma^{(i)}) + \langle \nabla f(\gamma^{(i)}) | \gamma - \gamma^{(i)} \rangle + \frac{1}{2\beta_U} \|\gamma - \gamma^{(i)}\|_{\text{U}}^2 + g(\gamma^{(i)}) \} \quad (14)$$

Thus, $U$ should be selected such that (iii) the majorizer $Q$ in (14) is close to $f$ ie. to the least squares term in (8).

### C. Hessian-like preconditioner design

This section focuses on designing a Hessian-like preconditioner based on criteria of subsection IV-B, aiming to produce tight upper bounds while leveraging local curvature information. The design involves balancing majorizer tightness with preconditioner conditioning. Consider a preconditioner of the form:

$$U = \mathcal{H}_w + \eta \mathbf{I}_G \quad (15)$$

with $\eta > 0$. This preconditioner choice satisfies criterion (ii). Then, to assess criterion (iii), the distance between the majorizer obtained using $U$ (15) and the optimal majorizer obtained using the Hessian matrix $\mathcal{H}_w$ (10) is:

$$\begin{aligned} \varepsilon(\gamma, U) &= \|\gamma - \gamma_0\|_{2,U}^2 - \|\gamma - \gamma_0\|_{2,\mathcal{H}_w}^2 \\ &= |(\gamma - \gamma_0)^H (U - \mathcal{H}_w)(\gamma - \gamma_0)| \quad (16) \\ &= \eta \times \|\gamma - \gamma_0\|_{2}^2 \end{aligned}$$

with:

$$\eta = \frac{\lambda_{\max}(\mathcal{H}_w)}{\rho(\mathbf{U}) - 1} \quad (17)$$

$$\rho(\mathbf{U}) = \frac{\lambda_{\max}(\mathbf{U})}{\lambda_{\min}(\mathbf{U})} = \frac{\eta + \lambda_{\max}(\mathcal{H}_w)}{\eta} \quad (18)$$

where $\rho(\mathbf{U})$ is the condition number of $U$ left to the user. Eq. (16) highlights the trade-off between the tightness of the majorizer and the condition number of the Hessian-like preconditioner, as $\epsilon(\gamma, U)$ is inversely proportional to $\rho(\mathbf{U})$. To construct the preconditioner, a condition number is thus chosen and the corresponding $\eta$ value is obtained via (18). As discussed in subsection IV-B, an effective preconditioner should yield both a low condition number (i) and a tight majorizer (iii). Several condition numbers $\rho(\mathbf{U})$ are evaluated in section V based on these criteria.

### D. Numerical evaluation of the proximal operator

Using a preconditioner $U$ in the PGA (9) modifies the proximal step (12). To circumvent the lack of a closed-form of (12), the proximal operator can be numerically evaluated using the CP [29]. Yet, the CP algorithm requires the convexity of $\phi_{\text{CEL0}}$ which is not the case [12]. To remediate the non-convexity of $\phi_{\text{CEL0}}$, we propose to employ the IRL1 algorithm [31] applied to $\phi_{\text{CEL0}}$ [41] whose convergence towards a

critical point of $\mathcal{J}_{\text{CEL0}}$ is proven in [31] as $\mathcal{J}_{\text{CEL0}}$ satisfies the Kurdyka-Łojasiewicz property [12]. At each iteration, IRL1 minimizes $\mathcal{J}_{\ell_1, w^{(i)}}$, a convex weighted-$\ell_1$ majorizer of $\mathcal{J}_{\text{CEL0}}$, equal to it at the current solution:

$$\mathcal{J}_{\ell_1, w^{(i)}}(\lambda, \gamma) = \frac{1}{2} \|\mathbf{y} - \mathbf{B}_w \gamma\|_2^2 + \|\gamma\|_{1, w^{(i)}} \quad (19)$$

$\mathcal{J}_{\ell_1, w^{(i)}}$ is optimized using preconditioned ISTA [15]:

$$\gamma^{(i+1)} = \operatorname{prox}_{\beta_U \ell_1, w^{(i)}} (\gamma^{(i)} - \beta_U U^{-1} \mathbf{B}_w^H (\mathbf{B}_w \gamma^{(i)} - \mathbf{y})) \quad (20)$$

where $w^{(i)}$ are the weights computed at each iteration using [41] and $\operatorname{prox}_{\beta\ell_1, w^{(i)}}$, the proximal operator of $\ell_1, w$ relatively to the norm induced by $U$ (12). Additionally, the efficiency of IRL1 for minimizing a CEL0-based objective has been demonstrated in various applications: [30], [36] (DoA), [41], [42] (microscopy). The proposed preconditioned IRL1-CEL0 is summarized in Alg.1.

### Algorithm 1 Preconditioned IRL1-CEL0

**Input:** $\mathbf{y}, \gamma^{(0)}, \mathbf{B}_w, \mathbf{U}, \lambda > 0$

**Output:** $\gamma^{(i)}$

1. **while** stopping criterion is not satisfied **do**

2. Update $\ell_1$-norm weights [41] with $\alpha = \|\mathbf{b}_w(\varphi_g)\|_2$:

$$w_g^{(i)} \leftarrow (\sqrt{2\lambda\alpha} - \alpha^2 |\gamma_g^{(i)}|) \mathbf{1}_{\{|\gamma_g^{(i)}| \le \sqrt{2\lambda\alpha}\}} (\gamma_g^{(i)})$$

3. **for** $j = 1, \dots, N_{\text{ISTA}}$ **do**

4. Preconditioned ISTA iteration: compute the preconditioned weighted $\ell_1$-norm proximal operator using the CP [20] algorithm:

$$\gamma^{(j)} \leftarrow \operatorname{prox}_{\beta_U \ell_1, w^{(i)}} (\gamma^{(j)} - \beta_U U^{-1} \mathbf{B}_w^H (\mathbf{B}_w \gamma^{(j)} - \mathbf{y}))$$

5. **end for**

6. Update $\gamma^{(i)}$: $\gamma^{(i)} \leftarrow \gamma^{(j)}$

7. **end while**

8. **return** $\gamma^{(i)}$

## V. NUMERICAL EXPERIMENTS

### A. Experimental setup

A circular array of $N = 4$ antennas is considered, with three elements uniformly positioned on a circle of radius $0.35\lambda_0$, with $\lambda_0$ the wavelength, and one at the center. $M = 2$ independent sources at $\theta_1 = 0^\circ$ and $\theta_2 = 20^\circ$ with equal SNR defined as SNR = $10 \log_{10} (\frac{\gamma_m}{\sigma^2})$ and $K = 200$ array snapshots are considered. An oracle grid with stepsize $2^\circ$ containing the ML estimates is employed to remove any grid bias (see [10] for details). To increase convergence speed of Alg.1, ISTA iterations within Alg.1 are replaced by FISTA [19] ones. The proposed method is compared against classical and variable stepsize PGA [21], while $\ell_1$-based convex relaxations of (8) are excluded due to regularization parameter selection issues [35]. All algorithms are applied with 1000 iterations, initialization $\gamma^{(0)} = 0_{G \times 1}$, and convergence is assessed using a joint criterion based on the relative difference between consecutive iterates, the difference of their corresponding criterion values and a maximum number of iterations.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Projections of majorizers $\mathcal{Q}$ onto directions $\{\varphi_1, \varphi_2\}$ in two cases with SNR = 10 dB. Top: projection onto the sources directions $\{\varphi_1, \varphi_2\} = \{\theta_1, \theta_2\} = \{0^\circ, 20^\circ\}$. Bottom: projection onto others directions $\{\varphi_1, \varphi_2\} = \{114^\circ, 115^\circ\} \neq \{\theta_1, \theta_2\}$. Black dots represent the global minimizers.</figcaption></figure>

## B. Preconditioner influence on the criterion

First, the preconditioner's influence on the criterion is visualized on Fig.1 for the following choices: $U_0 = I_G$ (PGA (9) without preconditioning) and $U = H_w + \eta I$ (15) with $H_w$ the Hessian matrix (10) for several condition numbers.

1) **In the sources directions:** All preconditioners match the true criterion circular geometry. Yet, $U$ with large condition numbers ($1.6 \times 10^4$, $1.6 \times 10^3$ and $1.6 \times 10^2$) yields tighter majorizers and so small distances $\epsilon(\gamma, U)$ (16).

2) **Others directions:** The high correlation between steering vectors yields elliptical contour lines [43]. The standard preconditioner $U_0$ fails to capture this geometry, resulting in a loose majorizer. In contrast, $U$ (with $\rho(U) \in \{1.6 \times 10^4, 1.6 \times 10^3, 1.6 \times 10^2\}$) yield tighter majorizers with matching elliptical contours. Yet, associated condition numbers remain high for $U$ with $1.6 \times 10^4$ and $1.6 \times 10^3$. The choice $U$ with $1.6 \times 10^2$ appears as a good compromise as it yields both a tight majorizer and a low condition number, thereby satisfying the desired criteria (i), (ii), and (iii), and is thus retained for subsequent analysis.

## C. Convergence analysis of the preconditioned Hessian PGA

Convergence speed is assessed using $\Delta\mathcal{J} = |\mathcal{J}_{LS}(\gamma^{(i)}) - \mathcal{J}_{ML}^*|$ the distance toward the least squares term in (8) optimal value is represented on Fig.2 for PGA (9) (ie. with $U_0$), IRL1 with $U_0$ and preconditioned IRL1 using $U$. Fig.2 shows the preconditioning benefits: it outperforms both PGA and standard IRL1. In both scenarios, it converges towards the criterion (8) global minimizer ie. the minimum of the ML criterion in few iterations (27 and 16 iterations). Conversely, standard IRL1 fails to reach the global minimizer in the low SNR scenario. Under more favourable conditions, it achieves the global minimizer yet requiring 800 iterations

while preconditioned IRL1 need only 27 iterations. Furthermore, preconditioned IRL1 outperforms the variable stepsize PGA in both scenarios. Finally, a performance analysis of the proposed method is conducted in [43].

<figure><img src="image_3.png" /><figcaption>Fig. 2. Convergence curves with $\mathcal{J}_{ML}^*$ the optimal value of the least squares term in the objective (8) for (top) SNR = 5 dB and (bottom) SNR = 20 dB.</figcaption></figure>

# VI. CONCLUSION

This paper introduces a Hessian-like preconditioned Proximal Gradient Algorithm for ML-DoA estimation. To this end, the initially non-convex penalty is minimized using the IRL1 scheme thus enabling the numerical evaluation of the proximal operator through the Chambolle-Pock algorithm. Numerical experiments show that exploiting the criterion local geometry through the Hessian-like preconditioner yields massive convergence speed improvements together with performance enhancements in low SNR scenarios.
REFERENCES

[1] Z. Yang, J. Li, P. Stoica, and L. Xie, "Chapter 11 - sparse methods for direction-of-arrival estimation," in *Academic Press Library in Signal Processing, Volume 7*, R. Chellappa and S. Theodoridis, Eds. Academic Press, 2018, pp. 509-581.

[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] R. Schmidt, "Multiple emitter location and signal parameter estimation," *IEEE Transactions on Antennas and Propagation*, vol. 34, no. 3, pp. 276-280, 1986.

[4] 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.

[5] 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, 2015.

[6] J. Liu, W. Zhou, F. H. Juwono, and D. D. Huang, "Reweighted smoothed l0-norm based DoA estimation for MIMO radar," *Signal Processing*, vol. 137, pp. 44-51, 2017.

[7] H. Zhang, X. Wang, J. Lin, and Z. Du, "Sparse optimization driven deep unfolding network for DoA estimation," in *2024 IEEE International Instrumentation and Measurement Technology Conference (I2MTC)*, 2024, pp. 1-6.

[8] 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.

[9] 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.

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

[11] 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.

[12] E. Soubies, L. Blanc-Féraud, and G. Aubert, "A Continuous Exact l0 penalty (CEL0) for least squares regularized problem," *SIAM Journal on Imaging Sciences*, vol. 8, no. 3, pp. pp. 1607-1639 (33 pages), Jul. 2015.

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

[14] T. Blumensath and M. E. Davies, "Iterative hard thresholding for compressed sensing," *Applied and Computational Harmonic Analysis*, vol. 27, no. 3, pp. 265-274, 2009.

[15] I. Daubechies, M. Defrise, and C. Mol, "An iterative thresholding algorithm for linear inverse problems with a sparsity constrains," *Communications on Pure and Applied Mathematics*, vol. 57, 11 2004.

[16] J. Bect, L. Blanc-Féraud, G. Aubert, and A. Chambolle, "A 11-unified variational framework for image restoration," in *Computer Vision - ECCV 2004. 8th European Conference on Computer Vision, Prague, Czech Republic, May 11-14, 2004. Proceedings, Part IV.*, ser. Lecture Notes in Computer Science, J. M. T. Pajdla, Ed., vol. 3024. Prague, Czech Republic: Springer, May 2004, pp. 1-13.

[17] G. H.-G. Chen and R. T. Rockafellar, "Convergence rates in forward-backward splitting," *SIAM J. Optim.*, vol. 7, pp. 421-444, 1997.

[18] Y. Nesterov, "A method for solving the convex programming problem with convergence rate $o(1/k^2)$," p. 543, 1983.

[19] A. Beck and M. Teboulle, "A fast iterative shrinkage-thresholding algorithm for linear inverse problems," *SIAM Journal on Imaging Sciences*, vol. 2, no. 1, pp. 183-202, 2009.

[20] A. Chambolle and C. H. Dossal, "On the convergence of the iterates of "FISTA"," *Journal of Optimization Theory and Applications*, vol. Volume 166, no. Issue 3, p. 25, Aug. 2015.

[21] T. Aussaguès, A. Ferréol, A. Delmer, and P. Larzabal, "Low Cost Quasi Gauss-Newton Algorithm Implementation for 2D ML-DoA Estimation using a Sparse Representation of Array Covariance," in *EUSIPCO 2025 - 33rd European Signal Processing Conference*, Sep. 2025.

[22] J. Bonnans, J. C. Gilbert, C. Lemaréchal, and C. Sagastizabal, "A family of variable metric proximal methods," *Mathematical Programming*, vol. 68, pp. 15-47, 05 1995.

[23] E. Chouzenoux, J.-C. Pesquet, and A. Repetti, "Variable metric forward-backward algorithm for minimizing the sum of a differentiable function and a convex function," *Journal of Optimization Theory and Applications*, vol. 162, pp. 107-132, 07 2014.

[24] P. L. Combettes and B. C. V. and, "Variable metric forward-backward splitting with applications to monotone inclusions in duality," *Optimization*, vol. 63, no. 9, pp. 1289-1318, 2014.

[25] M. Savanier, E. Chouzenoux, J.-C. Pesquet, and C. Riddell, "Unmatched preconditioning of the proximal gradient algorithm," *IEEE Signal Processing Letters*, vol. 29, pp. 1122-1126, 2022.

[26] A. Repetti, E. Chouzenoux, and J.-C. Pesquet, "A preconditioned forward-backward approach with application to large-scale nonconvex spectral unmixing problems," in *2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, 2014, pp. 1498-1502.

[27] H. Bester, A. Repetti, S. Perkins, O. Smirnov, and J. Kenyon, "A practical preconditioner for wide-field continuum imaging of radio interferometric data," *ADASS XXX 2020; Conference date: 08-11-2020 Through 12-11-2020*, Nov. 2020.

[28] X. Yao and W. Dai, "A preconditioned hessian proximal algorithm for spectral compressed sensing," in *2023 IEEE 9th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP)*, 2023, pp. 476-480.

[29] A. Chambolle and T. Pock, "A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging," *Journal of Mathematical Imaging and Vision*, vol. 40, no. 1, pp. 120-145, May 2011.

[30] E. Soubies, A. Chinatto, P. Larzabal, J. M. T. Romano, and L. Blanc-Féraud, "Direction-of-arrival estimation through exact continuous l0-norm relaxation," *IEEE Signal Processing Letters*, vol. 28, pp. 16-20, 2021.

[31] P. Ochs, A. Dosovitskiy, T. Brox, and T. Pock, "On iteratively reweighted algorithms for nonsmooth nonconvex optimization in computer vision," *SIAM Journal on Imaging Sciences*, vol. 8, no. 1, pp. 331-372, 2015.

[32] H. Krim, P. Forster, and J. Proakis, "Operator approach to performance analysis of root-music and root-min-norm," *IEEE Transactions on Signal Processing*, vol. 40, no. 7, pp. 1687-1696, 1992.

[33] R. J. Muirhead, "Aspects of multivariate statistical theory," in *Wiley Series in Probability and Statistics*, 1982.

[34] 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.

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

[36] —, "On the complementarity of sparse L0 and CEL0 regularized loss landscapes for DoA estimation," *Sensors*, vol. 21, no. 18, 2021.

[37] J. J. Moreau, "Proximité et dualité dans un espace hilbertien," *Bulletin de la société mathématique de France*, vol. 93, pp. 273-299, 1965.

[38] S. Becker and M.-J. Fadili, "A quasi-newton proximal splitting method," in *Neural Information Processing Systems*, 2012.

[39] S. Becker, J. Fadili, and P. Ochs, "On quasi-newton forward-backward splitting: Proximal calculus and convergence," *SIAM Journal on Optimization*, vol. 29, 01 2018.

[40] J. Hiriart-Urruty and C. Lemaréchal, *Convex Analysis and Minimization Algorithms*. Springer-Verlag, New-York.

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

[42] M. Lazzaretti, L. Calatroni, and C. Estatico, "Weighted-CEL0 sparse regularisation for molecule localisation in super-resolution microscopy with Poisson data," in *2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI)*, Nice, France, Apr. 2021.

[43] T. Aussaguès, A. Ferréol, A. Delmer, and P. Larzabal, "Whitening effects for ML-DoA estimation using a sparse representation of array
covariance,” in ICASSP 2025 - 2025 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2025.