Sequential Low-Rank Change Detection

Detecting emergence of a low-rank structure is a problem arising from many high-dimensional streaming data applications, such as video surveillance, financial time series, and sensor networks. The subspace structure may represent an anomaly or novelty that we would like to detect as soon as possible once it appears. One such example is swarm behavior monitoring. Biological swarms consist of many simple individuals following basic rules to form complex collective behaviors [1]. Examples include flocks of birds, schools of fish, and colonies of bacteria. The collective behavior and movement patterns of swarms have inspired much recent research into designing robotic swarms consisting of many agents that use simple algorithms to collectively accomplish complicated tasks, e.g., a swarm of UAVs [2]. Early detection of an emerging or transient behavior that leads to a specific form of behavior is very important for applications in swarm monitoring and control. One key observation, as shown in [1] for classification of behavior purposes, is that many forms of swarm behaviors are represented as low-dimensional linear subspaces. This leads to a low-rank change in the covariance structure of the observations.

In this paper, we propose a sequential change-point detection procedure based on the maximum eigenvalue of the sample covariance, which is a natural approach to detect the emergence of a low-rank signal. We characterize the performance of the detection procedure using two standard performance metrics, the average-run-length (ARL) that is related to the false-alarm-rate, and the expected detection delay. The distribution of the largest eigenvalue of the sample covariance matrix when data are Gaussian is presented in [3], which corresponds to the Tracey-Widom law of order one. However, it is intractable to directly analyze ARL using the Tracey-Widom law. Instead, we use a simple $\varepsilon$-set argument to decompose the detection statistic into a set of $\chi^{2}$-CUSUM procedures, whose performance is well understood in literature.

Moreover, we present a new approach to detect low-rank changes based on sketching of the observation vectors. Sketching corresponds to linear projection of the original observations: $y_{t}=Ax_{t}$, with $A\in\mathbb{R}^{m\times p}$. We may use sample covariance matrix of the sketches to detect emergence of a low-rank signal component, since it can be shown that for random projection, it can be shown that when $m$ is greater the rank of the signal covariance matrix, the sample covariance matrix of the linear sketches will have the same rank as the signal covariance matrix. We show that the minimum number of sketches $m$ is related to the property of the projection $A$ and the signal subspace. One key difference of the sketching for low-rank change detection from covariance sketching for low-rank covariance matrix recovery [4] is that, here we do not have to recover the covariance matrix which may require more number of sketches. Hence, sketching is a natural fit to detection problem.

The low-rank change is related to the so-called spiked covariance matrix [3], which assumes that a small number directions explain most of the variance. Such assumption is also made by sparse PCA, where the low-rank component is further assumed to be sparse. The goal of sparse PCA is to estimate such sparse subspaces (see, e.g., [5]). A fixed-sample hypothesis test to determine whether or not there exists a sparse and low-rank component, based on the largest eigenvalue statistic, is studied in [6], where it is shown to asymptotically minimax optimal. Another test statistic for such problems, the so-called Kac-Rice statistic, has been studied in [7]. The Kac-Rice statistic is the conditional survival function of the largest observed singular value conditioned on all other observed singular values, and it has a simple distribution form as uniformly distributed on [0, 1]. However, the statistic involves an infinite integral over the real line, which may not be easy to evaluate, and the test statistic needs to compute all eigenvalues of the sample covariance matrix instead of the largest eigenvalue.

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. Our goal is to detect such a change as quickly as possible. Formally, such a problem can be stated as the following hypothesis test:

Here $\sigma_{0}^{2}$ is the noise variance, which is assumed to be known or has been estimated from training data. Assume the signal covariance matrix $\Sigma$ is low-rank, meaning $\mbox{rank}(\Sigma)<p$. The signal covariance matrix is unknown.

One may construct a maximum likelihood ratio statistic. However, since the signal covariance matrix is 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 for detecting a low-rank signal. Assume a scanning window approach. We estimate the sample covariance matrix using samples in a time window of $[t-w,t]$, where $w$ is the window size. Assume $w$ is chosen sufficiently large and it is greater than the anticipated longest detection delay. Using the ratio of the largest eigenvalue of the sample covariance matrix relative to the noise variance, we form the maximum eigenvalue procedure which is a stopping time given by:

where $d$ is the pre-set drift parameter, $b>0$ is the threshold, and $\lambda_{1}(\Sigma)$ denotes the largest eigenvalue of a matrix $\Sigma$. Here the index $k$ represents the possible change-point location. Hence, samples between $[k+1,t]$ corresponds to post-change samples. The sample covariance matrix for post-change samples up to time $t$ is given by:

The maximization over $k$ corresponds to search for the unknown change-point location. An alarm is fired whenever the detection statistic exceeds the threshold $b$.

In this section, we characterize the performance of the procedure in terms of two standard performance metrics, the expected value of the stopping time when there is no change, called the average run length (ARL), and the expected detection (EDD), which is the expected time to stop in the extreme case when the change occurs immediately at $\kappa=0$.

Our approach is to decompose the maximum eigenvalue statistic to a set of $\chi^{2}$ CUSUM procedures. Our argument is based on the $\varepsilon$-net, which provides a convenient way to discretize unit sphere in our case. The number of such compact set is called the covering number, $\mathcal{C}(X,\varepsilon)$, which is the minimal cardinality of an $\varepsilon$-net. The covering number of a unit sphere is given by

The $\varepsilon$-net help to derive an upper bound of the largest eigenvalue, as stated in the following lemma

Our main theoretical results are the following, which are the lower bound on the ARL and the approximation to EDD of the largest eigenvalue procedure.

The lower bound above can be used to control false-alarm of the procedure. Given a target ARL, we may choose the corresponding threshold $b$.

Let $\|\Sigma\|$ denote the spectral norm of a matrix $\Sigma$, which corresponds to the largest eigenvalue. We have the following approximation:

In (1), $\rho^{2}/\sigma^{2}$ represents the signal-to-noise ratio. Note that the right-hand-side of (1) is a decreasing function of $\rho^{2}/\sigma^{2}$, which is expected, since the detection delay should be smaller when the signal-to-noise ratio is larger.

The following informal derivation justifies the theorems. First, from Lemma III.1, we have that the detection statistic

For each $q\in\mathcal{C}_{\varepsilon}(S)$, $\|q\|=1$, we have

Note that under $H_{0}$: $x_{i}\sim\mathcal{N}(0,\sigma_{0}^{2}I_{p})$, and hence $q^{\intercal}x_{i}/\sigma_{0}\sim\mathcal{N}(0,1)$, and

Alternatively, under $H_{1}$: $x_{i}\sim\mathcal{N}(0,\sigma_{0}^{2}I_{p}+\Sigma)$, and hence $q^{\intercal}x_{i}/\sigma_{0}\sim\mathcal{N}(0,1+q^{\intercal}\Sigma q/\sigma^{2})$, and hence

Hence, the distribution before the change is $\chi^{2}$ random variable, and after the change-point is a scaled $\chi^{2}$ random variable. Now define a set of procedures, for each $q\in\mathcal{C}_{\varepsilon}(S)$:

Note that each one of these procedures is a $\chi^{2}$ CUSUM procedure. Now if set

then due to the above relations, we may bound the average run length of the original procedure in terms of these $\chi^{2}$ procedures:

Since each $\widetilde{T}_{q}$ is a $\chi^{2}$ CUSUM procedure, whose properties are well understood, we may obtain a lower bound to the ARL of the maximum eigenvalue procedure.

When $p$ is large, a common practice is to use a linear projection $A$ to reduce data dimensionality. The linear projection maps the original $p$-dimensional vector into a lower dimension $m$-dimensional vector. We refer such linear projection as sketching. One implementation of sketching is illustrated in Fig. 3, where each signal vector $x_{t}$ is projected by $M$ vectors, $a_{1}$, …, $a_{M}$. The sketching corresponds to linear projection of the original vector:

Define a vector of observations $y_{t}=[y_{1,t},\ldots,y_{M,t}]^{\intercal}\in\mathbb{R}^{M}$, and

we have $y_{t}=A^{\intercal}x_{t}$, $t=1,2,\ldots$.

One intriguing question is whether we may perform detection of the low-rank change using the linear projections. The answer is yes, as we present in the following. We first show that each linear corresponds to a bi-linear projection of the original covariance matrix. Define the sample covariance matrix of the sketches

A key observation is that for certain choice of $A$, e.g., Gaussian random matrix, when $M>\mbox{rank}(\Sigma)$, the rank of $\widehat{\Sigma}_{y}$ is equal to the rank of $\widehat{\Sigma}_{x}$ with probability one. Hence, we may detect change using the largest eigenvalue of the sample covariance matrix of $y_{t}$.

A related line of work is covariance sketching, where the goal is to recover a low-rank covariance matrix using the quadratic sketches, which are square of each $(a_{i}^{\intercal}x_{t})^{2}$. This corresponds to only using diagonal entries of the sample covariance matrix (2) of the sketches.

The sketches $y_{t}$ are still multi-variate Gaussian distributed, and the projection changes the variance. For the sketches, we may consider the following hypothesis test based on the original test form:

Without loss of generality, to preserve noise property (or equivalently, to avoid amplifying noise), i.e., we choose the projection to be orthogonal, i.e., $A^{\intercal}A=I_{m}$. Moreover, due to the Gaussian form, the analysis for the ARL and EDD of the procedure in Section III still holds, with $\rho$ refined by $\|A^{\intercal}\Sigma A\|$. Hence, it can be seen from this analysis, that projection affects the signal-to-noise ratio, and hence, not surprisingly, although in principal we only need $M$ to be greater than the rank of $\Sigma$, in fact, we cannot choose $M$ to be arbitrarily small.

Fig. 4 illustrates this effect. In this example, we use a random subspace $A$ of dimension $M$ by $p$, and vary the number of sketches $M$, when the post-change covariance matrix of $x_{t}$ is generated as $\sigma_{0}^{2}I_{n}$ plus a Gaussian random matrix with rank 3. Then we plot the EDD of the largest eigenvalue procedure when sketches are used. We calibrate the thresholds in each setting so that the ARL is fixed to be 5000. Note that when the number of sketches is sufficiently large (and greater than $\mbox{rank}(\Sigma)$, the EDD approaches to a small number.

The following analysis may shed some light on how small we may choose $M$ to be. Let $s=\mbox{rank}(\Sigma)$, which is smaller than the ambient dimension $p$. Since the EDD (1) is a decreasing function of SNR, which is a proportion to $\|A^{\intercal}\Sigma A\|$. Denote the eigendecomposition of $\Sigma$ to be $U\Lambda U^{\intercal}$. Clearly, SNR for the sketching case will depend on $A^{\intercal}U$, which depends on the principal angle between two subspaces. Recall that we have required $A$ to be a subspace $A^{\intercal}A=I_{m}$, and hence the SNR will depend on the principal angle between the random projection and the signal subspace $U$. The principal angle between two random subspaces is studied in [10]. The eigenvalues of the sample covariance matrix are jointly Beta distribution. Based on this fact, the following lemma is obtained in [11], which helps us to understand the behavior of the procedure. It characterizes the $\ell_{2}$ norm of the fixed unit-norm vector projected by a random subspace.

Using this lemma in our setting, we may obtain a simple lower bound for the $\|A^{\intercal}U\|$:

with high probability, where the last inequality is due to the maximum of a set of Beta random variables. Hence, the number of sketches $M$ should scale as $p/\log s$.

There is a connection of the low-rank covariance and subspace model. Assume a rank $s$ post-change covariance matrix $\Sigma$, and its eigen-decomposion $U\Lambda U^{\intercal}$. Then we may expression each observation vector $x_{t}=U\beta_{t}+w_{t}$, where $w_{t}$ is a $p$-dimensional Gaussian random vector with covariance matrix $\sigma_{0}^{2}I_{p}$. Before the change, $\beta_{t}=0$. After the change, $\beta_{t}\neq 0$ and is a $s$-dimensional Gaussian random vector with covariance matrix being $\Lambda$. Hence, we may detect the low-rank change by detecting a non-zero $\beta_{t}$, if $U$ is known. When $U$ is unknown, we may perform an online subspace estimation from a sequence of data.

Based on such a connection, we may develop an efficient way to compute the detection statistic online via subspace tracking. This is related to the so-called matched-subspace detection [13], and here we further combine matched-subspace detection with online subspace estimation. Start with an initialization of the subspace $U_{0}$, using the sequence of observations $x_{1},x_{2},\ldots$, we may update the subspace using stochastic gradient descent on grassmannian manifolds, e.g., via the GROUSE algorithm [12]. Then we perform projection of the next vector $U_{t-1}^{\intercal}x_{t}$, which gives an estimate for $\beta_{t}$. We may claim a change when either $\max_{i}|[\beta_{t}]_{i}|$ (mimicking the largest eigenvalue) or the norm square $\|\beta_{t}\|^{2}$ (mimicking the sum of the eigenvalues) becomes large. Since the subspace tracking algorithm (e.g., GROUSE) can even deal with missing data, this approach allows us to compute the detection statistic even when we can only observe a subset of entries of $x_{t}$ at each time. Fig. 6 demonstrates the detection statistic computed via GROUSE subspace tracking when only about 70$\%$ of the entries can be observed at random locations. There is a change-point at time $500$. Clearly, the detection statistic computed via subspace tracking can detect such a change by raise a large peak.

We have presented a sequential change-point detection procedure based on the largest eigenvalue statistic for detecting a low-rank change to the signal covariance matrix. It is related to the so-called spiked covariance model. We present a lower-bound for the average-run-length (ARL) of the procedure when there is no change, which can be used to control false alarm rate. We also present an approximation to the expected detection delay (EDD) of the proposed procedure, which characterizes the dependence of EDD on the spectral norm of the post-change covariance matrix $\Sigma$. Our theoretical results are obtained using an $\varepsilon$-net argument, which leads to a decomposition of the proposed procedure into a set of $\chi^{2}$-CUSUM procedures. We further present a sketching procedure, which linearly projects the original observations into lower-dimensional sketches, and performs the low-rank change detection using these sketches. We demonstrate that when the number of sketches is sufficiently large (on the order of $p/\log s$ with $s$ being the rank of the post-change covariance), the sketching procedure can detect the change with little performance loss relative to the original procedure using full data. Finally, we present an approach to compute the detection statistic online via subspace tracking.

The author would like to thank Dake Zhang at Tsinghua University, Lauren Hue-Seversky and Matthew Berger at Air Force Research Lab, Information Directorate for stimulating discussions, and Shuang Li at Georgia Tech for help with numerical examples. This work is partially supported by a AFRI Visiting Faculty Fellowship.