A Subspace Estimator for Fixed Rank Perturbations of Large Random Matrices

Introduction

Parameter estimation algorithms based on the estimation of an eigenspace of the autocorrelation matrix of an observed multivariate time series are very popular in the areas of statistics and signal processing. Applications of such algorithms include the estimation of the angles of arrival of plane waves impinging on an array of antennas, the estimation of the frequencies of superimposed sine waves, or the resolution of multiple paths of a radio signal. Denoting by $N$ the signal dimension (e.g., the number of antennas) and by $n$ the length of the time observation window, the observed time series is represented by a $N\times n$ random matrix $\Sigma_{n}=X_{n}+P_{n}$ where $X_{n}$ and $P_{n}$ are respectively the so-called noise and signal matrices. In many applications, $P_{n}$ is represented as

We shall consider here “direction of arrival” vector functions $b(\varphi)$ that are typically met in the field of antenna processing. These functions are written

In practice, $\Pi$ is classically replaced with the orthogonal projection matrix $\widehat{\Pi}$ on the eigenspace associated with the $r$ largest eigenvalues of $\Sigma_{n}\Sigma_{n}^{*}$. Assuming $N$ is fixed and $n\to\infty$, and assuming furthermore that $S_{n}^{*}S_{n}$ converges to some matrix $O>0$ in this asymptotic regime, the $\Sigma\Sigma^{*}\xrightarrow{\text{a.s.}}BOB^{*}+I_{N}$ by the Law of Large Numbers (a.s. stands for almost surely). Hence, the random variable $\chi_{\text{classical}}(\varphi)=b(\varphi)^{*}\widehat{\Pi}b(\varphi)$ a.s. converges to $\chi(\varphi)$, and it is standard to estimate the arrival angles as local maxima of $\chi_{\text{classical}}(\varphi)$.

However, in many practical situations, the signal dimension $N$ and the window length $n$ are of the same order of magnitude in which case the spectral norm of $\widehat{\Pi}-\Pi$ is not small, as we shall see below. In these situations, it is often more relevant to assume that both $N$ and $n$ converge to infinity at the same pace, while the number of parameters $r$ is kept fixed. The subject of this paper is to develop a new estimator better suited to this asymptotic regime, and to study its first and second order behavior with the help of large random matrix theory.

In large random matrix theory, much has been said about the spectral behavior of $X_{n}X_{n}^{*}$ in this asymptotic regime, for a wide range of statistical models for $X_{n}$. In particular, it is frequent that the spectral measure of this matrix converge to a compactly supported limiting probability measure $\pi$, and that the extreme eigenvalues of $X_{n}X_{n}^{*}$ a.s. converge to the edges of this support. Considering that $\Sigma_{n}$ is the sum of $X_{n}$ and a fixed-rank perturbation, it is well-known that $\Sigma_{n}\Sigma_{n}^{*}$ also has the limiting spectral measure $\pi$ [2, Lemma 2.2]. However, the largest eigenvalues of $\Sigma_{n}\Sigma_{n}^{*}$ have a special behavior: Under some conditions, these eigenvalues leave the support of $\pi$, and in this case, their related eigenspaces give valuable information on the eigenspaces of $P_{n}$. This paper shows how the angles $\varphi_{k}$ can be estimated from these eigenspaces.

The problem of the behavior of the extreme eigenvalues of large random matrices subjected to additive or multiplicative low rank perturbations (often called “spiked models”) have received a great deal of interest in the recent years. In this regard, the authors of study the behavior of the extreme eigenvalues of a sample covariance matrix when the population covariance matrix has all but finitely many eigenvalues equal to one, a problem described in . Reference is devoted to the extreme eigenvalues of a Wigner matrix that incurs a fixed-rank additive perturbation. Fluctuations of these eigenvalues are studied in .

Recently, Benaych-Georges and Nadakuditi proposed in a powerful technique for characterizing the behavior of extreme eigenvalues and their associated eigenspaces for three generic spiked models: The models $X_{n}+P_{n}$ and $(I_{n}+P_{n})X_{n}$ when both $X_{n}$ and $P_{n}$ are Hermitian and $P_{n}$ is low-rank, and the model that encompasses ours $(X_{n}+P_{n})(X_{n}+P_{n})^{*}$ where $X_{n}$ and $P_{n}$ are rectangular. One feature of this approach is that it uncovers simple relations between the extreme eigenvalues and their associated eigenspaces on the one hand, and certain quadratic forms involving resolvents related with the non-perturbed matrix $X_{n}$ on the other. This makes the method particularly well-suited (but not limited to) the situation where $X_{n}$ is unitarily or bi-unitarily invariant, a situation that we shall consider in this paper. Indeed, in this situation, these quadratic forms exhibit a particularly simple behavior in the considered large dimensional asymptotic regime.

In this paper, we make use of the approach of to develop a new subspace estimator of the angles $\varphi_{k}$ based on the eigenspaces of the isolated eigenvalues of $\Sigma_{n}\Sigma_{n}^{*}$. We perform the first and second order analyses of this estimator that we call the “Spike MUSIC” estimator. Our mathematical developments differ somehow from those of and could have their own interest. They are based on two simple ingredients: The first is an analogue of the Poincaré-Nash inequality for the Haar distributed unitary matrices which has been recently discovered by Pastur and Vasilchuk , and the second is a contour integration method by means of which the first and second order analyses are done. The key step of the second order analysis of our estimator lies in the establishment of a Central Limit Theorem on the quadratic forms $b(\varphi_{i})^{*}\widehat{\Pi}_{i}b(\varphi_{i})$ where the $\widehat{\Pi}_{i}$ are the orthogonal projection matrices on certain eigenspaces of $\Sigma_{n}\Sigma_{n}^{*}$ associated with the isolated eigenvalues. The employed technique can easily be used to study the fluctuations of projections of other types of vectors on these eigenspaces.

We now state our general assumptions and introduce some notations.

We now state the general assumptions of the paper. Consider the sequence of $N\times n$ matrices $\Sigma_{n}=X_{n}+P_{n}$ where:

The dimensions $N,n$ satisfy: $N\leq n$, $n\to\infty$ and

(notation for this asymptotic regime: $n\to\infty$).

The following assumption on $X_{n}$ is widely used in the random matrix literature :

Matrices $X_{n}$ are random $N\times n$ bi-unitarily invariant matrices, i.e., each $X_{n}$ admits the singular value decomposition $X_{n}=L_{n}\Gamma_{n}R_{n}^{*}$ where $L_{n}$, the $N\times N$ matrix $\Gamma_{n}$ and $R_{n}$ are independent, $L_{n}$ is Haar distributed on the group ${\mathcal{U}}(N)$ of unitary $N\times N$ matrices, and $R_{n}$ is a $n\times N$ submatrix of a Haar distributed matrix on ${\mathcal{U}}(n)$.

We recall that the Stieltjes transform of a probability measure $\pi$ on the real line is the complex function

The quantity $\|X_{n}X_{n}^{*}\|$ a.s. converges to $\lambda_{+}$ as $n\to\infty$, where $\|\cdot\|$ denotes the spectral norm.

In the areas of signal processing and communication theory, the noise matrix $X_{n}$ satisfying Assumptions A2-A4 is such that $\sqrt{n}X_{n}$ is standard Gaussian - see for instance , .

We first make a general assumption on matrices $P_{n}$; it will be specified later, and adapted to the context of the MUSIC algorithm:

Matrices $P_{n}$ are deterministic with a fixed rank equal to $r$ for all $n$ large enough. Denoting by $P_{n}=U_{n}\Omega_{n}V_{n}^{*}$ a singular value decomposition of $P_{n}$, the matrix of singular values $\Omega_{n}=\operatorname*{diag}(\omega_{1,n},\ldots,\omega_{r,n})$ with $\omega_{1,n}\geq\omega_{2,n}\geq\cdots\geq\omega_{r,n}$ converges to

where $\omega_{1}>\cdots>\omega_{s}>0$ and $j_{1}+\cdots+j_{s}=r$.

The eigenvalues of $\Sigma_{n}\Sigma_{n}^{*}$ are $\hat{\lambda}_{1,n}\geq\hat{\lambda}_{2,n}\geq\cdots\geq\hat{\lambda}_{N,n}$. Associated eigenvectors will be denoted $\hat{u}_{1,n},\hat{u}_{2,n},\cdots,\hat{u}_{N,n}$. For $k\in\{1,\ldots,r\}$, we shall denote by $i(k)$ the index $i\in\{1,\ldots,s\}$ such that $j_{1}+\cdots+j_{i-1}<k\leq j_{1}+\cdots+j_{i}$. For $i=1,\ldots,s$, We shall denote by $\widehat{\Pi}_{i,n}$ the orthogonal projection matrix on the eigenspace of $\Sigma_{n}\Sigma_{n}^{*}$ associated with the eigenvalues $\hat{\lambda}_{k,n}$ such that $i(k)=i$, i.e., $\widehat{\Pi}_{i,n}=\sum_{k:i(k)=i}\hat{u}_{k,n}\hat{u}_{k,n}^{*}$ when this eigenspace is defined. Columns of $U_{n}$ (see A5) will be denoted $u_{1,n},\cdots,u_{r,n}$. Given $i$, the orthogonal projection matrix on the eigenspace of $P_{n}P_{n}^{*}$ associated with the eigenvalues $\omega_{k,n}^{2}$ such that $i(k)=i$ will be $\Pi_{i,n}=\sum_{k:i(k)=i}u_{k,n}u_{k,n}^{*}$. Indexes $n$ and $N$ will often be dropped for readability.

Paper organization

The paper is organized as follows. Section 2 is devoted to the mathematical preliminaries. The general approach is described in Section 3. The Spike MUSIC algorithm is presented in Section 4 along with a first order study of this algorithm. Fluctuations of the estimates of the $\varphi_{k}$ are studied in Section 5 under the form of a Central Limit Theorem.

Preliminary mathematical results

We shall need the two following results. The first one is well-known . The second result, due to Pastur and Vasilchuk, is the unitary analogue of the well-known Poincaré-Nash inequality.

Let $W=[w_{ij}]$ be a random matrix Haar distributed on ${\mathcal{U}}(n)$. Then

Given a small $\varepsilon_{1}>0$, let $O_{n}$ be the probability event

Let Assumption A2 holds true and let $u,v$ be two unit norm deterministic $N\times 1$ vectors such that $u^{*}v=0$. Then for any $z$ with $\Re(z)>\lambda_{+}+\varepsilon_{1}$,

Recall that $X=L\,\Gamma\,R^{*}$ by Assumption A2; let $D=(\Gamma^{2}-zI)^{-1}$; write:

We now proceed by induction; assume that the result is true until $p\geq 1$. Applying Lemma 2 to $\Phi_{i}^{(p+1)/2}$, we obtain:

Using again the induction hypothesis, we get:

Let Assumption A2 hold true; let $u,v$ be two unit norm deterministic vectors with respective dimensions $N\times 1$ and $n\times 1$. Then for any $z$ such as $\Re(z)>\lambda_{+}+\varepsilon_{1}$,

Recall the definition (3) of the set $O_{n}$ and assume that $\varepsilon_{1}$ is chosen such that $\min_{z\in{\mathcal{C}}}\Re(z)>\lambda_{+}+\varepsilon_{1}$; let

Fixed Rank Perturbations: First Order Behavior

We first recall a result on matrix analysis that can be found in [19, Th. 7.3.7]:

Given a $N\times n$ matrix $A$ with $N\leq n$, let ${\bf A}$ be the matrix:

Then $\sigma_{1},\cdots,\sigma_{N}$ are the singular values of $A$ if and only if $\sigma_{1},\cdots,\sigma_{N},-\sigma_{1},\cdots,-\sigma_{N}$ in addition to $n-N$ zeros are the eigenvalues of ${\bf A}$. Furthermore, a pair $(u,v)$ of unit norm vectors is a pair of (left,right) singular vectors of $A$ associated with the singular value $\sigma$ if and only if $\begin{bmatrix}u/\sqrt{2}\\ v/\sqrt{2}\end{bmatrix}$ is a unit norm eigenvector of ${\bf A}$ associated with the eigenvalue $\sigma$.

Along the ideas in , we now characterize the behavior of the largest eigenvalues of $\Sigma\Sigma^{*}$, and then focus on their eigenspaces.

We start with an informal description of the approach. By Lemma 6, $\lambda$ is an eigenvalue of $\Sigma\Sigma^{*}$ if and only if $\det(\boldsymbol{\Sigma}-\sqrt{\lambda}I)=0$ where ${\boldsymbol{\Sigma}}=\begin{bmatrix}0&\Sigma\\ \Sigma^{*}&0\end{bmatrix}$. Writing:

and assuming that $x>0$ is not a singular value of $X$, we have:

after noticing that $J=J^{-1}$. Using the formula for the inversion of a partitioned matrix (see )

whence for $n$ large enough, the isolated eigenvalues of $\Sigma\Sigma^{*}$ above $\lambda_{+}$ will coincide with the zeros of $\det\widehat{H}(\sqrt{x})$ that lie above $\lambda_{+}$. Under Assumptions A1-A5, Lemma 5 shows that $\widehat{H}(x)$ a.s. converges to

decreases from $g(\lambda_{+}^{+})=\lim_{x\downarrow\lambda_{+}}g(x)$ to zero on $(\lambda_{+},\infty)$. Let $\omega_{1}^{2}>\cdots>\omega_{q}^{2}$ be those among the diagonal elements of $O^{2}$ that satisfy $\omega_{i}^{2}>1/g(\lambda_{+}^{+})$. Equation $g(x)=\omega_{i}^{-2}$ will have a unique solution $x=\rho_{i}>\lambda_{+}$ for any $i=1,\cdots,q$, while it will have no solution larger than $\lambda_{+}$ for $i>q$. It is then expected that any eigenvalue $\hat{\lambda}_{k,n}$ of $\Sigma_{n}\Sigma_{n}^{*}$ for which $i(k)\leq q$ (remember the definition of $i(k)$ provided in the paragraph “Assumptions and Notations” in Section 1), will converge to $\rho_{i}$, while $\hat{\lambda}_{j_{1}+\cdots+j_{q}+1,n}\to\lambda_{+}$ almost surely.

These facts are formalized in the following theorem, shown in :

Let Assumptions A1-A5 hold true; let $q$ be the maximum index such that $\omega_{q}^{2}>1/g(\lambda_{+}^{+})$. Let $\rho_{i}$ be the unique real number $>\lambda_{+}$ satisfying $\omega_{i}^{2}g(\rho_{i})=1$ for $i=1,\cdots,q$. Then

In the case where $\sqrt{n}X$ is a standard Gaussian matrix, $\pi$ is the Marčenko-Pastur distribution with support $\operatorname*{supp}(\pi)=[\lambda_{-},\lambda_{+}]=[(1-\sqrt{c})^{2},(1+\sqrt{c})^{2}]$, and

for $x\in(\lambda_{+},\infty)$. After a few derivations, we obtain:

Assume $\sqrt{n}X$ is standard Gaussian. Let $q$ be the maximum index such that $\omega_{q}^{2}>\sqrt{c}$. Then

and $\hat{\lambda}_{j_{1}+\cdots+j_{q}+1,n}\xrightarrow[n\to\infty]{\text{a.s.}}(1+\sqrt{c})^{2}$.

We now turn our attention to the eigenspaces of the isolated eigenvalues.

Asymptotic behavior of certain bilinear forms.

Recall the definition of $s$ as provided in Assumption A5. Given $i\leq s$, assume that $\omega_{i}^{2}>1/g(\lambda_{+}^{+})$. Given two $N\times 1$ deterministic sequences of vectors $b_{1,n}$ and $b_{2,n}$ with bounded norms, we shall find here a simple asymptotic relation between $b_{1,n}^{*}\widehat{\Pi}_{i,n}b_{2,n}$ and $b_{1,n}^{*}\Pi_{i,n}b_{2,n}$, that will be at the basis of the Spike MUSIC algorithm. A close problem has been considered in . We consider here a different technique, based on a contour integration and on the use of Lemmas 3 and 4. This method lends itself easily to the first and second order analyses of the Spike MUSIC algorithm that we shall develop in the following sections.

Writing ${\bf b}_{i}=\begin{bmatrix}b_{i}\\ 0\end{bmatrix}$ with $i=1,2$, we have by virtue of Lemma 6:

where ${\mathcal{C}}_{i,n}$ is a positively oriented circle that encloses the only singular values $\sqrt{\hat{\lambda}_{k,n}}$ of $\Sigma_{n}$ for which $i(k)=i$. Recalling (4) and using Woodbury’s identity ([19, §0.7.4]) together with the fact that $J=J^{-1}$, we obtain:

Using (5), we obtain after a straightforward calculation:

Intuitively, the first integral is zero for $n$ large enough and the second is close to

The approximation $b_{1}^{*}\widehat{\Pi}_{i}b_{2}\simeq T_{i}$ will be justified rigorously below. For the moment, let us develop the expression of $T_{i}$. Defining the $r\times r$ matrices:

where the integers $j_{i}$ are defined in Assumption A5, we have

Let Assumptions A1-A5 hold true. For a given $i\leq s$, assume that $\omega_{i}^{2}>1/g(\lambda_{+}^{+})$. Let $(b_{1,n})$ and $(b_{2,n})$ be two sequences of deterministic vectors with bounded norms. Then

Then, with probability one, $b_{1}^{*}\widehat{\Pi}_{i}b_{2}=\widehat{T}_{i}$ for $n$ large enough. Indeed, on the set $O_{n}$ (as defined in (3)), the singular values of $\Sigma$ greater than $\sqrt{\lambda_{+}}+\varepsilon_{1}$ coincide with the poles of $\widehat{H}(z)$ which are greater than $\sqrt{\lambda_{+}}+\varepsilon_{1}$ by the argument preceding Theorem 1. On this set, the first integral on the right hand side (r.h.s.) of (9) is zero, and by Theorem 1, the second integral can be replaced with $\int_{\gamma_{i}}$ with probability one for $n$ large enough. By Lemma 5, the differences $\widehat{H}(z)-H(z)$, $\hat{a}_{1}(z)-a_{1}(z)$, and $\hat{a}_{2}(z)-a_{2}(z)$ a.s. converge to zero, uniformly on $\gamma_{i}$. Hence $\widehat{T}_{i}-T_{i}\xrightarrow{\text{a.s.}}0$. ∎

The Spike MUSIC Estimation Algorithm

In the area of signal processing, the positive real numbers $\omega_{i}^{2}$ are called the Signal to Noise Ratios (SNR) associated with the $r$ sources. Assumption A5 becomes:

Matrices $P_{n}$ of dimension $N\times n$ are deterministic and are written:

as $n\to\infty$, where $O$ is defined in Assumption A5, and ${\mathcal{O}}$ is the classical Landau notation.

The assumption over the speed of convergence of $S^{*}S$ will be needed only for the purpose of the second order analysis. It is satisfied by most practical systems met in the field of signal processing. We moreover observe that it is possible to relax the assumption that $O$ is diagonal at the expense of a more complicated second order analysis.

In order for the algorithm to be able to estimate the $r$ angles, it is necessary that the perturbation $P$ gives rise to $r$ isolated eigenvalues, a fact that is stated in the following assumption:

Recall the definition (6) of function $g$, let $\lambda_{+}$ as defined in A3 and let $g(\lambda_{+}^{+})=\lim_{x\downarrow\lambda_{+}}g(x)$. Let the $\omega_{i}$’s as defined in A5, then:

The Spike MUSIC algorithm goes like this. The localization function $\chi(\varphi)$ defined in the introduction is also written as $\chi(\varphi)=\sum_{i=1}^{s}b(\varphi)^{*}\Pi_{i}b(\varphi)$. Given $\varphi$, the results of the previous section (Theorems 1 and 2 with $b_{1}=b_{2}=b(\varphi)$) show us that:

is a consistent estimator of $\chi_{n}(\varphi)$ in the asymptotic regime described by A1. By searching for the maxima of $\hat{\chi}(\varphi)$, we infer that we obtain consistent estimates of the angles or arrival. Observe that this algorithm requires the knowledge of the Stieltjes Transform of the limit spectral measure of $XX^{*}$ (available if the statistical description of the noise is known) and the number $r$ of emitting sources. Notice that when this number is unknown, it can be estimated along the ideas described in e.g. . We now perform the first order analysis of this algorithm.

First order analysis of the Spike MUSIC algorithm

We now formalize the argument of the previous paragraph and we push it further to show the consistency “up to the order $n$” of the Spike MUSIC estimator. We shall need this speed to perform the second order analysis (Lemma 9 below).

Let Assumptions A1-A6 hold true. Then for all $k=1,\cdots,r$, there exists a local maximum $\hat{\varphi}_{k,n}$ of $\hat{\chi}_{n}(\varphi)$ such that

The proof of this theorem is performed in two steps. With an approach similar to the one used in Section 3, we first prove that $\hat{\chi}(\varphi)-\chi(\varphi)\xrightarrow{\text{a.s.}}0$, and the convergence is uniform on $\varphi\in[0,\pi/D]$ (Proposition 1 below). Next, following the technique of , we prove that this uniform a.s. convergence leads to Theorem 3.

Beware that $\hat{a}^{*}$ and $a^{*}$ are not the Hermitian adjoints of $\hat{a}$ and $a$ (see the footnote associated to Eq. (10)).

By Theorem 1 and the continuity of $\zeta$ on $(\lambda_{+},+\infty)$, the first term at the r.h.s. goes to zero a.s. and uniformly in $\varphi$. Consider the second term. Let $\gamma_{i}$ be a small enough positively oriented circle which does not meet $\operatorname*{supp}(\pi)\cup\{\sqrt{\rho_{1}},\cdots,\sqrt{\rho_{s}}\}$ and such that only $\sqrt{\rho_{i}}\in\operatorname*{Int}(\gamma_{i})$. Since $\hat{\lambda}_{k}\xrightarrow{\text{a.s.}}\rho_{i(k)}$,

Recalling Eq. (12), it will therefore be enough to prove that

where $R$ is the radius of $\gamma_{i}$ and where

Since $\|H^{-1}\|$, $\max_{\varphi}\|a\|$ and $\max_{\varphi}\|\hat{a}\|$ are bounded on $\gamma_{i}$, $e(z,\varphi)$ satisfies on this path

By Lemma 5 and the fact that $\|H^{-1}\|$ is bounded on $\gamma_{i}$, the term $\|\widehat{H}^{-1}-H^{-1}\|=\|\widehat{H}^{-1}(H-\widehat{H})H^{-1}\|$ converges to zero uniformly on $\gamma_{i}$ with probability one. To obtain the result, we prove that $\|\hat{a}-a\|\xrightarrow{\text{a.s.}}0$ and that this convergence is uniform on $(z,\varphi)\in\gamma_{i}\times[0,\pi/D]$. Let us focus on the first term $zu_{1}^{*}(Q(z^{2})-m(z^{2})I)b(\varphi)$ of $\hat{a}-a$, where we recall that $u_{1}$ is the first column of $U$. Since $\|b(\varphi)\|=\|u_{1}\|=1$,

With probability one, the second term converges to zero on $\gamma_{i}$, and the convergence is uniform (along the principle of the proof of Lemma 5). Since

for every $(z_{1},\varphi_{1})$, $(z_{2},\varphi_{2})$ in $\gamma_{i}\times[0,\pi/D]$. Therefore, it will be enough to prove that

where $A_{n}$ contains $n$ regularly spaced points in $\gamma_{i}$ and $B_{n}$ contains $n^{2}$ regularly spaced points in $[0,\pi/D]$. This can be obtained from Lemma 3 with $p=9$, Markov inequality and Borel Cantelli’s lemma. The other terms of $\hat{a}-a$ can be handled similarly. ∎

We now prove Theorem 3 by following the ideas of . To that end, we need the following lemma, proven in :

Let $(c_{N})$ be a sequence of real numbers belonging to a compact of $[-1/2,1/2]$ and converging to $c$. Let

where $\operatorname*{sinc}$ stands as usual for sine cardinal.

We start by observing that $\chi(\varphi)=d(\varphi)^{*}(B^{*}B)^{-1}d(\varphi)$ where $B$ is the matrix defined in A6 and where $d(\varphi)=\begin{bmatrix}b(\varphi_{k})^{*}b(\varphi)\end{bmatrix}_{k=1}^{r}$. By Lemma 7, $B^{*}B\to I_{r}$, hence $\chi(\varphi)-\|d(\varphi)\|^{2}\to 0$.

In the remainder of the proof, we shall stay in the probability one set where the uniform convergence in the statement of Proposition 1 holds true. Taking $k=1$ without loss of generality, we shall show that any sequence $\hat{\varphi}_{1,n}$ for which $\hat{\chi}(\hat{\varphi}_{1,n})$ attains its maximum in the closure of a small neighborhood of $\varphi_{1}$ satisfies $N(\hat{\varphi}_{1,n}-\varphi_{1})\to 0$. Given a sequence of such $\hat{\varphi}_{1,n}$, assume we can extract a subsequence $\hat{\varphi}_{1,n^{*}}$ such that $N|\hat{\varphi}_{1,n^{*}}-\varphi_{1}|\to\infty$. In this case, Lemma 7 and the observations made above on the structure of $\chi(\varphi)$ show that $\chi(\hat{\varphi}_{1,n^{*}})\to 0$. Since $\max_{\varphi}|\hat{\chi}(\varphi)-\chi(\varphi)|\to 0$, $\hat{\chi}(\hat{\varphi}_{1,n^{*}})\to 0$. But $\hat{\chi}(\varphi_{1})\to\chi(\varphi_{1})=1$, which contradicts the fact that $\hat{\varphi}_{1,n^{*}}$ maximizes $\hat{\chi}$. Hence the sequence $N(\hat{\varphi}_{1,n^{*}}-\varphi_{1})$ belongs to a compact. Assume $N(\hat{\varphi}_{1,n^{*}}-\varphi_{1})\not\to 0$. If we take a further subsequence of the latter that converges to a constant $d\neq 0$, then by Lemma 7, $\hat{\chi}$ converges to $\operatorname*{sinc}(d)^{2}<1$ along this subsequence, which also raises a contradiction. This proves the theorem.∎

Second Order Analysis of the Spike MUSIC Estimator

In order to perform the second order analysis, we also assume:

The main result of this section is the following:

Let Assumptions A1-A8 hold true. Then the estimates $\hat{\varphi}_{k,n}$ satisfy

When $\sqrt{n}X$ is standard Gaussian, plugging the r.h.s. of (7) into this expression leads after some derivations to:

If $\sqrt{n}X$ is standard Gaussian and if $\sqrt{n}(c_{n}-c)\to 0$, the convergence (16) holds true with

Recalling that $\omega_{i}^{2}>\sqrt{c}$ is the condition for the existence of a corresponding isolated eigenvalue (Corollary 1), we observe that the estimator variance for $\varphi_{k}$ goes to infinity as the corresponding $\omega_{i}^{2}$ decreases to $\sqrt{c}$. At the other extreme, this variance behaves like $6c^{-2}D^{-2}\omega_{i}^{-2}$ as $\omega_{i}^{2}\to\infty$. It is useful to notice that this asymptotic variance coincides with the Cramér-Rao bound for estimating $\varphi_{k}$ . In other words, the Spike MUSIC estimator is efficient at high SNR when the noise matrix is standard Gaussian.

In order to illustrate the convergence and the fluctuations of the Spike MUSIC algorithm, we simulate a radio signal transmission satisfying Assumptions A1-A8. We consider $r=2$ emitting sources located at the angles $0.5$ and $1$ radian, and a number of receiving antennas ranging from $N=5$ to $N=50$. The observation window length is set to $n=2N$ (hence $c=0.5$). The noise matrix $X_{n}$ is such that $\sqrt{n}X_{n}$ is standard Gaussian. The source powers are assumed equal, so that the matrix $O$ given by Equation (2) is written $O=\omega I_{2}$, and the Signal to Noise Ratio for any source is $\text{SNR}=10\log_{10}\omega^{2}$ decibels. In Figure 2, the SNR is set to $10$ dB, and the empirical variance of $\hat{\varphi}_{1,n}-\varphi_{1}$ (red curve) is computed over $2000$ runs. The variance provided by Corollary 2 is also plotted versus $N$. We observe a good fit between the variance predicted by Corollary 2 and the empirical variance after $N=15$ antennas.

In Figure 3, the variance is plotted as a function of the SNR, the number of antennas being fixed to $N=20$. The empirical variance is computed over $5000$ runs. The Cramér-Rao Bound is also plotted. The empirical variance fits the theoretical one from $\text{SNR}\approx 6$ dB upwards.

Proof of Theorem 4.

We start with some additional notations and definitions. Matrix $B=\begin{bmatrix}b(\varphi_{1}),\ldots,b(\varphi_{r})\end{bmatrix}$ will be often written as $B=[b_{1},\ldots,b_{r}]$ or in block form as $B=\begin{bmatrix}B_{1},\ldots,B_{s}\end{bmatrix}$ where $B_{i}$ has $j_{i}$ columns. We shall also write $B^{\prime}=\begin{bmatrix}b^{\prime}(\varphi_{1}),\ldots,b^{\prime}(\varphi_{r})\end{bmatrix}$ and $B^{\prime\prime}=\begin{bmatrix}b^{\prime\prime}(\varphi_{1}),\ldots,b^{\prime\prime}(\varphi_{r})\end{bmatrix}$ where $b^{\prime}(\varphi)$ and $b^{\prime\prime}(\varphi)$ are respectively the first and second derivatives of $b(\varphi)$. We shall also use the short hand notations $B^{\prime}=[b_{1}^{\prime},\ldots,b_{r}^{\prime}]$ and $B^{\prime\prime}=[b_{1}^{\prime\prime},\ldots,b_{r}^{\prime\prime}]$. Matrix $B^{\perp}=[b_{1}^{\perp},\ldots,b_{r}^{\perp}]$ will be defined by the equation

Finally, if $x_{n},y_{n}$ are random sequences, we denote by $x_{n}\asymp y_{n}$ the convergence $x_{n}-y_{n}\xrightarrow[]{\mathcal{P}}0$.

We now state some preliminary results. In the following, we say that the complex random vector $\eta$ is governed by the law ${\cal CN}(0,R)$ where $R$ is a nonnegative Hermitian matrix if the real vector $\begin{bmatrix}\Re(\eta)\\ \Im(\eta)\end{bmatrix}$ has the law {\cal N}\Bigl{(}0,\frac{1}{2}\begin{bmatrix}\Re(R)&-\Im(R)\\ \Im(R)&\Re(R)\end{bmatrix}\Bigr{)}. The following proposition, whose proof is postponed to A, is crucial:

Assume $t$ is even. Given real numbers $\rho_{1},\ldots,\rho_{t/2}$ all strictly greater than $\lambda_{+}$, the $t\times 1$ random vector

converges in distribution towards $\mathcal{CN}(0,R)$ with

Writing $Q-mI=(Q-\alpha I)+(\alpha-m)I$, and similarly for $\widetilde{Q}$, we obtain:

Assume in addition that Assumption A8 is satisfied. Then

Intuitively, tightness of $\xi_{n}$ leads to the tightness of the $\sqrt{n}(\hat{\lambda}_{k,n}-\rho_{i(k)})$. This is formalized by the following proposition, proven in B:

Assume the setting of Theorem 4. Then the sequences $\sqrt{n}(\hat{\lambda}_{k,n}-\rho_{i(k)})$ are tight for $1\leq k\leq r$.

Let Assumptions A5 and A6 hold true. Then the following convergences hold true:

where $\Pi_{B_{i}}$ is the orthogonal projection matrix on the column space of $B_{i}$.

Recall the definitions (13) and (14) of $\hat{\chi}$ and $\zeta$. In most of the proof, we shall focus on $\sqrt{n}(\hat{\varphi}_{1,n}-\varphi_{1})$. Recalling that $\hat{\chi}^{\prime}(\hat{\varphi}_{1})=0$ and performing a Taylor-Lagrange expansion of $\hat{\chi}^{\prime}$ around $\varphi_{1}$, we obtain

where $\hat{\chi}^{(3)}$ is the third derivative of $\hat{\chi}$ and where $\bar{\varphi}_{1}\in[\varphi_{1}\wedge\hat{\varphi}_{1},\varphi_{1}\vee\hat{\varphi}_{1}]$. Hence

We start by characterizing the asymptotic behavior of the denominator of this equation:

Assume that the setting of Theorem 4 holds true. Then,

Theorem 1 along with the continuity of $\zeta$ on $(\lambda_{+},\infty)$, and Theorem 2 show that

by the first, fourth and fifth assertions of Lemma 8. By the same lemma,

Hence $n^{-2}\hat{\chi}^{\prime\prime}(\varphi_{1})\to-c^{2}D^{2}/6$.

Furthermore, it is easily seen that $n^{-3}\hat{\chi}^{(3)}(\bar{\varphi}_{1})$ is bounded. Since $n(\hat{\varphi}_{1}-\varphi_{1})\xrightarrow{\text{a.s.}}0$ by Theorem 3, $n^{-2}(\hat{\varphi}_{1}-\varphi_{1})\hat{\chi}^{(3)}(\bar{\varphi}_{1})\xrightarrow{\text{a.s.}}0$, which establishes the result. ∎

We now turn to the numerator $n^{-1/2}\hat{\chi}^{\prime}(\varphi_{1})=2n^{-1/2}\sum_{k=1}^{r}\zeta(\hat{\lambda}_{k})\Re\left(b_{1}^{*}\hat{u}_{k}\hat{u}_{k}^{*}b^{\prime}_{1}\right)$, and start with the following lemma:

Assume that the setting of Theorem 4 holds true. Then

and where the deterministic circle $\gamma_{i}$ encloses $\rho_{i}^{1/2}$ only and:

Recall the definition of $\hat{\chi}$ as given in (13). A direct computation yields:

Recall that $r$ and $s$ are fixed and independent from $n$ by A5. We start by showing that

Since $\sqrt{n}(\zeta(\hat{\lambda}_{k,n})-\zeta(\rho_{i(k)}))$ is tight as a corollary of Proposition 3, it will be enough to prove that $n^{-1}\Re\left(b_{1}^{*}\hat{u}_{k}\hat{u}_{k}^{*}b^{\prime}_{1}\right)\to 0$ in probability for every $k$. By the definition (17) of $B^{\perp}$, we have

and by Lemma 8, $b_{1}^{*}\Pi_{i(k)}b_{1}\ (b_{1}^{\perp})^{*}\Pi_{i(k)}b_{1}^{\perp}\ \to\ 0$ (consider alternatively the cases $i(k)=1$ and $i(k)>1$) which proves (20).

Now, applying (9) and (15), and taking up an argument used in the proof of Theorem 2, we have

with probability one for $n$ large. On the other hand, recalling (12), we have

Write $\widehat{H}(z)=H(z)+E(z)$ and $\hat{a}(z,\varphi)=a(z,\varphi)+e(z,\varphi)$. To be more specific,

Write $e^{\prime}_{\varphi}(z,\varphi)=\partial e(z,\varphi)/\partial\varphi$. For a given $z\in\gamma_{i}$, $\widehat{H}^{-1}=H^{-1}-H^{-1}EH^{-1}+{\mathcal{O}}(\|E\|^{2})$. This suggests the following development

where the terms $q_{i}$ are “higher order terms” that appear when we expand the r.h.s. of (19). We first handle the terms $X_{k,i}$’s, then $q_{i}$.

Writing $U_{n}=\begin{bmatrix}U_{1,n}\cdots U_{s,n}\end{bmatrix}$ and $V_{n}=\begin{bmatrix}V_{1,n}\cdots V_{s,n}\end{bmatrix}$ where both $U_{i,n}$ and $V_{i,n}$ have $j_{i}$ columns, and recalling (11), we have

Due to the bounded character of $\|n^{-1}b^{\prime}\|$ and to Corollary 3, $X_{1,i}$ is tight for every $i$. By Lemma 8,

Thanks to Corollary 3 and Lemma 8, $\Re(G_{ii}(\rho_{i}))\xrightarrow{\mathcal{P}}0$. The same can be said about $G^{\prime}_{ii}(\rho_{i})$ after a simple modification of Proposition 2 and Corollary 3. In conclusion,

These are the higher order terms that appear when we expand the right hand side of (19). We shall work here on one of these terms, namely

and show that $\varepsilon\xrightarrow{\mathcal{P}}0$. The other higher order terms can be handled similarly. Writing $z=\sqrt{\rho_{i}}+R\exp(2\imath\pi\theta)$ on the circle $\gamma_{i}$, we have

where $K$ is a constant whose value can change from line to line, but which remains independent from $n$. Let $\phi$ be a function from $$to a normed vector space. If$\phi$is twice differentiable on$(0,1)$, then it is known that$\|\phi(1)-\phi(0)-\phi^{\prime}(0)\|\leq\sup_{t\in(0,1)}0.5\|\phi^{\prime\prime}(t)\|$.

Setting $\phi(t)=(H+tE)^{-1}$ and recalling that $\hat{H}=H+E$, we have $\phi(1)=\hat{H}$, $\phi(0)=H$ and $\phi^{\prime\prime}(t)=(H+tE)^{-1}E(H+tE)^{-1}E(H+tE)^{-1}$, hence

Consider any element of $E_{1}$, for instance $zu_{1}^{*}(Q(z^{2})-\alpha(z^{2})I)u_{1}$. By Lemma 3,

which shows that $\sqrt{n}\int_{0}^{1}\|E_{1}\|^{2}d\theta\xrightarrow{\mathcal{P}}0$.

Final derivations

Write $\hat{\boldsymbol{\chi}}^{\prime}=\begin{bmatrix}\hat{\chi}^{\prime}(\varphi_{1}),\ldots,\hat{\chi}^{\prime}(\varphi_{r})\end{bmatrix}$. Generalizing the previous argument to all the $\varphi_{k}$ and gathering the results, we obtain

By Lemma 8, matrix $A=\begin{bmatrix}V_{i(k)}U_{i(k)}^{*}b_{k}\end{bmatrix}_{k=1}^{r}$ satisfies $A^{*}A\to I_{r}$. Recall from the same lemma that $B^{*}B\to I_{r}$, $(B^{\perp})^{*}B^{\perp}\to I_{r}$ and $(B^{\perp})^{*}B\to 0$. Hence, Proposition 2 can be applied to the r.h.s. of this expression, and $n^{-1/2}\hat{\boldsymbol{\chi}^{\prime}}$ converges in law to

It remains to recall Lemmas 9 and 10 to terminate the proof of Theorem 4.

Appendix A Proof of Proposition 2

where $\widetilde{Z}[1;N]$ is $\widetilde{Z}$ truncated to its first $N$ rows. By the Law of Large Numbers, $N^{-1}Z^{*}Z\to I_{t}$ and $n^{-1}\widetilde{Z}^{*}\widetilde{Z}\to I_{t}$ almost surely. Hence, if we show that the multidimensional random variables $A_{k,n}=N^{-1/2}Z^{*}(D_{k}-N^{-1}\operatorname*{tr}D_{k})Z$ and $B_{k,n}=N^{-1/2}Z^{*}C_{k}\widetilde{Z}[1;N]$ are tight for $k=1,\ldots,t/2$, and

converges in law towards $\mathcal{CN}(0,R)$, the second result of Proposition 2 is proven. From A3 and A4,

Observe that covariance matrix of $\bar{\eta}_{n}$ conditional to $\Gamma_{n}$ converges almost surely to $R$. Moreover, thanks to A4, it is easy to see that the Lyapunov condition

is satisfied for any $a>0$, hence $\bar{\eta}_{n}\xrightarrow{\mathcal{L}}\mathcal{CN}(0,R)$ which completes the proof of Proposition 2.

Appendix B Sketch of the proof of Proposition 3.

For $k=1,\dots,r$, let $\bar{\rho}_{k,n}$ be the solutions of the equation $\omega_{k,n}^{2}g(\rho)=1$, where we recall that the $\omega_{k,n}^{2}$ are the diagonal elements of matrix $\Omega_{n}$. Then, by a simple extension to the case $r\geq 1$ of the proof of [9, Th. 2.15], one can show that the sequences $\sqrt{n}(\hat{\lambda}_{k,n}-\bar{\rho}_{k,n})$ are tight. To obtain the result, we show that $\sqrt{n}(\bar{\rho}_{k,n}-\rho_{i(k)})={\mathcal{O}}(1)$. Since $g$ is decreasing, this amounts to showing that $\sqrt{n}(\omega_{k,n}^{2}-\omega_{i(k)}^{2})={\mathcal{O}}(1)$. Since the non zero eigenvalues of $PP^{*}$ coincide with those of $B^{*}B\,S^{*}S$, it will be enough to prove that $\sqrt{n}(B^{*}B\,S^{*}S-O)={\mathcal{O}}(1)$. It is clear that $B^{*}B=I_{r}+n^{-1}A$ where $\sup_{n}\|A\|<\infty$, hence $\sqrt{n}(B^{*}BO-O)\to 0$. By the last item in Assumption A6, $\sqrt{n}B^{*}B(S^{*}S-O)={\mathcal{O}}(1)$, and the proposition is shown.

Appendix C Proof of Lemma 8.

References