Fluctuations of spiked random matrix models and failure diagnosis in sensor networks

I Introduction

In the field of fault detection and diagnosis, one of the elementary requests is the fast, reliable and computationally light identification of a system failure. In dynamical scenarios, these systems are composed of several fluctuating parameters whose evolutions are tracked by a mesh of sensors reporting successive correlated and noisy data measurements to a central decision unit. With the growth in size and complexity of such systems, it becomes increasingly difficult for decision units to process simultaneously and at a low computational cost the augmenting load of reported measurements. Examples of such systems are the recent cognitive radio networks and smart grid technologies . In the former, multiple cooperative wireless communication devices, referred to as the secondary network, exchange sensed data in order to decide collectively which communication bandwidths are left unused by the licensed, also called primary, network users. Fast detection of sudden changes, e.g. new primary user communications, is here demanded to minimize the interference generated by secondary users. In the smart-grid framework, a large dimensional graph of interconnected electricity producers, transportation systems, and consumers evolve in real-time, their behaviour being reported by diverse sensors such as voltage phasor measurements at the nodes of the electricity grid to regional controllers. Fast detection of link and node failures is requested in this scenario to minimize the risk of cascaded failures leading to regional blackouts . There exists a rich literature on failure detection, diagnosis and change-point estimation, ranging from off-line detection methods of uncorrelated data to fast change detection methods in time correlated signals . Subspace methods were in particular proposed to detect system changes from modifications in the eigenstructure of sampled covariance matrices for dynamical systems . In this article, we propose a novel subspace approach to solve the problem of off-line detection and identification of local failures from independent or linearly time-correlated samples.

The approach under consideration here follows the theory of large dimensional random matrices. Precisely, we consider the setting where both $N$ and $n$ grow large and such that $c_{N}=N/n\to c$, with $0<c<1$, as $N,n\to\infty$. Under this assumption, we develop asymptotic results on the extreme eigenvalues and associated eigenvectors of a certain family of random matrices to provide novel subspace methods for failure detection and localization. Our interest is on random matrices of the spiked model type, introduced by Johnstone , specifically here of matrices modeled as $\Sigma=(I_{N}+P)^{\frac{1}{2}}X$, where $X$ is a left-unitarily invariant random matrix and $P$ is a rank-$r$ Hermitian matrix with $r\ll N$. Such matrix models have been largely studied in the recent random matrix literature, very often in the special case where $X$ is a standard Gaussian matrix, which refers in this article to a random matrix with independent $\mathcal{CN}(0,1/n)$ entries. In , for $X$ a standard Gaussian matrix, it is first shown that there exists a natural mapping between the extreme (empirical) eigenvalues of $\Sigma\Sigma^{\ast}$ and the (population) eigenvalues of $P$. It is then proved that, almost surely, the extreme empirical eigenvalues converge to deterministic limits in the asymptotic setting, found either at the edges of the support of the Marc̆enko-Pastur law , i.e. the (almost sure) weak limit of the eigenvalue distribution of $XX^{\ast}$, or away from them, depending on the corresponding population eigenvalues. This induces a phase transition having important consequences on fault detectability in sensor networks . This observation is extended to the non-Gaussian case and generalized to other spiked models in . The fluctuations of the extreme eigenvalues are studied with different approaches depending on whether the limiting eigenvalues are found at the edge or outside the support of the Marc̆enko-Pastur law. When at the edge, it is proved successively in that the (centered and scaled) limiting eigenvalue has Tracy-Widom fluctuations. When outside the support, those fluctuations are linked to the distribution of the eigenvalues of GUE matrices, as shown in . In the specific case where the spiked eigenvalues of $P$ have unit multiplicity, the fluctuations are Gaussian.

In this article, the properties of the extreme eigenvalues in a spiked model will be used to provide failure detection tests, in the same line as . For failure localization, the information on the eigenvalue position can be used to reduce the number of hypotheses $K$. However, tests solely based on the limiting properties of the eigenvalues will turn out to be inefficient to discriminate the remaining hypotheses and we therefore develop novel results on the eigenspaces associated to these eigenvalues. In , it is shown, in the real Gaussian case, that the projection of the eigenvectors associated with the extreme empirical eigenvalues of $\Sigma\Sigma^{\ast}$ on the subspace of the corresponding population eigenvectors of $P$ has a positive limiting norm, which is close to $1$ for $c$ small. This remark is extended in to the non-Gaussian case. This property is the basis of our novel failure diagnosis method. However, the fluctuations of the eigenvector projections, fundamental here to derive test statistics for failure localization, have never been derived before for either the Gaussian or the non-Gaussian cases. The main mathematical result of this article, Theorem 4, provides the joint fluctuations of the eigenvalues and eigenspace projections for the eigenvalues found away from the limiting support of $XX^{\ast}$. Our proof technique is largely inspired by . We also use some tools from and . Based on these results, we suggest an original framework for local failure detection and identification in large sensor networks.

The remainder of this document unfolds as follows. Section II introduces elementary examples of sensor networks for which local failures translate into small rank perturbations of the identity matrix. Section III reminds important notions of random matrix theory and introduces the main mathematical results of this article. Practical application algorithms along with simulations are then carried out in Section IV. Finally, Section V concludes the article.

II Detection and localization of local failures

To motivate the study of the fluctuations of extreme eigenvalues and eigenvectors of sample covariance matrices in the context of local failures in large dimensional sensor networks, we introduce in the following two basic examples of sensor network failure scenarios, which can all be modeled as small rank perturbations of the identity matrix, as well as related engineering applications.

In case of failure of sensor $k$, $y(k)$, the $k^{th}$ entry of $y$, will start suddenly to return noisy outputs inconsistent with the model (1). Assuming this noise Gaussian with zero mean and variance $\sigma_{k}^{2}$ and denoting $y^{\prime}$ the observations of the network with failure at sensor $k$, we can write

Therefore, $y^{\prime}$ is Gaussian (as the sum of Gaussian variables) with zero mean and variance

Denoting $s=R^{-{\frac{1}{2}}}y^{\prime}$, we have

Therefore, the population covariance matrix ${\rm E}[ss^{\ast}]$ is a perturbation of the identity matrix by

Notice that the image of $P_{k}$ is included in the subspace ${\rm Span}(R^{-{\frac{1}{2}}}e_{k},R^{-{\frac{1}{2}}}HH^{\ast}e_{k})$ and is therefore at most of dimension two. Generalizing the above to $M$ node failures at nodes $k_{1},\ldots,k_{M}$, the vector $s$ is now such that

with $E=[e_{k_{1}},\ldots,e_{k_{M}}]$, $\Lambda=\operatorname{diag}(\sigma_{k_{1}},\ldots,\sigma_{k_{M}})$, where now (II-A) becomes

II-B Sudden parameter change

Consider again the elementary model (1) and now assume that, instead of a sensor failing, $\theta(k)$, the $k^{th}$ entry of $\theta$, experiences a sudden change in mean and variance. The resulting observation $y^{\prime}$ can be modeled as

Denoting $R=HH^{\ast}+\sigma^{2}I_{N}$ as in the previous example and taking $s=R^{-{\frac{1}{2}}}y^{\prime}$, we finally have

which is a rank-$1$ perturbation of the identity matrix by the matrix

with $\beta_{k}=\mu_{k}^{2}+(1+\alpha_{k})^{2}-1$. Note that, in this scenario, the eigenvector of $P_{k}$ associated with the non-zero eigenvalue is independent of $\mu_{k}$ and $\alpha_{k}$. For practical applications, this has the interesting advantage that simple localization can be performed even if $\mu_{k}$ and $\alpha_{k}$ are unknown. This is further discussed in Section IV.

The derivation above generalizes to sudden changes of multiple parameters. If the means and variances for the sensors $k_{1},\ldots,k_{M}$ are modified simultaneously with respective parameters $\mu_{k_{1}},\ldots,\mu_{k_{M}}$ and $\alpha_{k_{1}},\ldots,\alpha_{k_{M}}$, then

with $E=[e_{k_{1}},\ldots,e_{k_{M}}]$ and $\Lambda=\operatorname{diag}(\beta_{k_{1}},\ldots,\beta_{k_{M}})$, $\beta_{k_{i}}=\mu_{k_{i}}^{2}+(1+\alpha_{k_{i}})^{2}-1$, which is a rank-$M$ perturbation of the identity matrix by the matrix

Note that, contrary to the one-dimensional case, the eigenvectors of $P_{k_{1},\ldots,k_{M}}$ depend here explicitly on the parameters $\mu_{k_{i}}$ and $\alpha_{k_{i}}$.

In the following section, we introduce the novel detection and localization framework and we discuss engineering applications referencing the examples described in this section.

II-C Detection and localization in sensor networks

The natural approach to detect and identify a failure event in a sensor network upon the observations $s_{1},\ldots,s_{n}$ is to systematically perform a maximum likelihood test on the $K+1$ hypotheses $\mathcal{H}_{0},\ldots,\mathcal{H}_{K}$, with $\mathcal{H}_{k}$ defined as the event $s\sim\mathcal{CN}(0,I_{N}+P_{k})$. However, this optimal approach has some intrinsic limitations. From a computational aspect, evaluating the probability of each hypothesis $k$ requires to evaluate the term $\operatorname{tr}\Sigma^{\ast}(I_{N}+P_{k})^{-1}\Sigma$, an operation whose cost is of order $\mathcal{O}(N^{3})$ (which can be brought down to $\mathcal{O}(N^{2})$ using matrix inversion lemmas). When the number of hypotheses $K$ and the system size $N$ are large, these operations become extremely demanding. Pre-calculus of the inverses $(I_{N}+P_{k})^{-1}$ also requires possibly large memory storage.

Since the node failure information is entirely captured by the perturbation matrix $P_{k}$, we provide in the following a suboptimal test relying on the properties linking $P_{k}$ to the observation matrix $\Sigma$, for large system dimensions $(N,n)$. Precisely, based on recent advances in the field of large dimensional random matrix theory , we provide a two-step approach to successively (i) decide on the existence of a failure from the location of the extreme eigenvalues of $\Sigma\Sigma^{\ast}$ and (ii) identify the failure event from eigenspace projections. This diagnosis framework relies on the asymptotic statistics of these extreme eigenvalues and eigenspace projections. This subspace approach has multiple advantages compared to the optimal hypothesis testing method discussed above. From a computational aspect, step (i) requires to determine the eigenvalues of $\Sigma\Sigma^{\ast}$, hence a singular value decomposition. This step already provides sufficient information for step (ii) to become computationally cheap: on the one hand, the position of the extreme eigenvalues of $\Sigma\Sigma^{\ast}$ may be used to reduce the set $\mathcal{H}_{1},\ldots,\mathcal{H}_{K}$ to a possibly small subset of consistent hypotheses; on the other hand, for the remaining hypotheses, the localization test will merely consists in the characterization of eigenvector projections, an operation of computational cost $\mathcal{O}(N)$. No matrix inverse needs to be computed and only the eigenvectors and non-zero eigenvalues of $P_{k}$ need to be stored. The technique also has the advantage to be consistent in its usage of eigenvalues and eigenspace projections to perform hypothesis tests. Finally, as will be discussed in Section IV, the framework can be extended to account easily for unknown failure amplitudes, which would be much more involved from a maximum-likelihood approach.

The following section is dedicated to the study of the asymptotic eigenvalue and eigenspace projection statistics as the dimensions of the matrix $\Sigma$ grow large.

III Main results

The derivation arguments found in this section follow the ideas of , , and . In our proofs, we shall also borrow some of the arguments of whose context is close to ours.

We start by summarizing the major notations and facts needed here. We consider a generic small rank perturbation model and define

In the remainder of the paper, we shall consider the asymptotic regime where $n\to\infty$ and $N/n\to c\in(0,1)$. The notation $n\to\infty$ will henceforth refer to this asymptotic regime.

The probability law of $X$ is invariant by left multiplication by a deterministic unitary matrix.

Thanks to the left unitary invariance of $X$, $Q(z)$ writes as $Q(z)=W(\Lambda-zI_{n})^{-1}W^{*}$ where $\Lambda$ is the matrix of eigenvalues of $XX^{\ast}$, $W$ is a unitary random matrix Haar distributed on its unitary group, and $W$ and $\Lambda$ are independent.

We have $\|XX^{\ast}\|\overset{\rm a.s.}{\longrightarrow}b$ and $(\|(XX^{\ast})^{-1}\|)^{-1}\overset{\rm a.s.}{\longrightarrow}a$.

The most classical model of a matrix $X$ that satisfies A1-A3 is when $X$ is standard Gaussian, i.e. with independent ${\mathcal{C}N}(0,1/n)$ elements, as introduced in the system models of Section II-A and Section II-B. For this model, the limiting probability distribution $\pi$ is the well known Marc̆enko-Pastur distribution . Its Stieltjes transform $m(z)$ is given by

The unitary invariance of $X$ is the basis of the following important lemma, shown in using an inequality of which involves Haar unitary matrices:

We now start our analysis of the extreme eigenvalues and eigenspace projections of $\Sigma\Sigma^{\ast}$ by studying the first order behavior.

III-B First order behavior

after noticing that $I_{N}-(I_{N}+P)^{-1}=P(I_{N}+P)^{-1}$. Therefore, if $x$ is an eigenvalue of $\Sigma\Sigma^{\ast}$ but not of $XX^{\ast}$, it must cancel the rightmost determinant. This determinant can be further rewritten

From the identity $U^{\ast}(I_{N}+U\Omega U^{\ast})^{-1}=(I_{r}+\Omega)^{-1}U^{\ast}$, we then have

(take $u$ and $v$ in Lemma 1 as any couple of columns of $U$, take $p>2$ and use Borel-Cantelli’s lemma ). We therefore expect the solutions of the equation $\det H(x)=0$ which are outside $[a,b]$ to coincide with the limits of the isolated eigenvalues of $\Sigma\Sigma^{\ast}$.

having a unique real solution $\rho_{i}$ satisfying $\rho_{i}>b$ if and only if $h(b^{+})+({1+\omega_{i}})/{\omega_{i}}<0$.We denote by $x^{+}$ and $x^{-}$ any quantity infinitesimally greater and smaller than the real $x$, respectively. When $\omega_{i}<0$, (6) has a unique solution $0<\rho_{i}<a$ if and only if $h(a^{-})+({1+\omega_{i}})/{\omega_{i}}>0$. We therefore have the following result, for which a rigorous proof is found in :

Assume A1-A3. Let $p$ be zero or the maximum index such that $\omega_{p}>0$ and $h(b^{+})+({1+\omega_{p}})/{\omega_{p}}<0$. For $i=1,\ldots,p$, let $\rho_{i}$ be the unique solution of (6) such that $\rho_{i}>b$. Then,

Let $q$ be $t+1$ or the minimum index such that $\omega_{q}<0$ and $h(a^{-})+({1+\omega_{q}})/{\omega_{q}}>0$. For $i=q,\ldots,t$, let $\rho_{i}$ be the unique solution of (6) such that $\rho_{i}<a$. Then,

In the remainder of the article, the variables $\omega_{1},\ldots,\omega_{p}$ and $\omega_{q},\ldots,\omega_{t}$ satisfying the conditions of Theorem 1 will be said to satisfy the separation condition.

When $X$ is standard Gaussian, applying Theorem 1 shows after some simple derivations the following result:

Consider the setting of Theorem 1. Assume additionally that $X$ is standard Gaussian. Let $p$ be zero or the maximum index for which $\omega_{p}>\sqrt{c}$ and $q$ be $t+1$ or the minimum index such that $\omega_{q}<-\sqrt{c}$. Then

Corollary 1 implies that, for $\omega_{i}$ sufficiently far from zero (either positive or negative) or, equivalently, for $c$ sufficiently small, the spectrum of $\Sigma\Sigma^{\ast}$ exhibits $j_{i}$ eigenvalues outside the support $S_{\pi}$ of the Marc̆enko-Pastur law $\pi$ which all converge to $\rho_{i}$. For failure detection purposes, upon observation of $\Sigma$, we may then test the null hypothesis $\Sigma=X$ (call it hypothesis $\mathcal{H}_{0}$) against the hypothesis $\Sigma=(I_{N}+P)^{\frac{1}{2}}X$ (call it hypothesis $\bar{\mathcal{H}}_{0}$), depending on whether eigenvalues of $\Sigma\Sigma^{\ast}$ are found outside $S_{\pi}$. Depending on the scenario, for $c$ small enough, it may be that a mere evaluation of the number of eigenvalues outside the support suggests the number of simultaneous failures in the sensor network. This is the case of the two failure scenarios described in Section II-A and Section II-B. However, the information on the extreme eigenvalues of $\Sigma\Sigma^{\ast}$, if sufficient for failure detection purposes, is usually not good enough to perform accurate failure localization. This is because different failure scenarios, characterized by different perturbation matrices $P$, may exhibit very similar eigenvalues. Also, if the failure amplitude is a priori unknown, then eigenvalues are in general irrelevant to discriminate between failure hypotheses; see the application Section IV-C. In such scenarios, we then need to consider eigenspace properties of $P$. This is the target of the following section.

III-B2 Projections on eigenspaces

Using Woodbury’s matrix identity, we have

By Assumption A3 and Theorem 1, with probability one for all large $n$, the first term on the right-hand side is zero, while the second is equal to

where $\gamma_{i}$ is a deterministic positively oriented circle away from $[a,b]$ enclosing $\rho_{i}$ but none of the $\rho_{j}$, $j\neq i$. Using Lemma 1 in conjunction with the analyticity properties of the integrand, one can show that $\hat{a}_{1}(z)^{\ast}\widehat{H}(z)^{-1}\hat{a}_{2}(z)$ converges uniformly to $a_{1}(z)^{\ast}H(z)^{-1}a_{2}(z)$ on $\gamma_{i}$ in the almost sure sense, where

It results that $b_{1}^{\ast}\widehat{\Pi}_{i}b_{2}-T_{i}\overset{\rm a.s.}{\longrightarrow}0$, where

Details can be found in in a similar situation. Let us find the expression of $T_{i}$. Noticing that

In particular, we find after some derivations:

Under the assumptions of Theorem 2, let $X$ be standard Gaussian. Then

This result is consistent with derived in the real Gaussian case for eigenvalues with unit multiplicity.

Theorem 4 and Corollary 4 provide an interesting characterization of the eigenspaces of $P$ through limiting projections in the large dimensional setting. In the context of local failure in large sensor networks, it is therefore possible to detect and diagnose one or multiple failures by comparing eigenspace projection patterns associated with each failure type. Precisely, an appropriate diagnosis consists in determining the most likely failure type among all hypothetical failures, given the extreme eigenvalues and associated eigenspace projections of $\Sigma\Sigma^{\ast}$. To this end though, not only first order limits but also second order behaviour need be characterized precisely. This is the target of the following section.

III-C Second order behavior

Let $X$ be standard Gaussian, then if $0<\omega_{i}<\sqrt{c}$,

The tools used to derive Theorem 3 are much different from those exploited here and will not be discussed. Note that in , an extension to the case where $X$ may be correlated is provided but only considers the fluctuations of the largest eigenvalue. Similar to , Theorem 3 will be used to derive tests to decide on the presence of eigenvalues outside the support of the Marc̆enko-Pastur law. For failure detection purposes in sensor networks, this will be used to declare a failure prior to diagnose the fault. Then, to diagnose a failure, second order statistics of both eigenvalue and eigenspace projections when the separation property arises are needed. This is the aim of the remainder of the section.

We now turn to the second order analysis of the eigenspectrum of $\Sigma\Sigma^{\ast}$ when $\omega_{i}$ satisfies the separation condition, and when $X$ is only assumed to satisfy A1-A3. We first need the following additional assumption:

For practical purposes, we shall also assume:

Each $\omega_{i}$, $1\leq i\leq t$, satisfies the separation condition.

The main result of this section is the following theorem.

where $m^{(3)}$ is the third derivative of $m$. Consider the matrices

Theorem 4 provides a very general expression of the joint limiting fluctuations of both eigenvalues and eigenspace projections. It is particularly interesting to note that the fluctuations of $(V_{i,n},L_{i,n})$ are asymptotically independent across $i$.

Consider the setting of Theorem 4. Assume in addition that $j_{i}=1$ for all $i$. Then

After some calculus, in the standard Gaussian case, we further have:

Under the assumptions of Corollary 3, if $X$ is a standard Gaussian matrix, then $D(\rho_{i})R(\rho_{i})D(\rho_{i})^{\ast}=C(\rho_{i})$, where

Due to its simple expression, Corollary 4 is particularly handy to use in the context of failure diagnosis when hypothetical failures are characterized by distinct values of $\omega_{i}$, as will be shown in Section IV.

The remainder of this section is devoted to the proof of Theorem 4. We start with the following lemma, which deals with the asymptotic behavior of the $V_{i,n}$. This lemma will be proved in Appendix A-A.

We now consider the isolated extreme eigenvalues. In order to study the asymptotic behavior of these eigenvalues, we shall adapt to our situation the approach of . For $i=1,\ldots,t$, consider real numbers $x_{1}(i)>y_{1}(i)>x_{2}(i)>y_{2}(i)>\cdots>y_{j_{i}}(i)$. Since the separation condition is satisfied by assumption for each $i=1,\ldots,t$, the equation $\det\widehat{H}(x)=0$ has $r$ roots outside $[a,b]$ with probability one for all $n$ large. Therefore, we have the equivalence relation

where $\beta_{i}=\omega_{i}/(1+\omega_{i})$. Then

for every finite sequence $(x_{1},\ldots,x_{p})$.

for $i=1,\ldots,t$, similarly to the elements of $L_{n}$. From Lemma 2, Lemma 3, and the discussion preceding Lemma 3, we have

as $n\to\infty$, for arbitrary rectangles $B$ and arbitrary rectangles $C$ specified at the left hand side of (11). Observe that

In order to terminate the proof of Theorem 4, we shall make use of the following lemma, that we state in a slightly more general form than needed here.

Assume A1-A4. Let $f_{1},\ldots,f_{t}$ and $g_{1},\ldots,g_{t}$ be real functions analytical on a neighborhood of $[a,b]$. Let $S_{n}$ be the $t$-uple of random matrices

For $i=1,\ldots,t$, define the covariance matrices

Then $S_{n}$ converges in distribution towards

where the matrices $(M_{1,1},M_{2,1},\ldots,M_{1,t},M_{2,t})$ are independent GUE matrices such that $M_{1,i}$ and $M_{2,i}$ have dimensions $j_{i}\times j_{i}$.

Applying this lemma with $f_{i}(\lambda)=1/(\lambda-\rho_{i})$ and $g_{i}(\lambda)=1/(\lambda-\rho_{i})^{2}$, $R_{i}$ takes the value $R(\rho_{i})$ provided in the statement of Theorem 4. It results that

IV Application

In this section, we provide a general framework for local failure detection and diagnosis in large sensor networks, such as the examples proposed in Section II, based on the results of Section III. This framework is a two-step approach for successively (i) detecting failures within a given maximally acceptable false alarm rate and (ii) upon positive detection, diagnosing the failures with high probability. Simulations are then run to validate the proposed algorithms.

As stated in the introduction, the detection phase relies on existing results, and more specifically on the fluctuations of the largest eigenvalues given by Theorem 3. The detection algorithms proposed here parallel that introduced in in the context of collaborative signal sensing. The objective is to decide between hypothesis $\mathcal{H}_{0}$ and its complementary $\bar{\mathcal{H}}_{0}$.

First assume that all $P_{k}$ only have non-negative eigenvalues. From Theorem 1, the largest eigenvalue $\hat{\lambda}_{1}$ of $\Sigma\Sigma^{\ast}$ tends to the right edge of the support $S_{\pi}$ of the Marc̆enko-Pastur law for all large $n$ under $\mathcal{H}_{0}$, while $\hat{\lambda}_{1}$ is found away from this edge under $\bar{\mathcal{H}}_{0}$ if the largest eigenvalue $\omega_{k,1}$ of $P_{k}$ exceeds $\sqrt{c}$. We will therefore assume in the following that $\omega_{k,1}>\sqrt{c}$ is verified for all $k$. That is, we assume that $c_{N}\triangleq N/n<c_{+}$ where $c_{+}$ is defined as

This condition allows for a theoretically almost sure error detection, as $N,n\to\infty$. We then rely on Theorem 3 to design an appropriate hypothesis test. Our test consists in rejecting hypothesis $\mathcal{H}_{0}$ if the probability in favor of $\mathcal{H}_{0}$ is sufficiently low. That is, for a given acceptable false alarm rate $\eta$,We recall that the false alarm rate is the probability of declaring $\bar{\mathcal{H}}_{0}$ under true hypothesis $\mathcal{H}_{0}$. the statistical test is defined as

where $\hat{\lambda}^{\prime}_{1}$ is given by

That is, the test verifies whether $\hat{\lambda}_{1}$ exceeds some threshold above which the probability for $\mathcal{H}_{0}$ is less than $\eta$.

If the matrices $P_{k}$ are now all non-positive definite, then, symmetrically, we need to set $c_{N}$ such that the smallest eigenvalue $\hat{\lambda}_{N}$ of $\Sigma\Sigma^{\ast}$ is visible on the left-hand side of $S_{\pi}$. That is, we take $N,n$ to be such that $c_{N}<c_{-}$, with $c_{-}$ defined as

The decision test is in that case given by

where $\hat{\lambda}^{\prime}_{N}$ is defined as

The above test is particularly suited to the model of Section II-B in which the matrices $P_{k}$ are non-positive definite when $\mu_{k}=0$ and $\alpha_{k}=-1$ for all $k$, corresponding to a sudden drop of a zero mean random parameter $\theta(k)$ to zero.

The choice of $a(\eta),b(\eta)$ depends primarily on the structure of $P_{k}$ and will impact the correct detection rate for fixed false alarm rates.

In , the asymptotic independence of the fluctuations of the largest and the smallest eigenvalues of GUE matrices is proved, while the same result for the eigenvalues $\hat{\lambda}_{1}$ and $\hat{\lambda}_{N}$ of $\Sigma\Sigma^{\ast}$ under $\mathcal{H}_{0}$ is conjectured. Following this conjecture, (15) would become asymptotically

For any fixed $b$, taking $b(\eta)=b$, the hypothesis test now becomes

In particular, for $b=\infty$, $T_{2}(b)=1$ and then the test reduces to

which is the same test as proposed in (13). Taking instead $a(\eta)=-\infty$, we obtain the test (14).

For rather symmetrical distributions of the eigenvalues $\omega_{k,i}$ of $P_{k}$ around zero, it may be interesting to set $b(\eta)=a(\eta)$, in which case

In this setting, the decision test is now

In the following section, we assume that the procedure of failure detection was achieved successfully and that we are now interested in localizing the failures.

IV-B Localization algorithm

We now wish to detect all possible failure events from a set of failures indexed by $k\in\{1,\ldots,K\}$. The index set $\{1,\ldots,K\}$ may gather all events accounting for a single, as well as multiple, local failures. Similar to the previous sections, we denote $\rho_{k,i}=1+\omega_{k,i}+c\frac{1+\omega_{k,i}}{\omega_{k,i}}$ and $\zeta_{k,i}=\frac{1-c\omega_{k,i}^{-2}}{1+c\omega_{k,i}^{-1}}$, we define the mapping $\mathcal{K}_{k}$ to be such that $\mathcal{K}_{k}(i)=j_{k,1}+\ldots+j_{k,i-1}$ if $1\leq i\leq s_{k}$ and $\mathcal{K}_{k}(i)=N-(j_{k,i}+\ldots+j_{k,t_{k}})$ if $s_{k}+1\leq i\leq t_{k}$. Finally, we denote $\widehat{\Pi}_{k,i}$ any projector on the subspace generated by the eigenvalues $\hat{\lambda}_{\mathcal{K}_{k}(i)+1},\ldots,\hat{\lambda}_{\mathcal{K}_{k}(i)+j_{k,i}}$.

An initial hypothesis rejection may be performed at this stage to select only those hypotheses $\mathcal{H}_{k}$ such that the $\rho_{k,i}$ are consistent with the observations $\hat{\lambda}_{1},\ldots,\hat{\lambda}_{N}$. For instance, if the largest eigenvalue $\hat{\lambda}_{1}$ is significant in the system model, one may preselect from the set $\{1,\ldots,K\}$ the $L$ hypothesis indexes defined by

However, if $K$ is large, many hypotheses may have very close parameters $\rho_{k,i}$, so that it is hazardous to conclude on the most likely hypothesis based only on the eigenvalues of $\Sigma\Sigma^{\ast}$. However, since different $P_{k}$ matrices have in general very distinct eigenspaces, we propose the following subspace localization test, which decides on the hypothesis $\mathcal{H}_{k^{\star}}$ for which $k^{\star}$ is given by

with $g_{k}$ the actual density of the vector $\left(V^{k}_{i,n},L^{k}_{i,n}\right)_{i\in\mathcal{L}(p_{k},q_{k})}$, $\mathcal{L}(p_{k},q_{k})=\{1,\ldots,p_{k},q_{k},\ldots,r_{k}\}$, $S$ the set of (remaining) indexes $k$ such that $\mathcal{L}(p_{k},q_{k})$ is non-empty, and where

Note that we need here to specify the indexation $i\in\mathcal{L}(p_{k},q_{k})$ since we do not assume A5.

From Theorem 4, this probability can be approximated for large $n$, which provides immediately a maximum likelihood test for the most asymptotically likely $\mathcal{H}_{k}$ hypothesis. In the particular case where the $\omega_{k,i}$ all have multiplicity one, according to Corollary 4, as $N,n$ grow large, the vectors in the test (16) are asymptotically independent and Gaussian. We therefore substitute the test (16) by the following test, leading to the estimator $\hat{k}$ defined as

We provide below some remarks and discuss the advantages of the detection tests proposed in Section IV-A and the localization algorithms (17)–(IV-B) compared to the optimum maximum likelihood approach:

the detection algorithms proposed in Section IV-A are very versatile, as they adapt to multiple failure scenarios showing small rank perturbations in the population covariance, and provide a theoretical expression of the minimum ratio $c_{N}=N/n$ necessary for detectability;

unlike the traditional maximum-likelihood approach which tests the joint distribution of $\Sigma$ for all hypotheses $1\leq k\leq K$, and therefore leads to calculus of the order $\mathcal{O}(N^{3})$ (or $\mathcal{O}(N^{2})$ with some simplification methods) for each $k$, the proposed localization algorithm (17) is based on a test requiring for each $k$ (taken from a possibly reduced subset of $\{1,\ldots,K\}$) eigenvector projections of computational load of order $\mathcal{O}(N)$;

we may decide not to consider the joint fluctuations of all eigenvalues found outside $S_{\pi}$, but only some of them. This leads to an asymptotically less efficient, although much faster, algorithm, where $\mathcal{L}(p_{k},q_{k})$ in (17) is replaced by $\mathcal{L}(p^{\prime},q^{\prime})$ for given $p^{\prime}\leq p_{k}$, $q^{\prime}\leq q_{k}$, for all $k$. For $N$ not too large, it is in fact preferable to consider only a few eigenvalues and eigenspace projections simultaneously, due to convergence speed limitations of the limiting normal distributions;

the entries $L^{k}_{i,n}$ of the vector $(V^{k}_{i,n},L^{k}_{i,n})$ may also be discarded, especially in scenarios where eigenvalues of $P_{k}$ are very similar for each hypothesis $\mathcal{H}_{k}$. This may again increase the convergence speed of the asymptotic approximation for not-too-large $N$, while it is expected to perform worse for large $N$.

So far, we have performed failure detection under the important assumption that the failure scenarios form a discrete set $\{1,\ldots,K\}$. This assumes in particular that the failure amplitudes are known prior to detection and localization. In the next section, we use Corollary 4 to improve this approach in the particularly simple example of Section II-B, when the failure amplitude is a priori unknown.

IV-C Extension to unknown failure amplitude

In this section, we assume the scenario where the eigenvectors of the perturbation matrix $P_{k}$ are independent of the amplitude of the failure parameters, in the sense that a change in magnitude of the failure of type $k$ does not affect the eigenspaces of $P_{k}$. This is for instance the case of the single-failure scenario of Section II-B, for which we recall that $P_{k}$ expresses as $P_{k}=\beta_{k}R^{-{\frac{1}{2}}}He_{k}e_{k}^{\ast}H^{\ast}R^{-{\frac{1}{2}}}$ with $\beta_{k}$ the failure parameter. We now assume $\beta_{k}$ unknown, which is a more realistic assumption than assuming it perfectly known in advance. We also suppose that $X$ has i.i.d. Gaussian entries. Based on a simple extension of the algorithm presented in Section IV-B to unknown $\omega_{k}$, we provide hereafter a second localization algorithm.

For notational convenience, we assume $P_{k}=\omega_{k}u_{k}u_{k}^{\ast}$ for each $k$ and that $\omega_{k}>\sqrt{c}$, unknown. We then denote $\hat{\lambda}$ the largest eigenvalue of $\Sigma\Sigma^{\ast}$ and $\hat{u}$ its associated eigenvector.

Obviously, since $\omega_{k}$ is not known, neither is $\rho_{k}$. Therefore, we cannot proceed here to localization based on the fluctuations of $\hat{\lambda}$. Instead, we will use $\hat{\lambda}$ precisely as an estimate of $\rho_{k}$, which we know is consistent with growing $N,n$. From $\hat{\lambda}$, assumed larger than $(1+\sqrt{c})^{2}$, we want to derive an estimate $\hat{\omega}$ of $\omega_{k}$ ($k$ is the effective failure index). This is obtained from an inversion of the relation (7). Precisely, we obtain

From this estimate, we then obtain an estimate $\hat{\zeta}$ of $\zeta_{k}$ as follows

A natural object to consider for the failure localization is now $|u_{k}^{\ast}\hat{u}|^{2}-\hat{\zeta}$. To provide a diagnosis test, we need to derive the fluctuations of this random variable. From Theorem 4, the fluctuations of $\sqrt{N}(|u_{k}^{\ast}\hat{u}|^{2}-\zeta_{k})$ depend on $\omega_{k}$ but not on $u_{k}$. From the expression of $\hat{\zeta}$, it is immediate that the fluctuations of $\sqrt{N}(\hat{\zeta}-\zeta_{k})$ also depend on $\omega_{k}$ only. But since $\omega_{k}$ is estimated by $\hat{\omega}$, irrespective of the failure index $k$, the diagnosis test leads to finding the most likely argument $\hat{k}^{\prime}$ among $K$ variables with same Gaussian statistics. This therefore simplifies the estimator $\hat{k}^{\prime}$ of the most likely index $k$ to the following minimum-distance estimator

Note importantly that, contrary to our proposed scheme, the optimal maximum-likelihood localization method cannot be easily extended to the scenario of unknown failure amplitude, therefore bringing another significant advantage of the subspace approach.

In the next section, we provide simulation results for single failure localization for the detection and localization algorithms assuming the failure amplitude known or unknown, applied to the scenarios of Section II-A and Section II-B, respectively.

IV-D Simulations

In this section, we focus on the application of the algorithms designed in Sections IV-A and IV-B for single node failure in the scenario of Section II-A and single parameter change in the scenario of Section II-B.

Our first application example relates to the sensor network model $y=H\theta+\sigma w$ of Section II-A for $N=10$ nodes, $p=N$, and $\sigma^{2}=-20~{}{\rm dB}$. This is depicted in Figure 2, where the entries of $HH^{\ast}+\sigma^{2}I_{N}$ are presented. We also take $\sigma_{k}^{2}=\sum_{i=1}^{N}(HH^{\ast})_{ki}$, which is a natural assumption to avoid that a mere energy detector on $y(k)$ provide a simpler solution to our problem. This failure amplitude is assumed known by the experimenter. In practical scenarios, this may arise if a sensor starts returning time delayed data, supposedly uncorrelated with real-time data but with same variance. We assume a single failure scenario. In this context, it appears that, for all $k$, $\omega_{k,1}>0$, $\omega_{k,2}<0$ and $\omega_{k,1}$ is much larger than $|\omega_{k,2}|$. It is therefore more interesting only to consider the largest eigenvalue of $\Sigma\Sigma^{\ast}$ to detect and locate an hypothetical node failure. Under these conditions, the theoretical threshold for $c_{N}=N/n$ (if $N,n$ were large) is $0.8$ with the worst-case failure corresponding to a failure of node $10$. We therefore carry out $100\,000$ Monte Carlo simulations of node $10$ failures for $n$ varying from $8$ to $140$ and under false alarm rates varying from $10^{-2}$ to $10^{-4}$. This is depicted in Figure 3, where it can be observed that, for $n=8$, detection and localization are barely possible, although it is clearly the starting point where detection becomes feasible. For not too large $n$, while detection rates increase, we observe that localization capabilities are still unsatisfying. This is mainly due to the inappropriate fit of the large dimensional model with $N=10$ and with the eigenvectors corresponding to the extreme eigenvalues of $\Sigma\Sigma^{\ast}$ being too loosely correlated to their associated population eigenvectors. Larger values of $n$ show much better performance with miss localization probability going to zero as $n\to\infty$. In particular, about five times the asymptotically optimal ratio $n/N$ is required for localization to be very efficient. In this case, the large dimensional model for the fluctuations of the eigenvalues and eigenvectors is more adapted.

The same conditions are simulated for a system with $N=100$ nodes in which each node has eight neighbors and with correlation values of the same order of magnitude as in Figure 2. The detectability threshold for $N/n$ is here $0.85$ and we still consider the worst case failure scenario. This is depicted in Figure 4, where one can see that smaller ratios $n/N$ over the asymptotically optimal threshold are demanded for high detectability and localization ability to appear, when compared to the scenario $N=10$.

IV-D2 Sudden unknown parameter change

In this section, we consider the parameter change scenario of Section II-B. We still consider the network of Figure 2, and $\sigma^{2}=-20~{}{\rm dB}$ as above. We now assume a sudden change of parameter $\theta(10)$ with $\beta_{10}=2$, being the worst case scenario for failure identification if $\beta_{k}=2$ for all $k$. We depict the performance of the failure detection and localization algorithms and compare the settings where $\beta_{k}$ is known or unknown in advance to the experimenter. In the former scenario, we apply the localization algorithm of Section IV-B based on the joint fluctuations of the extreme eigenvalues and eigenspace projections, while in the latter, we apply the localization algorithm of Section IV-C, where a prior step of eigenvalue inference is performed before the study of the fluctuations of the eigenspace projections. The results are presented in Figure 5.

It appears from Figure 5 that the suboptimal algorithm of Section IV-C performs only slightly worse than the algorithm of Section IV-B for large $n$, and that it even performs better for small $n$. This last observation is explained by the inadequacy of the theoretical value of $\zeta_{k}$ for too small values of $n$. It is therefore interesting to see that, for practical purposes, the absence of prior knowledge on the amplitude of the failure does not severely reduce the efficiency of the localization algorithm.

V Conclusion

In this article, a characterization of the joint fluctuations of the extreme eigenvalues and corresponding eigenspace projections of a certain class of random matrices is provided. This characterization was used to perform fast and computationally reasonable detection and localization of multiple failures in large sensor networks through a general hypothesis testing framework. The main practical outcomes of this article lie first in a characterization of the minimum number of observations necessary to ensure failure detectability in large networks and second in the design of flexible but simple algorithms that can be adapted to multiple types of failure scenarios consistent with the small rank perturbation random matrix model. We also extend the detection and diagnosis approach to scenarios where the amplitudes of the hypothetical failures are not a priori known. Practical simulations suggest that the proposed algorithms allow for high failure detection and localization performance even for networks of small sizes, although for those much more observations than theoretically predicted are in general demanded.

Appendix A Proofs of results of Section III

We shall assume without loss of generality that $i=1$. In Section III-B2, we saw that

(take $b_{1}$ and $b_{2}$ as any two columns of $U_{1}$ in (8)). Similarly,

where $\mathcal{E}$ contains all the higher order terms that appear when we develop the integrand at the right hand side of the first equality.

In what follows, we successively study each of the terms at the right hand side of this equation. Recalling (9), the term $Z_{1}$ writes

The denominator has one simple zero in ${\rm{Int}}(\gamma_{1})$. With probability one, the numerator has no zero in ${\rm{Int}}(\gamma_{1})$. Using the residue theorem and the identity $(1+\omega_{1})^{-1}=(1+h(\rho_{1}))/h(\rho_{1})$, we obtain

which shows that we have a pole with degree $2$ in ${\rm{Int}}(\gamma_{1})$. Write the integrand as $G(z)/g(z)$ and recall that the residue of a meromorphic function $f(z)$ associated with a degree $2$ pole at $z_{0}$ is $\lim_{z\to z_{0}}d\left((z-z_{0})^{2}f(z)\right)/dz$. After some simple calculations, this results in

We now show briefly that the last term $\mathcal{E}$ in the expression of $V_{1,n}$ converges to zero in probability. A more detailed argument is given in . Recall that $\mathcal{E}$ accounts for all the higher order terms that show up when we expand the integrand $\widehat{A}_{1}\widehat{H}^{-1}\widehat{A}_{2}-A_{1}H^{-1}A_{2}$. Let us focus on one of these terms, namely

and show that it converges in probability to zero. The other terms can be treated similarly. First, we can show that

on $\gamma_{1}$ where $K$ is some constant. Now we write

Noticing that $A_{1}(z)$ and $A_{2}(z)$ are bounded on $\gamma_{1}$, and writing $z=\rho_{1}+R\exp(2\imath\pi\theta)$ on $\gamma_{1}$, the result is shown if we show that

for $i=1,2$. Lemma 1 shows that ${\rm E}\|E_{1}(z)\|^{2}\leq K^{\prime}/n$ on $\gamma_{1}$ where the constant $K^{\prime}$ is independent of $z$. By Markov’s inequality, (19) is true for $E_{1}$. Convergence for $E_{2}$ is obtained from Assumption A4 in conjunction with the analyticity of $\alpha(z)-m(z)$, as shown in .

Taking the sum $Z_{1}+Z_{2}+Z_{3}$, we obtain the desired result.

A-B Proof of Lemma 3

Let $y_{n}=\rho_{1}+x/\sqrt{N}$. Write $U=\begin{bmatrix}U_{1}&\widetilde{U}_{1}\end{bmatrix}$ and

where $K$ and $K^{\prime}$ are some constants, hence

Turning to the second term, we can show using A4 that

Again by A4, $Z_{3}\overset{\rm a.s.}{\longrightarrow}0$. This results in

The same argument for $i>1$ leads to the result.

A-C Proof of Lemma 4

Recall that $XX^{*}$ admits the spectral factorization $XX^{*}=W\Lambda W^{*}$ where $W$ and $\Lambda=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{N})$ are independent, and where $W$ is Haar distributed on the group of $N\times N$ unitary matrices.

From Assumption A4 and the analyticity of $f_{1}$, we can show as in that

hence the lemma is shown if we show the result on

where $[Z(Z^{*}Z)^{-\frac{1}{2}}]_{i}$ is the matrix formed by the columns $j_{1}+\ldots+j_{i-1}$ to $j_{1}+\ldots+j_{i}$ of $Z(Z^{*}Z)^{-\frac{1}{2}}$. By the law of large numbers, $N^{-1}Z^{*}Z\overset{\rm a.s.}{\longrightarrow}I_{r}$, hence it will be enough to show the result on

We observe that the summands of (20) are centered and are independent conditionally to $\Lambda$. Observe also that for every $k$, the vectors $(v_{i,k})_{i=1}^{t}$ are independent and that the elements of each of these vectors are decorrelated. Based on A2 and A3, we have

which coincides with the covariance matrix of (12) after the rearrangement.

Furthermore, thanks to A3, it is easy to see that the Lyapunov condition

is valid for any $\eta>0$, hence (20) satisfies the conditions of the central limit theorem, which proves the lemma.

References