Source Enumeration via Toeplitz Matrix Completion

This paper addresses the problem of source enumeration by an array of sensors in the presence of noise whose spatial covariance structure is a diagonal matrix with possibly different variances, referred to non-iid noise hereafter, when the sources are uncorrelated. The diagonal terms of the sample covariance matrix are removed and, after applying Toeplitz rectification as a denoising step, the signal covariance matrix is reconstructed by using a low-rank matrix completion method adapted to enforce the Toeplitz structure of the sought solution. The proposed source enumeration criterion is based on the Frobenius norm of the reconstructed signal covariance matrix obtained for increasing rank values. As illustrated by simulation examples, the proposed method performs robustly for both small and large-scale arrays with few snapshots, i.e. small-sample regime.


I. INTRODUCTION
Source enumeration is a classical problem in array signal processing which consists in estimating the number of signals received by an array of sensors. It has applications in numerous fields such as wireless communications, radar, biomedical and geophysical signal processing [1], [2]. Classical approaches to solve this problem are based on information-theoretic criteria [3]- [6], which use order fitting rules based on functions of the eigenvalues of the sample covariance matrix (SCM) penalized by the model complexity. These criteria are derived under a white noise model assumption and their performance degrades when the number of snapshots is relatively small in comparison to the number of antennas, the so-called smallsample regime. Different methods have been proposed for order estimation in the small-sample regime, such as those based on random matrix theory [7]- [9], an exponential fitting test of eigenvalues [10], and bootstrapping [11]. These methods, however, are designed for white noise, and usually provide poor results under non-iid noises. Recently, a subspace averaging (SA) order estimation method has been proposed for white noise [12], [13], and then extended to account for noniid noise models in [14]. The SA method provides good results with low sample support, but it requires high-dimensional scenarios, i.e. large arrays to work satisfactorily since otherwise the averaging procedure is not effective. When the noise covariance matrix is diagonal with unknown elements (a model that becomes relevant in situations with uncalibrated receivers or due to hardware nonidealities), estimating the number of sources is equivalent to estimating the number of common factors in a multivariate factor analysis problem, for which several methods have been proposed in the statistics literature [15], [16]. Algorithms to maximize the likelihood function for this problem can be found in [17], [18]. But again, these methods perform poorly in the small-sample regime.
In this paper we propose a method for enumerating sources in the presence of non-iid noises that performs satisfactorily in a wide number of scenarios, including the small-sample regime both with large and short arrays, and low signal-to-noise-ratio (SNR) environments. The method obtains an estimate of the signal covariance matrix that is i) Toeplitz, and ii) low-rank. To this end, we use a low-rank matrix completion (MC) method that includes an additional regularization term to enforce a Toeplitz structure in the solution. The MC algorithm takes as input a denoised version of the SCM where the elements along the main diagonal have been removed. An additional denoising process known as Toeplitz rectification [19], [20] is used to find a better estimate of the signal covariance matrix. Finally, we propose an order estimation criterion which is based on the Frobenius norm of the matrices reconstructed for increasing rank values.
Low-rank matrix completion methods have been used before in array signal processing problems. In [21], for instance, MC is used for direction of arrival (DOA) estimation, when the number of sources exceeds the number of sensors. An iterative reweighted nuclear norm minimization method is used in [22] for DOA estimation with nested arrays. The case of non-iid noises is considered in [23], where matrix completion algorithms are used to reconstruct the zeroed entries of the SCM along its diagonal. All these methods, however, address the DOA estimation problem and assume that the number of sources is known. Differently from these works, in this paper we exploit the eigenvalue sparsity provided by MC algorithms to devise a novel order estimation criterion.
Notation. The superscripts (·) T and (·) H denote transpose and Hermitian, respectively. The trace and Frobenius norm of a matrix B will be denoted, respectively, as tr(B) and ||B|| F , whereas | · | denotes the norm of a complex number. Furthermore, x ∼ CN n (0, R) denotes a complex Gaussian vector in C n with zero mean and covariance R.

II. PROBLEM STATEMENT
Let us consider K narrowband signals impinging on a uniform linear array (ULA) with M antennas. The received signal is (1) where a(θ k ) = 1, e −jθk , ... ,e −jθk(M −1) T is the M × 1 complex array response for the kth source, with unknown DOA θ k . The signals are assumed to be uncorrelated and are modelled as s[n] ∼ CN K (0, Ψ), where Ψ is a diagonal matrix. The noises e[n] ∼ CN M (0, Σ), are assumed to be temporarily white, uncorrelated across antennas with different variances at each sensor (non-iid). Therefore, the noise covariance matrix is Σ = diag σ 2 1 , σ 2 2 , . . . , σ 2 M , where σ 2 m is the noise variance at the m-th sensor. The noise variances are modeled as uniformly distributed independent random variables: where σ 2 is a common noise variance and 0 ≤ ≤ 1 allows us to control the spatial non-whiteness of the noise. Notice that for = 0 the noise is spatially white with covariance matrix Σ = σ 2 I. From the signal model (1), the covariance matrix is where , then the sample covariance matrix iŝ The source enumeration problem consists in estimating K from X orR.

III. ORDER ESTIMATION VIA TOEPLITZ MATRIX COMPLETION
The proposed order determination criterion builds upon the reconstruction of the signal covariance matrix R s in (2) for increasing values of its rank. To this end, we apply a matrix completion (MC) approach, which is described in the following subsection. Then, we introduce the proposed order fitting rule, which is based on the Frobenius norm of the reconstructed matrices.

A. Toeplitz Matrix Completion
The signal covariance matrix R s is Toeplitz and low-rank. These two properties are not fulfilled by the sample covariancê R, which is symmetric (but non-Toeplitz) and full-rank. In addition, since Σ in (2) is diagonal, the off-diagonal terms of R are unaffected by the noise covariance matrix. Therefore, the diagonal entries ofR will be more affected by the noise. Matrix completion algorithms can then be used to reconstruct the low-rank signal covariance matrix from the off-diagonal terms ofR.
Let Ω = {(i, j) : i = j, i = 1, . . . , M, j = 1, . . . , M} be the set of indices for the off-diagonal entries ofR. According to [24], we can recover R s by solving where P Ω denotes the projection operator that sets to zero the entries with indices not belonging to Ω and leaves the rest unchanged, and ||R s || * denotes the nuclear norm of R s . Due to the fact that we are dealing with a limited number of snapshots, the nonzero entries in P Ω (R) are still noisy [19] and thus recovering the full R s via (4) might yield unreliable estimates in the low-sample regime. Therefore, we propose to use a denoising step called Toeplitz rectification [19], [20], [25] before applying matrix completion, and enforce the Toeplitz structure in the reconstruction process.
An unbiased estimator of sample covariance with Toeplitz structure is obtained by averaging its entries along each subdiagonal as [25] This process reduces the noise in the off-diagonal terms and provides a better reconstruction of the signal covariance matrix.
After Toeplitz correction, we proceed to perform matrix completion by means of matrix factorization [26]. Since the signal covariance matrix is symmetric, we can factorize it as R s = WW H , where W ∈ C M ×p and p is a fixed value that limits the rank of the reconstructed matrix. Then, using the identity ||R s || * = min Rs=WW H ||W|| 2 F [26], R s can be estimated by solving the following optimization problem where μ is a regularization parameter and the constraint restricts the solution to a set T , which we define as the set of Toeplitz matrices in C M ×M . Note that the term weighted by μ regulates the nuclear norm of the solution, hence leading to a sparse eigenvalue distribution as μ grows. An approximate solution to (6) can be obtained by introducing the Toeplitz constraint in the form of an additional regularization term as followŝ where w H i is the i-th row vector of W, and α is a regularization scalar. The solutionR mc =ŴŴ H can be obtained by iteratively optimizing over each w i until convergence, at a cost of O(p 3 ) per iteration.

B. Order Estimation Criterion
The main insight for the proposed order estimation criterion is that, due to the eigenvalue sparsity enforced by the MC algorithm, as long as p (the rank used in the factorization) is larger than K, the reconstructed signal covariance matrixR mc should not change significantly. This intuition is corroborated in Fig. 1a which shows how the Frobenius norm ofR mc changes with p. It can be observed that the norm grows until p = K, but once p exceeds K it is almost constant. This suggests to use the difference function D(p) defined as whereR mc (p) denotes the reconstructed signal covariance matrix at a particular value of p. Since D(p) will take very small values for p > K, finding the position at which this sharp decline change occurs will yield an order estimate. To this end, we propose the following criterion where p max is an overestimation of K and δ is a small constant to avoid numerical issues. Fig. 1b shows for different K, where each line has been normalized by its maximum value to enhance visibility. We observe that the peaks are positioned at p = K, thus evidencing that the chosen criterion is adequate. A summary of the order estimation via Toeplitz matrix completion (TMC) algorithm is shown in Algorithm 1.

IV. SIMULATION RESULTS
We compare the performance of the proposed TMC with some representative methods for order estimation in the presence of non-iid noises and/or with low sample support, which are briefly reviewed in the following.
• Subspace Averaging (SA). By exploiting the shift invariance property of ULAs, the SA method estimates the low-dimensional signal subspace (and its dimension) as the average of a collection of subspaces extracted from consecutive sub-arrays [12], [13], [27]. The SA method, which was originally proposed for white noise and low sample support, was generalized to account for non-iid noises in [14] through a majority-vote approach.
whereR ML denotes the Maximum Likelihood (ML) estimate for a covariance matrix with the required structure (low-rank plus diagonal) for a fixed order k, and ν k = M + k(2M − k) is the number of free-adjusted real parameters [28]. Although a closed-form ML estimate under non-iid noises is not possible, it can be obtained by using iterative algorithms [17], [18], [29]. • Linear Shrinkage-MDL (LSMDL). The standard MDL method proposed by Wax and Kailath [5] under the assumption of white noise iŝ where a(k) and g(k) are the geometric and the arithmetic mean, respectively, of the M − k smallest eigenvalues of R. The linear shrinkage MDL (LS-MDL) [8], replaces the noise eigenvalues in the MDL criterion by a linear shrinkage. For all simulations we assume that K uncorrelated narrowband signals with a separation of Δ θ are impinging on an uniform linear array with M half-wavelength separated antennas. For TMC, we use μ = M 2 and α = M 10 . The signalto-noise-ratio is SNR= 10 log tr(Rs) tr(Σ) . Sources have unequal powers so that Ψ = σ 2 s diag(1, 0.8, 0.6), where σ 2 s is the common signal variance. We observe that TMC provides better performance than competing methods. Thanks to the Toeplitz rectification and the denoising ofR, the TMC method reliably detects the rank of the signal covariance matrix at lower SNRs than the rest of methods and hence yields a more robust source enumeration method.   In the last experiment the robustness of the TMC method against the non-whiteness parameter, , is examined in different scenarios. Since the diagonal terms ofR are eliminated as a pre-processing step, it is expected that the TMC results will not be affected by changes in . This behavior can be observed for two different scenarios: i) K = 2 and N = 200 and ii) K = 6 and N = 50, in Fig. 4, which shows P D with respect to SNR for M = 50, Δ θ = 10 • and equal power sources. It can be noticed that for both the scenarios, results are almost unaffected by . Let us recall that = 0 represents the white noise case; therefore, TMC is robust under both uncorrelated non-iid and iid noises.

V. CONCLUSIONS
We addressed the problem of source enumeration when the noise covariance matrix is diagonal with unknown entries and with relatively few snapshots. The diagonal terms of the sample covariance matrix, which are most affected by this noise model, are eliminated and the off-diagonal terms are further denoised by Toeplitz rectification. The low-rank and Toeplitz signal covariance matrix is then reconstructed with matrix completion techniques. We have shown that the Frobenius norm of the signal covariance matrix reconstructed by the proposed denoising+MC technique provides an accurate order determination criterion. The method performs robustly for iid and non-iid noises, as well as for small and large arrays, in the small-sample regime.