Quickest Hub Discovery in Correlation Graphs

TL;DR

This paper introduces a sequential testing method for quickly detecting and isolating hubs in high-dimensional correlation graphs, using novel summary statistics to improve accuracy and speed.

Contribution

It proposes a new sequential test leveraging local and global statistics for efficient hub detection in high-dimensional Gaussian correlation graphs.

Findings

The test is consistent in identifying hubs as false alarm rate approaches zero.

Delay and false alarm performance are analytically characterized.

Numerical results demonstrate the method's effectiveness.

Abstract

A sequential test is proposed for detection and isolation of hubs in a correlation graph. Hubs in a correlation graph of a random vector are variables (nodes) that have a strong correlation edge. It is assumed that the random vectors are high-dimensional and are multivariate Gaussian distributed. The test employs a family of novel local and global summary statistics generated from small samples of the random vectors. Delay and false alarm analysis of the test is obtained and numerical results are provided to show that the test is consistent in identifying hubs, as the false alarm rate goes to zero.

Equations

$$f_{\mathbf{X}(k)}(\mathbf{x})=\frac{1}{{\sqrt{2\pi }}^p|{\Sigma }_k{|}^{p/2}}{e}^{\left(-\frac{1}{2}(\mathbf{x}-{\mathbf{\mu }}_k{)}^T{\mathbf{\Sigma }}_k^{-1}(\mathbf{x}-{\mathbf{\mu }}_k)\right)}.$$

$$\begin{array}{rl}\mathbf{X}(k)& \sim f_{\mathbf{X}}^0(\mathbf{x})=\frac{1}{{\sqrt{2\pi }}^p|{\Sigma }_0{|}^{p/2}}{e}^{-\frac{1}{2}(\mathbf{x}-{\mathbf{\mu }}_k{)}^T{\mathbf{\Sigma }}_0^{-1}(\mathbf{x}-{\mathbf{\mu }}_k)},\\ & k<\gamma ,\\ & \sim f_{\mathbf{X}}^1(\mathbf{x})=\frac{1}{{\sqrt{2\pi }}^p|\Sigma {|}^{p/2}}{e}^{-\frac{1}{2}(\mathbf{x}-{\mathbf{\mu }}_k{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{x}-{\mathbf{\mu }}_k)},\\ & k\ge \gamma \end{array}$$

$$V_k(\mathbf{\Sigma })=\underset{i\ne k}{\max }|{\rho }_{ki}|,\text{mbox}fork\in \{1,\cdots ,p\},$$

$$\begin{array}{rl}\mathsf{H}& =\{k:V_k=\underset{1\le j\le p}{\max }{V}_j\}.\end{array}$$

$$\mathcal{D}\left(\mathbf{X}(1),\cdots ,\mathbf{X}(\tau )\right)\in 2^{\{1,\cdots ,p\}}$$

$$\begin{array}{rl}\underset{\{\tau ,\mathcal{D}\}}{\min }& \underset{\gamma \ge 1}{\sup }{\mathsf{E}}_{\gamma }[\tau -\gamma |\tau \ge \gamma ]\\ \text{mbox}subj.to& {\mathsf{E}}_{\infty }[\tau ]\ge \beta ,\\ \text{mbox}and& \underset{\gamma }{\sup }{\mathsf{P}}_{\gamma }[\mathcal{D}\ne \mathsf{H}]\le \zeta ,\end{array}$$

$$\begin{array}{rl}\mathbb{X}(1)& =[\mathbf{X}(1),\cdots ,\mathbf{X}(n){]}^T\\ \mathbb{X}(2)& =[\mathbf{X}(n+1),\cdots ,\mathbf{X}(2n){]}^T,\text{mbox}etc.\end{array}$$

$$\mathbb{X}=[{\mathbf{X}}_1,\cdots ,{\mathbf{X}}_p]=[{\mathbf{X}}_{(1)}^T,\cdots ,{\mathbf{X}}_{(n)}^T{]}^T,$$

$$\mathbf{S}=\frac{1}{n-1}\sum _{i=1}^n({\mathbf{X}}_{(i)}-\bar{\mathbf{X}}{)}^T({\mathbf{X}}_{(i)}-\bar{\mathbf{X}}),$$

$$\mathbf{R}={{\mathbf{D}}_{\mathbf{S}}}^{-1/2}\mathbf{S}{{\mathbf{D}}_{\mathbf{S}}}^{-1/2},$$

$$V_k(\mathbb{X}):=V_k(\mathbf{R}):=\underset{i\ne k}{\max }|{\mathbf{R}}_{ki}|,k\in \{1,\cdots ,p\}.$$

$$\begin{array}{rl}{d}_k(\mathbf{\Sigma },\rho )& :=\#\{i\ne k:|{\rho }_{ki}|\ge \rho \},\\ & k\in \{1,\cdots ,p\},\rho \in [0,1],\end{array}$$

$$\mathsf{P}(V_k(\mathbb{X})\le \rho )\to e^{-{\Lambda }_k},$$

$${\Lambda }_k=\underset{p\to \infty ,\rho \to 1}{\lim }{\Lambda }_{k,\rho },$$

$${\Lambda }_{k,\rho }=(p-1)P_0(\rho )J_k,$$

$$P_0(\rho )=I_{1-{\rho }^2}((n-2)/2,1/2),$$

$$J_k=1+O((j/p{)}^2).$$

$$J_k=1.$$

$$\mathsf{P}(V_k(\mathbb{X})\ge \rho )=\mathsf{P}(d_k(\mathbb{X},\rho )>0).$$

$$\begin{array}{rl}\mathsf{P}(V_k(\mathbb{X})\le \rho )& \approx \exp (-{\Lambda }_{k,\rho }),\\ & =\exp (-(p-1)P_0(\rho )J_k)\rho \in [0,1].\end{array}$$

$$f_{V_k}(\rho )\approx -(p-1)P_0^{\prime }(\rho )J_k\exp (-(p-1)P_0(\rho )J_k),\rho \in (0,1].$$

$$\begin{array}{rl}{f}_{V_k}& (\rho ;J_k)\\ =& \frac{2(p-1)J_k(1-{\rho }^2{)}^{\frac{n-4}{2}}}{B((n-2)/2,1/2)}\exp (-(p-1)J_k{P}_0(\rho )).\end{array}$$

$$\begin{array}{rl}{f}_V& (y;J)=JC(p,n)(1-y^2{)}^{\frac{n-4}{2}}{e}^{-(p-1)J{P}_0(y)},\end{array}$$

$$\begin{array}{rl}C(p,n)& =\frac{2(p-1)}{B((n-2)/2,1/2)},\\ P_0(y)& =I_{1-y^2}((n-2)/2,1/2),\\ I_x(a,b)& =\frac{B(x;a,b)}{B(a,b)},\\ B(x,a,b)& ={\int }_0^x{t}^{a-1}(1-t{)}^{b-1}dt,\\ B(a,b)& =B(1,a,b).\end{array}$$

$$U(\mathbb{X})=\underset{k}{\max }{V}_k(\mathbb{X}),$$

$$g(u;\theta )=\frac{D\theta }{2}(1-u^2{)}^{\frac{n-4}{2}}\exp \left(-\frac{D}{2}\theta T(u)\right),u\in (0,1],$$

$$\begin{array}{rl}D=D_{p,n}& =2p(p-1)B((n-2)/2,1/2),\\ T(u)& ={\int }_u^1(1-s^2{)}^{\frac{n-4}{2}}ds.\end{array}$$

$$\begin{array}{rl}U(m)& \sim g(\cdot ;1),m<\gamma \\ & \sim g(\cdot ;\theta ),\theta \ne 1,m\ge \gamma .\end{array}$$

$$\begin{array}{rl}& {\tau }_{\text{U}}\\ & =\min \left\{m:\underset{1\le \ell \le m}{\max }\underset{\theta :|\theta -1|\ge {ϵ}_u}{\sup }\sum _{i=\ell }^m\log \frac{g(U(i);\theta )}{g(U(i);1)}>A_u\right\},\end{array}$$

$$\begin{array}{rl}{V}_k(m)& \sim f_V(\cdot ;1),m<\gamma \\ & \sim f_V(\cdot ;J_k),J_k\ne 1,m\ge \gamma .\end{array}$$

Full text

I Introduction

In this paper we consider the problem of quickest detection of a change in correlation between variables in a sequence of high-dimensional vectors. We are also interested in isolating or identifying the variables whose magnitude correlation with other variables has changed the most. We refer to such variables as hubs. A precise definition will be given below.

The problem of correlation detection and estimation is an important one in science and engineering. The problem is particularly difficult when the data involved is high-dimensional; see [1], [2], and [3] for a survey of the literature. We are interested in the sequential version of the problem where the variables are initially uncorrelated or independent, and due to some event, the correlation between the variables abruptly changes. This has applications in problems of fault detection, anomaly detection, detection of changes in time-series data or financial data, etc. For simplicity, we restrict our attention to multivariate Gaussian random vectors. However, the results in this paper are application to the broader class of elliptically contoured distributions [4].

If the pre-change parameter are known (means and variances of the uncorrelated variables), the classical approach to sequential detection of a change in the covariance matrix of a sequence of Gaussian random vectors is to use the generalized likelihood ratio (GLR) based Cumulative Sum (CUSUM) method of Lorden [5]. In this algorithm one sequentially estimates the covariance of the data and applies classical quickest change detection (QCD) test using the estimated covariance. But, this algorithm works only for low dimensional data. One approach to handle high-dimensional data, is to make assumptions about the structure of the covariance of the matrix, like sparsity, and then estimate the high-dimensional covariance matrix. Although this approach is appealing, it is generally hard to design sequential tests based on such an approach. When the pre-change parameters are not known, as is the case in this paper, the problem is even harder. For low dimensional data, one can use a double-GLR approach [6], but optimality of such tests is still under investigation.

In this paper we take a random matrix based approach to solving this problem. We collect the sequence of high-dimensional vectors in batches of successive samples from which we form a random matrix. We then define summary statistics for the random matrix and obtain approximations for the asymptotic distribution for the summary statistics as the dimension of the vector goes to infinity. The asymptotic distribution we obtain belongs to a one-parameter exponential family, and a change in correlation of the random vectors induces a change in the parameter of the random matrix. We use this fact to detect the change in correlation. This works differs from our effort in [4] in that we also define a family of local summary statistics, one for each variable, and obtain their approximate asymptotic distributions. We use these local summary statistics to isolate or identify variables that have seen the most amount of change in correlation with other variables. This fault isolation process is what we call hub discovery. We obtain asymptotic expressions for delay and false alarm for our proposed procedure, and show via numerical results that the test is consistent in hub discovery.

II Problem Description

We have a sequence of samples of independent and identically distributed (i.i.d.) Gaussian $p$-variate random vectors $\{\mathbf{X}(k)\}$ with densities

[TABLE]

For some time parameter $\gamma$ the samples are assumed to have common covariance parameter $\mathbf{\Sigma}_{0}=D_{p}$ for $k<\gamma$, where $D_{p}$ is a digonal matrix, and common covariance parameter $\mathbf{\Sigma}\neq D_{p}$ and for $k\geq\gamma$. $\gamma$ is called the change point and the pre-change and post-change distributions of $\mathbf{X}(k)$ are denoted by $f_{\mathbf{X}}^{0}$ and $f_{\mathbf{X}}^{1}$, respectively. No assumptions are made about the mean parameter $\boldsymbol{\mu}_{k}$, and can take different values for different $k$. More specifically, the change-point model is described by:

[TABLE]

Thus, in the above change point model, the variables are uncorrelated before change, and correlated after change. The objective is to detect this change in correlation as quickly as possible, while avoiding false alarms, and also correctly isolating the variables that have experienced the most change in their correlation level (to be made precise below). The mean parameters $\{\boldsymbol{\mu}_{k}\}$, the diagonal entries in $\boldsymbol{\Sigma}_{0}$ and the post-change covariance matrix $\boldsymbol{\Sigma}$, are unknown.

A decision-maker sequentially acquires the samples $\{\mathbf{X}(k)\}$. At each time point $k$ the decision-maker decides to either stop sampling, declaring that the change has occurred, i.e., $k\geq\gamma$, or to continue sampling. The decision to stop at time $k$ is only a function of $(\mathbf{X}(1),\cdots,\mathbf{X}(k))$. Thus, the time at which the decision-maker decides to stop sampling, say $\tau$, is a stopping time for the vector sequence $\{\mathbf{X}(k)\}$. At the stopping time $\tau$, the decision maker incurs a delay of $\tau-\gamma$, if $\tau\geq\gamma$. If $\tau<\gamma$, we say that we have an event of false alarm. At time $\tau$, using the available data $\left(\mathbf{X}(1),\cdots,\mathbf{X}(\tau)\right)$, the decision maker also has to identify or isolate variables that have undergone the most amount of change in correlation. We refer to such variables as hubs. We now formalize the notion of a hub.

For a non-diagonal covariance matrix $\boldsymbol{\Sigma}$ with correlation coefficients $\{\rho_{ki}\}$ define

[TABLE]

to be the maximum magnitude correlation coefficient for the $k$th variable. Hubs are defined as

[TABLE]

Define the correlation graph $\mathcal{G}(\boldsymbol{\Sigma})$ associated with the post-change correlation matrix $\boldsymbol{\Sigma}$ as an undirected graph with $p$ vertices, each representing a variable in the vector $\mathbf{X}$. An edge is present between vertices $i$ and $j$ if the magnitude of the correlation coefficient $\rho_{ij}$ between the $i^{th}$ and $j^{th}$ components of the random vector $\mathbf{X}$ is nonzero, i.e., if $|\rho_{ij}|>0$, $i\neq j$. The correlation graph $\mathcal{G}(\boldsymbol{\Sigma})$ can be treated as a weighted graph with weight on edge between node $i$ and $j$ equal to $|\rho_{ij}|$, provided the latter is nonzero. Then, hubs are the nodes with highest maximum weights.

There are also other interesting ways to define hubs. For example, one can select a threshold $\rho$ and define hubs as those variables for which $V_{k}>\rho$. Another way to define hubs is to redefine correlation graphs with edges present only of $|\rho_{ij}|\geq\rho$, and then define hubs as the nodes with the highest degree. Results on hub discovery for these other type of hubs will be reported elsewhere. In this paper we only discuss hub discovery for hubs defined in (4).

At time $\tau$, let

[TABLE]

be the decision function that selects a subset of the $p$ variables as the hubs. The decision-maker’s objective is to choose a test, i.e. $\{\tau,\mathcal{D}\}$, to detect this change in correlation of the random vectors, as quickly as possible, subject to a constraint on the rate of false alarms, and on the rate of false isolation.

A problem formulation that captures the above performance trade-offs is the following. We seek $\{\tau,\mathcal{D}\}$ to solve

[TABLE]

where $\mathsf{E}_{\gamma}$ is the expectation with respect to the probability measure under which the change occurs at $\gamma$, $\mathsf{E}_{\infty}$ is the expectation when the change never occurs, and $\beta$ and $\zeta$ are user-specified constraints on the mean time to false alarm and the probability of false isolation, respectively.

This problem formulation is also studied in [7] in the context of detecting and isolating a single affected stream out of a finite number of independent streams of observations, when the exact pre- and post-change distributions are known. Even with independent streams, and known distributions, the author in [7] obtained a solution for only a relaxed version of the above problem, by replacing the false isolation probability metric by $\mathsf{P}_{1}[\mathcal{D}\neq\mathsf{H}]$.

In our case there are three major challenges towards developing a test or extending/applying the results from [7]:

The means $\{\mu_{k}\}$, covariance matrices $\boldsymbol{\Sigma}_{0}$ and $\boldsymbol{\Sigma}$ are unknown 2. 2.

The high dimension $p$ that is much larger than the total number of samples we will possibly accumulate. 3. 3.

The streams (one for each dimension) in our observation models are not independent: In fact, we are trying to detect a change in the level of dependence between the streams.

In the sequential analysis literature, for such intractable problems, an alternative approach often taken is to propose a sub-optimal test, and provide a detailed performance analysis of the proposed test. The performance analysis is often used to design the test (to choose thresholds). The effectiveness of the proposed test is then verified by simulations. See for example [8], [9], and [10].

In [4], we proposed a test for the non-parametric family of vector elliptically contoured densities, of which the Gaussian case considered here is a special case, with the objective of minimizing detection delay while avoiding false alarms. In our test in [4], we used a global summary statistic to detect the change in correlation. However, the global nature of our summary statistics in [4] does not allow for classification of where the changes occurred in the correlation matrix, i.e., correlation hub discovery. In this paper we propose a family of $p$ local summary statistics, one for each variable, and obtain asymptotic densities for these local statistics. We use these and the global statistic from [4] to perform sequential classification of the hubs in addition to detection of the existence of hubs.

III A Random Matrix Solution Approach

We propose a random matrix based solution to our detection problem. The decision maker analyzes the $p$-dimensional vector sequence $\{\mathbf{X}(k)\}$ in batches of size $n$, where $p\gg n$. This leads to a $n\times p$ random matrix sequence $\{\mathbb{X}(m)\}_{m\geq 1}$. Specifically,

[TABLE]

Thus, for each $m$, each of the $n$ rows of the random matrix $\mathbb{X}(m)$ is an i.i.d. samples of the multivariate Gaussian $p$-variate random vector $\mathbf{X}$ $=[X_{1},\cdots,X_{p}]^{T}$ (1).

A change in the covariance of the random vector sequence $\{\mathbf{X}(k)\}$ as per the description in (2) will change the covariance of the random matrix sequence $\{\mathbb{X}(m)\}$. We propose to detect the change by using a stopping rule for the sequence $\{\mathbb{X}(m)\}$. Note that:

A stopping rule for the sequence $\{\mathbb{X}(m)\}$ is also a stopping rule for the sequence $\{\mathbf{X}(k)\}$, but it takes values in multiples of $n$ (with respect to index $k$). The precise choice of $n$ will be discussed below. 2. 2.

The law of the random matrix $\mathbb{X}$ at the change point may consist of random vectors $\mathbf{X}$ from both pre- and post-change distributions. We ignore this issue in the rest of the paper and assume that all the vectors in the matrix either have density $f_{\mathbf{X}}^{0}$ or have density $f_{\mathbf{X}}^{1}$ (2). However, see [11] and [12], where such transient behavior is addressed.

In random matrix theory (RMT), asymptotic distributions are obtained for various functions of the random matrix, e.g., a semicircle law for the empirical eigen value distribution of a Wigner matrix, etc [13]. Such results are in general valid when both the dimensions of the matrix are taken to infinity. For our random matrix $\mathbb{X}$, however, the batch size $n$ is fixed, and the dimension $p$ is large. This is called the purely high-dimensional regime [14]. Thus, traditional RMT theorems cannot be applied in this regime. In this paper, we obtain approximate asymptotic distribution for a proposed summary statistic (in fact, a family of them) for the random matrix $\mathbb{X}$ in this regime. We then use this approximate distribution(s) to detect a change in the law of the random matrix. We also propose techniques to isolate the hubs.

IV Local Summary Statistics for the Random Data Matrix

In this section we define a family of local summary statistics to detect hubs (4). The summary statistics are a function of the sample correlation matrix obtained from the random data matrix $\mathbb{X}$. We then obtain their asymptotic distribution in the purely high dimensional regime of $p\rightarrow\infty$, $n$ fixed. The asymptotic distributions are members of a one-parameter exponential family to be specified below. In Section VI, we will use these summary statistics and their asymptotic distributions for sequential hub discovery.

For random data matrix $\mathbb{X}$ we write

[TABLE]

where $\mathbf{X}_{i}=[X_{1i},\cdots,X_{ni}]^{T}$ is the $i^{th}$ column and $\mathbf{X}_{(i)}=[X_{i1},\cdots,X_{ip}]$ is the $i^{th}$ row. Define the sample covariance matrix as

[TABLE]

where $\bar{\mathbf{X}}$ is the sample mean of the $n$ rows of $\mathbb{X}$. Also define the sample correlation matrix as

[TABLE]

where $\mathbf{D_{A}}$ denotes the matrix obtained by zeroing out all but the diagonal elements of the matrix $\mathbf{A}$. Note that, under our assumption that the covariance matrix $\mathbf{\Sigma}$ of the rows of $\mathbb{X}$ is positive definite, $\mathbf{D_{S}}$ is invertible with probability one. Thus $\mathbf{R}_{ij}$, the element in the $i^{th}$ row and the $j^{th}$ column of the matrix $\mathbf{R}$, is the sample correlation coefficient between the $i^{th}$ and $j^{th}$ columns of $\mathbb{X}$.

Define the summary statistics (see also (3))

[TABLE]

This defines a set of $p$ local statistics $\{V_{k}\}_{k\leq p}$, one for each random variable in the random vector. We will use these statistics for detection of correlation strength hubs.

Note that the summary statistics are the empirical estimates of the variables $\{V_{k}\}$ defined for the population correlation matrix $\boldsymbol{\Sigma}$ in (3). If we could collect enough samples such that $n>p$, then we can rely on these estimates to check if they are close to zero, because the population version of these quantities are identically zero when the variables are uncorrelated. But, since $p\gg n$, we may have to collect a large number of samples before we can rely on these estimates. Such a delay may not be acceptable in practice, or those many samples may not even be available. In fact, for large $p$, $V_{k}(\mathbb{X})$ can be arbitrarily close to $1$; see [15] and Fig. 2 in this paper.

However, as we show below, one can obtain asymptotic law of these estimates as $p\to\infty$. The precise statement is provided in the theorems below. A change in correlation between the variables may be reflected in the asymptotic law of these estimates. This phenomena can be utilized for detecting a change in correlation. The rest of this section is devoted to obtaining the asymptotic laws. The algorithms are discussed in Section VI.

We say that a matrix is row sparse of degree $j$ if there are no more than $j$ nonzero entries in any row. We say that a matrix is block sparse of degree $j$ if the matrix can be reduced to block diagonal form having a single $j\times j$ block, via row-column permutations. For a threshold parameter $\rho\in[0,1]$ define the correlation graph $\mathcal{G}_{\rho}(\boldsymbol{\Sigma})$ associated with the post-change correlation matrix $\boldsymbol{\Sigma}$ as an undirected graph with $p$ vertices, each representing a variable in the vector $\mathbf{X}$. An edge is present between vertices $i$ and $j$ if the magnitude of the correlation coefficient $\rho_{ij}$ between the $i^{th}$ and $j^{th}$ components of the random vector $\mathbf{X}$ is greater than $\rho$, i.e., if $|\rho_{ij}|\geq\rho$, $i\neq j$. We define $d_{k}$ to be the degree of vertex $k$ in the graph $\mathcal{G}_{\rho}(\boldsymbol{\Sigma})$:

[TABLE]

where $\#A$ is the size of set $A$. Finally, throughout the index $k$ will be associated with the $k$th variable.

Theorem IV.1

Let $\mathbf{\Sigma}$ be the covariance matrix of $\mathbf{X}$ or the rows of $\mathbb{X}$, $\mathbf{\Sigma}$ could be diagonal or non-diagonal depending on whether it is from the pre- or post-change data, respectively. Let $\mathbf{\Sigma}$ be row sparse of degree $j=o(p)$. Also let $p\to\infty$ and $\rho=\rho_{p}\to 1$ such that $p(1-\rho^{2})^{n/2}\to 0$. Then,

[TABLE]

where

[TABLE]

with

[TABLE]

[TABLE]

with $I_{x}(a,b)$ being the regularized incomplete beta function with parameter $a$ and $b$. The parameter $J_{k}$ is a positive real number that is a function of the asymptotic joint density of $\mathbf{X}$. Further, if the covariance matrix $\mathbf{\Sigma}$ of the p-variate vector $\mathbf{X}$ is block sparse of degree $j$, then

[TABLE]

In particular, if $\mathbf{\Sigma}$ is diagonal then

[TABLE]

Proof:

The events $\{V_{k}(\mathbb{X})\geq\rho\}$ and $\{d_{k}>0\}$ are equivalent. Hence

[TABLE]

The result now follows from Theorem IV.2 below. ∎

Theorem IV.2** ([16])**

Let $\mathbf{\Sigma}$ be row sparse of degree $j=o(p)$. Also let $p\to\infty$ and $\rho=\rho_{p}\to 1$ such that $p(1-\rho^{2})^{n/2}\to 0$. Then $d_{k}(\mathbb{X},\rho)$ converges in distribution to a Poisson random variable with rate parameter $\Lambda_{k}$.

Based on the asymptotic distribution obtained in Theorem IV.1, the large $p$ distribution of $V_{k}$ defined in (8) can be approximated by

[TABLE]

where $\Lambda_{k,\rho}$ is as defined in Theorem IV.1. Using [16, Thm. 3.1] it can easily be shown that the approximation error associated with (14) decays to zero at least as fast as $p(1-\rho^{2})^{n/2}\Delta$ where $\Delta$ is a dependency coefficient associated with the set of U-scores..

The distribution (12) is differentiable everywhere except at $\rho=0$ since $P(V_{k}(\mathbb{X})=0)>0$. For $\rho>0$ and large $p$, $V_{k}$ has density

[TABLE]

Note that $f_{V}$ in (13) is the density of the Lebesgue continuous component of the distribution (12) and that it integrates to $1-O(e^{-p^{2}})$ over $\rho\in(0,1]$.

For each $k$, the density $f_{V_{k}}$ is a member of a one-parameter exponential family with $J_{k}$ as the unknown parameter. Then, since $P_{0}(\rho)^{\prime}=-\frac{2(1-\rho^{2})^{(n-4)/2}}{B((n-2)/2,1/2)}$, where $B(a,b)$ is the beta function, we have for $\rho\in[0,1]$, the exponential family form of the density $f_{V}$ with parameter $J_{k}$:

[TABLE]

Note that the density $f_{V_{k}}$ corresponds to the $k$th variable, and the whole family $\{f_{V_{k}}\}$ depends on the index $k$ only through the parameter $J_{k}$, which can be different for different $k$. We thus rewrite the family in (14) (with complete specification for easy accessibility) as

[TABLE]

where

[TABLE]

In Fig. 1 we have plotted the density $f_{V}$ for various values of $J$ for $n=10$, and $p=100$. We note that for the chosen values of $n$ and $p$, the density is concentrated close to $1$, consistent with large values of $\rho$ arising in the purely high dimensional regime assumed in Theorem IV.1 [14]. Note that the variable $V$ tends to take on higher values as the parameter $J$ increases. This is the fundamental idea used for hub discovery in this paper.

In Fig. 2 we have plotted the normalized histogram of the samples $V_{1}(\mathbb{X})$ for multivariate Gaussian independent random vectors and compare it with the approximation $f_{V}(\cdot\;;1)$. Note that the approximation is quite accurate. Also, note that for this case the population $V_{1}=0$, but the samples are concentrated in the interval $(0.55,0.95)$. This is due to the fact that $p\gg n$.

V Global Summary Statistics

In [4] we proposed a global summary statistic and obtained its asymptotic distribution. Specifically, we showed that the statistic

[TABLE]

in the limit $p\to\infty$ (taken in a particular sense) and under sparsity assumptions similar to that in Theorem IV.1, has approximate density in an exponential family given by

[TABLE]

where

[TABLE]

The parameter $\theta$ is a positive parameter and has similar interpretation and nature as that of the parameter $J$ in (15). Specifically, if the rows of $\mathbb{X}$ have diagonal covariance $\mathbf{\Sigma}$ then $\theta=1$, otherwise $\theta\neq 1$. Compared to $V_{k}(\mathbb{X})$, $U(\mathbb{X})$ has a global view of the covariance matrix. If any of the $V_{k}$s undergo a change, this change will be reflected through a change in the distribution of $U$. However, it is not possible to isolate hubs based only on the statistic $U$. Sequential tests for detecting a change in correlation using the statistic $U$ is studied in detail in [4]. The sequential test from [4] will play a fundamental role in what follows.

VI Multivariate Quickest Hub Discovery

We now discuss how to use the results from the previous sections to detect and isolate correlation hubs.

Recall that we have a sequence of $p$-variate Gaussian distributed vectors $\{\mathbf{X}(k)\}$ with a change point model given in (2) and $p\gg 1$. Thus, the vectors are uncorrelated before change and correlated with covariance matrix $\mathbf{\Sigma}$ after change. We collect these vectors in batches of size $n$ and obtain a sequence of $n\times p$ random matrices $\{\mathbb{X}(m)\}$; see (7). For quick detection, we should choose a batch size $n$. Thus, we also have $p\gg n$.

Assuming that the post-change covariance matrix $\mathbf{\Sigma}$ is row sparse, a change in the law of $\mathbf{X}$ changes the law of the random matrix $\mathbb{X}$, which in turn affects the law of summary statistics $U(\mathbb{X})$ and $\{V_{k}(\mathbb{X})\}_{k=1}^{p}$. Since, the distributions of these statistic each belong to a one-parameter exponential family (15) and (18), a change in law within this family is a simple change in the value of the parameters $\{J_{k}\}$ and $\theta$. We propose to detect the change in correlation in $\mathbf{X}$ through a change in the parameters $\{J_{k}\}$ and $\theta$

From (15) and (18), it follows that we can choose $n>4$, because the densities are well defined for these values of $n$. Note also that here we assume that the change in correlation in $\mathbf{X}$ affects the parameters $\{J_{k}\}$ and $\theta$, otherwise our test is insensitive to the change in correlation.

VI-A Correlation Change Detection Using Global Statistic

Assuming that the post-change covariance matrix $\Sigma$ is row sparse, we map the random matrix sequence $\{\mathbb{X}(m)\}$ to the summary statistic sequence $\{U(\mathbb{X}(m))\}$. The random variables $\{U(\mathbb{X}(m))\}$ have approximate density given by (18) which is in an exponential family with parameter $\theta$ which is $1$ before the change and some value $\theta\neq 1$ after change. For simplicity, we refer to $U(\mathbb{X}(m))$ simply by $U(m)$. The QCD problem on the density $f_{\mathbf{X}}$, depicted in (2), is reduced to the following QCD problem

[TABLE]

Consider the following generalized likelihood ratio (GLR) based QCD tests (see [5], [4]) defined by the stopping time:

[TABLE]

where we assign $\tau_{\scriptscriptstyle\text{U}}=\infty$ if the right hand side of (23) is the empty set. Here $A_{u}>0$ and $\epsilon_{u}>0$ are user-defined parameters, and $\epsilon_{u}$ represents the minimum magnitude of change, away from $\theta=1$, that the user wishes to detect.

VI-B Correlation Change Detection Using Local Statistics

We can also detect a change in correlation using the local sample correlation statistics $\{V_{k}(\mathbb{X})\}_{k=1}^{p}$ defined in (8). For simplicity, $V_{k}(\mathbb{X}(m))$ is denoted simply by $V_{k}(m)$. For large $p$ and row sparse covariance, the random variables $\{V_{k}(m)\}$ have approximate density given by (15) which is in an exponential family with parameter $J_{k}$. This parameter is $1$ before the change and takes on some value $J_{k}\neq 1$ after the change. A natural locally pooled sqequential change detection test is to implement $p$ local sequential change detection tests in parallel, one for each variable $\{V_{k}(\mathbb{X})\}_{k=1}^{p}$ (equivalently, each variable in $\mathbf{X}$). An alarm is raised when any of the local tests raises an alarm.

Mathematically, the QCD problem on the density $f_{\mathbf{X}}$, depicted in (2), is reduced to the following family of QCD problems (one for each $k$):

[TABLE]

Consider the following GLR-based QCD tests defined by the family of local stopping times [5]:

[TABLE]

where, again we assign $\tau_{\scriptscriptstyle\text{V}}^{(k)}=\infty$ if the right hand side of (23) is the empty set, and $A_{v}>0$ and $\epsilon_{v}>0$ are user-defined parameters. This stopping time can be used to detect a change in the parameter $J_{k}$ for the $k$th variable. Define the overall local stopping time to be

[TABLE]

where

[TABLE]

is the GLR statistic. The test $\tau_{\scriptscriptstyle\text{V}}$ can be used to detect if the parameter of any of the $p$ variables is affected by the change. $\epsilon_{v}$ represents the minimum magnitude of change, away from $J_{k}=1$, that the user wishes to detect.

VI-C Joint Detection and Hub Discovery

For joint detection and hub discovery, we combine the local and global tests and propose to use the following maximum of the stopping times:

[TABLE]

Let $D_{k}(\mathbb{X}(1),\cdots,\mathbb{X}(\tau_{\scriptscriptstyle\text{HB}}))$ be the binary decision variable with value $1$ if variable $k$ is declared a hub. We use the following rule for hub discovery: fix a positive integer $q$ and set

[TABLE]

To appreciate the above rule, recall that a correlation strength hub is defined as the set of those variables that achieve the maximum in (17); see (4). Thus, there are at least two variables that are correlation strength hubs in the correlation graph. These hub variables in general have value of parameter $J_{k}$ that is higher than the parameter values for other variables. However, there may be many variables that have experienced a change in correlation, and their $V_{k}$ value, and hence the value of parameter $J_{k}$, might be close to the values attained by the hubs. To illustrate this idea, we have plotted in Fig. 4 estimated $\{J_{k}\}$ values from samples of a multivariate Gaussian vector with a row sparse covariance matrix (see Section VII on how the covariance matrix is generated). Here the hubs are the two variables with the highest values of $J_{k}$. But, as seen in the figure, there are many other variables that have experiences significant changes in correlation, and hence also have high values of $J_{k}$.

In Fig. 4, we have plotted the values of the GLR statistic $\{G_{k}\}$ at the time of stopping $\tau_{\scriptscriptstyle\text{HB}}$. As seen in the figure, in addition to the $G_{k}$ values of hubs, the values corresponding to other variables are also significant. Thus, a simple strategy like choosing the variable with the highest $G_{k}$ value as hub, will not work well.

However, what is clear from Fig. 4 and Fig. 4 is that high $G_{k}$ values at the time of stopping correspond to high $J_{k}$ values of the summary statistics, and these values together have significantly higher magnitude than the values for the rest of the variables. This motivates the pick-the-top approach. In fact, in practice, one can plot the $G_{k}$ values this way and identify all the variables that have experienced a change in correlation.

Picking the variables corresponding to the top $q$ values where $q>2$ ensures that the hubs are discovered with high probability, provided a slightly modified definition of hub discovery is used: the variables in the hub $\mathsf{H}$ are considered discovered if

[TABLE]

Indeed, this fact is used in Section VII for numerical computations, where it is shown that the proposed rule’s probability of false hub discovery goes to zero, as the false alarm rate goes to zero; see Fig. 5.

VI-D Delay and False Alarm Analysis

The approximate delay and false alarm performance of the stopping rule $\tau_{\scriptscriptstyle\text{HB}}$ is given in the theorem below.

Theorem VI.1

Fix any $\epsilon_{u}>0$ and $\epsilon_{v}>0$. If the densities (15) and (18) are the true distributions of the samples then we have the following:

For the stopping rule $\tau_{\scriptscriptstyle\text{HB}}$, the supremum in (6) is achieved at $\gamma=1$, i.e.,

[TABLE] 2. 2.

Setting $A_{u}=A_{v}=\log\beta$ ensures that as $\beta\to\infty$,

[TABLE] 3. 3.

For each possible true post-change parameters $J$ and $\theta$, with $|J-1|\geq\epsilon_{v}$, and $|\theta-1|\geq\epsilon_{u}$, as $\beta\to\infty$

[TABLE]

where $I(J)$ is the Kullback-Leibler divergence between the densities $f_{V}(\cdot;J)$ and $f_{V}(\cdot;1)$, and $I(\theta)$ is the Kullback-Leibler divergence between the densities $g(\cdot;\theta)$ and $g(\cdot;1)$.

Proof:

The proof follows from standard QCD arguments in the literature (see [17]) and because

[TABLE]

and

[TABLE]

∎

VII Numerical Results

In Fig. 5 we report the hub discovery performance of the proposed algorithm. We set $\mathbf{\Sigma}_{0}=I_{p}$, the identity matrix. We assume that the post-change covariance matrix $\mathbf{\Sigma}$ is a row-sparse matrix of degree $j$, obtained as follows. A $p\times p$ sample from the Wishart distribution is generated and some of the entries are forced to be zero in such a way that no row has more than $j$ non-zero elements. Specifically, we retain the top left $j\times j$ block of the matrix, and for each row $k$, $j+1\leq k\leq(p+j)/2$, all but the diagonal and the $(p+j+1-k)$th element is forced to zero. Each time an entry $(k,i)$ is set to zero, the entry $(i,k)$ is also set to zero, to maintain symmetry. Finally, a scaled diagonal matrix is also added to $\mathbf{\Sigma}$ to restore its positive definiteness. We set $n=10$, $p=100$, and $j=5$.

To implement $\tau_{\scriptscriptstyle\text{V}}$ we have chosen $\epsilon_{v}=1$, and we use the the maximum likelihood estimator for $J_{k}$ which, as a function of $m$ samples $(V_{k}(1),\cdots,V_{k}(m))$ from $f_{V}(\cdot,J)$, is given by

[TABLE]

Specifically,

[TABLE]

Similar techniques were used to implement $\tau_{\scriptscriptstyle\text{U}}$; see [4].

We consider three scenarios, each with two hubs:

[TABLE]

The parameter $J_{k}$s are estimated using $1000$ samples. For each case, the total number of false isolations is calculated using $1000$ sample paths by assuming $\gamma=1$. The parameter $q$ was chosen to be $10$. The definition (28) was used as a criterion calculating false isolation or false hub discovery. The results for all the three cases is plotted in Fig. 5 as a function of the threshold $A=A_{v}=A_{u}$. As shown in Fig. 5, the number of false discoveries goes to zero as the threshold is increased. Thus, empirical evidence suggests that the test is consistent.

In Fig. 6 we plot the delay ($\mathsf{E}_{1}[\tau_{\scriptscriptstyle\text{HB}}]$) vs the log of mean time to false alarm ($\log\mathsf{E}_{\infty}[\tau_{\scriptscriptstyle\text{HB}}]$) for various cases specified in the figure. The values in the figure are obtained by choosing different values of the threshold $A$ and estimating the delay by choosing the change point $\gamma=1$ and simulating the test for $500$ sample paths. The mean time to false alarm values are estimated by simulating the test for $1500$ sample paths.

VIII Conclusions and Future Work

We have introduced the quickest hub detection and localization problem for high dimensional correlation graphs. Local and global summary statistics were proposed based on purely high dimensional asymptotic theory for the presence of hub nodes, defined as vertices having either high degree or large correlations. The asymptotic densities of these summary statistics were shown to be in the one dimensional exponential family of densities. These densities were used to define a hybrid local and global quickest detection test for hubs and the asymptotic performance of this test was evaluated in terms of mean time to detect a change and mean time to false alarm.

IX Acknowledgments

This work was partially supported by the Consortium for Verification Technology under Department of Energy National Nuclear Security Administration award number DOE-NA0002534.