\subsubsection{Exponential tilting}Let the log moment generating function of $Z'_t$ be
\begin{equation*}
\psi(a) = \log\mathbb{E}_\infty[e^{a Z'_t}].
\end{equation*}
Define a family of new measures
\begin{equation*}
\frac{d\mathbb{P}_t}{d\mathbb{P}_\infty} = \exp\{a Z'_t - \psi(a)\},
\end{equation*}
where $\mathbb{P}_t$ denotes the new measure after the transformation.
The new measure takes the form of the exponential family, and $a$ can be viewed as the natural parameter. It can be verified that $\mathbb{P}_t$ is indeed a probability measure since
\begin{equation*}
\int d\mathbb{P}_t = \int \exp\{a Z'_t-\psi(a)\}d\mathbb{P} = 1.
\end{equation*}
It can also be shown that $\dot{\psi}(a)$ is the expected value of $Z'_t$ under $\mathbb{P}_t$, indeed
$$
\dot{\psi}(a)=\frac{\mathbb{E}_\infty[Z'_te^{a Z'_t}]}{\mathbb{E}_\infty[e^{a Z'_t}]} = \mathbb{E}_\infty[Z'_te^{a Z'_t-\psi(a)}] = \mathbb{E}_t[Z'_t],
$$
and similarly that $\ddot{\psi}(a)$ is equal to the variance under the tilted measure.

Under the assumption that $Z'_t$ is a standard Gaussian random variable, its log moment generating function satisfies $\psi(a) = a^2/2$.
Since we have the freedom to select $a$, we can set its value so that the mean under the tilted measure is equal to a given threshold $b$.
This can be done by choosing $a$ such that $\dot{\psi}(a)= b$ which yields $a = b$ a value independent from $t$. Given this choice, the transformed measure is given by $d\mathbb{P}_t = \exp(bZ'_t - b^2/2)d\mathbb{P}_\infty$. We also define, for each $t$, the log-likelihood ratio $\log(d\mathbb{P}_t/d\mathbb{P}_\infty)$ of the form
\begin{equation*}
\ell_t = b Z'_t - \frac{1}{2}b^2.
\end{equation*}

\subsubsection{Change-of-measure by the likelihood ratio identity} Now we convert the original problem of finding the small probability that the maximum of a random field exceeds a large threshold, to another problem: finding an alternative measure under which the event happens with a much higher probability.
Here, the alternative measure will be a mixture of simple exponential tilted measures. Define the maximum and the sum for likelihood ratio differences relative for a particular $t$:
\begin{equation*}
M_t = \max_{m \in \{1, \ldots, M\}} e^{\ell_m -\ell_t}, \quad S_t = \sum_{m \in \{1, \ldots, M\}} e^{\ell_m -\ell_t}.
\end{equation*}
Also define a re-centered likelihood ratio, which we call the global term
\begin{equation*}
\tilde{\ell}_t = b(Z'_t -b).
\end{equation*}
Then we have the following likelihood ratio identity:
\begin{align}
& \mathbb{P}_\infty(\max\limits_{1\le m \le M} Z'_m \geq b) = \mathbb{E}_\infty [\mathbbm{1}_{\{ \max\limits_{1\le m \le M} Z'_m \geq b\}}]\nonumber\\
&= \mathbb{E}_\infty \left[\frac{\sum_{t=1}^{M}e^{\ell_t}}{\sum_{n=1}^{M}e^{\ell_n}} \cdot \mathbbm{1}_{\{\max\limits_{1\le m \le M} Z'_m \geq b\}}\right]\nonumber\\
&=  \sum_{t=1}^M \mathbb{E}_\infty \left[\frac{e^{\ell_t}}{\sum_n e^{\ell_n}} \cdot \mathbbm{1}_{\{ \max\limits_{1\le m \le M} Z'_m \geq b\}} \right]\label{eq:tail_prob}\\
&=   \sum_{t=1}^M \mathbb{E}_t \left[\frac{1}{\sum_n e^{\ell_n}} \cdot \mathbbm{1}_{\{ \max\limits_{1\le m \le M} Z'_m \geq b\} }\right]\nonumber\\
& =  e^{-b^2/2} \sum_{t=1}^M \mathbb{E}_t \left[\frac{M_t}{S_t}e^{-(\tilde{\ell}_t +\log M_t)} \cdot \mathbbm{1}_{\{ \tilde{\ell}_t +\log M_t \geq 0 \} }\right].\nonumber
\end{align}
The last equation in (\ref{eq:tail_prob}) converts the tail probability to a product of two terms: a deterministic term $e^{-b^2/2}$ associated with the large deviation rate, and a sum of conditional expectations under the transformed measures. The conditional expectation involves a product of the ratio $M_t/S_t$, and an exponential function that depends on $\tilde{\ell}_t$, which plays the role of a weight. Under the new measure $\mathbb{P}_t$, $\tilde{\ell}_t$ has zero mean and variance equal to $b^2$ and it dominates the other term $\log M_t$, hence, the probability of exceeding zero is much higher. Next, we characterize the limiting ratio and the other factors precisely, by the localization theorem.

\subsubsection{Establish properties of local field and global term} In (\ref{eq:tail_prob}), our target probability has been decomposed into terms that only depend on (i) the local field $\{\ell_m - \ell_t\}, 1 \le m \le M$, which are the differences between the log-likelihood ratios with parameter $t$ and $m$, and (ii) the global term $\tilde{\ell}_t$, which is the centered and scaled likelihood ratio with parameter $t$. We need to first establish some useful properties of the local field and global term before applying the localization theorem. We will eventually show that the local field and the global term are asymptotically independent.

It is easy to show that under the tilted measure,
\begin{equation*}
\mathbb{E}_t[\tilde{\ell}_t] = 0, \quad \text{Var}_t[\tilde{\ell}_t] = b^2,
\end{equation*}
i.e., the global term $\tilde{\ell}_t$ is zero mean for any $t$, with variance diverging with $b$.

For the local field $\{\ell_m - \ell_t\}$, let $r_{m,t}$ denote the correlation between $Z'_m$ and $Z_t^\prime$ (given by Lemma\,\ref{CovZ}), then we have
\[
\begin{aligned}
\mathbb{E}_t(\ell_m - \ell_t) &= -b^2(1-r_{m,t}),\\
\mbox{Var}_t(\ell_m - \ell_t) &= 2b^2(1-r_{m,t}),\\
\mbox{Cov}_t(\ell_{m_1} - \ell_t,\ell_{m_2} - \ell_t) &= b^2(1+r_{m_1,m_2} - r_{m_1,t} - r_{m_2,t}).
\end{aligned}
\]

Since we assume $Z'_t$ is approximately Gaussian, the local field $\ell_m - \ell_t$ (or equivalently $b(Z'_m - Z'_t)$ ) and the global term $\tilde{\ell}_t$ (or equivalently $b(Z'_t - b)$) are also approximately Gaussian.
Therefore, when $|\delta|$ is small (i.e., in the neighborhood of zero), we can approximate the local field using a two-sided Gaussian random walk with drift $\mu^2/2$ and variance of the increment equal to $\mu^2$:
\begin{equation}
\ell_{t+\delta} - \ell_t \stackrel{\Delta}{=} \mu \sum\limits_{i=1}^{|\delta|} \vartheta_i - \frac{1}{2}\mu^2\delta, \delta = \pm1,\pm2,\ldots,
\label{eq:local_decomp}
\end{equation}
where $\vartheta_i$ are i.i.d. standard normal random variables and $\mu = b\sqrt{2\beta}$.

We have Lemma\,\ref{CovZ} to characterize the local correlation, which offers reasonably good approximation for ${\rm corr}(Z_t,Z_{t+\delta})$ and leads to some insights for our analysis.

\subsubsection{Approximation using localization theorem} Then we use the localization theorem (Theorem 5.1 in \cite{siegmund2010tail}) given below.
\begin{lemma}[Localization theorem]
Under mild assumptions, we have
\begin{equation*}
\begin{aligned}
&\lim_{\kappa\rightarrow\infty}\kappa^{1/2}\mathbb{E} \left[ \frac{M_\kappa}{S_\kappa}e^{-[\tilde{\ell}_\kappa+\log M_\kappa]}\cdot \mathbbm{1}_{\{\tilde{\ell}_\kappa+\log M_\kappa\ge 0\}} \right] \\
= & \sigma^{-1}\phi\left(\frac\mu\sigma \right)\mathbb{E} \left[\frac{M}{S}\right],
\end{aligned}
\end{equation*}
where $\phi(\cdot)$ is the density of the standard normal distribution.
\end{lemma}

Intuitively, the localization theorem says the following. To find the desired limit of $\mathbb{E}[\frac{M_\kappa}{S_\kappa}e^{-[\tilde{\ell}_\kappa+\log M_\kappa]}\mathbbm{1}_{\{ \tilde{\ell}_\kappa+\log M_\kappa\ge 0 \}} ]$ as $\kappa \rightarrow \infty$, one first approximates $M_\kappa$ and $S_\kappa$ by their localized versions, which are obtained by restricting the maximization and summation in a neighborhood of parameter values. Then one can show that the localized ratio $M_\kappa/S_\kappa$ is asymptotically independent of the global term $\tilde{\ell}_\kappa$ as $\kappa\rightarrow\infty$. The asymptotic analysis is then performed on the local field and the global term separately. The expected value of the localized ratio $\mathbb{E}[M_\kappa/S_\kappa]$ converges to a constant independent of $\kappa$, and the limiting conditional distribution of $\tilde{\ell}_\kappa$ can be found using the local central limit theorem. Thus, one can calculate the remaining conditional expectation involving $\tilde{\ell}_\kappa$.

After verifing the validity of the conditions in our setting, we have
\begin{align*}
&\mathbb{P}_\infty(T\le M)= \mathbb{P}_\infty\left( \max_{1\le t\le M} Z'_t \geq b \right) \\
&= e^{-b^2/2}\sum\limits_{t=1}^M\mathbb{E}_t \left[ \frac {M_t}{S_t}e^{-[\tilde{\ell}_t+\log M_t]}\cdot \mathbbm{1}_{\{ \tilde{\ell}_t+\log M_t\ge0\}} \right]\\
&= Mb\phi(b)\beta v(b\sqrt{2\beta})(1+o(1)).
\end{align*}
Once the cumulative distribution function of $T$ is available, we can approximate it as an exponential distribution with parameter $b\phi(b)\beta v(b\sqrt{2\beta})$ yielding as mean the value $1/[b\phi(b)\beta v(b\sqrt{2\beta})]$. Since $Z'_t$ is standardized, here the threshold $b$ need to be converted to the original threshold using a simple formula \[b' = \frac{b - [\sigma_{w,k}\mathbb{E}(W_1) + \mu_{w,k}]}{\sigma_{w,k}\sqrt{\rm{Var}(W_1)}}.\] This completes the proof.
\end{proof}

\begin{proof}[Proof of Proposition \ref{EDD}]
We first relate the largest eigenvalue procedure to a CUSUM procedure, note that
\begin{equation}
\lambda_{\max}(\hat{\Sigma}_{t,w}) = \max_{\|q\|=1}
q^\intercal \hat{\Sigma}_{t,w} q.
\label{eq:q}
\end{equation}
For each $q$, we have
\[
q^\intercal \hat{\Sigma}_{t,w} q =  \sum_{i=t-w+1}^t (q^\intercal x_i)^2.
\]
According to the Grothendieck's Inequality \cite{guedon2016community}, the $q$ that attains the maximum in equation (\ref{eq:q}) is very close to $u$ under the alternative. Therefore, assuming the optimal $q$ always equals to $u$ will only cause a small error but will bring great convenience to our analysis.

Now we have under $\mathbb{P}_\infty$, $q^\intercal x_i \sim \mathcal{N}(0, \sigma^2)$ and under $\mathbb{P}_0$, $q^\intercal x_i \sim \mathcal{N}(0, \sigma^2+\theta)$. Let $f_\infty$ denote the pdf of $\mathcal{N}(0, \sigma^2)$ and $f_0$ the pdf of $\mathcal{N}(0, \sigma^2+\theta)$. For each observation $y$, we can derive the one-sample log-likelihood ratio:
\[
\log\frac{f_0(y)}{f_\infty(y)} = -\frac12\log(1+\rho)+\frac1{2\sigma^2}\left(1-\frac1{1+\rho}\right)y^2.
\]
Define the CUSUM procedure
\begin{multline*}
\widetilde{T}= \inf\bigg\{ t : \max_{0\leq k<t}\sum_{i=k+1}^t \Big[\frac1{2\sigma^2}\left(1-\frac1{1+\rho}\right)(q^\intercal x_i)^2\\
-\frac{\log(1\!+\!\rho)}{2}\Big] \! \geq \! b'\bigg\},
\end{multline*}
where $b'=\frac1{2\sigma^2(1+\rho)}\left[b\rho-(1+\rho)\sigma^2\log(1+\rho)\right]$, we then have
\[
\mathbb{E}_0[T_{\rm E}] \geq \mathbb{E}_0[\widetilde{T}].
\]
Since $\widetilde{T}$ is a CUSUM procedure with
\[
\int\log\left[\frac{f_0(y)}{f_\infty(y)}\right]f_0(y)dy=-\frac12\log(1+\rho)+\frac{\rho}{2},
\]
from \cite{siegmund2013sequential} we have:
\[
\mathbb{E}_0[\widetilde{T}]=\frac{e^{-b^\prime}+b^\prime-1}{-\log(1+\rho)/2+\rho/2}.
\]
This completes the proof.
\end{proof}

\section{Introduction} \label{sec:intro}

\IEEEPARstart{D}{etecting} the changepoint from high-dimensional streaming data is a fundamental problem in various applications such as video surveillance, sensor networks, and seismic events detection. In various scenarios, the change happens to the covariance structure and can be represented as a linear subspace. For example, the covariance matrix may shift from an identity matrix to a spiked covariance model \cite{berthet2013optimal,johnstone2001distribution}.

Given a sequence of observed vectors $x_1, x_2, \ldots,$ where $x_t \in \mathbb{R}^k$ and $k$ is the signal dimension, there may be a changepoint time $\tau$ when the distribution of the data stream  changes. Our goal is to detect this change as quickly as possible from streaming (sequentially obtained) data using online techniques.
We are particularly interested in the structured change occurring in the signal covariance. We study two related settings: the {\it emerging subspace}, meaning that the change is a subspace emerging from a noisy background; and the {\it switching subspace}, meaning that the change is a switch in the direction of the subspace. The emerging subspace problem can arise, for instance, from coherent weak signal detection from seismic sensor arrays, and the switching subspace detection can be used for principal component analysis for streaming data. In these settings, the changes can be shown to be equivalent to a low-rank component added to the original covariance matrix. On the other hand, the switching subspace problem, as we will see, can be reduced to the emerging subspace problem, if we are willing to tolerate a performance loss. Therefore, we will focus on the analysis of the emerging subspace problem.

In this paper, we consider three detection procedures. We start with the exact CUSUM which is known to be optimum when we have complete knowledge of the pre- and post-change statistics and parameters. Since the post-change parameters are usually unknown we propose two alternatives to deal with the case of unknown parameters. Specifically we consider the largest eigenvalue procedure where we use as test statistic the largest eigenvalue of the sample covariance matrix of the data contained within a sliding time-window. This can be regarded as a straightforward extension of its offline counterpart \cite{berthet2013optimal}. The second method which we call the Subspace-CUSUM uses the structure of the exact CUSUM but in place of the parameters that are known it uses their estimates which are computed using data within, again, a sliding window. We perform a theoretical analysis of the largest eigenvalue procedure. A similar analysis, which is far more complicated and extended, is postponed for a future article. The three algorithms are compared using simulations

The rest of the paper is organized as follows. Section\,\ref{sec:formulation} details on the two problems of emerging and switching subspace and how they can be related. Section\,\ref{sec:detect} presents the three sequential change detection procedures. In Section\,\ref{sec:theo} we develop theoretical bounds for the average run length and the expected detection delay of the largest eigenvalue procedure. In Section\,\ref{sec:numerical} we presents numerical results and comparisons for the competing algorithms. Finally Section\,\ref{sec:conclusion} contains our concluding remarks.

\subsection{Related Work}

Classical approaches to covariance change detection usually consider generic settings without assuming any structure.
The CUSUM statistics can be derived if the pre-change and post-change distributions are known.
For the multivariate case, the Hotelling $T^2$ control chart is the traditional way to detect the covariance change. The determinant of the sample covariance matrix was also used in \cite{alt2004multivariate} to detect change of the determinant of the covariance matrix. A multivariate CUSUM based on likelihood functions of multivariate Gaussian is studied in \cite{healy1987note}  but it only considers the covariance change from $\Sigma$ to $c\Sigma$ for a constant $c$.
Offline change detection of covariance change from $\Sigma_1$ to $\Sigma_2$ is studied in \cite{chen2004statistical} using the Schwarz information criterion \cite{schwarz1978estimating}, where the changepoint location must satisfy certain regularity condition to ensure the existence of the maximum likelihood estimator. In \cite{arias2012detection} we find a hypothesis testing approach to detect a shift in an off-diagonal sub-matrix of the covariance matrix using likelihood ratios. Recently, \cite{jiao2018subspace} studies a CUSUM-like procedure for detection of switching subspaces, when the distributions (as well as the subspaces) before and after the changepoint are exactly known; this is different from our work since we assume the subspace after the change is unknown.

The most related work to our present effort is the hypothesis testing methods developed in \cite{berthet2013optimal}, which uses the largest eigenvalue of the sample covariance matrix to detect a sparse spiked covariance model given a fixed number of samples. The largest eigenvalue statistic is shown to be asymptotically minimax optimal for determining whether there exists a sparse and low-rank component in the offline setting. A natural sequential version of this idea is to use a sliding window and estimate the largest eigenvalue of the corresponding sample covariance matrix. However, this approach, under a sequential setting does not enjoy any form of (asymptotic) optimality.

A different test statistic, the so-called Kac-Rice statistic \cite{TaylorRice}, has been considered for testing spiked covariance model. The Kac-Rice statistic is the conditional survival function of the largest observed singular value conditioned on all other observed singular values, and is characterized by a simple asymptotic distribution (uniform in $[0, 1]$). However, the statistic involves the computation of an integral over the whole real line, and it is not clear how this can be carried over to the sequential formulation.

\section{Problem Setup} \label{sec:formulation}
We first introduce the spiked covariance model \cite{johnstone2001distribution}, which assumes that a small number of directions explain most of the variance. In particular, we consider the rank-one spiked covariance matrix, which is given by
\[
\Sigma = \sigma^2I_k + \theta uu^\intercal,
\]
where $I_k$ denotes an identity matrix of size $k$;  $\theta$ is the signal strength; $u \in \mathbb{R}^k$ represents a basis for the subspace with unit norm $\left\Vert u \right\Vert=1$; $\sigma^2$ is the noise variance, which will be considered known since it can be estimated from training data\footnote{In fact it is possible to consider $\sigma^2$ unknown as well and provide estimates of this parameter along with the necessary estimates of $u$. However, in order to simplify our presentation, we decided to consider $\sigma^2$ known.}. The Signal-to-Noise Ratio (SNR) is defined as
$\rho = \theta/\sigma^2$.

Formally, the \textit{emerging} subspace problem can be cast as follows:
\begin{equation}
\begin{array}{ll}
x_t  \stackrel{\text{iid}}{\sim} \mathcal{N}(0, \sigma^2I_k), &t = 1,2,\ldots,\tau, \\
x_t \stackrel{\text{iid}}{\sim} \mathcal{N}(0,  \sigma^2I_k+ \theta uu^\intercal), &t = \tau+1,\tau+2,\ldots
\end{array}
\label{eq:hypothesis}
\end{equation}
where $\tau$ is the unknown changepoint that we would like to detect from data that are acquired sequentially.

Similarly, the \textit{switching} subspace problem can be formulated as follows:
\begin{equation}
\begin{array}{ll}
x_t  \stackrel{\text{iid}}{\sim} \mathcal{N}(0, \sigma^2I_k + \theta u_1u_1^\intercal), &\!\!\!t = 1,2,\ldots,\tau, \\
x_t \stackrel{\text{iid}}{\sim} \mathcal{N}(0,  \sigma^2I_k  + \theta u_2 u_2^\intercal),&\!\!\!t = \tau+1,\tau+2,\ldots
\end{array}
\label{eq:hypothesis_switch}
\end{equation}
where $u_1, u_2 \in \mathbb R^{k}$ represent bases for the subspaces before and after the change, with $\|u_1\| = \|u_2\| = 1$ and $u_1$ is considered known. In both settings, our goal is to detect the change as quickly as possible.

The switching subspace problem \eqref{eq:hypothesis_switch} can be reduced into the emerging subspace problem \eqref{eq:hypothesis} by a simple data projection. Specifically, select any orthonormal matrix $Q \in \mathbb{R}^{(k-1)\times k}$ such that
\[ Q u_1 = 0, \quad QQ^\intercal = I_{k-1},\]
which means that all rows of $Q$ are orthogonal to $u_1$, they are orthogonal to each other and have unit norm. Then, using the matrix $Q$, we project each observation $x_t$ onto a $k-1$ dimensional space and obtain a new sequence
\[ y_t = Qx_t \in \mathbb{R}^{k-1}, t = 1,2,\ldots . \]
Then $y_t$ is a zero-mean random vector with covariance matrix $\sigma^2I_{k-1}$ before the change and $\sigma^2I_{k-1} + \theta Qu_2u_2^\intercal Q^\intercal$ after the change.
Let $u = Qu_2/\left\Vert Qu_2 \right\Vert$, and
\[ \tilde{\theta} = \theta \left\Vert Qu_2 \right\Vert^2 = \theta [ 1- (u_1^\intercal u_2)^2]. \]
Thus, problem \eqref{eq:hypothesis_switch} can be reduced to the following
\begin{equation}
\begin{array}{ll}
y_t  \stackrel{\text{iid}}{\sim} \mathcal{N}(0, \sigma^2I_{k-1}), &\!\!\!t = 1,2,\ldots,\tau, \\
y_t \stackrel{\text{iid}}{\sim} \mathcal{N}(0,  \sigma^2I_{k-1}  +\tilde{ \theta} u u^\intercal),&\!\!\!t = \tau+1,\tau+2,\ldots
\end{array}
\label{eq:hypothesis_switch_transform}
\end{equation}
Note that this way the switching subspace problem is reduced into the emerging subspace problem, where the new signal power $\tilde \theta$ depends on the angle between $u_1$ and $u_2$, which is consistent with our intuition.

We would like to emphasize that by projecting the observations onto a lower dimensional space we lose information, suggesting that the two versions of the problem are not equivalent. Indeed, the optimum detector for the transformed data in \eqref{eq:hypothesis_switch_transform} and the one for the original data in \eqref{eq:hypothesis_switch} do not coincide. This can be easily verified by computing the corresponding CUSUM tests and their optimum performance. Despite this difference, it is clear that with the proposed approach we put both problems under the same framework, offering, as we will see, computationally simple methods to solve the original problem in \eqref{eq:hypothesis_switch}. Consequently, in the following analysis, we focus solely on problem \eqref{eq:hypothesis}.

\section{Detection Procedures} \label{sec:detect}

As we mentioned before, we consider three methods: The exact CUSUM procedure where all parameters are considered known, the largest eigenvalue procedure, and the Subspace-CUSUM procedure. It is clear that since CUSUM is optimum it will be regarded as a point of reference for the other two approaches. We first introduce some necessary notation.
Denote with $\mathbb{P}_\tau$ and $\mathbb{E}_\tau$ the probability and expectation induced when there is a changepoint at the \textit{deterministic} time $\tau$. Under this definition $\mathbb{P}_\infty$ and $\mathbb{E}_\infty$ is the probability and the expectation under the nominal regime (change never happens) while $\mathbb{P}_0$ and $\mathbb{E}_0$ the probability and expectation under the alternative regime (change happens before we take any data).

\subsection{Optimal CUSUM Procedure}\label{detect:exactcusum}

The CUSUM procedure \cite{page1954continuous,siegmund2013sequential} is the most popular sequential test for change detection. When the observations are i.i.d. before and after the change, CUSUM is known to be exactly optimum \cite{moustakides1986optimal} in the sense that it solves a very well defined constrained optimization problem introduced in \cite{lorden1971procedures}. However, the CUSUM procedure can be applied only when we have exact knowledge of the pre- and post-change distributions. Thus, for our problem, it requires complete specification of all parameters namely the subspace $u$, noise power $\sigma^2$ and SNR $\rho$.

To derive the CUSUM procedure, let $f_\infty(\cdot),f_0(\cdot)$ denote the pre- and post-change probability density function (pdf) of the observations. Then the CUSUM statistics is defined by maximizing the log-likelihood ratio statistic over all possible changepoint locations
\begin{equation}
S_ t =  \max_{1 \leq j  \leq t } \sum_{i=j}^t \log \frac{f_0(x_i)}{f_{\infty}(x_i)},
\label{cusum_exact}
\end{equation}
which has the recursive implementation
\begin{equation}
S_t = (S_{t-1})^+ + \log \frac{f_0(x_t)}{f_{\infty}(x_t)},
\label{cusum_recursive}
\end{equation}
that enables its efficient calculation \cite{moustakides1986optimal}. The CUSUM stopping time in turn is defined as
\begin{equation}
T_{\rm C} = \inf\{t>0: S_t \geq b\},
\label{cusum_procedure}
\end{equation}
where $b$ is a threshold selected to meet a suitable false alarm constraint.
For our problem of interest we can derive that
\begin{multline}
\log \frac{f_0(x_t)}{f_\infty(x_t)}  = \log\bigg[ \frac{[(2\pi)^p|\sigma^2I_k + \theta u u^\intercal|]^{-1/2}}{[(2\pi)^p\sigma^{2p}]^{-1/2}} \\
 \times
 \frac{\exp\{ - x_t^\intercal(\sigma^2I_k+\theta u u^\intercal)^{-1}x_t / 2  \} }{\exp\{- x_t^\intercal x_t /(2\sigma^2) \} } \bigg] \\
= \log\left[
\left|I_k + \rho u u^\intercal \right|^{-\frac 1 2}
\exp\left\{ \frac 1 2 \frac{\theta}{\theta+\sigma^2} \frac{(u^\intercal x_t)^2}{\sigma^2}  \right\} \right]  \\
 =\frac{\rho}{2\sigma^2(1+\rho)} \Big\{(u^\intercal x_t)^2- \sigma^2 \left(1+ \frac1{\rho}\right)\log(1+\rho)\Big\}.
\label{expansion}
\end{multline}
The second equality is due to the matrix inversion lemma \cite{woodbury1950inverting} that allows us to write
\[
(\sigma^2I_k+\theta u u^\intercal)^{-1} =
 \frac1{\sigma^2}I_k - \frac{\theta}{\theta+\sigma^2}\frac{uu^\intercal}{\sigma^2},
\]
which, after substitution into the equation, yields the desired result. Note that the multiplicative factor $\rho/[2\sigma^2(1+\rho)]$ is positive, so we can omit it from the log-likelihood ratio when forming the CUSUM statistic. This leads to
\begin{equation}
S_t = (S_{t-1})^+ + (u^\intercal x_t)^2 - \sigma^2\left(1+\frac{1}{\rho}\right)\log(1+\rho).
\label{cusum_recur}
\end{equation}

\begin{remark}\label{drift_property_1}
We can show that the increment in \eqref{cusum_recur}, i.e.,
\[(u^\intercal x_t)^2 - \sigma^2\left(1+\frac{1}{\rho}\right)\log(1+\rho),
\] has the following property: its expected value is negative under the pre-change and positive under the post-change probability measure. The proof relies on a simple argument based on Jensen's inequality. Due to this property, before the change, the CUSUM statistics $S_t$ will oscillate near $0$ while it will exhibit, on average, a positive linear drift after the occurrence of the change forcing it, eventually, to hit or exceed the threshold.
\end{remark}

\subsection{Largest Eigenvalue Procedure} \label{detect:eigen}

Motivated by the off-line test in \cite{berthet2013optimal}, a natural strategy to detect the change is to use the largest eigenvalue of the sample covariance matrix. Under the sequential setting, we adopt a sliding window approach and form the sample covariance matrix using observations that lie within a time window of length $w$. For each time $t>0$, the \textit{un-normalized} sample covariance matrix using the available samples is given by
\begin{equation}
\hat{\Sigma}_{t,\min\{t,w\}}=\left\{
\begin{array}{l}
x_1 x_1^\intercal+\cdots+x_t x_t^\intercal,~\text{for}~t<w\\
x_{t-w+1} x_{t-w+1}^\intercal+\cdots+x_t x_t^\intercal,~\text{for}~t\geq w.
\end{array}\right.
\label{sampleCov}
\end{equation}
We note that for $t=1$ the matrix contains a single outer product and as time progresses the number of outer products increases linearly until it reaches $w$. After this point, namely for $t\geq w$, the number of outer products remains equal to $w$.

Let $\lambda_{\max}(X)$ denote the largest eigenvalue of a symmetric matrix $X$. We define the {\it largest eigenvalue procedure}, as the one that stops according to the following rule:
\begin{equation}
\label{eigenvalue_procedure}
T_{\rm E} = \inf\{t>0: \lambda_{\max}(\hat{\Sigma}_{t,\min\{t,w\}}) \geq b \},
\end{equation}
where $b > 0$ is a constant threshold selected to meet a suitable false alarm constraint.
We need to emphasize that we do not divide by $\min\{t,w\}$ when forming the un-normalized sample covariance matrix. As we explain in Section\,\ref{sec:metric}, it is better for $T_{\rm E}$ to always divide by $w$ instead of $\min\{t,w\}$. Consequently, we can omit the normalization with $w$ from our detection statistics by absorbing it into the threshold.

\subsection{Subspace-CUSUM Procedure}\label{sec:S-CUSUM}

Usually the subspace $u$ and SNR $\rho$ are unknown. In this case it is impossible to form the exact CUSUM statistic depicted in \eqref{cusum_recur}. One option is to estimate the unknown parameters and substitute them back into the likelihood function. Here we propose to estimate only $u$ and call $d=\sigma^2(1+1/\rho)\log(1+\rho)$ which leads to the following Subspace-CUSUM update
\begin{equation}
\mathcal S_t = (\mathcal S_{t-1})^{+} + (\hat{u}_{t+w}^\intercal x_t)^2 - d,
\label{sscusum_update}
\end{equation}
We denote the estimator of $u$ as $\hat{u}_{t+w}$. This is because at time $t$ the estimate will rely on the data $x_{t+w},\ldots,x_{t+1}$ that are in the ``future'' of $t$. Practically, this is always possible by properly delaying our data by $w$ samples. Stopping occurs similarly to CUSUM, that is
\[
T_{\rm SC} = \inf\{t>0: \mathcal S_t \geq b\}.
\]
Of course, in order to be fair, at the time of stopping we must make the appropriate correction, namely if $\mathcal{S}_t$ exceeds the threshold at $t$ for the first time, then the actual stopping takes place at $t+w$. The reason we use estimates based on ``future'' data is to make $x_t$ and $\hat{u}_{t+w}$ independent which in turn will help us decide what is the appropriate choice for the drift constant $d$.

For the drift parameter $d$ we need the following double inequality to be true
\begin{equation}
\mathbb{E}_\infty[(\hat{u}_{t+w}^\intercal x_t)^2] < d < \mathbb{E}_0[(\hat{u}_{t+w}^\intercal x_t)^2].
\label{expected_d}
\end{equation}
With \eqref{expected_d} we can guarantee that $\mathcal{S}_t$ mimics the behavior of the exact CUSUM statistic $S_t$ mentioned in Remark\,\ref{drift_property_1}, namely, it exhibits a negative drift before and a positive after the change. To apply \eqref{sscusum_update}, we need to specify $d$ and of course provide the estimate $\hat u_{t+w}$. Regarding the latter we simply use the \textit{unit-norm} eigenvector corresponding to the largest eigenvalue of $\hat \Sigma_{t+w,w}$ depicted in \eqref{sampleCov}. As we mentioned, the main advantage of using $\hat \Sigma_{t+w,w}$ is that it provides estimates $\hat u_{t+w}$ that are \textit{independent} from $x_t$. This independence property allows for the straightforward computation of the two expectations in \eqref{expected_d} and contributes towards the proper selection of $d$. Note that under the pre-change distribution we can write
\begin{multline}
\mathbb{E}_\infty[(\hat{u}_{t+w}^\intercal x_t)^2] =  \mathbb{E}_\infty[ \hat{u}_{t+w}^\intercal\mathbb{E}_\infty[ x_t x_t^\intercal] \hat{u}_{t+w} ]\\
 = \sigma^2\mathbb{E}_\infty [\hat{u}_{t+w}^\intercal\hat{u}_{t+w}]= \sigma^2,
\label{drift_infty}
\end{multline}
where the first equation is due to the independence of $x_t$ and $\hat{u}_{t+w}$, the next one due to $x_t$ having covariance $\sigma^2I_k$ and the last equality due to $\hat{u}_{t+w}$ being of unit norm.

Under the post-change regime, we need to specify the statistical behavior of $\hat u_{t+w}$ for the computation of $\mathbb{E}_0[(\hat{u}_{t+w}^\intercal x_t)^2]$. We will assume that the window size $w$ is sufficiently large so that Central Limit Theorem (CLT) approximations \cite{anderson1963asymptotic, paul2007asymptotics} are possible for $\hat u_{t+w}$. The required result appears in the next lemma.
\begin{lemma
\label{lem:eigenvector}
Suppose vectors $x_1,\ldots,x_w$ are of dimension $k$ and follow the distribution  $\mathcal{N}(0,  \sigma^2I_k + \theta u u^\intercal)$. Let $\hat{\varphi}_w$ be the eigenvector corresponding to the largest eigenvalue of the sample covariance matrix $(1/w)(x_1x_1^\intercal+\cdots+x_wx_w^\intercal)$, then, as $w \rightarrow \infty$, we have the following CLT version for $\hat{\varphi}_w$
\[
\sqrt{w}(\hat{\varphi}_w - u)  \rightarrow \mathcal N\left(0, \frac{ 1+\rho}{\rho^2}  ( I_k - uu^\intercal)   \right).
\]
\end{lemma}

\begin{proof} The proof is detailed in the Appendix.
\end{proof}
Lemma\,\ref{lem:eigenvector} provides an asymptotic statistical description of the \textit{un-normalized} estimate of $u$. More precisely it characterizes the estimation error $v_w=\hat{\varphi}_w-u$. In our case we estimate the eigenvector from the matrix $\hat \Sigma_{t+w,w}$ but, as mentioned before, we adopt a \textit{normalized} (unit norm) version $\hat{u}_{t}$. Therefore if we fix $w$ at a sufficiently large value and $v_t$ denotes the estimation error of the un-normalized estimate at time $t$ then, from Lemma\,\ref{lem:eigenvector}, we can deduce
\begin{multline}\label{hat_u_separate}
 \hat{u}_{t+w} = \frac{\hat{\varphi}_{t+w}}{\|\hat{\varphi}_{t+w}\|}=\frac{u + v_{t+w}}{\left\Vert u + v_{t+w} \right\Vert},\\
   v_{t+w} \sim \mathcal{N}\left(0, \frac{1+\rho}{w\rho^2}(I_k-uu^\intercal)\right).
 \end{multline}
Note that $v_{t+w}$ is also independent from $x_t$ and orthogonal to $u$, the latter being true because the covariance matrix of $v_{t+w}$ is $I_k -uu^\intercal$. This implies $ \left\Vert u + v_{t+w}\right\Vert^2 = 1+\left\Vert v_{t+w}\right\Vert^2 $. Combining the above results, we have
\begin{multline*}
\mathbb{E}_0\left[(\hat{u}_{t+w}^\intercal x_t)^2 \right] =\\
\mathbb{E}_0\left[\frac{(u + v_{t+w})^\intercal (\sigma^2I_k + \theta uu^\intercal)(u + v_{t+w})} {1 + \left\Vert v_{t+w}\right\Vert^2} \right].
\end{multline*}
Because $v_{t+w}^\intercal u = 0, \|u\|=1$, the above expression simplifies to
\begin{multline}\label{expect_0}
\mathbb{E}_0\left[ \frac{\sigma^2 +\theta + \sigma^2 \left\Vert v_{t+w}\right\Vert^2 } { 1+\left\Vert v_{t+w}\right\Vert^2  } \right] = \sigma^2 +\theta\,\mathbb{E}_0\left[\frac{1} { 1+\left\Vert v_{t+w}\right\Vert^2 } \right]\\
=\sigma^2 +\theta\,\left\{1-\mathbb{E}_0[\|v_{t+w}\|^2]+\mathbb{E}_0\left[\frac{\|v_{t+w}\|^4}{1+\|v_{t+w}\|^2}\right]\right\}.
\end{multline}
For the two expectations in \eqref{expect_0}, using the Gaussian approximation from Lemma\,\ref{lem:eigenvector}, we have
\[
 \mathbb{E}_0 [\left\Vert v_{t+w}\right\Vert^2] = \frac{1+\rho}{w\rho^2}(k-1),
\]
and
$$
\mathbb{E}_0\left[\frac{\|v_{t+w}\|^4}{1+\|v_{t+w}\|^2}\right]\leq\mathbb{E}_0[\|v_{t+w}\|^4]=\frac{(1+\rho)^2}{w^2\rho^4}(k^2-1).
$$
Consequently
 \begin{equation}\label{lower_bound}
\mathbb{E}_0\left[(\hat{u}_{t+w}^\intercal x_t)^2\right]  =  \sigma^2(1+\rho)\left( 1 - \frac{k-1}{w\rho} +o\left(\frac{1}{w}\right)\right),
 \end{equation}
with the $o(\cdot)$ term being negligible compared to the other two when $\frac{k}{w}\ll1$.

Consider now the case where $\rho$ is unknown but exceeds some pre-set minimal SNR $\rho_{\min}$. From the above derivation, given the worst-case SNR and an estimation for the noise variance $\hat{\sigma}^2$, we can give a lower bound for $\mathbb{E}_0[(\hat{u}_{t+w}^\intercal x_t)^2]$. Consequently, the drift $d$ can be anything between
$\hat{\sigma}^2$
  and
  $\hat{\sigma}^2(1+\rho_{\min}) (1-  (k-1)/(w\rho_{\min}) )$ where, we observe, that the latter quantity exceeds $\hat{\sigma}^2$ when $w > (k-1)(1+\rho_{\min})/\rho_{\min}^2$.
 Below, for simplicity, for $d$ we use the average of the two bounds.

Alternatively, and in particular when $w$ does not satisfy $w\gg k$, we can estimate $\mathbb{E}_0[(\hat{u}_{t+w}^\intercal x_t)^2]$ by Monto Carlo simulation.
This method requires: (i) estimating the noise level $\hat{\sigma}^2$, which can be obtained from training data without a changepoint; (ii) the pre-set worst-case SNR $\rho_{\min}$; (iii) a unit norm vector $u_0$ that is generated randomly. Under the nominal regime we have $\mathbb{E}_\infty[(\hat{u}_{t+w}^\intercal x_t)^2]=\hat{\sigma}^2$. Under the alternative $\mathbb{E}_0[(\hat{u}_{t+w}^\intercal x_t)^2]$ depends only on the SNR $\rho$ as shown in \eqref{lower_bound}. We can therefore simulate the worst-case scenario $\rho_{\min}$ using the randomly generated vector $u_0$ by generating samples from the distribution $\mathcal{N}(0,  \hat{\sigma}^2I_k + \rho_{\min}u_0u_0^\intercal )$.

Even though the average of the update in \eqref{sscusum_update} does not depend on $u$, the computation of the test statistic $\mathcal{S}_t$ requires the estimate $\hat{u}_{t+w}$ of the eigenvector. This can be accomplished by applying the singular value decomposition (SVD) (or the power method \cite{mises1929praktische}) on the un-normalized sample covariance matrix $\hat \Sigma_{t+w,w}$.
\begin{remark}
An alternative possibility is to use the generalized likelihood ratio (GLR) statistic, where both $\rho$ and $u$ are estimated for each possible change location $\kappa$. The GLR statistic is
\[
\max_{\kappa<t} \left\{-\frac{t-\kappa}{2} \log(1+\hat{\rho}_{\kappa,t}) + \frac1{2\hat{\sigma}^2}\frac{\hat{\rho}_{\kappa,t}}{1+\hat{\rho}_{\kappa,t}}\sum_{i = \kappa+1}^t (\hat{u}_{\kappa,t}^\intercal x_i )^2 \right\},
\]
where $\hat{\rho}_{\kappa,t}$, $\hat{u}_{\kappa,t}$ are estimated from samples $\{x_i\}_{i = \kappa+1}^t$. However, this computation is more intensive since there is no recursive implementation for the GLR statistic, furthermore it requires growing memory\footnote{There are finite memory versions \cite{lai1999efficient} which, unfortunately, are equally complicated in their implementation.}. Therefore, we do not consider the GLR statistic in this paper.
\end{remark}

\subsection{Calibration} \label{sec:metric}

To fairly compare the detection procedures discussed in the previous section we need to properly calibrate them. Clearly the calibration process must be consistent with the performance measure we are interested in. It is exactly this point we are discussing next.

For a given stopping time $T$ we measure false alarms through the Average Run Length (ARL) expressed with $\mathbb{E}_\infty[T]$. For the detection capability of $T$ we use the (worst-case) Expected Detection Delay (EDD) proposed by Lorden \cite{lorden1971procedures}
\begin{equation}
\sup\limits_{\tau \geq 0} \supess \mathbb{E}_\tau[ (T - \tau)^{+} | T>\tau, x_1, \ldots, x_\tau],
\label{eq:lorden}
\end{equation}
which considers the worst possible data before the change (expressed through the $\supess$) and the worst possible change-time $\tau$.

We now consider scenarios that lead to the worst-case detection delay. For the largest eigenvalue procedure, assume $1\leq t-w+1\leq\tau< t$. Since for the detection we use $\lambda_{\max} (\hat \Sigma_{t,w})$ and compare it to a threshold, it is clear that the worst-case data before $\tau$ are the ones that will make $\lambda_{\max}$ as small as possible. We observe that
\begin{multline}
\lambda_{\max} (\hat{\Sigma}_{t,w})\\
=\lambda_{\max} ( x_{t-w+1} x_{t-w+1}^\intercal + \cdots +
x_{\tau} x_{\tau}^\intercal+\cdots+ x_{t} x_t^\intercal)\\
\geq\lambda_{\max} ( x_{\tau+1} x_{\tau+1}^\intercal +
\cdots+ x_{t} x_{t}^\intercal)
=\lambda_{\max} (\hat{\Sigma}_{t,t-\tau}),
\label{eq:worst}
\end{multline}
which corresponds to the data $x_{t-w+1},\ldots,x_\tau$, before the change, being all equal to zero. In fact, the worst-case scenario at any time instant $\tau$ is equivalent to forgetting all data before and including $\tau$ and restarting the procedure from $\tau+1$ using, initially, one, then two, etc.~outer products in the un-normalized sample covariance matrix, exactly as we do when we start at time 0. Due to stationarity, this suggests that we can limit ourselves to the case $\tau=0$ and compute $\mathbb{E}_0[T_{\rm E}]$ and this will constitute the worst-case EDD. Furthermore, the fact that in the beginning we do not normalize with the number of outer products, is beneficial for $T_{\rm E}$ since it improves its ARL.

We should emphasize that if we do not force the data before the change to become zero and use simulations to evaluate the detector with a change occurring at some time different from 0, then it is possible to arrive at misleading conclusions. Indeed, it is not uncommon this test to appear outperforming the exact CUSUM test for low ARL values. Of course this is impossible since the exact CUSUM is optimum for \textit{any} ARL in the sense that it minimizes the \textit{worst-case} EDD depicted in \eqref{eq:lorden}.

Let us now consider the worst-case scenario for Subspace-CUSUM. We observe that
\[
\mathcal S_t = (\mathcal S_{t-1})^{+} + (\hat{u}_{t+w}^\intercal x_t)^2 - d \geq 0 +(\hat{u}_{t+w}^\intercal x_t)^2 - d,
\]
suggesting that when $\mathcal S_t$ restarts this is the worst it can happen for the detection delay. We therefore understand that the well-known property of the worst-case scenario in the exact CUSUM carries over to Subspace-CUSUM. Again, because of stationarity, this allows us to fix the changetime at $\tau=0$. Of course, as mentioned before, because $\hat{u}_{t+w}$ uses data coming from the future of $t$, if our detector stops at some time $t$ (namely when for the first time we experience $\mathcal{S}_t\geq b$) then the \textit{actual} time of stopping must be corrected to $t+w$. A similar correction is not necessary for CUSUM because this test has the exact information for all parameters.

Threshold $b$ is chosen so that the ARL meets a pre-specified value. In practice, $b$ is determined by simulation. A very convenient tool in accelerating the estimation of ARL (which is usually large) is the usage of the following formula that connects the ARL of CUSUM to the average of the SPRT stopping time \cite{siegmund2013sequential}
\begin{equation}
\mathbb{E}_\infty [T_{\rm C}] = \frac{\mathbb{E}_\infty[T_{\rm SPRT}]}{\mathbb{P}_\infty (S_{T_{\rm SPRT}} \ge b)},
\label{eq:CUSUM-SPRT}
\end{equation}
where the SPRT stopping time is defined as
\[
T_{\rm SPRT} = \inf\{t>0: S_t \notin [0, b]\}.
\]
The validity of this formula relies on the CUSUM property that after each restart, $S_t$ is independent from the data before the time of the restart. Unfortunately this key characteristic is no longer true in the proposed Subspace-CUSUM scheme due to the fact that $\hat{u}_{t+w}$ uses data from the future of $t$. We could, however, argue that this dependence is weak. Indeed, as we have seen in Lemma\,\ref{lem:eigenvector}, each $\hat{u}_t$ is basically equal to $u$ plus some small random perturbation (estimation error with power of the order of $1/w$), with these perturbations being practically independent in time. As we observed with numerous simulations, estimating the ARL directly and through \eqref{eq:CUSUM-SPRT} (with $S_t$ replaced by $\mathcal{S}_t$), results in almost indistinguishable values even for moderate window sizes $w$. This suggests that we can use \eqref{eq:CUSUM-SPRT} to estimate the ARL of the Subspace-CUSUM as well. As we mentioned, in the final result we need to add $w$ to account for the future data used by the estimate $\hat{u}_{t+w}$.

\section{Analysis of the Largest Eigenvalue Procedure} \label{sec:theo}

It is clear that, in this work we are interested in promoting the Subspace-CUSUM detection procedure for the change detection problem of interest. Therefore, it would have been very supportive to this method to offer a theoretical analysis and derive formulas for the corresponding ARL and EDD. Even though such an analysis is possible it is unfortunately overly lengthy, for this reason, we postpone its presentation for a future publication. In this section we intend to characterize the ARL and EDD of the largest eigenvalue procedure which turns out to be simpler. In doing so we will also introduce some of the mathematical tools we are going to use in the (future) analysis of the Subspace-CUSUM.

\subsection{Link with Random Matrix Theory}
Since the study of ARL requires the understanding of the property of the largest eigenvalue under the null, i.e., the samples are i.i.d. Gaussian random vectors with zero-mean and identity covariance matrix, we first review some related results from random matrix theory.

There has been an extensive literature on the distribution of the largest eigenvalue of the sample covariance matrix, see, e.g., \cite{johnstone2001distribution, yin1988limit, baik2006eigenvalues}. There are two kinds of results typically available for eigenvalue distributions: the so-called {\it bulk} \cite{edelman2013random}, which treats a continuum of eigenvalues, and the {\it extremes}, which are the (first few) largest and smallest eigenvalues. Assume there are $w$ samples which are $k$-dimensional Gaussian random vectors with zero-mean and identity covariance matrix. Let $\hat{\Sigma}_w=\sum_{i=1}^w x_i x_i^\intercal$ denote the un-normalized sample covariance matrix. If $k/w \rightarrow \gamma >0$, the largest eigenvalue of the sample covariance matrix converges to $(1+\sqrt{\gamma})^2$ almost surely \cite{geman1980limit}.
To characterize the distribution of the largest eigenvalue, \cite{johnstone2001distribution} uses the  Tracy-Widom law \cite{TracyWidom96}. Define the center and scaling constants
\begin{equation}
\begin{aligned}
 \mu_{w,k}&=\big(\sqrt{w-1}+\sqrt{k}\big)^2, \\
 \sigma_{w,k}&=\big(\sqrt{w-1}+\sqrt{k}\big)\left(\frac1{\sqrt{w-1}}+\frac1{\sqrt{k}}\right)^{1/3}.
 \end{aligned}
 \label{eq:parameter}
 \end{equation}
 If $k/w \rightarrow \gamma <1$, then the centered and scaled largest eigenvalue converges in distribution to a random variable $W_1$ with the so-called Tracy-Widom law of order one \cite{johnstone2001distribution}:
\begin{equation}
 \frac{\lambda_{\max}(\hat{\Sigma}_{w})-\mu_{w,k}}{\sigma_{w,k}} \rightarrow W_1.
 \label{eq:tw}
 \end{equation}
The Tracy-Widom law can be described in terms of a partial differential equation and the Airy function, and its tail can be computed numerically (using for example the R-package {\small \textsf{RMTstat}}).

\subsection{Approximation of ARL Ignoring Temporal Correlation}
If we ignore the temporal correlation of the largest eigenvalues produced by the sliding window, we can obtain a simple approximation for the ARL. If we call $p = \mathbb{P}_\infty(\lambda_{\max}(\hat{\Sigma}_{t,w} ) > b)$ for $t \geq w$ then the probability to stop at $t$ is geometric and it is easy to see that the ARL can be expressed as $1/p$. Clearly, to obtain this result we must assume that $\mathbb{P}_\infty(\lambda_{\max}(\hat{\Sigma}_{t,w} ) > b)=p$ for $t<w$ as well, which is clearly not true. Since for $t<w$ the un-normalized sample covariance has less than $w$ terms, the corresponding probability is smaller than $p$. This suggests that $1/p$ is actually a lower bound to the ARL while $w+1/p$ an upper bound. If $w\ll 1/p$ then approximating the ARL with $1/p$ is quite acceptable. We can use the Tracy-Widom law to obtain an asymptotic expression relating the ARL with the threshold $b$. The desired formula is depicted in the following proposition.

\begin{proposition}[Approximation of ARL by ignoring temporal correlation]
For any $0<p\ll1$ we have $\mathbb{E}_\infty[T_{\rm E}]\approx 1/p$, if we select
\begin{equation}
b = \sigma_{w,k}b_p + \mu_{w,k},
\label{eq:threshold1}
\end{equation}
where $b_p$ denotes the $p$-upper-percentage point of $W_1$ namely $\mathbb{P}(W_1\geq b_p)=p$.
\label{ARL1}
\end{proposition}
\begin{proof}
The proof is straightforward and therefore ommitted.
\end{proof}

\subsection{Approximation of ARL Including Temporal Correlation}
Now we aim to capture the temporal correlation between detection statistics due to overlapping time windows. We leverage a proof technique developed in \cite{SiegmundYakirZhang2010}, which can obtain satisfactory approximation for the tail probability of the maximum of a random field.
For each $\hat{\Sigma}_{t,w}$,  define
\begin{equation}
Z_{t}=\lambda_{\max}(\hat{\Sigma}_{t,w}).
\end{equation}
Fig. \ref{overlapping} illustrates the overlap of two sample covariance matrices and provides necessary notation.
\begin{figure}[H]
\centerline{
\includegraphics[width = 0.4\textwidth]{fig1}

\caption{Illustration of two sample covariance matrices.}
\label{overlapping}
\end{figure}
\noindent We note that for any given $M > 0$,
\[
\mathbb{P}_\infty(T\leq M) = \mathbb{P}_\infty\left(\max_{1\leq t\leq M} Z_t  \geq b \right),
\]
which is the max over a set of correlated variables $\{ Z_t\}_{t=1}^M$.
Capturing the temporal dependence of $\{Z_t\}$ is challenging. For our analysis we recall Pearson's correlation between two random variables as
\[
{\rm corr}(X, Y) = \frac{\mathbb{E}[X Y] - \mathbb E[X] \mathbb E[Y]}{\sqrt{{\rm Var}(X)}\sqrt{{\rm Var}(Y)}}.
\]
We then have the following lemma that addresses the problem of interest.
\begin{lemma}[Approximation of local correlation]\label{CovZ}
Let
\begin{equation}
\beta = 1+\frac{2k^{\frac13}+3k^{\frac16}\frac{c_1}{\sqrt{w}}+\frac{c_1^2}{w}}{c_2^2},
\label{corr_approx}
\end{equation}
where $c_1 = \mathbb{E}[W_1]=-1.21$ and $c_2=\sqrt{{\rm Var}(W_1)}= 1.27$.
Then when $k/w\to\gamma\ll1$ and $\delta \ll w$,
\begin{equation}
{\rm corr}(Z_t,Z_{t+\delta}) \leq 1 - \beta\delta + o(\delta).
\label{corr}
\end{equation}
\end{lemma}

\begin{proof}
The proof is given in the Appendix.
\end{proof}

By leveraging the properties of the local approximation in \eqref{corr}, we can obtain an asymptotic approximation using the localization theorem \cite{siegmund2010tail,SiegmundYakirZhang2010}.  Define a special function $v(\cdot)$ which is closely related to the Laplace transform of the overshoot over the boundary of a random walk \cite{siegmund2007statistics}:
\[
v(x) \approx \frac{\frac{2}{x}[\Phi\big(\frac{x}{2}\big)-0.5]}{\frac{x}{2}\Phi\big(\frac{x}{2}\big)+\phi\big(\frac{x}{2}\big)},
\]
where $\phi(x)$ and $\Phi(x)$ are the pdf and cdf of the standard normal distribution $\mathcal{N}(0,1)$.
\begin{proposition}[ARL with temporal correlation]
For large values of $b$ we can write
\begin{equation}
\mathbb{E}_\infty[T_{\rm E}] = \left[\beta b' \phi(b') v\big(b'\sqrt{2\beta}\big)\right]^{-1}\big(1+o(1)\big),
\label{eq:ARL}
\end{equation}
where \[b' = \frac{b - [\sigma_{w,k}\mathbb{E}(W_1) + \mu_{w,k}]}{\sigma_{w,k}\sqrt{\rm{Var}(W_1)}}.\]
\label{ARL2}
\end{proposition}
\begin{proof}
The proof is detailed in the Appendix.
\end{proof}

We perform simulations to verify the accuracy of the threshold value obtained without and with considering the temporal correlation (Proposition\,\ref{ARL1} and Proposition\,\ref{ARL2}, respectively). The results are shown in Table\,\ref{eq:threshold}. We find that, indeed, the threshold, when temporal correlation \eqref{eq:ARL} is taken into account, is more accurate than its counterpart obtained by using the Tracy-Widom law \eqref{eq:threshold1}.
\begin{table}[!h]
\begin{center}
\caption{Comparison of the threshold $b$ obtained from simulations and using the approximations. Window length $w=200$, dimension $k = 10$. The numbers shown are $b/w$.}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
Target ARL & 5k  &  10k  &  20k  & 30k &  40k & 50k\\
\hline
Simulation   &    1.633  &  1.661 &   1.688   & 1.702   & 1.713  &  1.722   \\
\hline
Approx \eqref{eq:threshold1} & 1.738 & 1.763 & 1.787 & 1.800 &  1.809 & 1.816\\
\hline
Approx \eqref{eq:ARL} & 1.699  &  1.713  &  1.727   & 1.735 &   1.740  &  1.744 \\
\hline
\end{tabular}
\label{eq:threshold}
\end{center}
\end{table}

\subsection{Lower Bound on EDD using Marginal Power}
We now focus on the detection performance and present a tight lower bound for the EDD of the largest eigenvalue procedure.
\begin{proposition}
For large values of $b$ we have
\begin{equation}
\mathbb{E}_0[T_{\rm E}] \geq 2\frac{b^\prime + e^{-b^\prime}-1}{-\log(1+\rho)+\rho}\big(1+ o(1)\big).
\label{eq:EDD}
\end{equation}
where
\[
b' = \frac1{2\sigma^2(1+\rho)}\left[b\rho-(1+\rho)\sigma^2\log(1+\rho)\right].
\]
\label{EDD}
\end{proposition}
\begin{proof}
The proof is based on a known result for CUSUM \cite{siegmund2013sequential} and requires the derivation of the Kullback-Leibler divergence for our problem. Details are given in the Appendix.
\end{proof}
Consistent with intuition, in Proposition \ref{EDD}, the right-hand-side of (\ref{eq:EDD}) is a decreasing function of the SNR $\rho$,

\begin{figure}[H]
\centerline{\includegraphics[width = 0.4\textwidth]{fig2.pdf}}
\caption{Simulated EDD and lower bound as a function of the threshold $b$.}
\label{fig:EDD}
\end{figure}

Comparing the lower bound in Proposition \ref{EDD} with simulated average delay, as shown in Fig. \ref{fig:EDD}, we can show that in the regime of small detection delay (which is the main regime of interest), the lower bound serves as a reasonably good approximation.

\section{Numerical examples} \label{sec:numerical}
In this section, numerical results are presented to compare the three detection procedures. The tests are first
applied to synthetic data and the performance of the Subspace-CUSUM and largest eigenvalue test are compared against the CUSUM optimum performance. Then the performance of Subspace-CUSUM is optimized by selecting the most appropriate window size.

\subsection{Performance Comparison}
We perform simulations to compare the largest eigenvalue procedure, the Subspace-CUSUM procedure, and the exact CUSUM procedure. The threshold for each procedure is determined by Monte-Carlo simulation, as discussed in Section \ref{sec:metric}. Fig.\,\ref{fig:eddcompare} depicts EDD versus log-ARL for parameter values $k=5$, $\theta=1$, $\sigma^2 = 1$ and window length $w = 20$. The black line corresponds to the exact CUSUM procedure, which is clearly the best and it lies below the other curves. Subspace-CUSUM has always smaller EDD than the largest eigenvalue procedure and the difference increases with increasing ARL.

\begin{figure}[h!]
\centerline{\includegraphics[width = 0.4\textwidth]{fig3.pdf}}
\caption{Comparison of CUSUM, Subspace-CUSUM and the Largest eigenvalue procedure.}
\label{fig:eddcompare}
\end{figure}

\subsection{Optimal window size.}
\begin{figure}[h!]
\centerline{\includegraphics[width = 0.4\textwidth]{fig4}}
\caption{Minimal EDD vs ARL among window sizes from 1 to 50.}
\label{fig:minimalEDD}
\end{figure}
\begin{figure}[h!]
\centerline{\includegraphics[width = 0.4\textwidth]{fig5}}
\caption{Optimal window size resulting in smallest EDD.}
\label{fig:optimalw}
\end{figure}
We also consider the EDD/ARL curve where $w$ is optimized to minimize the detection delay at every ARL.
We first compute the EDD for window sizes $w=1,2,\ldots,50$ given each ARL value. Then we plot in Fig.\,\ref{fig:minimalEDD}  the lower envelope of EDDs corresponding to the optimal EDD achieved by varying $w$. We also plot the optimal value of $w$ as a function of ARL in Fig.\,\ref{fig:optimalw}. Even though the best EDD of the Subspace-CUSUM is diverging from the performance enjoyed by CUSUM this divergence we believe is slower than the increase of the optimum CUSUM EDD. One of the goals in the future publication regarding the analysis of Subspace-CUSUM is to show that this is indeed the case, which in turn will demonstrate that this detection structure is first-order asymptotically optimum.

\section{Conclusion}\label{sec:conclusion}

In this paper, we considered three detection procedures for the rank-one change in the covariance matrix: the largest eigenvalue procedure, the exact CUSUM procedure, and the Subspace-CUSUM procedure. For Subspace-CUSUM we perform a simultaneous estimate of the required subspace in parallel with its sequential detection. We avoid estimating all unknown parameters by following a worst-case analysis with respect to the subspace power. We were able to derive theoretical expressions for the ARL and an interesting lower bound for the EDD of the largest eigenvalue procedure. In particular we were able to handle the correlations resulting from the usage of a sliding window which is an issue that is not present in the off-line version of the same procedure. For the comparisons of the three competing detectors we discuss how it is necessary to calibrate each detector so that comparisons are fair. Comparisons were performed using simulated data and Subspace-CUSUM was found to exhibit a significantly better performance than the largest eigenvalue procedure. Ongoing work involves establishing first-order asymptotic optimality of the Subspace-CUSUM procedure by determining the optimal drift parameter $d$ and by relating the sliding window length $w$ to the desired ARL.

\section*{Acknowledgment}
The work of Liyan Xie and Yao Xie was supported by the US National Science Foundation under Grants CCF\,1442635, CMMI\,1538746, DMS\,1830210, and the Career Award CCF\,1650913. The work of George Moustakides was supported by the US National Science Foundation under Grant CIF\,1513373, through Rutgers University.

\bibliographystyle{IEEEtran}

\subsubsection{Exponential tilting}Let the log moment generating function of $Z'_t$ be
\begin{equation*}
\psi(a) = \log\mathbb{E}_\infty[e^{a Z'_t}].
\end{equation*}
Define a family of new measures
\begin{equation*}
\frac{d\mathbb{P}_t}{d\mathbb{P}_\infty} = \exp\{a Z'_t - \psi(a)\},
\end{equation*}
where $\mathbb{P}_t$ denotes the new measure after the transformation.
The new measure takes the form of the exponential family, and $a$ can be viewed as the natural parameter. It can be verified that $\mathbb{P}_t$ is indeed a probability measure since
\begin{equation*}
\int d\mathbb{P}_t = \int \exp\{a Z'_t-\psi(a)\}d\mathbb{P} = 1.
\end{equation*}
It can also be shown that $\dot{\psi}(a)$ is the expected value of $Z'_t$ under $\mathbb{P}_t$, indeed
$$
\dot{\psi}(a)=\frac{\mathbb{E}_\infty[Z'_te^{a Z'_t}]}{\mathbb{E}_\infty[e^{a Z'_t}]} = \mathbb{E}_\infty[Z'_te^{a Z'_t-\psi(a)}] = \mathbb{E}_t[Z'_t],
$$
and similarly that $\ddot{\psi}(a)$ is equal to the variance under the tilted measure.

Under the assumption that $Z'_t$ is a standard Gaussian random variable, its log moment generating function satisfies $\psi(a) = a^2/2$.
Since we have the freedom to select $a$, we can set its value so that the mean under the tilted measure is equal to a given threshold $b$.
This can be done by choosing $a$ such that $\dot{\psi}(a)= b$ which yields $a = b$ a value independent from $t$. Given this choice, the transformed measure is given by $d\mathbb{P}_t = \exp(bZ'_t - b^2/2)d\mathbb{P}_\infty$. We also define, for each $t$, the log-likelihood ratio $\log(d\mathbb{P}_t/d\mathbb{P}_\infty)$ of the form
\begin{equation*}
\ell_t = b Z'_t - \frac{1}{2}b^2.
\end{equation*}

\subsubsection{Change-of-measure by the likelihood ratio identity} Now we convert the original problem of finding the small probability that the maximum of a random field exceeds a large threshold, to another problem: finding an alternative measure under which the event happens with a much higher probability.
Here, the alternative measure will be a mixture of simple exponential tilted measures. Define the maximum and the sum for likelihood ratio differences relative for a particular $t$:
\begin{equation*}
M_t = \max_{m \in \{1, \ldots, M\}} e^{\ell_m -\ell_t}, \quad S_t = \sum_{m \in \{1, \ldots, M\}} e^{\ell_m -\ell_t}.
\end{equation*}
Also define a re-centered likelihood ratio, which we call the global term
\begin{equation*}
\tilde{\ell}_t = b(Z'_t -b).
\end{equation*}
Then we have the following likelihood ratio identity:
\begin{align}
& \mathbb{P}_\infty(\max\limits_{1\le m \le M} Z'_m \geq b) = \mathbb{E}_\infty [\mathbbm{1}_{\{ \max\limits_{1\le m \le M} Z'_m \geq b\}}]\nonumber\\
&= \mathbb{E}_\infty \left[\frac{\sum_{t=1}^{M}e^{\ell_t}}{\sum_{n=1}^{M}e^{\ell_n}} \cdot \mathbbm{1}_{\{\max\limits_{1\le m \le M} Z'_m \geq b\}}\right]\nonumber\\
&=  \sum_{t=1}^M \mathbb{E}_\infty \left[\frac{e^{\ell_t}}{\sum_n e^{\ell_n}} \cdot \mathbbm{1}_{\{ \max\limits_{1\le m \le M} Z'_m \geq b\}} \right]\label{eq:tail_prob}\\
&=   \sum_{t=1}^M \mathbb{E}_t \left[\frac{1}{\sum_n e^{\ell_n}} \cdot \mathbbm{1}_{\{ \max\limits_{1\le m \le M} Z'_m \geq b\} }\right]\nonumber\\
& =  e^{-b^2/2} \sum_{t=1}^M \mathbb{E}_t \left[\frac{M_t}{S_t}e^{-(\tilde{\ell}_t +\log M_t)} \cdot \mathbbm{1}_{\{ \tilde{\ell}_t +\log M_t \geq 0 \} }\right].\nonumber
\end{align}
The last equation in (\ref{eq:tail_prob}) converts the tail probability to a product of two terms: a deterministic term $e^{-b^2/2}$ associated with the large deviation rate, and a sum of conditional expectations under the transformed measures. The conditional expectation involves a product of the ratio $M_t/S_t$, and an exponential function that depends on $\tilde{\ell}_t$, which plays the role of a weight. Under the new measure $\mathbb{P}_t$, $\tilde{\ell}_t$ has zero mean and variance equal to $b^2$ and it dominates the other term $\log M_t$, hence, the probability of exceeding zero is much higher. Next, we characterize the limiting ratio and the other factors precisely, by the localization theorem.

\subsubsection{Establish properties of local field and global term} In (\ref{eq:tail_prob}), our target probability has been decomposed into terms that only depend on (i) the local field $\{\ell_m - \ell_t\}, 1 \le m \le M$, which are the differences between the log-likelihood ratios with parameter $t$ and $m$, and (ii) the global term $\tilde{\ell}_t$, which is the centered and scaled likelihood ratio with parameter $t$. We need to first establish some useful properties of the local field and global term before applying the localization theorem. We will eventually show that the local field and the global term are asymptotically independent.

It is easy to show that under the tilted measure,
\begin{equation*}
\mathbb{E}_t[\tilde{\ell}_t] = 0, \quad \text{Var}_t[\tilde{\ell}_t] = b^2,
\end{equation*}
i.e., the global term $\tilde{\ell}_t$ is zero mean for any $t$, with variance diverging with $b$.

For the local field $\{\ell_m - \ell_t\}$, let $r_{m,t}$ denote the correlation between $Z'_m$ and $Z_t^\prime$ (given by Lemma\,\ref{CovZ}), then we have
\[
\begin{aligned}
\mathbb{E}_t(\ell_m - \ell_t) &= -b^2(1-r_{m,t}),\\
\mbox{Var}_t(\ell_m - \ell_t) &= 2b^2(1-r_{m,t}),\\
\mbox{Cov}_t(\ell_{m_1} - \ell_t,\ell_{m_2} - \ell_t) &= b^2(1+r_{m_1,m_2} - r_{m_1,t} - r_{m_2,t}).
\end{aligned}
\]

Since we assume $Z'_t$ is approximately Gaussian, the local field $\ell_m - \ell_t$ (or equivalently $b(Z'_m - Z'_t)$ ) and the global term $\tilde{\ell}_t$ (or equivalently $b(Z'_t - b)$) are also approximately Gaussian.
Therefore, when $|\delta|$ is small (i.e., in the neighborhood of zero), we can approximate the local field using a two-sided Gaussian random walk with drift $\mu^2/2$ and variance of the increment equal to $\mu^2$:
\begin{equation}
\ell_{t+\delta} - \ell_t \stackrel{\Delta}{=} \mu \sum\limits_{i=1}^{|\delta|} \vartheta_i - \frac{1}{2}\mu^2\delta, \delta = \pm1,\pm2,\ldots,
\label{eq:local_decomp}
\end{equation}
where $\vartheta_i$ are i.i.d. standard normal random variables and $\mu = b\sqrt{2\beta}$.

We have Lemma\,\ref{CovZ} to characterize the local correlation, which offers reasonably good approximation for ${\rm corr}(Z_t,Z_{t+\delta})$ and leads to some insights for our analysis.

\subsubsection{Approximation using localization theorem} Then we use the localization theorem (Theorem 5.1 in \cite{siegmund2010tail}) given below.
\begin{lemma}[Localization theorem]
Under mild assumptions, we have
\begin{equation*}
\begin{aligned}
&\lim_{\kappa\rightarrow\infty}\kappa^{1/2}\mathbb{E} \left[ \frac{M_\kappa}{S_\kappa}e^{-[\tilde{\ell}_\kappa+\log M_\kappa]}\cdot \mathbbm{1}_{\{\tilde{\ell}_\kappa+\log M_\kappa\ge 0\}} \right] \\
= & \sigma^{-1}\phi\left(\frac\mu\sigma \right)\mathbb{E} \left[\frac{M}{S}\right],
\end{aligned}
\end{equation*}
where $\phi(\cdot)$ is the density of the standard normal distribution.
\end{lemma}

Intuitively, the localization theorem says the following. To find the desired limit of $\mathbb{E}[\frac{M_\kappa}{S_\kappa}e^{-[\tilde{\ell}_\kappa+\log M_\kappa]}\mathbbm{1}_{\{ \tilde{\ell}_\kappa+\log M_\kappa\ge 0 \}} ]$ as $\kappa \rightarrow \infty$, one first approximates $M_\kappa$ and $S_\kappa$ by their localized versions, which are obtained by restricting the maximization and summation in a neighborhood of parameter values. Then one can show that the localized ratio $M_\kappa/S_\kappa$ is asymptotically independent of the global term $\tilde{\ell}_\kappa$ as $\kappa\rightarrow\infty$. The asymptotic analysis is then performed on the local field and the global term separately. The expected value of the localized ratio $\mathbb{E}[M_\kappa/S_\kappa]$ converges to a constant independent of $\kappa$, and the limiting conditional distribution of $\tilde{\ell}_\kappa$ can be found using the local central limit theorem. Thus, one can calculate the remaining conditional expectation involving $\tilde{\ell}_\kappa$.

After verifing the validity of the conditions in our setting, we have
\begin{align*}
&\mathbb{P}_\infty(T\le M)= \mathbb{P}_\infty\left( \max_{1\le t\le M} Z'_t \geq b \right) \\
&= e^{-b^2/2}\sum\limits_{t=1}^M\mathbb{E}_t \left[ \frac {M_t}{S_t}e^{-[\tilde{\ell}_t+\log M_t]}\cdot \mathbbm{1}_{\{ \tilde{\ell}_t+\log M_t\ge0\}} \right]\\
&= Mb\phi(b)\beta v(b\sqrt{2\beta})(1+o(1)).
\end{align*}
Once the cumulative distribution function of $T$ is available, we can approximate it as an exponential distribution with parameter $b\phi(b)\beta v(b\sqrt{2\beta})$ yielding as mean the value $1/[b\phi(b)\beta v(b\sqrt{2\beta})]$. Since $Z'_t$ is standardized, here the threshold $b$ need to be converted to the original threshold using a simple formula \[b' = \frac{b - [\sigma_{w,k}\mathbb{E}(W_1) + \mu_{w,k}]}{\sigma_{w,k}\sqrt{\rm{Var}(W_1)}}.\] This completes the proof.
\end{proof}

\begin{proof}[Proof of Proposition \ref{EDD}]
We first relate the largest eigenvalue procedure to a CUSUM procedure, note that
\begin{equation}
\lambda_{\max}(\hat{\Sigma}_{t,w}) = \max_{\|q\|=1}
q^\intercal \hat{\Sigma}_{t,w} q.
\label{eq:q}
\end{equation}
For each $q$, we have
\[
q^\intercal \hat{\Sigma}_{t,w} q =  \sum_{i=t-w+1}^t (q^\intercal x_i)^2.
\]
According to the Grothendieck's Inequality \cite{guedon2016community}, the $q$ that attains the maximum in equation (\ref{eq:q}) is very close to $u$ under the alternative. Therefore, assuming the optimal $q$ always equals to $u$ will only cause a small error but will bring great convenience to our analysis.

Now we have under $\mathbb{P}_\infty$, $q^\intercal x_i \sim \mathcal{N}(0, \sigma^2)$ and under $\mathbb{P}_0$, $q^\intercal x_i \sim \mathcal{N}(0, \sigma^2+\theta)$. Let $f_\infty$ denote the pdf of $\mathcal{N}(0, \sigma^2)$ and $f_0$ the pdf of $\mathcal{N}(0, \sigma^2+\theta)$. For each observation $y$, we can derive the one-sample log-likelihood ratio:
\[
\log\frac{f_0(y)}{f_\infty(y)} = -\frac12\log(1+\rho)+\frac1{2\sigma^2}\left(1-\frac1{1+\rho}\right)y^2.
\]
Define the CUSUM procedure
\begin{multline*}
\widetilde{T}= \inf\bigg\{ t : \max_{0\leq k<t}\sum_{i=k+1}^t \Big[\frac1{2\sigma^2}\left(1-\frac1{1+\rho}\right)(q^\intercal x_i)^2\\
-\frac{\log(1\!+\!\rho)}{2}\Big] \! \geq \! b'\bigg\},
\end{multline*}
where $b'=\frac1{2\sigma^2(1+\rho)}\left[b\rho-(1+\rho)\sigma^2\log(1+\rho)\right]$, we then have
\[
\mathbb{E}_0[T_{\rm E}] \geq \mathbb{E}_0[\widetilde{T}].
\]
Since $\widetilde{T}$ is a CUSUM procedure with
\[
\int\log\left[\frac{f_0(y)}{f_\infty(y)}\right]f_0(y)dy=-\frac12\log(1+\rho)+\frac{\rho}{2},
\]
from \cite{siegmund2013sequential} we have:
\[
\mathbb{E}_0[\widetilde{T}]=\frac{e^{-b^\prime}+b^\prime-1}{-\log(1+\rho)/2+\rho/2}.
\]
This completes the proof.
\end{proof}

\section{Introduction} \label{sec:intro}

\IEEEPARstart{D}{etecting} the changepoint from high-dimensional streaming data is a fundamental problem in various applications such as video surveillance, sensor networks, and seismic events detection. In various scenarios, the change happens to the covariance structure and can be represented as a linear subspace. For example, the covariance matrix may shift from an identity matrix to a spiked covariance model \cite{berthet2013optimal,johnstone2001distribution}.

Given a sequence of observed vectors $x_1, x_2, \ldots,$ where $x_t \in \mathbb{R}^k$ and $k$ is the signal dimension, there may be a changepoint time $\tau$ when the distribution of the data stream  changes. Our goal is to detect this change as quickly as possible from streaming (sequentially obtained) data using online techniques.
We are particularly interested in the structured change occurring in the signal covariance. We study two related settings: the {\it emerging subspace}, meaning that the change is a subspace emerging from a noisy background; and the {\it switching subspace}, meaning that the change is a switch in the direction of the subspace. The emerging subspace problem can arise, for instance, from coherent weak signal detection from seismic sensor arrays, and the switching subspace detection can be used for principal component analysis for streaming data. In these settings, the changes can be shown to be equivalent to a low-rank component added to the original covariance matrix. On the other hand, the switching subspace problem, as we will see, can be reduced to the emerging subspace problem, if we are willing to tolerate a performance loss. Therefore, we will focus on the analysis of the emerging subspace problem.

In this paper, we consider three detection procedures. We start with the exact CUSUM which is known to be optimum when we have complete knowledge of the pre- and post-change statistics and parameters. Since the post-change parameters are usually unknown we propose two alternatives to deal with the case of unknown parameters. Specifically we consider the largest eigenvalue procedure where we use as test statistic the largest eigenvalue of the sample covariance matrix of the data contained within a sliding time-window. This can be regarded as a straightforward extension of its offline counterpart \cite{berthet2013optimal}. The second method which we call the Subspace-CUSUM uses the structure of the exact CUSUM but in place of the parameters that are known it uses their estimates which are computed using data within, again, a sliding window. We perform a theoretical analysis of the largest eigenvalue procedure. A similar analysis, which is far more complicated and extended, is postponed for a future article. The three algorithms are compared using simulations

The rest of the paper is organized as follows. Section\,\ref{sec:formulation} details on the two problems of emerging and switching subspace and how they can be related. Section\,\ref{sec:detect} presents the three sequential change detection procedures. In Section\,\ref{sec:theo} we develop theoretical bounds for the average run length and the expected detection delay of the largest eigenvalue procedure. In Section\,\ref{sec:numerical} we presents numerical results and comparisons for the competing algorithms. Finally Section\,\ref{sec:conclusion} contains our concluding remarks.

\subsection{Related Work}

Classical approaches to covariance change detection usually consider generic settings without assuming any structure.
The CUSUM statistics can be derived if the pre-change and post-change distributions are known.
For the multivariate case, the Hotelling $T^2$ control chart is the traditional way to detect the covariance change. The determinant of the sample covariance matrix was also used in \cite{alt2004multivariate} to detect change of the determinant of the covariance matrix. A multivariate CUSUM based on likelihood functions of multivariate Gaussian is studied in \cite{healy1987note}  but it only considers the covariance change from $\Sigma$ to $c\Sigma$ for a constant $c$.
Offline change detection of covariance change from $\Sigma_1$ to $\Sigma_2$ is studied in \cite{chen2004statistical} using the Schwarz information criterion \cite{schwarz1978estimating}, where the changepoint location must satisfy certain regularity condition to ensure the existence of the maximum likelihood estimator. In \cite{arias2012detection} we find a hypothesis testing approach to detect a shift in an off-diagonal sub-matrix of the covariance matrix using likelihood ratios. Recently, \cite{jiao2018subspace} studies a CUSUM-like procedure for detection of switching subspaces, when the distributions (as well as the subspaces) before and after the changepoint are exactly known; this is different from our work since we assume the subspace after the change is unknown.

The most related work to our present effort is the hypothesis testing methods developed in \cite{berthet2013optimal}, which uses the largest eigenvalue of the sample covariance matrix to detect a sparse spiked covariance model given a fixed number of samples. The largest eigenvalue statistic is shown to be asymptotically minimax optimal for determining whether there exists a sparse and low-rank component in the offline setting. A natural sequential version of this idea is to use a sliding window and estimate the largest eigenvalue of the corresponding sample covariance matrix. However, this approach, under a sequential setting does not enjoy any form of (asymptotic) optimality.

A different test statistic, the so-called Kac-Rice statistic \cite{TaylorRice}, has been considered for testing spiked covariance model. The Kac-Rice statistic is the conditional survival function of the largest observed singular value conditioned on all other observed singular values, and is characterized by a simple asymptotic distribution (uniform in $[0, 1]$). However, the statistic involves the computation of an integral over the whole real line, and it is not clear how this can be carried over to the sequential formulation.

\section{Problem Setup} \label{sec:formulation}
We first introduce the spiked covariance model \cite{johnstone2001distribution}, which assumes that a small number of directions explain most of the variance. In particular, we consider the rank-one spiked covariance matrix, which is given by
\[
\Sigma = \sigma^2I_k + \theta uu^\intercal,
\]
where $I_k$ denotes an identity matrix of size $k$;  $\theta$ is the signal strength; $u \in \mathbb{R}^k$ represents a basis for the subspace with unit norm $\left\Vert u \right\Vert=1$; $\sigma^2$ is the noise variance, which will be considered known since it can be estimated from training data\footnote{In fact it is possible to consider $\sigma^2$ unknown as well and provide estimates of this parameter along with the necessary estimates of $u$. However, in order to simplify our presentation, we decided to consider $\sigma^2$ known.}. The Signal-to-Noise Ratio (SNR) is defined as
$\rho = \theta/\sigma^2$.

Formally, the \textit{emerging} subspace problem can be cast as follows:
\begin{equation}
\begin{array}{ll}
x_t  \stackrel{\text{iid}}{\sim} \mathcal{N}(0, \sigma^2I_k), &t = 1,2,\ldots,\tau, \\
x_t \stackrel{\text{iid}}{\sim} \mathcal{N}(0,  \sigma^2I_k+ \theta uu^\intercal), &t = \tau+1,\tau+2,\ldots
\end{array}
\label{eq:hypothesis}
\end{equation}
where $\tau$ is the unknown changepoint that we would like to detect from data that are acquired sequentially.

Similarly, the \textit{switching} subspace problem can be formulated as follows:
\begin{equation}
\begin{array}{ll}
x_t  \stackrel{\text{iid}}{\sim} \mathcal{N}(0, \sigma^2I_k + \theta u_1u_1^\intercal), &\!\!\!t = 1,2,\ldots,\tau, \\
x_t \stackrel{\text{iid}}{\sim} \mathcal{N}(0,  \sigma^2I_k  + \theta u_2 u_2^\intercal),&\!\!\!t = \tau+1,\tau+2,\ldots
\end{array}
\label{eq:hypothesis_switch}
\end{equation}
where $u_1, u_2 \in \mathbb R^{k}$ represent bases for the subspaces before and after the change, with $\|u_1\| = \|u_2\| = 1$ and $u_1$ is considered known. In both settings, our goal is to detect the change as quickly as possible.

The switching subspace problem \eqref{eq:hypothesis_switch} can be reduced into the emerging subspace problem \eqref{eq:hypothesis} by a simple data projection. Specifically, select any orthonormal matrix $Q \in \mathbb{R}^{(k-1)\times k}$ such that
\[ Q u_1 = 0, \quad QQ^\intercal = I_{k-1},\]
which means that all rows of $Q$ are orthogonal to $u_1$, they are orthogonal to each other and have unit norm. Then, using the matrix $Q$, we project each observation $x_t$ onto a $k-1$ dimensional space and obtain a new sequence
\[ y_t = Qx_t \in \mathbb{R}^{k-1}, t = 1,2,\ldots . \]
Then $y_t$ is a zero-mean random vector with covariance matrix $\sigma^2I_{k-1}$ before the change and $\sigma^2I_{k-1} + \theta Qu_2u_2^\intercal Q^\intercal$ after the change.
Let $u = Qu_2/\left\Vert Qu_2 \right\Vert$, and
\[ \tilde{\theta} = \theta \left\Vert Qu_2 \right\Vert^2 = \theta [ 1- (u_1^\intercal u_2)^2]. \]
Thus, problem \eqref{eq:hypothesis_switch} can be reduced to the following
\begin{equation}
\begin{array}{ll}
y_t  \stackrel{\text{iid}}{\sim} \mathcal{N}(0, \sigma^2I_{k-1}), &\!\!\!t = 1,2,\ldots,\tau, \\
y_t \stackrel{\text{iid}}{\sim} \mathcal{N}(0,  \sigma^2I_{k-1}  +\tilde{ \theta} u u^\intercal),&\!\!\!t = \tau+1,\tau+2,\ldots
\end{array}
\label{eq:hypothesis_switch_transform}
\end{equation}
Note that this way the switching subspace problem is reduced into the emerging subspace problem, where the new signal power $\tilde \theta$ depends on the angle between $u_1$ and $u_2$, which is consistent with our intuition.

We would like to emphasize that by projecting the observations onto a lower dimensional space we lose information, suggesting that the two versions of the problem are not equivalent. Indeed, the optimum detector for the transformed data in \eqref{eq:hypothesis_switch_transform} and the one for the original data in \eqref{eq:hypothesis_switch} do not coincide. This can be easily verified by computing the corresponding CUSUM tests and their optimum performance. Despite this difference, it is clear that with the proposed approach we put both problems under the same framework, offering, as we will see, computationally simple methods to solve the original problem in \eqref{eq:hypothesis_switch}. Consequently, in the following analysis, we focus solely on problem \eqref{eq:hypothesis}.

\section{Detection Procedures} \label{sec:detect}

As we mentioned before, we consider three methods: The exact CUSUM procedure where all parameters are considered known, the largest eigenvalue procedure, and the Subspace-CUSUM procedure. It is clear that since CUSUM is optimum it will be regarded as a point of reference for the other two approaches. We first introduce some necessary notation.
Denote with $\mathbb{P}_\tau$ and $\mathbb{E}_\tau$ the probability and expectation induced when there is a changepoint at the \textit{deterministic} time $\tau$. Under this definition $\mathbb{P}_\infty$ and $\mathbb{E}_\infty$ is the probability and the expectation under the nominal regime (change never happens) while $\mathbb{P}_0$ and $\mathbb{E}_0$ the probability and expectation under the alternative regime (change happens before we take any data).

\subsection{Optimal CUSUM Procedure}\label{detect:exactcusum}

The CUSUM procedure \cite{page1954continuous,siegmund2013sequential} is the most popular sequential test for change detection. When the observations are i.i.d. before and after the change, CUSUM is known to be exactly optimum \cite{moustakides1986optimal} in the sense that it solves a very well defined constrained optimization problem introduced in \cite{lorden1971procedures}. However, the CUSUM procedure can be applied only when we have exact knowledge of the pre- and post-change distributions. Thus, for our problem, it requires complete specification of all parameters namely the subspace $u$, noise power $\sigma^2$ and SNR $\rho$.

To derive the CUSUM procedure, let $f_\infty(\cdot),f_0(\cdot)$ denote the pre- and post-change probability density function (pdf) of the observations. Then the CUSUM statistics is defined by maximizing the log-likelihood ratio statistic over all possible changepoint locations
\begin{equation}
S_ t =  \max_{1 \leq j  \leq t } \sum_{i=j}^t \log \frac{f_0(x_i)}{f_{\infty}(x_i)},
\label{cusum_exact}
\end{equation}
which has the recursive implementation
\begin{equation}
S_t = (S_{t-1})^+ + \log \frac{f_0(x_t)}{f_{\infty}(x_t)},
\label{cusum_recursive}
\end{equation}
that enables its efficient calculation \cite{moustakides1986optimal}. The CUSUM stopping time in turn is defined as
\begin{equation}
T_{\rm C} = \inf\{t>0: S_t \geq b\},
\label{cusum_procedure}
\end{equation}
where $b$ is a threshold selected to meet a suitable false alarm constraint.
For our problem of interest we can derive that
\begin{multline}
\log \frac{f_0(x_t)}{f_\infty(x_t)}  = \log\bigg[ \frac{[(2\pi)^p|\sigma^2I_k + \theta u u^\intercal|]^{-1/2}}{[(2\pi)^p\sigma^{2p}]^{-1/2}} \\
 \times
 \frac{\exp\{ - x_t^\intercal(\sigma^2I_k+\theta u u^\intercal)^{-1}x_t / 2  \} }{\exp\{- x_t^\intercal x_t /(2\sigma^2) \} } \bigg] \\
= \log\left[
\left|I_k + \rho u u^\intercal \right|^{-\frac 1 2}
\exp\left\{ \frac 1 2 \frac{\theta}{\theta+\sigma^2} \frac{(u^\intercal x_t)^2}{\sigma^2}  \right\} \right]  \\
 =\frac{\rho}{2\sigma^2(1+\rho)} \Big\{(u^\intercal x_t)^2- \sigma^2 \left(1+ \frac1{\rho}\right)\log(1+\rho)\Big\}.
\label{expansion}
\end{multline}
The second equality is due to the matrix inversion lemma \cite{woodbury1950inverting} that allows us to write
\[
(\sigma^2I_k+\theta u u^\intercal)^{-1} =
 \frac1{\sigma^2}I_k - \frac{\theta}{\theta+\sigma^2}\frac{uu^\intercal}{\sigma^2},
\]
which, after substitution into the equation, yields the desired result. Note that the multiplicative factor $\rho/[2\sigma^2(1+\rho)]$ is positive, so we can omit it from the log-likelihood ratio when forming the CUSUM statistic. This leads to
\begin{equation}
S_t = (S_{t-1})^+ + (u^\intercal x_t)^2 - \sigma^2\left(1+\frac{1}{\rho}\right)\log(1+\rho).
\label{cusum_recur}
\end{equation}

\begin{remark}\label{drift_property_1}
We can show that the increment in \eqref{cusum_recur}, i.e.,
\[(u^\intercal x_t)^2 - \sigma^2\left(1+\frac{1}{\rho}\right)\log(1+\rho),
\] has the following property: its expected value is negative under the pre-change and positive under the post-change probability measure. The proof relies on a simple argument based on Jensen's inequality. Due to this property, before the change, the CUSUM statistics $S_t$ will oscillate near $0$ while it will exhibit, on average, a positive linear drift after the occurrence of the change forcing it, eventually, to hit or exceed the threshold.
\end{remark}

\subsection{Largest Eigenvalue Procedure} \label{detect:eigen}

Motivated by the off-line test in \cite{berthet2013optimal}, a natural strategy to detect the change is to use the largest eigenvalue of the sample covariance matrix. Under the sequential setting, we adopt a sliding window approach and form the sample covariance matrix using observations that lie within a time window of length $w$. For each time $t>0$, the \textit{un-normalized} sample covariance matrix using the available samples is given by
\begin{equation}
\hat{\Sigma}_{t,\min\{t,w\}}=\left\{
\begin{array}{l}
x_1 x_1^\intercal+\cdots+x_t x_t^\intercal,~\text{for}~t<w\\
x_{t-w+1} x_{t-w+1}^\intercal+\cdots+x_t x_t^\intercal,~\text{for}~t\geq w.
\end{array}\right.
\label{sampleCov}
\end{equation}
We note that for $t=1$ the matrix contains a single outer product and as time progresses the number of outer products increases linearly until it reaches $w$. After this point, namely for $t\geq w$, the number of outer products remains equal to $w$.

Let $\lambda_{\max}(X)$ denote the largest eigenvalue of a symmetric matrix $X$. We define the {\it largest eigenvalue procedure}, as the one that stops according to the following rule:
\begin{equation}
\label{eigenvalue_procedure}
T_{\rm E} = \inf\{t>0: \lambda_{\max}(\hat{\Sigma}_{t,\min\{t,w\}}) \geq b \},
\end{equation}
where $b > 0$ is a constant threshold selected to meet a suitable false alarm constraint.
We need to emphasize that we do not divide by $\min\{t,w\}$ when forming the un-normalized sample covariance matrix. As we explain in Section\,\ref{sec:metric}, it is better for $T_{\rm E}$ to always divide by $w$ instead of $\min\{t,w\}$. Consequently, we can omit the normalization with $w$ from our detection statistics by absorbing it into the threshold.

\subsection{Subspace-CUSUM Procedure}\label{sec:S-CUSUM}

Usually the subspace $u$ and SNR $\rho$ are unknown. In this case it is impossible to form the exact CUSUM statistic depicted in \eqref{cusum_recur}. One option is to estimate the unknown parameters and substitute them back into the likelihood function. Here we propose to estimate only $u$ and call $d=\sigma^2(1+1/\rho)\log(1+\rho)$ which leads to the following Subspace-CUSUM update
\begin{equation}
\mathcal S_t = (\mathcal S_{t-1})^{+} + (\hat{u}_{t+w}^\intercal x_t)^2 - d,
\label{sscusum_update}
\end{equation}
We denote the estimator of $u$ as $\hat{u}_{t+w}$. This is because at time $t$ the estimate will rely on the data $x_{t+w},\ldots,x_{t+1}$ that are in the ``future'' of $t$. Practically, this is always possible by properly delaying our data by $w$ samples. Stopping occurs similarly to CUSUM, that is
\[
T_{\rm SC} = \inf\{t>0: \mathcal S_t \geq b\}.
\]
Of course, in order to be fair, at the time of stopping we must make the appropriate correction, namely if $\mathcal{S}_t$ exceeds the threshold at $t$ for the first time, then the actual stopping takes place at $t+w$. The reason we use estimates based on ``future'' data is to make $x_t$ and $\hat{u}_{t+w}$ independent which in turn will help us decide what is the appropriate choice for the drift constant $d$.

For the drift parameter $d$ we need the following double inequality to be true
\begin{equation}
\mathbb{E}_\infty[(\hat{u}_{t+w}^\intercal x_t)^2] < d < \mathbb{E}_0[(\hat{u}_{t+w}^\intercal x_t)^2].
\label{expected_d}
\end{equation}
With \eqref{expected_d} we can guarantee that $\mathcal{S}_t$ mimics the behavior of the exact CUSUM statistic $S_t$ mentioned in Remark\,\ref{drift_property_1}, namely, it exhibits a negative drift before and a positive after the change. To apply \eqref{sscusum_update}, we need to specify $d$ and of course provide the estimate $\hat u_{t+w}$. Regarding the latter we simply use the \textit{unit-norm} eigenvector corresponding to the largest eigenvalue of $\hat \Sigma_{t+w,w}$ depicted in \eqref{sampleCov}. As we mentioned, the main advantage of using $\hat \Sigma_{t+w,w}$ is that it provides estimates $\hat u_{t+w}$ that are \textit{independent} from $x_t$. This independence property allows for the straightforward computation of the two expectations in \eqref{expected_d} and contributes towards the proper selection of $d$. Note that under the pre-change distribution we can write
\begin{multline}
\mathbb{E}_\infty[(\hat{u}_{t+w}^\intercal x_t)^2] =  \mathbb{E}_\infty[ \hat{u}_{t+w}^\intercal\mathbb{E}_\infty[ x_t x_t^\intercal] \hat{u}_{t+w} ]\\
 = \sigma^2\mathbb{E}_\infty [\hat{u}_{t+w}^\intercal\hat{u}_{t+w}]= \sigma^2,
\label{drift_infty}
\end{multline}
where the first equation is due to the independence of $x_t$ and $\hat{u}_{t+w}$, the next one due to $x_t$ having covariance $\sigma^2I_k$ and the last equality due to $\hat{u}_{t+w}$ being of unit norm.

Under the post-change regime, we need to specify the statistical behavior of $\hat u_{t+w}$ for the computation of $\mathbb{E}_0[(\hat{u}_{t+w}^\intercal x_t)^2]$. We will assume that the window size $w$ is sufficiently large so that Central Limit Theorem (CLT) approximations \cite{anderson1963asymptotic, paul2007asymptotics} are possible for $\hat u_{t+w}$. The required result appears in the next lemma.
\begin{lemma
\label{lem:eigenvector}
Suppose vectors $x_1,\ldots,x_w$ are of dimension $k$ and follow the distribution  $\mathcal{N}(0,  \sigma^2I_k + \theta u u^\intercal)$. Let $\hat{\varphi}_w$ be the eigenvector corresponding to the largest eigenvalue of the sample covariance matrix $(1/w)(x_1x_1^\intercal+\cdots+x_wx_w^\intercal)$, then, as $w \rightarrow \infty$, we have the following CLT version for $\hat{\varphi}_w$
\[
\sqrt{w}(\hat{\varphi}_w - u)  \rightarrow \mathcal N\left(0, \frac{ 1+\rho}{\rho^2}  ( I_k - uu^\intercal)   \right).
\]
\end{lemma}

\begin{proof} The proof is detailed in the Appendix.
\end{proof}
Lemma\,\ref{lem:eigenvector} provides an asymptotic statistical description of the \textit{un-normalized} estimate of $u$. More precisely it characterizes the estimation error $v_w=\hat{\varphi}_w-u$. In our case we estimate the eigenvector from the matrix $\hat \Sigma_{t+w,w}$ but, as mentioned before, we adopt a \textit{normalized} (unit norm) version $\hat{u}_{t}$. Therefore if we fix $w$ at a sufficiently large value and $v_t$ denotes the estimation error of the un-normalized estimate at time $t$ then, from Lemma\,\ref{lem:eigenvector}, we can deduce
\begin{multline}\label{hat_u_separate}
 \hat{u}_{t+w} = \frac{\hat{\varphi}_{t+w}}{\|\hat{\varphi}_{t+w}\|}=\frac{u + v_{t+w}}{\left\Vert u + v_{t+w} \right\Vert},\\
   v_{t+w} \sim \mathcal{N}\left(0, \frac{1+\rho}{w\rho^2}(I_k-uu^\intercal)\right).
 \end{multline}
Note that $v_{t+w}$ is also independent from $x_t$ and orthogonal to $u$, the latter being true because the covariance matrix of $v_{t+w}$ is $I_k -uu^\intercal$. This implies $ \left\Vert u + v_{t+w}\right\Vert^2 = 1+\left\Vert v_{t+w}\right\Vert^2 $. Combining the above results, we have
\begin{multline*}
\mathbb{E}_0\left[(\hat{u}_{t+w}^\intercal x_t)^2 \right] =\\
\mathbb{E}_0\left[\frac{(u + v_{t+w})^\intercal (\sigma^2I_k + \theta uu^\intercal)(u + v_{t+w})} {1 + \left\Vert v_{t+w}\right\Vert^2} \right].
\end{multline*}
Because $v_{t+w}^\intercal u = 0, \|u\|=1$, the above expression simplifies to
\begin{multline}\label{expect_0}
\mathbb{E}_0\left[ \frac{\sigma^2 +\theta + \sigma^2 \left\Vert v_{t+w}\right\Vert^2 } { 1+\left\Vert v_{t+w}\right\Vert^2  } \right] = \sigma^2 +\theta\,\mathbb{E}_0\left[\frac{1} { 1+\left\Vert v_{t+w}\right\Vert^2 } \right]\\
=\sigma^2 +\theta\,\left\{1-\mathbb{E}_0[\|v_{t+w}\|^2]+\mathbb{E}_0\left[\frac{\|v_{t+w}\|^4}{1+\|v_{t+w}\|^2}\right]\right\}.
\end{multline}
For the two expectations in \eqref{expect_0}, using the Gaussian approximation from Lemma\,\ref{lem:eigenvector}, we have
\[
 \mathbb{E}_0 [\left\Vert v_{t+w}\right\Vert^2] = \frac{1+\rho}{w\rho^2}(k-1),
\]
and
$$
\mathbb{E}_0\left[\frac{\|v_{t+w}\|^4}{1+\|v_{t+w}\|^2}\right]\leq\mathbb{E}_0[\|v_{t+w}\|^4]=\frac{(1+\rho)^2}{w^2\rho^4}(k^2-1).
$$
Consequently
 \begin{equation}\label{lower_bound}
\mathbb{E}_0\left[(\hat{u}_{t+w}^\intercal x_t)^2\right]  =  \sigma^2(1+\rho)\left( 1 - \frac{k-1}{w\rho} +o\left(\frac{1}{w}\right)\right),
 \end{equation}
with the $o(\cdot)$ term being negligible compared to the other two when $\frac{k}{w}\ll1$.

Consider now the case where $\rho$ is unknown but exceeds some pre-set minimal SNR $\rho_{\min}$. From the above derivation, given the worst-case SNR and an estimation for the noise variance $\hat{\sigma}^2$, we can give a lower bound for $\mathbb{E}_0[(\hat{u}_{t+w}^\intercal x_t)^2]$. Consequently, the drift $d$ can be anything between
$\hat{\sigma}^2$
  and
  $\hat{\sigma}^2(1+\rho_{\min}) (1-  (k-1)/(w\rho_{\min}) )$ where, we observe, that the latter quantity exceeds $\hat{\sigma}^2$ when $w > (k-1)(1+\rho_{\min})/\rho_{\min}^2$.
 Below, for simplicity, for $d$ we use the average of the two bounds.

Alternatively, and in particular when $w$ does not satisfy $w\gg k$, we can estimate $\mathbb{E}_0[(\hat{u}_{t+w}^\intercal x_t)^2]$ by Monto Carlo simulation.
This method requires: (i) estimating the noise level $\hat{\sigma}^2$, which can be obtained from training data without a changepoint; (ii) the pre-set worst-case SNR $\rho_{\min}$; (iii) a unit norm vector $u_0$ that is generated randomly. Under the nominal regime we have $\mathbb{E}_\infty[(\hat{u}_{t+w}^\intercal x_t)^2]=\hat{\sigma}^2$. Under the alternative $\mathbb{E}_0[(\hat{u}_{t+w}^\intercal x_t)^2]$ depends only on the SNR $\rho$ as shown in \eqref{lower_bound}. We can therefore simulate the worst-case scenario $\rho_{\min}$ using the randomly generated vector $u_0$ by generating samples from the distribution $\mathcal{N}(0,  \hat{\sigma}^2I_k + \rho_{\min}u_0u_0^\intercal )$.

Even though the average of the update in \eqref{sscusum_update} does not depend on $u$, the computation of the test statistic $\mathcal{S}_t$ requires the estimate $\hat{u}_{t+w}$ of the eigenvector. This can be accomplished by applying the singular value decomposition (SVD) (or the power method \cite{mises1929praktische}) on the un-normalized sample covariance matrix $\hat \Sigma_{t+w,w}$.
\begin{remark}
An alternative possibility is to use the generalized likelihood ratio (GLR) statistic, where both $\rho$ and $u$ are estimated for each possible change location $\kappa$. The GLR statistic is
\[
\max_{\kappa<t} \left\{-\frac{t-\kappa}{2} \log(1+\hat{\rho}_{\kappa,t}) + \frac1{2\hat{\sigma}^2}\frac{\hat{\rho}_{\kappa,t}}{1+\hat{\rho}_{\kappa,t}}\sum_{i = \kappa+1}^t (\hat{u}_{\kappa,t}^\intercal x_i )^2 \right\},
\]
where $\hat{\rho}_{\kappa,t}$, $\hat{u}_{\kappa,t}$ are estimated from samples $\{x_i\}_{i = \kappa+1}^t$. However, this computation is more intensive since there is no recursive implementation for the GLR statistic, furthermore it requires growing memory\footnote{There are finite memory versions \cite{lai1999efficient} which, unfortunately, are equally complicated in their implementation.}. Therefore, we do not consider the GLR statistic in this paper.
\end{remark}

\subsection{Calibration} \label{sec:metric}

To fairly compare the detection procedures discussed in the previous section we need to properly calibrate them. Clearly the calibration process must be consistent with the performance measure we are interested in. It is exactly this point we are discussing next.

For a given stopping time $T$ we measure false alarms through the Average Run Length (ARL) expressed with $\mathbb{E}_\infty[T]$. For the detection capability of $T$ we use the (worst-case) Expected Detection Delay (EDD) proposed by Lorden \cite{lorden1971procedures}
\begin{equation}
\sup\limits_{\tau \geq 0} \supess \mathbb{E}_\tau[ (T - \tau)^{+} | T>\tau, x_1, \ldots, x_\tau],
\label{eq:lorden}
\end{equation}
which considers the worst possible data before the change (expressed through the $\supess$) and the worst possible change-time $\tau$.

We now consider scenarios that lead to the worst-case detection delay. For the largest eigenvalue procedure, assume $1\leq t-w+1\leq\tau< t$. Since for the detection we use $\lambda_{\max} (\hat \Sigma_{t,w})$ and compare it to a threshold, it is clear that the worst-case data before $\tau$ are the ones that will make $\lambda_{\max}$ as small as possible. We observe that
\begin{multline}
\lambda_{\max} (\hat{\Sigma}_{t,w})\\
=\lambda_{\max} ( x_{t-w+1} x_{t-w+1}^\intercal + \cdots +
x_{\tau} x_{\tau}^\intercal+\cdots+ x_{t} x_t^\intercal)\\
\geq\lambda_{\max} ( x_{\tau+1} x_{\tau+1}^\intercal +
\cdots+ x_{t} x_{t}^\intercal)
=\lambda_{\max} (\hat{\Sigma}_{t,t-\tau}),
\label{eq:worst}
\end{multline}
which corresponds to the data $x_{t-w+1},\ldots,x_\tau$, before the change, being all equal to zero. In fact, the worst-case scenario at any time instant $\tau$ is equivalent to forgetting all data before and including $\tau$ and restarting the procedure from $\tau+1$ using, initially, one, then two, etc.~outer products in the un-normalized sample covariance matrix, exactly as we do when we start at time 0. Due to stationarity, this suggests that we can limit ourselves to the case $\tau=0$ and compute $\mathbb{E}_0[T_{\rm E}]$ and this will constitute the worst-case EDD. Furthermore, the fact that in the beginning we do not normalize with the number of outer products, is beneficial for $T_{\rm E}$ since it improves its ARL.

We should emphasize that if we do not force the data before the change to become zero and use simulations to evaluate the detector with a change occurring at some time different from 0, then it is possible to arrive at misleading conclusions. Indeed, it is not uncommon this test to appear outperforming the exact CUSUM test for low ARL values. Of course this is impossible since the exact CUSUM is optimum for \textit{any} ARL in the sense that it minimizes the \textit{worst-case} EDD depicted in \eqref{eq:lorden}.

Let us now consider the worst-case scenario for Subspace-CUSUM. We observe that
\[
\mathcal S_t = (\mathcal S_{t-1})^{+} + (\hat{u}_{t+w}^\intercal x_t)^2 - d \geq 0 +(\hat{u}_{t+w}^\intercal x_t)^2 - d,
\]
suggesting that when $\mathcal S_t$ restarts this is the worst it can happen for the detection delay. We therefore understand that the well-known property of the worst-case scenario in the exact CUSUM carries over to Subspace-CUSUM. Again, because of stationarity, this allows us to fix the changetime at $\tau=0$. Of course, as mentioned before, because $\hat{u}_{t+w}$ uses data coming from the future of $t$, if our detector stops at some time $t$ (namely when for the first time we experience $\mathcal{S}_t\geq b$) then the \textit{actual} time of stopping must be corrected to $t+w$. A similar correction is not necessary for CUSUM because this test has the exact information for all parameters.

Threshold $b$ is chosen so that the ARL meets a pre-specified value. In practice, $b$ is determined by simulation. A very convenient tool in accelerating the estimation of ARL (which is usually large) is the usage of the following formula that connects the ARL of CUSUM to the average of the SPRT stopping time \cite{siegmund2013sequential}
\begin{equation}
\mathbb{E}_\infty [T_{\rm C}] = \frac{\mathbb{E}_\infty[T_{\rm SPRT}]}{\mathbb{P}_\infty (S_{T_{\rm SPRT}} \ge b)},
\label{eq:CUSUM-SPRT}
\end{equation}
where the SPRT stopping time is defined as
\[
T_{\rm SPRT} = \inf\{t>0: S_t \notin [0, b]\}.
\]
The validity of this formula relies on the CUSUM property that after each restart, $S_t$ is independent from the data before the time of the restart. Unfortunately this key characteristic is no longer true in the proposed Subspace-CUSUM scheme due to the fact that $\hat{u}_{t+w}$ uses data from the future of $t$. We could, however, argue that this dependence is weak. Indeed, as we have seen in Lemma\,\ref{lem:eigenvector}, each $\hat{u}_t$ is basically equal to $u$ plus some small random perturbation (estimation error with power of the order of $1/w$), with these perturbations being practically independent in time. As we observed with numerous simulations, estimating the ARL directly and through \eqref{eq:CUSUM-SPRT} (with $S_t$ replaced by $\mathcal{S}_t$), results in almost indistinguishable values even for moderate window sizes $w$. This suggests that we can use \eqref{eq:CUSUM-SPRT} to estimate the ARL of the Subspace-CUSUM as well. As we mentioned, in the final result we need to add $w$ to account for the future data used by the estimate $\hat{u}_{t+w}$.

\section{Analysis of the Largest Eigenvalue Procedure} \label{sec:theo}

It is clear that, in this work we are interested in promoting the Subspace-CUSUM detection procedure for the change detection problem of interest. Therefore, it would have been very supportive to this method to offer a theoretical analysis and derive formulas for the corresponding ARL and EDD. Even though such an analysis is possible it is unfortunately overly lengthy, for this reason, we postpone its presentation for a future publication. In this section we intend to characterize the ARL and EDD of the largest eigenvalue procedure which turns out to be simpler. In doing so we will also introduce some of the mathematical tools we are going to use in the (future) analysis of the Subspace-CUSUM.

\subsection{Link with Random Matrix Theory}
Since the study of ARL requires the understanding of the property of the largest eigenvalue under the null, i.e., the samples are i.i.d. Gaussian random vectors with zero-mean and identity covariance matrix, we first review some related results from random matrix theory.

There has been an extensive literature on the distribution of the largest eigenvalue of the sample covariance matrix, see, e.g., \cite{johnstone2001distribution, yin1988limit, baik2006eigenvalues}. There are two kinds of results typically available for eigenvalue distributions: the so-called {\it bulk} \cite{edelman2013random}, which treats a continuum of eigenvalues, and the {\it extremes}, which are the (first few) largest and smallest eigenvalues. Assume there are $w$ samples which are $k$-dimensional Gaussian random vectors with zero-mean and identity covariance matrix. Let $\hat{\Sigma}_w=\sum_{i=1}^w x_i x_i^\intercal$ denote the un-normalized sample covariance matrix. If $k/w \rightarrow \gamma >0$, the largest eigenvalue of the sample covariance matrix converges to $(1+\sqrt{\gamma})^2$ almost surely \cite{geman1980limit}.
To characterize the distribution of the largest eigenvalue, \cite{johnstone2001distribution} uses the  Tracy-Widom law \cite{TracyWidom96}. Define the center and scaling constants
\begin{equation}
\begin{aligned}
 \mu_{w,k}&=\big(\sqrt{w-1}+\sqrt{k}\big)^2, \\
 \sigma_{w,k}&=\big(\sqrt{w-1}+\sqrt{k}\big)\left(\frac1{\sqrt{w-1}}+\frac1{\sqrt{k}}\right)^{1/3}.
 \end{aligned}
 \label{eq:parameter}
 \end{equation}
 If $k/w \rightarrow \gamma <1$, then the centered and scaled largest eigenvalue converges in distribution to a random variable $W_1$ with the so-called Tracy-Widom law of order one \cite{johnstone2001distribution}:
\begin{equation}
 \frac{\lambda_{\max}(\hat{\Sigma}_{w})-\mu_{w,k}}{\sigma_{w,k}} \rightarrow W_1.
 \label{eq:tw}
 \end{equation}
The Tracy-Widom law can be described in terms of a partial differential equation and the Airy function, and its tail can be computed numerically (using for example the R-package {\small \textsf{RMTstat}}).

\subsection{Approximation of ARL Ignoring Temporal Correlation}
If we ignore the temporal correlation of the largest eigenvalues produced by the sliding window, we can obtain a simple approximation for the ARL. If we call $p = \mathbb{P}_\infty(\lambda_{\max}(\hat{\Sigma}_{t,w} ) > b)$ for $t \geq w$ then the probability to stop at $t$ is geometric and it is easy to see that the ARL can be expressed as $1/p$. Clearly, to obtain this result we must assume that $\mathbb{P}_\infty(\lambda_{\max}(\hat{\Sigma}_{t,w} ) > b)=p$ for $t<w$ as well, which is clearly not true. Since for $t<w$ the un-normalized sample covariance has less than $w$ terms, the corresponding probability is smaller than $p$. This suggests that $1/p$ is actually a lower bound to the ARL while $w+1/p$ an upper bound. If $w\ll 1/p$ then approximating the ARL with $1/p$ is quite acceptable. We can use the Tracy-Widom law to obtain an asymptotic expression relating the ARL with the threshold $b$. The desired formula is depicted in the following proposition.

\begin{proposition}[Approximation of ARL by ignoring temporal correlation]
For any $0<p\ll1$ we have $\mathbb{E}_\infty[T_{\rm E}]\approx 1/p$, if we select
\begin{equation}
b = \sigma_{w,k}b_p + \mu_{w,k},
\label{eq:threshold1}
\end{equation}
where $b_p$ denotes the $p$-upper-percentage point of $W_1$ namely $\mathbb{P}(W_1\geq b_p)=p$.
\label{ARL1}
\end{proposition}
\begin{proof}
The proof is straightforward and therefore ommitted.
\end{proof}

\subsection{Approximation of ARL Including Temporal Correlation}
Now we aim to capture the temporal correlation between detection statistics due to overlapping time windows. We leverage a proof technique developed in \cite{SiegmundYakirZhang2010}, which can obtain satisfactory approximation for the tail probability of the maximum of a random field.
For each $\hat{\Sigma}_{t,w}$,  define
\begin{equation}
Z_{t}=\lambda_{\max}(\hat{\Sigma}_{t,w}).
\end{equation}
Fig. \ref{overlapping} illustrates the overlap of two sample covariance matrices and provides necessary notation.
\begin{figure}[H]
\centerline{
\includegraphics[width = 0.4\textwidth]{fig1}

\caption{Illustration of two sample covariance matrices.}
\label{overlapping}
\end{figure}
\noindent We note that for any given $M > 0$,
\[
\mathbb{P}_\infty(T\leq M) = \mathbb{P}_\infty\left(\max_{1\leq t\leq M} Z_t  \geq b \right),
\]
which is the max over a set of correlated variables $\{ Z_t\}_{t=1}^M$.
Capturing the temporal dependence of $\{Z_t\}$ is challenging. For our analysis we recall Pearson's correlation between two random variables as
\[
{\rm corr}(X, Y) = \frac{\mathbb{E}[X Y] - \mathbb E[X] \mathbb E[Y]}{\sqrt{{\rm Var}(X)}\sqrt{{\rm Var}(Y)}}.
\]
We then have the following lemma that addresses the problem of interest.
\begin{lemma}[Approximation of local correlation]\label{CovZ}
Let
\begin{equation}
\beta = 1+\frac{2k^{\frac13}+3k^{\frac16}\frac{c_1}{\sqrt{w}}+\frac{c_1^2}{w}}{c_2^2},
\label{corr_approx}
\end{equation}
where $c_1 = \mathbb{E}[W_1]=-1.21$ and $c_2=\sqrt{{\rm Var}(W_1)}= 1.27$.
Then when $k/w\to\gamma\ll1$ and $\delta \ll w$,
\begin{equation}
{\rm corr}(Z_t,Z_{t+\delta}) \leq 1 - \beta\delta + o(\delta).
\label{corr}
\end{equation}
\end{lemma}

\begin{proof}
The proof is given in the Appendix.
\end{proof}

By leveraging the properties of the local approximation in \eqref{corr}, we can obtain an asymptotic approximation using the localization theorem \cite{siegmund2010tail,SiegmundYakirZhang2010}.  Define a special function $v(\cdot)$ which is closely related to the Laplace transform of the overshoot over the boundary of a random walk \cite{siegmund2007statistics}:
\[
v(x) \approx \frac{\frac{2}{x}[\Phi\big(\frac{x}{2}\big)-0.5]}{\frac{x}{2}\Phi\big(\frac{x}{2}\big)+\phi\big(\frac{x}{2}\big)},
\]
where $\phi(x)$ and $\Phi(x)$ are the pdf and cdf of the standard normal distribution $\mathcal{N}(0,1)$.
\begin{proposition}[ARL with temporal correlation]
For large values of $b$ we can write
\begin{equation}
\mathbb{E}_\infty[T_{\rm E}] = \left[\beta b' \phi(b') v\big(b'\sqrt{2\beta}\big)\right]^{-1}\big(1+o(1)\big),
\label{eq:ARL}
\end{equation}
where \[b' = \frac{b - [\sigma_{w,k}\mathbb{E}(W_1) + \mu_{w,k}]}{\sigma_{w,k}\sqrt{\rm{Var}(W_1)}}.\]
\label{ARL2}
\end{proposition}
\begin{proof}
The proof is detailed in the Appendix.
\end{proof}

We perform simulations to verify the accuracy of the threshold value obtained without and with considering the temporal correlation (Proposition\,\ref{ARL1} and Proposition\,\ref{ARL2}, respectively). The results are shown in Table\,\ref{eq:threshold}. We find that, indeed, the threshold, when temporal correlation \eqref{eq:ARL} is taken into account, is more accurate than its counterpart obtained by using the Tracy-Widom law \eqref{eq:threshold1}.
\begin{table}[!h]
\begin{center}
\caption{Comparison of the threshold $b$ obtained from simulations and using the approximations. Window length $w=200$, dimension $k = 10$. The numbers shown are $b/w$.}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
Target ARL & 5k  &  10k  &  20k  & 30k &  40k & 50k\\
\hline
Simulation   &    1.633  &  1.661 &   1.688   & 1.702   & 1.713  &  1.722   \\
\hline
Approx \eqref{eq:threshold1} & 1.738 & 1.763 & 1.787 & 1.800 &  1.809 & 1.816\\
\hline
Approx \eqref{eq:ARL} & 1.699  &  1.713  &  1.727   & 1.735 &   1.740  &  1.744 \\
\hline
\end{tabular}
\label{eq:threshold}
\end{center}
\end{table}

\subsection{Lower Bound on EDD using Marginal Power}
We now focus on the detection performance and present a tight lower bound for the EDD of the largest eigenvalue procedure.
\begin{proposition}
For large values of $b$ we have
\begin{equation}
\mathbb{E}_0[T_{\rm E}] \geq 2\frac{b^\prime + e^{-b^\prime}-1}{-\log(1+\rho)+\rho}\big(1+ o(1)\big).
\label{eq:EDD}
\end{equation}
where
\[
b' = \frac1{2\sigma^2(1+\rho)}\left[b\rho-(1+\rho)\sigma^2\log(1+\rho)\right].
\]
\label{EDD}
\end{proposition}
\begin{proof}
The proof is based on a known result for CUSUM \cite{siegmund2013sequential} and requires the derivation of the Kullback-Leibler divergence for our problem. Details are given in the Appendix.
\end{proof}
Consistent with intuition, in Proposition \ref{EDD}, the right-hand-side of (\ref{eq:EDD}) is a decreasing function of the SNR $\rho$,

\begin{figure}[H]
\centerline{\includegraphics[width = 0.4\textwidth]{fig2.pdf}}
\caption{Simulated EDD and lower bound as a function of the threshold $b$.}
\label{fig:EDD}
\end{figure}

Comparing the lower bound in Proposition \ref{EDD} with simulated average delay, as shown in Fig. \ref{fig:EDD}, we can show that in the regime of small detection delay (which is the main regime of interest), the lower bound serves as a reasonably good approximation.

\section{Numerical examples} \label{sec:numerical}
In this section, numerical results are presented to compare the three detection procedures. The tests are first
applied to synthetic data and the performance of the Subspace-CUSUM and largest eigenvalue test are compared against the CUSUM optimum performance. Then the performance of Subspace-CUSUM is optimized by selecting the most appropriate window size.

\subsection{Performance Comparison}
We perform simulations to compare the largest eigenvalue procedure, the Subspace-CUSUM procedure, and the exact CUSUM procedure. The threshold for each procedure is determined by Monte-Carlo simulation, as discussed in Section \ref{sec:metric}. Fig.\,\ref{fig:eddcompare} depicts EDD versus log-ARL for parameter values $k=5$, $\theta=1$, $\sigma^2 = 1$ and window length $w = 20$. The black line corresponds to the exact CUSUM procedure, which is clearly the best and it lies below the other curves. Subspace-CUSUM has always smaller EDD than the largest eigenvalue procedure and the difference increases with increasing ARL.

\begin{figure}[h!]
\centerline{\includegraphics[width = 0.4\textwidth]{fig3.pdf}}
\caption{Comparison of CUSUM, Subspace-CUSUM and the Largest eigenvalue procedure.}
\label{fig:eddcompare}
\end{figure}

\subsection{Optimal window size.}
\begin{figure}[h!]
\centerline{\includegraphics[width = 0.4\textwidth]{fig4}}
\caption{Minimal EDD vs ARL among window sizes from 1 to 50.}
\label{fig:minimalEDD}
\end{figure}
\begin{figure}[h!]
\centerline{\includegraphics[width = 0.4\textwidth]{fig5}}
\caption{Optimal window size resulting in smallest EDD.}
\label{fig:optimalw}
\end{figure}
We also consider the EDD/ARL curve where $w$ is optimized to minimize the detection delay at every ARL.
We first compute the EDD for window sizes $w=1,2,\ldots,50$ given each ARL value. Then we plot in Fig.\,\ref{fig:minimalEDD}  the lower envelope of EDDs corresponding to the optimal EDD achieved by varying $w$. We also plot the optimal value of $w$ as a function of ARL in Fig.\,\ref{fig:optimalw}. Even though the best EDD of the Subspace-CUSUM is diverging from the performance enjoyed by CUSUM this divergence we believe is slower than the increase of the optimum CUSUM EDD. One of the goals in the future publication regarding the analysis of Subspace-CUSUM is to show that this is indeed the case, which in turn will demonstrate that this detection structure is first-order asymptotically optimum.

\section{Conclusion}\label{sec:conclusion}

In this paper, we considered three detection procedures for the rank-one change in the covariance matrix: the largest eigenvalue procedure, the exact CUSUM procedure, and the Subspace-CUSUM procedure. For Subspace-CUSUM we perform a simultaneous estimate of the required subspace in parallel with its sequential detection. We avoid estimating all unknown parameters by following a worst-case analysis with respect to the subspace power. We were able to derive theoretical expressions for the ARL and an interesting lower bound for the EDD of the largest eigenvalue procedure. In particular we were able to handle the correlations resulting from the usage of a sliding window which is an issue that is not present in the off-line version of the same procedure. For the comparisons of the three competing detectors we discuss how it is necessary to calibrate each detector so that comparisons are fair. Comparisons were performed using simulated data and Subspace-CUSUM was found to exhibit a significantly better performance than the largest eigenvalue procedure. Ongoing work involves establishing first-order asymptotic optimality of the Subspace-CUSUM procedure by determining the optimal drift parameter $d$ and by relating the sliding window length $w$ to the desired ARL.

\section*{Acknowledgment}
The work of Liyan Xie and Yao Xie was supported by the US National Science Foundation under Grants CCF\,1442635, CMMI\,1538746, DMS\,1830210, and the Career Award CCF\,1650913. The work of George Moustakides was supported by the US National Science Foundation under Grant CIF\,1513373, through Rutgers University.

\bibliographystyle{IEEEtran}