Sequential detection of low-rank changes using extreme eigenvalues

TL;DR

This paper introduces sequential change detection methods for identifying abrupt shifts in covariance matrices from identity to low-rank structures using extreme eigenvalues, with theoretical analysis and real data validation.

Contribution

It proposes novel eigenvalue-based procedures for sequential change detection, incorporating accurate ARL and EDD approximations considering temporal correlations.

Findings

Methods effectively detect covariance changes in real data.

Theoretical approximations match empirical performance.

Approach is robust for swarm behavior analysis.

Abstract

We study the problem of detecting an abrupt change to the signal covariance matrix. In particular, the covariance changes from a "white" identity matrix to an unknown spiked or low-rank matrix. Two sequential change-point detection procedures are presented, based on the largest and the smallest eigenvalues of the sample covariance matrix. To control false-alarm-rate, we present an accurate theoretical approximation to the average-run-length (ARL) and expected detection delay (EDD) of the detection, leveraging the extreme eigenvalue distributions from random matrix theory and by capturing a non-negligible temporal correlation in the sequence of scan statistics due to the sliding window approach. Real data examples demonstrate the good performance of our method for detecting behavior change of a swarm.

Tables

Table 1. TABLE I: Comparison of upper bound for 𝔼 ​ [ Z t ​ Z t + δ ] 𝔼 delimited-[] subscript 𝑍 𝑡 subscript 𝑍 𝑡 𝛿 \mathbb{E}[Z_{t}Z_{t+\delta}] in Lemma 1 with simulation, for w = 200 𝑤 200 w=200 and p = 10 𝑝 10 p=10 .

	$\delta =2$	$\delta =6$	$\delta =10$	$\delta =15$	$\delta =20$
simulation	1.980	1.980	1.979	1.979	1.978
upper bound	2.144	2.241	2.294	2.342	2.379

Table 2. TABLE II: Comparison of the threshold b 𝑏 b obtained from simulation and the approximation. Window length w = 200 𝑤 200 w=200 , dimension p = 10 𝑝 10 p=10 . Simulation is the reference threshold.

Target ARL	5k	10k	20k	30k	40k	50k
Simulation	1.633	1.661	1.688	1.702	1.713	1.722
TW approx (8)	1.738	1.763	1.787	1.800	1.809	1.816
(12), $\beta$ in (11)	1.699	1.713	1.727	1.735	1.740	1.744

Equations

$$\left\{\begin{array}[]{ll}\textsf{H}_{0}:&x_{1},x_{2},\ldots,x_{t}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,I_{p})\\ \textsf{H}_{1}:&x_{1},x_{2},\ldots,x_{\tau}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,I_{p}),\\ &~{}~{}~{}x_{\tau+1},\ldots,x_{t}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,I_{p}+\theta uu^{\intercal}).\end{array}\right.$$

$${\hat{\Sigma }}_{t-w,t}=\frac{1}{w}\sum _{i=t-w+1}^t{x}_i{x}_i^{⊺}.$$

$$T_1=\inf \{t:{\lambda }_{\max }({\hat{\Sigma }}_{t-w,t})\ge b\},$$

$$\left\{\begin{array}[]{ll}\textsf{H}_{0}:&x_{1},x_{2},\ldots,x_{t}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,I_{p})\\ \textsf{H}_{1}:&x_{1},x_{2},\ldots,x_{\tau}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,I_{p}),\\ &~{}~{}~{}x_{\tau+1},\ldots,x_{t}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,\theta uu^{\intercal}).\end{array}\right.$$

$$T_2=\inf \{t:[{\lambda }_{\min }({\hat{\Sigma }}_{t-w,t}){]}^{-1}\ge b\},$$

$${\mu }_{np}$$

$${\sigma }_{np}$$

$$\frac{{\lambda }_{\max }({\hat{\Sigma }}_n)-{\mu }_{np}/n}{{\sigma }_{np}/n}\to W_1\sim F_1,$$

$$\frac{{\lambda }_{\min }({\hat{\Sigma }}_n)-{\mu }_{np}^{\prime }/n}{{\sigma }_{np}^{\prime }/n}\to W_1\sim F_1,$$

$${\mu }_{np}^{\prime }=(\sqrt{p}-\sqrt{n}{)}^2,$$

$${\sigma }_{np}^{\prime }=(\sqrt{p}-\sqrt{n})(\frac{1}{\sqrt{p}}-\frac{1}{\sqrt{n}}{)}^{1/3}.$$

$$b=\frac{{\sigma }_{wp}}{w}{\mathcal{T}}_{\alpha }+\frac{{\mu }_{wp}}{w},$$

$$b={\left[\frac{{\sigma }_{wp}^{\prime }}{w}{\mathcal{T}}_{1-\alpha }+\frac{{\mu }_{wp}^{\prime }}{w}\right]}^{-1},$$

$$Z_t={\lambda }_{\max }({\hat{\Sigma }}_{t-w,t}).$$

$$\mathrm{corr}(X,Y)=\frac{\mathbb{E}[XY]-\mathbb{E}(X)\mathbb{E}(Y)}{\sqrt{\mathrm{Var}(X)}\sqrt{\mathrm{Var}(Y)}}.$$

$$P=\sum _{i=t-w+1}^{t-w+\delta }{x}_i{x}_i^{\mathrm{T}},Q=\sum _{i=t-w+\delta +1}^t{x}_i{x}_i^{\mathrm{T}},W=\sum _{i=t+1}^{t+\delta }{x}_i{x}_i^{\mathrm{T}},$$

$$\mathbb{E}[Z_t{Z}_{t+\delta }]$$

$$+\mathbb{E}[{\lambda }_{\max }(Q)](\mathbb{E}[{\lambda }_{\max }(P)]+E[{\lambda }_{\max }(W)])$$

$$+\mathbb{E}[{\lambda }_{\max }(P)]\mathbb{E}[{\lambda }_{\max }(W)]\}$$

$$\mathrm{corr}(Z_t,Z_{t+\delta })\lesssim 1-\left(1+\frac{2p^{\frac{1}{3}}+3p^{\frac{1}{6}}\frac{c_1}{\sqrt{w}}+\frac{c_1^2}{w}}{c_2^2}\right)\delta +o(\delta )$$

$${\mathbb{E}}^{\infty }[T_1]\approx [\beta b\varphi (b)v(b\sqrt{2\beta }){]}^{-1},$$

$$v(x)\approx \frac{\frac{2}{x}\left[\Phi (\frac{x}{2})-\frac{1}{2}\right]}{\frac{x}{2}\Phi (\frac{x}{2})+\varphi (\frac{x}{2})},$$

$${\mathbb{E}}^1[T]\gtrsim \frac{b^{\prime }+e^{-b^{\prime }}-1}{\frac{\theta }{2}-\frac{1}{2}\log (1+\theta )}.$$

$$\mathbb{P}\left\{{\lambda }_{\max }({\hat{\Sigma }}_{t-w,t})\ge b\right\}$$

$$=$$

$$\approx$$

$$=$$

$$\frac{b-{\mu }_{wp}/w}{{\sigma }_{wp}/w}={\mathcal{T}}_{\alpha }\Rightarrow b=\frac{{\sigma }_{wp}}{w}{\mathcal{T}}_{\alpha }+\frac{{\mu }_{wp}}{w}$$

$$\mathbb{P}\left\{[{\lambda }_{\min }({\hat{\Sigma }}_{t-w,t}){]}^{-1}\ge b\right\}$$

$$=$$

$$\approx$$

$$=$$

$$\frac{1/b-{\mu }_{wp}^{\prime }/w}{{\sigma }_{wp}^{\prime }/w}={\mathcal{T}}_{1-\alpha }\Rightarrow b=\frac{1}{\frac{{\sigma }_{wp}^{\prime }}{w}{\mathcal{T}}_{1-\alpha }+\frac{{\mu }_{wp}^{\prime }}{w}}$$

$$Z_t$$

Full text

I Introduction

Detecting change-points from high-dimensional streaming data is a fundamental problem in various applications such as video surveillance, sensor networks, and biological swarm behavior study. In various scenarios, the change happens to the covariance. For instance, before the change, the signal covariance matrix is an identity matrix, and after the change, the signal covariance matrix becomes a spike low-rank matrix [1], or becomes a low-rank matrix due to the emergence of a subspace structure [2].

To detect such changes, it is natural to consider extreme (largest or smallest) eigenvalues of the sample covariance matrices over a time sliding window. Our prior work [2] presents an initial study of such method using the largest eigenvalue. In this paper, we present a more thorough theoretical analysis of these procedures with largest and smallest eigenvalues. In controlling the false alarm rate, i.e., deriving the average run length (ARL), we take a different approach than [2] by leveraging random matrix theory. In particular, we use the Tracy-Widom law for extreme eigenvalues to derive an ARL approximation. We further refined the precision our ARL approximation by capturing the temporal dependence of the detection statistic which is an inherent problem that needs to be addressed for change-point detection procedure due to its sliding window approach. Our technique is based on change-of-measure for Gaussian random field developed in [3]. Real data examples demonstrate the good performance of our algorithm.

The rest of this paper is organized as follows. Section II sets up the formalism for low-rank change detection. Section III presents the low-rank detection performance after ignoring temporal correlation of scan-statistics, and show an approximation to the temporal correlation of scan-statistic and give a more accurate theoretical analysis. Section IV presents the numerical example to demonstrate the performance of our method. Proofs can be found in the arXiv version this paper.

II Setup and detection procedures

Assuming a sequence of $p$-dimensional vectors $x_{1},x_{2},\ldots,x_{t}$, $t=1,2,\ldots$. There may be a change-point at time $\tau$ such that the distribution of the data stream changes. Assume the change happens at the covariance structure of the matrix. We will consider two related problems arising from specific applications: (i) the post-change covariance matrix model is a spiked covariance matrix; (ii) the covariance matrix before the change is full rank, and after the change becomes low-rank. Our goal is to detect such a change as quickly as possible.

II-A Largest eigenvalue procedure

Formally, problem (i) can be stated as the following hypothesis test:

[TABLE]

where $u\in\mathbb{R}^{p\times 1}$ represents a basis for a subspace and $u^{\intercal}u=1$.

One may construct a maximum likelihood ratio statistic. However, since the signal covariance matrix or $u$ are unknown, we may have to form the generalized likelihood ratio statistic, which replaces the covariance matrix with the sample covariance. This may cause an issue since the statistic involves inversion of the sample covariance matrix, whose numerical property (such as condition number) is usually poor when $p$ is large.

Alternatively, we consider the largest eigenvalue of the sample covariance matrix which is a natural detection statistic here. We adopt a scanning window approach to sequential methods. Consider samples in a time window of $[t-w,t]$ at each time $t$, where $w$ is the window size. Form the sample covariance matrix

[TABLE]

Using the largest eigenvalue of the sample covariance matrix, we form the maximum eigenvalue procedure, which is a stopping time given by:

[TABLE]

where $b>0$ is the threshold, and $\lambda_{\max}(\Sigma)$ denotes the largest eigenvalue of a matrix $\Sigma$. An alarm is fired whenever the detection statistic exceeds the threshold $b$.

II-B Smallest eigenvalue procedure

Problem (ii) can be formulated as

[TABLE]

where $u\in\mathbb{R}^{p\times 1}$ represents the basis for a one-dimensional subspace, $u^{\intercal}u=1$, and $\theta>0$ represents the signal power.

For this problem, we consider the smallest eigenvalue of the sample covariance matrix as the testing statistic, since the change is from full-rank to low-rank. We also adopt a scanning window approach assuming a window size $w$. Using the smallest eigenvalue of the sample covariance matrix, we form the minimum eigenvalue procedure, which is a stopping time given by:

[TABLE]

where $b$ is a pre-specified threshold.

II-C Choice of window-size $w$

There are two considerations when choosing $w$. First, $w$ has to be larger than the anticipated detection delay to not losing performance. Intuitively, we know that the signal strength $\theta$ need to be large enough in order to tell apart $\textsf{H}_{0}$ and $\textsf{H}_{1}$. The result in [4] shows that there is a critical value of signal strength $\theta$, i.e., when $\theta$ is above $\sqrt{p/n}\sigma$ the largest eigenvalue can be separated from the background noise. In our setting, this translate to that we need the window length needs to satisfy $\theta\geq\sigma\sqrt{p/w}$, which poses a requirement on the minimum window length $w\geq p/(\theta^{2}/\sigma^{2})$, where $\theta^{2}/\sigma^{2}$ can be interpreted as the signal-to-noise ratio. Second, $w$ should be chosen greater than the anticipated longest detection delay.

III Theoretical analysis

Two commonly used key performance metrics include the (i) the average run length (ARL), denoted as $\mathbb{E}^{\infty}[T_{i}]$, $i=1,2$, which is the expected duration in between two false alarms when there is no change, and (ii) the expected detection delay (EDD), denoted as $\mathbb{E}^{1}[T_{i}]$, $i=1,2$, which is expected number of samples before the procedure stops when the change-point happens at the first time $t=1$. In the following, we characterize these two metrics theoretically. Usually, the threshold $b$ is chosen such the ARL meets a certain large targeted value (e.g., 5,000 or 10,000).

III-A Background on extreme eigenvalue distributions

Since our detection procedures are based on extreme eigenvalues of the sample covariance matrix, it is essential to review their distributions from random matrix theory. There are two kinds of results typically available for eigenvalue distributions: one for the so-called bulk, which refers to the properties of the full set of eigenvalues, and one for the extremes, which are the (first few) largest and smallest eigenvalues.

For bulk spectrum, two well-known results are the the Wigner’s semicircle law (see, e.g., [5]), which describes the limiting density of eigenvalues of square symmetric random matrices, and the Marchenko-Pastur law for covariance matrices.

Based on the bulk spectrum distribution, various results have been obtained for extreme eigenvalue distributions. Assume there are $n$ samples are i.i.d. $p$-dimensional Gaussian random vectors with zero-mean and identity covariance matrix. Let the sample covariance matrix be $\widehat{\Sigma}_{n}=\frac{1}{n}\sum_{i=1}^{n}x_{i}x_{i}^{\intercal}$. [6] shows that if $p/n\rightarrow\gamma>0$, the largest eigenvalue of the sample covariance matrix converges to $(1+\sqrt{\gamma})^{2}$ almost surely. Although this is a very good approximation, it does not describe the variability of the largest eigenvalue. To characterize the distribution of the largest eigenvalues, [7] uses the Tracy-Widom law [8]. Define the center and scaling constants

[TABLE]

If $p/n\rightarrow\gamma<1$, then the centered and scaled largest eigenvalue converges in distribution to the so-called Tracy-Widom law of order one $W_{1}\sim F_{1}$ [7]:

[TABLE]

The Tracy-Widom law can be described in terms of partial differential equation and the Airy function, and its tail can be computed conveniently using an R-package RMTstat.

Similar result exists for the smallest eigenvalue of the sample covariance matrix [9]:

[TABLE]

where

[TABLE]

[TABLE]

III-B Approximation to Average-Run-Length (ARL)

Note that the scan-statistics in the Procedures (3) over time has temporal correlation due to overlapping data. For instance, $\widehat{\Sigma}_{t-w-1,t-1}$ and $\widehat{\Sigma}_{t-w,t}$ both involve samples $\{x_{t-w+1},\ldots,x_{t-1}\}$. To simplify the analysis, we first ignore the temporal correlation of scan-statistics and suppose different sample covariance matrixes $\widehat{\Sigma}_{t-w,t}$ are independent even if they may have some overlaps.

The Tracy-Widom law above gives an asymptotic distribution of the largest eigenvalue, which is very useful for us to analyze the distribution of our scan-statistics and choose the threshold of the detection procedure.

Proposition 1** (Approximation to ARL based on Tracy-Widom, ignore temporal correlation).**

For any $\alpha\in(0,1)$, when choosing for $T_{1}$

[TABLE]

$\mathbb{E}^{\infty}[T_{1}]\approx 1/\alpha$. Above, $\mathcal{T}_{\alpha}$ is the $\alpha$-upper-percentage point for Tracy-Widom law of order one. When choosing for $T_{2}$,

[TABLE]

$\mathbb{E}^{\infty}[T_{2}]\approx 1/\alpha$.

Remark: Note that the threshold choice here is different from that in [1], since here we used the tail probability of the precise distribution of the extreme eigenvalues based on Tracy-Widom law, whereas the later uses an approximate threshold to characterization the rate of change for extreme eigenvalues, which may be less accurate.

III-C Correlation of temporal scan-statistics

Now we aim to capture the temporal dependence in the scan statistics, which may be significant due to overlapping of adjacent time windows. To start, we consider $T_{1}$ which uses the largest eigenvalue. This is a very challenging task, as it is an open question what is the correlation between the largest eigenvalues of two sample covariance matrices share partially common data.

For sample covariance matrix defined over a time window $[t-w,t]$ in (2), define

[TABLE]

Define the correlation between two random variables as usual

[TABLE]

We have the following lemma, which leads to some insights for our analysis.

Lemma 1** (Approximation to local correlation).**

Under the null hypothesis, $x_{t}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,I_{p})$, For $\delta\in\mathbb{Z}^{+}$, let

[TABLE]

then $P$, $Q$ and $W$ are mutually independent random matrices. We have that

[TABLE]

where the mean and second-order moments can be computed using Tracy-Widom law shown in (7). Furthermore, an approximation to the local correlation for $w\gg p$ and small $\delta\ll w$, is given by

[TABLE]

where $c_{1}=\mathbb{E}(W_{1})=-1.21$ and $c_{2}=\sqrt{{\rm Var}(W_{1})}=1.27$.

Lemma 1 offers reasonably good approximation the the $\mathbb{E}[Z_{t}Z_{t+\delta}]$, as shown in Table I. We also run simulations to verify correlation between the largest eigenvalues between two overlapping data blocks ${\rm corr}(Z_{t},Z_{t+\delta})$, shown in Fig. 2. Note that the correlation increases as the “overlapping ratio” $1-\delta/w$ increases. Indeed, the correlation function can be linearly approximated in a local region when $\delta/w\ll 1$.

Proposition 2** (Approximation to ARL, considering temporal correlation).**

When $b\rightarrow\infty$, if ${\rm corr}(Z_{t},Z_{t+\delta})=1-\beta\delta+o(\delta)$, then the average-run-length (ARL) of $T_{1}$ in (3) can be approximated as

[TABLE]

Above, $v(\cdot)$ is a special function closely related to the Laplace transform of the overshoot over the boundary of a random walk ([10]):

[TABLE]

and $\phi(x)$ and $\Phi(x)$ are the probability density function and the cumulative density function of the standard normal distribution.

The proof uses the change-of-measure techniques in [3].

Now we verify the accuracy of the threshold obtained without and with considering the temporal correlation (Proposition 1 and Proposition 2, respectively). We find that indeed the threshold with temporal correlation correlation (12) with $\beta$ in (11) is more accurate than the threshold obtained from Tracy-Widom law ignoring temporal correlation (which seems to perform reasonable well).

III-D Approximation to EDD

A lower bound to EDD is computed using the same strategy as in [2] but here the approximation is more accurate.

Proposition 3** (Lower bound to EDD).**

When $b\rightarrow\infty$

[TABLE]

where $b^{\prime}=\frac{1}{2}[1-\frac{1}{1+\theta}][b-\frac{\log(1+\theta)}{1-1/(1+\theta)}]w$.

The proof is based on an result for the CUSUM procedure in [11]. Consistent with intuition, Proposition 3, the right-hand-side of (13) is a decreasing function of $\theta$, representing the signal-to-noise ratio. Moreover, we compare the lower bound in Proposition 3 with simulated average delay, as shown in Fig. 3. In the regime of small detection delay (which is the main regime of interest), the lower bound serves as a reasonably good approximation. There is still room for improvement, which we are working on leveraging Grothendieck’s Inequality [12].

IV Real-data examples

In this section we present two examples for detecting swarm behavior changes [13]. We are interested in detecting the transition of a swarm (represented by dots in Fig. 4 by dots) transitioning from random behavior to some organized behavior: such as “flock” or “torus”. The swarm dataset contains the coordinates of each individual at each unit time between a specific time interval. The swarm dataset contains a sequence of the coordinates and velocity of each individual at each time. We verify that in these cases the change can be well represented by the setting in (4), and applying $T_{2}$ in (5) can quickly detect the transition from random behavior to “flock” or “torus”.