Probability of Resolution of Partially Relaxed DML an Asymptotic Approach

The Partial Relaxation approach has recently been proposed to solve the Direction-of-Arrival estimation problem [1], [2]. In this paper, we investigate the outlier production mechanism of the Partially Relaxed Deterministic Maximum Likelihood (PR-DML) Direction-of-Arrival estimator using tools from Random Matrix Theory. Instead of applying a single source approximation to multi-source estimation criteria, which is the case for the MUSIC algorithm, the conventional beamformer, or the Capon beamformer, the Partial Relaxation approach accounts for the existence of multiple sources using a non-parametric modification of the signal model. In this paper, an accurate description of the probability of resolution for the PR-DML estimator is provided by analyzing the asymptotic behavior of the PR-DML cost function, assuming that both the number of antennas and the number of snapshots increase without bound at the same rate. The finite dimensional distribution of the PR-DML cost function is shown to be Gaussian in this asymptotic regime and this result is used to compute the probability of resolution.


I. INTRODUCTION
Direction-of-Arrival (DoA) estimation has been a major area of research in the past, mainly due to its wide spread applications in radar, sonar, seismology or electronic surveillance and mobile communication [3]- [6]. Several high resolution algorithms, such as Multiple Signal Classification (MUSIC) [7], the minimum variance method of Capon [8], Estimation of Signal Parameters via Rotational Invariance Technique (ESPRIT) [9] have been proposed [10], [11]. However, the performance of conventional "low-cost" methods strongly degrades when two or multiple sources are closely spaced [12], [13]. This is due to the fact that conventional spectral search based approaches ignore the presence of interfering sources and therefore treat multi-source scenarios as single source scenarios.
The Partial Relaxation (PR) framework was introduced to overcome the aforementioned disadvantages of the conventional spectral-based DoA methods [1], [2]. Instead of ignoring the presence of multiple sources, the PR approach considers both the signal impinging from the current direction of interest as well as the interfering ones. To reduce the computational demand, the manifold structure of the undesired signal components is relaxed, whereas the manifold structure of the desired signal component is kept unchanged. The multi-dimensional optimization problem reduces to a one-dimensional problem This work was supported in part by the Catalan Government under grant 2017-SGR-01479, the Spanish Government under grant RTI2018-099722-B-I00 that admits a simple spectral based grid search applicable to any array geometry.
The main objective of this paper is the performance characterization of the recently introduced PR-DML technique in the threshold region, whereby both the number of samples per antenna and the Signal-to-Noise Ratio (SNR) take moderate values. This region is typically characterized by a systematic appearance of outliers in the DoA estimates, which are mainly caused by the incapability of resolving closely spaced sources.
The performance characterization of PR DoA estimators has received little attention in the literature so far. Up to now, researchers have focused on the performance analysis of traditional Deterministic Maximum Likelihood (DML) and Stochastic Maximum Likelihood (SML), see e.g. [14], where the probability of resolving closely spaced sources for finite values of snapshots and SNR is studied. Some additional insights were given in [15], where the outlier production mechanism was related to the asymptotic eigendecomposition of the observation covariance matrix. In [16] the probability of resolution of DML and SML was investigated by studying the stochastic behavior of the corresponding cost functions under the assumption that the received signals are Gaussian random variables. The approach in [16] is asymptotic in both the number of antennas and the number of snapshots. In this paper, we follow a similar approach as in [14], [16] and compute the probability of resolution of the PR-DML technique by deriving the stochastic behavior of the PR-DML cost function for Gaussian distributed observations. Tools from Random Matrix Theory (RMT) are used to compute the asymptotic equivalent of the PR-DML cost function as well as the asymptotic covariance where both the number of snapshots and the number of antennas are large quantities but their quotient converges to a fixed finite value.
The paper is organized as follows. The signal model is introduced in Section II followed by the PR-DML DoA estimation technique in Section III. The asymptotic stochastic behavior of the PR-DML cost function is given in Section IV and the corresponding sketch of the proof is given in Section V. In Section VI an expression for the probability of resolution is provided. Simulation results are presented in Section VII. Finally, Section VIII concludes this paper.
II. SIGNAL MODEL Let us consider an antenna array equipped with M sensors and K impinging narrowband signals that satisfy M > K. The source signal at time instant n is denoted by s(n) = [s 1 (n),...,s K (n)] T ∈ C K×1 . The corresponding DoAs of the signals are denoted by θ = [θ 1 ,...,θ K ] T . Furthermore, the fullrank steering matrix is given by A(θ) = [a(θ 1 ),...,a(θ K )] ∈ C M ×K where a(θ i ) ∈ C M denotes the sensor array response for the i-th impinging signal. The number of sources K is assumed to be known. The received baseband signal x(n) ∈ C M at the n-th time instant is given by where N denotes the number of snapshots and n(n) ∈ C M the sensor noise. Assuming that signal and noise variables are statistically independent zero-mean circularly symmetric Gaussian distributed, the covariance matrix of the received signal R ∈ C M ×M is given by where R s = E s(n)s H (n) denotes the covariance of the transmitted signal and σ 2 I M is the noise covariance matrix.
Since the true covariance matrix is unavailable in practice, the sample covariance matrixR = 1 III. PARTIALLY RELAXED DETERMINISTIC MAXIMUM LIKELIHOOD In the framework of PR, not only the signals from the desired directions but also from the interfering directions are considered [1]. However, the structure of the interfering signals is relaxed and consequently the computational complexity of the multi-source criteria is greatly reduced. Unlike in conventional DML and SML DoA estimation criteria, the steering matrix A is not described by the fully structured array manifold Instead, A describes the partially relaxed array manifold which still retains some geometric structure of the sensor array [2]. Applying the PR approach to the conventional DML method yields where P ⊥ [a(θ),B] denotes the orthogonal projection matrix and K argmin a(θ)∈A1 the K arguments in θ that minimize the cost. At first, the objective function in (2) is minimized with respect to B. A closed-form solution for B is obtained and substituted back into the objective function. The concentrated objective function yields [1], [2] where the eigenvalues of the modified sample covariance matrixR(θ) = P ⊥ a (θ)RP ⊥ a (θ) are sorted in non-descending order 0 =λ 1 (θ) ≤ ··· ≤λ M (θ). The estimatesθ are given by K arguments that correspond to the K deepest minima of the concentrated objective function in (3). An efficient implementation of the PR-DML method is provided in [1].
The distinct true eigenvalues of the modified true covariance matrix and their corresponding multiplicities are given by K m (θ), for m = 0,...,M (θ). The number of distinct true eigenvalues of the modified true covariance matrix R(θ) amounts toM (θ)+1 and the sum of the multiplicities satisfies The non-necessarily Hermitian M ×M positive definite square root of the modified true covariance matrix R(θ) can also be expressed using the singular value decomposition where U r (θ) ∈ C M ×Kr(θ) and V r (θ) ∈ C M ×Kr(θ) generate the left and right orthonormal basis of R(θ) 1/2 .

IV. ASYMPTOTIC BEHAVIOR OF THE PARTIALLY RELAXED DETERMINISTIC MAXIMUM LIKELIHOOD COST FUNCTION
In the following, the first and second order moments of the cost function in (3)  In RMT it is well known that under all the previously mentioned assumptions, the empirical eigenvalue distribution of R(θ) is almost surely close to an asymptotic non-random distribution which is absolutely continuous with density q M (x,θ) [17]. With increasing number of snapshots and therefore decreasing c = M/N , q M (x,θ) tends to concentrate around the true eigenvalues forming different eigenvalue clusters. The number of eigenvalue clusters increases with decreasing c as clusters begin to split [18]. Assuming there are S distinct clusters, the support of the cluster is given by the set of . Furthermore, it can be observed that each modified true and distinct eigenvalue γ m (θ) is associated to only one cluster. However, one cluster may be associated to multiple true eigenvalues which results in a non-bijective correspondence [18]. For sufficiently small c, there always exist exactly as many clusters as distinct true eigenvalues M (θ) [18].
In order to distinguish between the eigenvalues that are considered by the PR-DML cost function in (3) and the remaining ones it is crucial that the (M −K+1)-th modified sample eigenvalue asymptotically splits from the (M −K+2)-th one, which can be formalized as follows. We assume that there exists an integer m(θ) such that M −K+1 = m(θ) r=0 K r (θ), and the cluster [x − s (θ),x + s (θ)] associated to the eigenvalue γ m(θ) (θ) separates from the one associated to γ m(θ)+1 (θ) in the asymptotic eigenvalue distribution ofR(θ). Under these conditions, the first and second order asymptotic behavior of the PR-DML cost-function are given below.

A. First Order Asymptotic Behavior
It can be shown that the PR-DML cost functionη(θ) in (3) becomes asymptotically close to its deterministic counterpart in the sense that |η(θ)−η(θ)| → 0 almost surely pointwise in θ as M,N → ∞ at the same rate.

B. Second Order Asymptotic Behavior
In the following, the nature of the fluctuation of the PR-DML cost function is studied in the asymptotic regime. Let us consider the two real-valued L×1 vectorŝ whereθ = θ 1 ,...,θ L T denotes a set of L points within the field of view Θ. The asymptotic covariance matrix of the random vectorη(θ) can be computed by solving the realvalued integral Im I p,q ω x,θ p dx, (7) for p,q = 1,...,L using numerical integration techniques such as the Riemann sum or the Simpsons' rule. Taking z ∈ C + ≡ {z ∈ C : Im(z) > 0}, we define ω(z,θ) as the unique solution on C + of Furthermore, the integrand in (7) is given by where ω 1 = ω z,θ p . The complex-valued scalars ϕ r (ω 1 ), r = 1,...,M θ q are sorted according to their real-part in ascending order and given by the solutions to the polynomial equation

C. Final Result: Asymptotic Distribution
Under the previously mentioned assumptions and as M,N → ∞, M/N → c and 0 < c < ∞, the random vector converges in distribution to a multivariate standardized Gaussian distribution. The asymptotic mean of the random vector η(θ) in (5) is denoted byη(θ) in (6) and the corresponding L×L asymptotic covariance matrix is given by Γ θ /M 2 in (7).

V. SKETCH OF THE PROOF
Let us introduce the complex random function for z ∈ C + and its asymptotic equivalent [17] m(z,θ) = ω(z,θ) z Consider the limit of ω(z,θ) andm(z,θ) as z goes to the real axis and let S(θ) be the positive support 1 of Im(m(z,θ)) on R. The PR-DML cost function in (3) can be equivalently expressed asη where C(θ) denotes a negatively oriented contour enclosing the M −K+1 smallest sample eigenvalues ofR(θ) only. By replacingm(z,θ) with its asymptotic equivalentm(z,θ) one obtains the asymptotic first order behaviorη(θ) in (4). Furthermore, the asymptotic second order behavior can be computed according to the following Theorem.
Theorem 1. According to the above definitions, consider an L×L matrix Γ θ with entries
Proof. The proof can be obtained using the approach in [19]. Also see Theorem 2 of [16].
Alternatively, the complex double contour integral in (10) can be expressed as where Using conventional residue calculus the closed-form expression for I p,q (ω 1 ) in (8) is obtained which can be substituted back into (12). By carefully choosing a parametrization of the contour C ω θ p the simplified expression in (7) is obtained.
VI. PROBABILITY OF RESOLUTION To analyze the probability of resolution we study a scenario with two sources K = 2 located at θ 1 and θ 2 . Considering the minimization problem in (3), the cost functionη(θ) is evaluated at the true DoAs θ 1 and θ 2 and the mid-angle (θ 1 +θ 2 )/2. We declare loss of resolution if the cost function evaluated at the mid-angle is lower than evaluated at both true DoAs [20]. The constraint to declare resolution can compactly be written as u Hη (θ) < 0 where By utilizing the previously obtained stochastic behavior of the PR-DML cost function in (9), the asymptotic probability density function (pdf) of the test quantity u Hη (θ) can easily be computed for M,N → ∞ at the same rate M u H Γ θ u −1/2 u Hη θ −u Hη θ → N (0,1) [21]. The probability of resolution is therefore obtained by computing the cumulative distribution function [22] P res = Pr u Hη where f u Hη (θ) (x) denotes the asymptotic pdf of the test quantity u Hη (θ).

VII. SIMULATION RESULTS
In this Section, the predicted probability of resolution is compared to the simulated one. All simulations are carried out for 10000 Monte-Carlo runs. We consider two uncorrelated and closely spaced sources at θ = [45°,50°] T and a Uniform Linear Array (ULA) equipped with M = 10 sensors. The transmitted signals are zero-mean and statistically independent with unit power. The SNR is given by SNR = 1/σ 2 n . The separation boundary is defined as the smallest SNR that provides separation between the eigenvalue clusters associated to the m θ l -th true eigenvalue and larger adjacent true eigenvalues smaller than the separation boundary it is not possible to distinguish between the eigenvalues that are considered by the PR-DML cost-function and the remaining ones. Figures 1 and 2 depict the probability of resolution versus the SNR for N = 10 and N = 100 snapshots. As expected, the probability of declaring resolution increases with increasing number of snapshots. Clearly, we are able to predict the probability of resolution remarkably well in both scenarios.

VIII. CONCLUSION
In this paper we have investigated the asymptotic behavior of the PR-DML DoA estimator under the setting of RMT where both the number of snapshots N and the number of sensors M go to infinity at the same rate. The finite dimensional distribution of the PR-DML cost function has been derived and shown to be asymptotically Gaussian. Furthermore, the probability distribution of the PR-DML method was used to characterize the probability of resolution in the threshold region, where the generation of outliers causes a total performance breakdown.