Sketching for Sequential Change-Point Detection
Yang Cao\thanks{Yang Cao (Email: caoyang@gatech.edu) and Yao Xie (Email: yao.xie@isye.gatech.edu) are with the H. Milton Stewart School of Industrial and Systems Engineering
Georgia Institute of Technology
Atlanta
GA
USA. Andrew Thompson (Email: thompson@maths.ox.ac.uk) is with the Mathematical Institute
University of Oxford
Oxford
UK. Meng Wang (Email: wangm7@rpi.edu) is with the Department of Electrical
Computer \& Systems Engineering
Rensselaer Polytechnic Institute
Troy
NY
USA. Authors are listed alphabetically. Part of the paper has appeared in our GlobalSIP 2015 paper \cite{xie2015sketching
We study sequential change-point detection using sketches (linear projections) of high-dimensional signal vectors, by presenting the sketching procedures that are derived based on the generalized likelihood ratio statistic. We consider both fixed and time-varying projections, and derive theoretical approximations to two fundamental performance metrics: the average run length (ARL) and the expected detection delay (EDD); these approximations are shown to be highly accurate by numerical simulations. We also characterize the performance of the procedure when the projection is a Gaussian random projection or a sparse 0-1 matrix (in particular, an expander graph). Finally, we demonstrate the good performance of the sketching performance using simulation and real-data examples on solar flare detection and failure detection in power networks.
changepoint-sketching-journal_3rd_arxiv.tex
% Template for ICASSP-2013 paper; to be used with:
%          spconf.sty  - ICASSP/ICIP LaTeX style file, and
%          IEEEbib.bst - IEEE bibliography style file.
% --------------------------------------------------------------------------
%\documentclass[letterpaper,12pt,onecolumn]{IEEEtran}
\documentclass[letterpaper, journal]{IEEEtran}

\usepackage{amsmath,amsthm,graphicx,amssymb,fullpage,wrapfig}

\usepackage{setspace}

\usepackage{amsmath,amsthm,amssymb,algorithmic,algorithm,hyperref}
\usepackage[margin=1in]{geometry}
\usepackage{setspace}
%\doublespacing
\usepackage{color}
\usepackage{soul}
\usepackage{epstopdf}

\theoremstyle{remark}

\addtolength{\topmargin}{-0.2in}         % moves top margin up .2 in
\addtolength{\textheight}{.4in}
\addtolength{\oddsidemargin}{-0.25in}    % moves left margin
\addtolength{\textwidth}{0.55in}         % moves right margin


\hypersetup{colorlinks=true}

%\pagestyle{myheadings}

% Example definitions.
% --------------------
%\def\x{{\mathbf x}}
%\def\L{{\cal L}}
%\renewcommand{\baselinestretch}{1.4}

%\pagestyle{myheadings}

\newcommand{\ben}{\begin{eqnarray}}
\newcommand{\no}{\nonumber \\}
\newcommand{\een}{\end{eqnarray}}
\newcommand{\Expect}{\mathbb{E}}
\newcommand{\sign}{\mathop{\bf sgn}}
\newcommand{\transpose}{^{\intercal}}
\newtheorem{thm}{Theorem}
\newtheorem{corollary}{Corollary}
\newtheorem{lemma}{Lemma}
\newtheorem{mydef}{Definition}
\newcommand{\trace}{\mathop{\bf tr}}
\newcommand{\ka}{\kappa}
\newcommand{\vect}{\boldsymbol} %alternative use \vec
\newcommand{\SNR}{\mbox{SNR}}
\newcommand{\lbb}{\left\{}
\newcommand{\sgn}{{sgn}}
\newcommand{\rbb}{\right\}}
\newcommand{\lsb}{\left(}
\newcommand{\rsb}{\right)}
\newcommand{\lmb}{\left[}
\newcommand{\rmb}{\right]}
\newcommand{\y}{{y}}
\newcommand{\z}{{z}}
\newcommand{\x}{{x}}
\newcommand{\Z}{{Z}}
\newcommand{\w}{{w}}
\newcommand{\ab}{{a}}
\newcommand{\bb}{{b}}
\newcommand{\bm}{{b}}
\newcommand{\cb}{{c}}
\newcommand{\eb}{{e}}
\newcommand{\gb}{{g}}
\newcommand{\xb}{{x}}
\newcommand{\rb}{{r}}
\newcommand{\nb}{{n}}
\newcommand{\zb}{{z}}
\newcommand{\A}{{A}}
\newcommand{\B}{{B}}
\newcommand{\G}{{G}}
\newcommand{\R}{{R}}
\newcommand{\T}{{T}}
\newcommand{\X}{{X}}
\newcommand{\I}{{I}}
\newcommand{\D}{{D}}
\newcommand{\Pb}{{P}}
\newcommand{\E}{{E}}
\newcommand{\M}{{M}}
\newcommand{\V}{{V}}
\newcommand{\U}{{U}}
\newcommand{\W}{{W}}
\newcommand{\yb}{{y}}
\newcommand{\s}{{s}}
\newcommand{\e}{{e}}
\newcommand{\n}{{n}}
\newcommand{\vb}{{v}}
\newcommand{\wb}{{w}}
\newcommand{\la}{\langle}
\newcommand{\ra}{\rangle}
\newcommand{\iset}{\mathcal{I}}
\newcommand{\ihat}{\hat{\mathcal{I}}}
\newcommand{\kN}{|\mathcal{N}|}
\def\bfzero{{\mathbf{0}}}
\def\bfone{{\mathbf{1}}}
\def\bfa{{\mathbf a}}
\def\bfb{{\mathbf b}}
\def\bfc{{\mathbf c}}
\def\bfd{{\mathbf d}}
\def\bfe{{\mathbf e}}
\def\bff{{\mathbf f}}
\def\bfg{{\mathbf g}}
\def\bfh{{\mathbf h}}
\def\bfi{{\mathbf i}}
\def\bfj{{\mathbf j}}
\def\bfk{{\mathbf k}}
\def\bfl{{\mathbf l}}
\def\bfm{{\mathbf m}}
\def\bfn{{\mathbf n}}
\def\bfo{{\mathbf o}}
\def\bfp{{\mathbf p}}
\def\bfq{{\mathbf q}}
\def\bfr{{\mathbf r}}
\def\bfs{{\mathbf s}}
\def\bft{{\mathbf t}}
\def\bfu{{\mathbf u}}
\def\bfv{{\mathbf v}}
\def\bfw{{\mathbf w}}
\def\bfx{{\mathbf x}}
\def\bfy{{\mathbf y}}
\def\bfz{{\mathbf z}}
%bold caps
\def\bfA{{\mathbf A}}
\def\bfB{{\mathbf B}}
\def\bfC{{\mathbf C}}
\def\bfD{{\mathbf D}}
\def\bfE{{\mathbf E}}
\def\bfF{{\mathbf F}}
\def\bfG{{\mathbf G}}
\def\bfH{{\mathbf H}}
\def\bfI{{\mathbf I}}
\def\bfJ{{\mathbf J}}
\def\bfK{{\mathbf K}}
\def\bfL{{\mathbf L}}
\def\bfM{{\mathbf M}}
\def\bfN{{\mathbf N}}
\def\bfO{{\mathbf O}}
\def\bfP{{\mathbf P}}
\def\bfQ{{\mathbf Q}}
\def\bfR{{\mathbf R}}
\def\bfS{{\mathbf S}}
\def\bfT{{\mathbf T}}
\def\bfU{{\mathbf U}}
\def\bfV{{\mathbf V}}
\def\bfW{{\mathbf W}}
\def\bfX{{\mathbf X}}
\def\bfY{{\mathbf Y}}
\def\bfZ{{\mathbf Z}}

\def\yc{\textcolor{black}}

%	Department C-D\\
%	Address C-D}
%

%\setlength{\intercalmargin}{+0.5in}
\newtheorem{theorem}{Theorem}
\newtheorem{cor}{Corollary}
\newtheorem{remark}{Remark}
\newtheorem{defi}{Definition}
\newtheorem{prop}{Proposition}
\bibliographystyle{ieee}

%\usepackage{MnSymbol}

% Title.
% ------
%\title{Compressed Sensing Meets \\Sequential Change-Point Detection}

\title{Sketching for Sequential Change-Point Detection}
%
% Single address.
% ---------------
%\name{Author(s) Name(s)\thanks{Thanks to XYZ agency for funding.}}
%\address{Author Affiliation(s)}
%
% For example:
% ------------
%\address{School\\
%	Department\\
%	Address}
%
% Two addresses (uncomment and modify for two-address case).
% ----------------------------------------------------------
%\twoauthors
%  {A. Author-one, B. Author-two\sthanks{Thanks to XYZ agency for funding.}}
%	{School A-B\\
%	Department A-B\\
%	Address A-B}
%  {C. Author-three, D. Author-four\sthanks{The fourth author performed the work
%	while at ...}}
%	{School C-D\\
%	Department C-D\\
%	Address C-D}
%
\author{Yang Cao\thanks{Yang Cao
    (Email: caoyang@gatech.edu) and Yao Xie (Email: yao.xie@isye.gatech.edu)
   are with the H. Milton Stewart School of
    Industrial and Systems Engineering, Georgia Institute of
    Technology, Atlanta, GA, USA. Andrew Thompson (Email: thompson@maths.ox.ac.uk) is with the Mathematical Institute, University of Oxford, Oxford, UK. Meng Wang (Email: wangm7@rpi.edu) is with the Department of Electrical, Computer \& Systems Engineering, Rensselaer Polytechnic Institute, Troy, NY, USA.
    Authors are listed alphabetically. Part of the paper has appeared in our GlobalSIP 2015 paper \cite{xie2015sketching}.}, Andrew Thompson, Meng Wang, and Yao Xie
%and Robert Calderbank
}
\begin{document}
%\ninept
%
\maketitle
%
\begin{abstract}
We study sequential change-point detection using sketches (linear projections) of high-dimensional signal vectors, by presenting the sketching procedures that are derived based on the generalized likelihood ratio statistic. We consider both fixed and time-varying projections, and derive theoretical approximations to two fundamental performance metrics: the average run length (ARL) and the expected detection delay (EDD); these approximations are shown to be highly accurate by numerical simulations. We also characterize the performance of the procedure when the projection is a Gaussian random projection or a sparse 0-1 matrix (in particular, an expander graph). Finally, we demonstrate the good performance of the sketching performance using simulation and real-data examples on solar flare detection and failure detection in power networks.
\end{abstract}
%
%\begin{keywords}
%One, two, three, four, five
%\end{keywords}
%%

\section{Introduction}

Detecting change-points online from high-dimensional streaming data is a problem arising frequently from applications such as real-time monitoring for sensor networks, computer network anomaly detection and computer vision (e.g. \cite{changepoint_new_book2014, PoorHadjiliadis2008}). To reduce data dimensionality, a common approach is \emph{sketching} \cite{sketchingWoodruff2014}, which essentially performs random projection of the high-dimensional data vectors into lower-dimensional ones.

In this paper, we consider performing change-point detection using sketches of high-dimensional data vectors, and present a new sequential \emph{sketching procedure} utilizing the generalized likelihood ratio (GLR) statistics. In particular, suppose we may choose an $M\times N$ matrix $A$ with $M\ll N$, and consider projections of the vectors $y_t = Ax_t$, $t = 1, 2, \ldots$. Assume the pre-change vector is zero-mean Gaussian distributed, and the post-change vector is Gaussian distributed with an unknown mean vector $\mu$, and with the covariance matrix unchanged. Here we assume the mean vector is unknown since it typically represents an anomaly. The sketching procedure detects such a change by forming a GLR statistic, where the unknown $\mu$ is replaced with its maximum likelihood ratio estimator. Due to random projection of the data into a lower dimensional space, we are not able to find a unique estimator for the mean. Instead, there will be a family of estimators and we show that picking one such estimator from the family suffices. In this sense, we find that the task of detecting the presence of a signal from its linear projection is much easier than the corresponding tasks of recovering the signal or its support.

We obtain analytic expressions for two fundamental performance metrics: the average run length (ARL) when there is no change and the expected detection delay (EDD) when there is a change-point by extending the results in \cite{XieSiegmund2012}. Our approximations are shown to be highly accurate using simulations. These approximations are quite useful in determining the threshold of the detection procedure to control false-alarms, without having to resort to the onerous numerical simulations.

Our theoretical results also enable us to choose the sketching matrix. We examine the EDD ratio between the sketching procedure and a procedure using the original data, when the sketching matrix $A$ is either a random Gaussian matrix or a sparse 0-1 matrix (in particular, an expander graph). We find that the sketching procedure may achieve a performance very similar to that using the original data when the signal-to-noise ratio is sufficiently large, even when the post-change mean vector is not sparse. This result is consistent with our numerical results. Furthermore, we present an algorithm for time-varying projection and analyze its performance. Finally, we demonstrate the good performance of our procedures on a real data example of solar flare detection and a synthetic example of power failure detection in a network with real-world power network topology.


\yc{Change-point detection problems are closely related to industrial quality control and multivariate statistical control charts (SPC), where an observed process is assumed initially to be in control and at a change-point becomes out of control. The idea of using random projections for change detection has been explored for SPC in a pioneering work \cite{Runger96} based on $U^2$ multivariate control chart, the follow-up work \cite{Bodnar05} for cumulative sum (CUSUM) control chart and the exponential weighting moving average (EWMA) schemes, and in \cite{Skubalska2013,change-point-bookchapt15} based on the Hotelling statistic. These works provide a complementary perspective from SPC design, while our method takes a different approach and is based on sequential hypothesis testing, treating both the change-point location and the post-change mean vector as unknown parameters. By treating the change-point location as an unknown parameter when deriving the detection statistic, the sequential hypothesis testing approach overcomes the drawback of some SPC methods due to a lack-of-memory, such as the Shewhart chart and the Hotelling chart, since they cannot utilize the information embedded in the entire sequence \cite{SQC06}. Moreover, our sequential GLR statistic may be preferred over the CUSUM procedure in the setting when it is difficult to specify the post-change mean vector.}

%It is known that the GLR approach is very efficient when a change in the mean is relatively large in comparison to the total variance of the noise (CUSUM and EWMA methods based on the whole sequence of observations are usually chosen when a change in the mean is expected to be relatively small).
%

\yc{Besides the above distinctions from the SPC methods, other novelty of our methods also include: (1) we developed new theoretical results for the sequential GLR statistic, combining analytical techniques for sequential change-point detection and compressed sensing; (2) we consider the sparse 0-1 and time-varying projections (the sparse 0-1 projection corresponds to sampling the observations and is easily usable in practice); the analysis of these two cases are completely new and not taken from literature; (3) we study the amount of dimensionality reduction can be performed (i.e., the minimum $M$) with little performance loss.}

This paper \yc{extends} on our preliminary work reported in \cite{xie2015sketching} but we provide several important extensions: (1) we introduce time-varying sketching projections, which is more widely applicable, and provide theoretical performance analysis; (2) we include more extensive numerical examples to verify the accuracy of our theoretical results analysis; (3) we include new real-data examples to show the good performance of our sketching procedures for solar flare detection and power failure detection.

Our work is related to compressive signal processing \cite{compressiveSP2010}, where the problem of interest is to estimate or detect (in the fixed-sample setting) a sparse signal using compressive measurements. In \cite{arias2012detecting}, an offline test for a non-zero vector buried in Gaussian noise using noise linear measurements is studied; interestingly, a conclusion similar to ours is drawn that the task of detection within this setting is much easier than the tasks of estimation and support recovery. Another related work is \cite{GengXuLai2013}, which considers a problem of identifying a subset of data streams within a larger set, where the data streams in the subset follow a distribution (representing anomaly) that is different from the original distribution; the problem considered therein is not a sequential change-point detection problem as the ``change-point''  happens at the onset ($t = 1$). %Their approach is the classic dynamic programming approach (in the book by Poor and Hadjiliadis \cite{PoorHadjiliadis2008}), which is different from our log-likelihood test approach.
In \cite{XuLai2013}, the authors also consider an offline setting with one sequence of samples and the goal is to identify $k$ out of $n$ samples whose distributions are different from the normal distribution $f_0$. They use a ``temporal'' mixing of the samples over this finite horizon. This is different from our setting, where we project over the signal dimension at each time.  %change-point from a single data sequence where the number of sensors affected by the change-point is much smaller than the pre-change samples. They mix the sequences over time. Our setup is multi-sensor change-point detection, where the number of sensors affected by the change-point is small. So we ``mix'' over space instead of ``mixing'' over time. Moreover, their current version is not on-line algorithm -- it's basically a fixed sample change-point detection; whereas ours is a online change-point detection procedure.
%
In \cite{compressive_change2014}, change-point detection using ``compressive'' measurements of a high-dimensional vector is studied in a Bayesian framework (i.e., Shiryaev's procedure); theoretical analysis of average detection delay is presented. Other related methods include applying kernel methods \cite{harchaoui2009kernel} and \cite{arlot2012kernel} but they focus on offline change-point detection. Finally, using change-point detection for detecting a transient change in power systems has been studied in \cite{chen2016quickest}, and the method is tailored to a physical dynamic model of the real-world power system.

%Here, our sketches correspond to their compressive measurements; however, we consider a non-Bayesian setting, and find that the signal sparsity is not essential to achieve a performance similar to the one without sketching. %Moreover, we consider the non-Bayesian setting.
%We also give a precise approximation to ARL (as well as EDD), which is very useful for setting up the threshold.


Another possible approach is to obtain compressed representations of the data streams by applying principal component analysis (PCA) \cite{mishin2014real}. The dimensionality reduction by PCA is achieved by projecting the signal along a fixed subspace, which is defined by the eigenspace of the signal associated with the dominant eigenvalues. Our theoretical approximation for ARL and EDD can also be applied in those settings. \yc{However, when the dimension of the data is really large, it is hard to implement PCA-base method \cite{PCA_changepoint}.}


Our notations are standard: $\chi^2_k$ denotes the Chi-square distribution with degree-of-freedom $k$, $I_n$ denotes an identity matrix of size $n$; $X^\dag$ denotes the pseudoinverse of a matrix $X$; $[x]_i$ denotes the $i$th coordinate of a vector $x$; $[X]_{ij}$ denotes the $ij$th element of a matrix $X$; $x\transpose$ denotes the transpose of a vector or matrix $x$.
%\footnote{Due to space limit, the proofs for Theorems are omitted and they can be found at \href{http://www.isye.gatech.edu/~yxie77/sketching-changepoint.pdf}{http://www.isye.gatech.edu/$\sim$yxie77/sketching-changepoint.pdf}}.
%are delegated to appendix.

The rest of the sections are organized as follows. Section \ref{sec:formulation} sets up the formulation of the sketching problem for sequential change-point detection. Section \ref{sec:sketching} presents the sketching procedures for solving the problem. Section \ref{sec:analysis} presents the performance analysis for the sketching procedures. Section \ref{sec:numerical} demonstrates the accuracy of our performance analysis and the good performance of our sketching procedures on both simulated and real data. Section \ref{sec:summary} concludes the paper. All proofs are delegated to the appendix.

%\textcolor{red}{add a section about pair-wise comparison. It's a special case of this sketching. }

\section{Formulation}
\label{sec:formulation}

Consider a sequence of observations over a possibly infinite time horizon $x_1, x_2, \ldots, x_t$, $t = 1, 2, \ldots$, where $x_t \in \mathbb{R}^N$ and $N$ is the signal dimension. Initially the observations are due to noise. There can be a time $\kappa$ such that an unknown change-point occurs and it changes the mean of the signal vector. Such a problem can be formulated as the following hypothesis testing problem:
\begin{equation}
\begin{array}{ll}
\textsf{H}_0: & x_t  \sim \mathcal{N}(0, I_N), \quad t = 1, 2, \ldots\\
\textsf{H}_1: & x_t \sim \mathcal{N}(0, I_N), \quad t = 1, 2, \ldots, \ka, \\
& x_t \sim \mathcal{N}(\mu, I_N), \quad t = \ka + 1, \ka+2, \ldots
\end{array}
\label{hypothese}
\end{equation}
where the {\it unknown} mean vector is defined as
\[\mu \triangleq [\mu_1, \ldots, \mu_N]^\intercal \in \mathbb{R}^N.\]
%
Our goal is to detect the change-point as soon as possible after it occurs, while subject to the false alarm constraints. Here, we assume the covariance of the data to be an identity matrix, and the change only occurs to the mean. This models a large class of problems, such as environment monitoring and genomic copy number variation detection \cite{siegmund2011detecting} \cite{basseville1993detection}.


To reduce data dimensionality, we may project each signal $x_t$ into a lower dimensional space using a linear projection, which we refer to as \emph{sketching}. We aim to develop procedures that can detect a change-point using the resulted \emph{sketches}. In the following, we consider two approaches to sketching: the fixed projection and the time-varying projection.

\vspace{0.05in}

\noindent{\it Fixed projection.} Choose an $M\times N$ projection matrix $A$ with $M\ll N$ and obtain low dimensional sketches via:
\begin{equation}
y_t\triangleq Ax_t, \quad t = 1, 2, \ldots
\label{sketching}
\end{equation}

From the earlier model (\ref{hypothese}), the hypothesis test based on the sketches in (\ref{sketching}) can be reformulated as
\begin{equation}
\begin{array}{ll}
\textsf{H}_0: & y_t  \sim \mathcal{N}(0, AA\transpose), \quad t = 1, 2, \ldots\\
\textsf{H}_1: & y_t \sim \mathcal{N}(0, AA\transpose), \quad t = 1, 2, \ldots, \ka, \\
& y_t \sim \mathcal{N}(A\mu, AA\transpose), \quad t = \ka + 1, \ka+2, \ldots
\end{array}\label{new_hypothesis}
\end{equation}
Note that both mean and covariance structures are changed by the projections.

\vspace{0.05in}

\noindent{\it Time-varying projection.}  In certain applications, we may allow the sketching matrix to vary at each time. This may offer, for instance, a more robust data collection process.  Assume at each time, the projection $A_t \in \mathbb{R}^{M_t\times N}$, where the number of sketches $M_t$ can be time-varying as well. Under the time-varying model, the hypothesis test based on the sketches becomes:
\begin{equation}
\begin{array}{ll}
\textsf{H}_0: & y_t  \sim \mathcal{N}(0, A_tA_t\transpose), \quad t = 1, 2, \ldots\\
\textsf{H}_1: & y_t \sim \mathcal{N}(0, A_tA_t\transpose), \quad t = 1, 2, \ldots, \ka, \\
& y_t \sim \mathcal{N}(A_t\mu, A_tA_t\transpose), \quad t = \ka + 1, \ka+2, \ldots
\end{array}
\end{equation}

Note that we may allow the number of sketches to be different at each time as well. The above models also capture several important cases:
\begin{enumerate}
\item[(i)] (Pairwise comparison.)
In certain problems such as social network study and computer vision, we are interested in pairwise comparison of  variables \cite{chen2015information} \cite{massiminoone}. This will entail a total of $N^2$ possible comparison results, and we may select $M$ out of $N^2$ such comparisons to reduce complexity. Such pairwise comparison problems can be modeled as structured fixed projections $A$ with only $\{0, 1, -1\}$ entries, for instance:
\[
A =
\begin{bmatrix}
    1       & -1 & 0 & \dots & 0 & 0 \\
    1       & 0 & 0 & \dots & -1 & 0\\
    \hdotsfor{6} \\
    0       & 1 & 0 & \dots & 0& -1
\end{bmatrix} \in \mathbb{R}^{M\times N}.
\]

\item[(ii)] (Missing data.) In many scenarios, each time we are only able to observe a subset of entries of the signal, and the location of the entries that we can observe can vary with time \cite{balzano2015local}. Suppose at each time we observe a subset of entries $\Omega_t \in \{1,\ldots,N\}$. Assume $|\Omega_t| = M_t$, then $A_t \in \mathbb{R}^{M_t \times N}$ is a submatrix of an identity matrix with rows selected according to $\Omega_t$.
%then this corresponds to having a $A_t = \mathcal{R}(I_{\Omega_t}) \in \mathbb{R}^{M_t\times N}$, where $I_{\Omega_t}$ denotes the $N$-by-$N$ diagonal matrix whose $(i,i)$-entry is $1$ if $i \in \Omega_t$ and is $0$ otherwise, and $\mathcal{R}(X)$ denotes a operator that removes the whole-zero row of matrix $X$.

\item[(iii)] (PCA.) \yc{There are also methods using principal component analysis (PCA) of data streams for change-point detection (e.g., \cite{mishin2014real, Jackson91}). The main idea is to extract the several principle components of the data streams in a sliding time window, and performing change detection (e.g., using SPC control chart a CUSUM procedure) with these principle components.}  Note that this method can be viewed as using fixed projections $A$, with $A$ being the eigenspace of the signal associated with the dominant eigenvalues.
\end{enumerate}


\section{Sketching procedure}
\label{sec:sketching}

\subsection{Fixed projection}
We may derive a likelihood ratio based detection procedure for the hypothesis test (\ref{new_hypothesis}) using sketches.
%
%The statistics will involve the following two quantities.
%Define the running partial sum of the data up to time $t$ as \[S_{t} = \sum_{i=1}^t y_i, \quad t = 1, 2, \ldots, \] and
The statistic  involves an average of samples within a window $[k, t]$, which is defined as
\ben
\bar{y}_{k, t} =\frac{\sum_{i=k+1}^t y_i}{t-k}.
\label{ykt}
\een
%
Since the observations are i.i.d. over time, for an assumed value of the change-point $\ka = k$, \yc{for our hypothesis testing model in} (\ref{new_hypothesis}), the log-likelihood of observations accumulated up to time $t > k$ is given by
\begin{eqnarray}
&&\ell(t, k, \mu) \nonumber = \sum_{i=k+1}^t [y_i\transpose(AA\transpose)^{-1}A\mu - \frac{1}{2} \mu\transpose A\transpose(AA\transpose)^{-1}A\mu] \nonumber \\
&&= (t-k)[\bar{y}_{k, t} \transpose(AA\transpose)^{-1}A\mu - \frac{1}{2} \mu\transpose A\transpose(AA\transpose)^{-1}A\mu],
\label{stats}
\end{eqnarray}
\yc{where (\ref{stats}) is the log likelihood ratio (log-LR) statistic for detection.}


%\textcolor{red}{in this case we do not need to specify the estimation exactly.}

Since $\mu$ is unknown, we replace it with a maximum likelihood estimator for fixed $k$ and $t$ in the likelihood ratio (\ref{stats}), to obtain \yc{the log generalized likelihood ratio (log-GLR) statistic (one example of log-LR and log-GLR statistics is given in \cite{logGLR06}).}
%
Taking the derivative of $\ell(t, k, \mu)$ in (\ref{stats}) with respect to $\mu$ and setting it to zero, we obtain an equation for the maximum likelihood estimate $\mu^*$ of the post-change mean vector. As a function of the current number of observations $t$ and putative change-point location $k$, $\mu^*$ needs to satisfy
\begin{equation}
 A\transpose(AA\transpose)^{-1}A\mu^*=A\transpose(AA\transpose)^{-1}\bar{y}_{t, k}.\label{normal_eqn}
\end{equation}
Note that we may also write this as
\begin{equation*}
 A\transpose[(AA\transpose)^{-1}A\mu^*-(AA\transpose)^{-1}\bar{y}_{t, k}]=0.%\label{normal_eqn2}
\end{equation*}
Hence, a maximum likelihood estimator for the post-change mean vector is any $\mu^*$ that satisfies
\[(AA\transpose)^{-1}A\mu^* = (AA\transpose)^{-1}\bar{y}_{t, k} + c,\]
for a vector $c\in \mathbb{R}^N$ that lies in the null space of $A$: $A\transpose c = 0$. In particular, we may choose the zero vector $c = 0$, and use the estimator such that
\begin{equation}
(AA\transpose)^{-1}A\mu^* = (AA\transpose)^{-1}\bar{y}_{t, k}. \label{mu_star}
\end{equation}
%
Substituting such a $\mu^*$ into (\ref{stats}), we form the log generalized ratio statistic (GLR). Using (\ref{mu_star}), the first and second terms in (\ref{stats}) become, respectively,
\begin{align*}
\bar{y}_{k, t} \transpose(AA\transpose)^{-1}A\mu^*
&= \bar{y}_{k, t} \transpose  (AA\transpose)^{-1}\bar{y}_{t, k}, \\
\frac{1}{2} {\mu^*}\transpose A\transpose(AA\transpose)^{-1}A\mu^* &=
%\frac{1}{2} {\mu^*}\transpose A\transpose(AA\transpose)^{-1}\bar{y}_{t,k}
%=
 \frac{1}{2} \bar{y}_{k,t} \transpose(AA\transpose)^{-1} \bar{y}_{t,k}.
 \end{align*}
%
Finally, (\ref{stats}) gives \yc{the log-GLR statistic}:
\begin{equation}
\begin{split}
%&(t-k)[\bar{y}_{k,t} \transpose(AA\transpose)^{-1}A\mu^* - \frac{1}{2} \bar{y}_{k,t} \transpose(AA\transpose)^{-1}A\mu^*] \\
%= &
\ell(t, k, \mu^*) = \frac{t-k}{2}\bar{y}_{k,t} \transpose (AA\transpose)^{-1}\bar{y}_{k,t}.
\end{split}
\label{12}
\end{equation}
%which is the log generalized likelihood ratio (GLR) statistics.
Finally, we define a sketching procedure that stops and raises an alarm whenever the log-GLR statistic raises above a threshold $b > 0$:
\ben
T = \inf\{t\geq 1:
\max_{t-w\leq k < t}\frac{t-k}{2}\bar{y}_{k,t} \transpose (AA\transpose)^{-1}\bar{y}_{k,t} > b\}. \label{proc}
\een
Here $w>0$ is a window-size, i.e., we only consider the past $w$ samples from the current time $t$. The role of the window-size is two-fold. On the one hand, it reduces the memory requirements to implement the detection procedure, and on the other hand, it effectively establishes a minimum level of change that we want to detect.

We may further simplify the log-GLR statistic in (\ref{12}) by using the singular value decomposition (SVD) of $A$ and derive an equivalent procedure. This facilitates the performance analysis in Section \ref{sec:analysis}. Let the SVD of $A$ be given by
\begin{equation}
A = U\Sigma V\transpose, \label{A_SVD}
\end{equation}
where $U\in \mathbb{R}^{M\times M}$, $\Sigma \in \mathbb{R}^{M\times M}$ is a diagonal matrix containing all non-zero singular values, and $V\in \mathbb{R}^{N\times M}$. Using the SVD of $A$, $(AA^\intercal)^{-1} =  U\Sigma^{-2} U^\intercal$, we can write the log-GLR statistic (\ref{12}) as
\ben
\ell(t, k, \mu^*) = \frac{t-k}{2}\bar{y}_{k,t}\transpose U\Sigma^{-2}U\transpose \bar{y}_{k,t}.
\label{final3}\een
Substitution of (\ref{ykt}) into (\ref{final3}) results in a detection statistic
\[ \frac{\left\|\Sigma^{-1}U\transpose \left(\sum_{i=k+1}^t y_i\right)\right\|^2}{2(t-k)}.\] Define $z_i \triangleq \Sigma^{-1} U\transpose y_i$. Since under the null hypothesis $y_i \sim \mathcal{N}(0, AA\transpose)$, and hence, $z_i \sim \mathcal{N}(0, I_M)$. Similarly, under the alternative hypothesis $y_i \sim \mathcal{N}(A\mu, AA\transpose)$, and hence,  $z_i \sim \mathcal{N}(V\transpose \mu, I_M)$. Then we have the following equivalent form for the sketching procedure $T$ in (\ref{proc}):
\ben
T' = \inf\{t\geq 1:
\max_{t-w\leq k < t} \frac{\left\|\sum_{i=k+1}^t z_i\right\|^2}{2(t-k)}> b\}. \label{proc2}
\een
The form of $T'$ offers one intuitive explaination: the sketching procedure essentially projects the data into $M$ (less than $N$) independent data streams and detects the change by forming a log-GLR statistic using these independent data streams.

%where $w>0$ is a window-size (i.e., we only consider the past $w$ samples), and $b>0$ is a chosen threshold.%, and $w>0$ specifies the window size.
%\begin{remark}
%
%\end{remark}
%
%\begin{remark}
%Note that the procedure $T_2$ based on sketches essentially performs a partial sum of the data vector within a time window $[k+1, t]$ and scales by $\sqrt{t-k}$ (interestingly, instead of performing an average by dividing by $t-k$; a similar pattern has been observed in other Gaussian distribution based detection procedures), perform a transform by a $m$-by-$m$ matrix $\Sigma^{-1}U^\intercal$ and then compute the total energy.
%\end{remark}
%
%\begin{remark} This detection procedure can be implemented efficiently through a simple recursion for streaming data.
%\end{remark}

\subsection{Time-varying projection}

Similarly, for an assumed value of the change-point $\ka = k$, the \yc{log-LR} statistic for observations accumulated up to time $t$ is given by
\begin{equation}
\begin{split}
&\ell(t, k, \mu) \\
=& \sum_{i=k+1}^t [y_i\transpose(A_iA_i\transpose)^{-1}A_i\mu - \frac{1}{2} \mu\transpose A_i\transpose(A_iA_i\transpose)^{-1}A_i\mu].
\end{split}
\label{stats2}
\end{equation}
Again, \yc{to obtain the log-GLR statistic}, we replace the unknown $\mu$ with a maximum likelihood estimator using data between the hypothetical change location $k$ and current time $t$.
%
Taking the derivative of $\ell(t, k, \mu)$ in (\ref{stats2}) with respect to $\mu$ and setting it to zero, we obtain an equation for the maximum likelihood estimate $\mu^*$. As a function of fixed $k$ and $t$, $\mu^*$ needs to satisfy
\begin{equation}
\left[\sum_{i=k+1}^t A_i\transpose(A_iA_i\transpose)^{-1}A_i\right]\mu^*=\sum_{i=k+1}^t A_i\transpose(A_iA_i\transpose)^{-1}y_i.\label{normal_eqn2}
\end{equation}

To solve for $\mu^*$, one needs to discuss the rank of the matrix $\sum_{i=k+1}^t A_i\transpose(A_iA_i\transpose)^{-1}A_i$ in the bracket in (\ref{normal_eqn2}). Considering the SVD of $A_i = U_i D_i V_i^{\transpose}$ with $U_i \in \mathbb{R}^{M_i \times \M_i}$, $D_i \in \mathbb{R}^{M_i \times \M_i}$ and $V_i \in \mathbb{R}^{N \times M_i}$. We have that
\[
\sum_{i=k+1}^t A_i\transpose(A_iA_i\transpose)^{-1}A_i = \sum_{i=k+1}^t V_i V_i^{\transpose} = QQ^{\transpose},
\]
where $Q=[V_{k+1}, \ldots, V_t] \in \mathbb{R}^{N \times S}$ and $S = \sum_{i=k+1}^t M_i$. Consider two cases:
\begin{enumerate}
\item[(i)] $A_i'$s are independent Gaussian random matrices.  Note that the columns within each $V_i$ are linearly independent since they are orthonormal. Moreover, the columns in different $U_i$ blocks are independent since $U_i'$s are independent, and their entries are drawn as independent continuous random variables. Therefore, the columns of $Q$ are linearly independent with probability $1$ and rank($QQ^{\transpose}$) = $\min(N, (t-k)M)$ with probability one. We can claim that if $t-k<N/M$, $QQ^{\transpose}$ is not full rank and if $t-k\geq N/M$, $QQ^{\transpose}$ is full rank with probability one,
\item[(ii)] $A_i'$s are independent random matrices with entries drawn from a discrete distribution. In this case, we can also claim that if $t-k<N/M$, $QQ^{\transpose}$ is not full rank and if $t-k\geq N/M$, $QQ^{\transpose}$ is full rank with high probability, however, not with probability one.
\end{enumerate}
   %Since rank of $A_i\transpose(A_iA_i\transpose)^{-1}A_i \in \mathbb{R}^{M\times N}$ is less than or equal to $M$ for any $i\in [k+1,t]$, the matrix $\sum_{i=k+1}^t A_i\transpose(A_iA_i\transpose)^{-1}A_i$ is not full rank and thus is not invertible whenever $t-k < N/M$. If

To avoid the issue of invertibility, we use the pseudo-inverse instead. Define
\[B_{k,t} \triangleq \left(\sum_{i=k+1}^t A_i\transpose(A_iA_i\transpose)^{-1}A_i\right)^{\dag} \in \mathbb{R}^{N\times N},\] then from (\ref{normal_eqn2}), $\mu^*$ is given by
\[
\mu^* = B_{k,t}\sum_{i=k+1}^t A_i\transpose(A_iA_i\transpose)^{-1}y_i.
\]
Substituting such a $\mu^*$ into (\ref{stats2}), we obtain the \yc{log-GLR statistic}:
\begin{equation}
\begin{split}
&\ell(t, k, \mu^*) \\
=&\left(\sum_{i=k+1}^t A_i\transpose(A_iA_i\transpose)^{-1}y_i\right)\transpose B_{k,t} \left(\sum_{i=k+1}^t A_i\transpose(A_iA_i\transpose)^{-1}y_i\right).
\label{GLR2}
\end{split}
\end{equation}
We cannot further simplify the expression of GLR in (\ref{GLR2}) without specifying the form of $A_t, t=1,2,\ldots$. In such cases, computation of the pseudo-inverse may occur a high cost on the order of $\mathcal{O}(N^3)$.

In the special case when $A_t$ has $0$-$1$ entries only, which corresponds to the case for missing data and the expander graph, we can obtain a simpler expression for the statistic in (\ref{GLR2}) that can be computed efficiently. In this case, $A_t A_t \transpose$ is an $M_t$-by-$M_t$ identity matrix, and $A_t \transpose A_t$ is a diagonal matrix, which simplifies our computation. For a diagonal matrix $D\in \mathbb{R}^{N\times N}$ with diagonal entries $\lambda_1,\ldots, \lambda_N$, the pseudo-inverse of $D$ is also a diagonal matrix with diagonal entries $\lambda_i^{-1}$ if $\lambda_i \neq 0$ and with diagonal entries beginning with $0$ if $\lambda_i =0$. Therefore, the statistic in (\ref{GLR2}) can be simplified to
\begin{equation}
\ell(t, k, \mu^*) = \sum_{n=1}^N \frac{\left[\sum_{i=k+1}^t \mathcal{P}_{A_i}(x_i)\right]_n^2}{V_n(k,t)},
\end{equation}
where $\mathcal{P}_{A_i}(x_i)$ denotes an $N$ dimensional vector with the unobserved entries at time $i$ filled by $0$, and $V_n(k,t)$ denotes the number of times that the $n$th entry is observed during the time interval $[k+1,t]$.
Hence, for this special case, the sketching procedure with time-varying projection becomes:
\begin{equation}
\begin{split}
&T_2= \\
& \inf\{t\geq 1:
\max_{t-w\leq k < t} \frac{1}{2}\sum_{n=1}^N \frac{\left[\sum_{i=k+1}^t \mathcal{P}_{A_i}(x_i)\right]_n^2}{V_n(k,t)} > b\},
\end{split}
\label{proc3}
\end{equation}
where $b>0$ is the prescribed threshold. Note that the statistic essentially ``averages'' the observed vector within the time window $[t-w,t)$, which performs a kind of ``missing data imputation''.


\section{Performance analysis}
\label{sec:analysis}

\subsection{Approximations to ARL and EDD, fixed projection}

First, we introduce some notation for theoretical analysis. Under the null hypothesis $\textsf{H}_0$ (\ref{hypothese}), the observations have zero mean. Probability and expectation in this case are denoted by $\mathbb{P}^{\infty}$ and $\mathbb{E}^{\infty}$, respectively. Under the hypothesis $\textsf{H}_1$, alternatively, there exists a change-point $\kappa$, $0\leq \kappa < \infty$ such that the observations have mean $\mu$ for all $t>\kappa$. Probability and expectation in this case are denoted by $\mathbb{P}^{\kappa}$ and $\mathbb{E}^{\kappa}$, respectively.

The choice of the threshold $b$ involves a tradeoff between two standard performance metrics that are commonly used for analyzing change-point detection procedures \cite{XieSiegmund2012}: (i) the expected value of the stopping time when there is
no change, the average run length (ARL);
(ii) the expected detection delay (EDD), defined to be
the expected stopping time in the extreme case where
a change occurs immediately at $\ka = 0$, which is denoted as $\mathbb{E}^0\{T\}$. \yc{The following argument from \cite{SiegmundVenkat1996} explains why we consider $\mathbb{E}^0\{T\}$. When there is a change at $\kappa$, we are interested in the expected delay until its detection, i.e., the conditional expectation $\mathbb{E}^{\kappa}\{T-\kappa|T > \kappa\}$, which is a function of $\kappa$. When the shift in the mean only occurs in the positive direction $[\mu]_i \geq 0$, it can be shown that $\sup_\kappa \mathbb{E}^{\kappa}\{T-\kappa|T > \kappa\} = \mathbb{E}^0\{T\}$. It is not obvious that this remains true when $[\mu]_i$ can be either positive or negative. However, since $\mathbb{E}^0\{T\}$ is certainly of interest and reasonably easy to analyze, it is common to consider $\mathbb{E}^0\{T\}$.}





%
By extending the results in \cite{XieSiegmund2012} based on the observations in (\ref{proc2}), we can obtain  approximations to the ARL and the EDD as follows.
%
\begin{thm}[ARL]
\label{thmARL}
Assume that $1\leq M\leq N$, $b \rightarrow \infty$ with $M \rightarrow \infty$ and $b/M$ fixed.
Then for $w = o(b^r)$ for some
positive integer $r$, we have that the ARL of the sketching procedure defined in (\ref{proc}) is given by
%
\begin{equation}
\begin{split}
&\mathbb{E}^\infty\{T\} = \\
&\frac{2\sqrt{\pi}}{c(M, b, w) }  \frac{1}{1-\frac{M}{2b}} \frac{1}{\sqrt{M}}
\left(\frac{M}{2b}\right)^{\frac{M}{2}} e^{b-\frac{M}{2}} + o(1), \label{ET}
\end{split}
\end{equation}
where
\begin{equation}
c(M, b, w)  = \int_{\sqrt{\frac{2b}{w}}(1-\frac{M}{2b})}^{\sqrt{2b}(1-\frac{M}{2b})} u \nu^2(u) du, \label{c_def}
\end{equation}
and the special function \[\nu(u) = 2u^{-2} \exp[-2\sum_{i=1}^\infty i^{-1} \Phi(-|u|i^{1/2}/2)],\] (cf. \cite{Siegmund1985}, page 82). For numerical purposes, a simple and accurate approximation is given by (cf. \cite{SiegmundYakir2007})
\[
\nu(u) \approx \frac{2/u[\Phi(u/2) - 0.5]}{(u/2)\Phi(u/2) + \phi(u/2)}.
\]
\end{thm}



\begin{thm}[EDD]
\label{thmEDD}
Suppose $b \rightarrow \infty$ with other parameters held fixed.
Then for a given matrix $A$ with right singular vectors $V$, the EDD of the sketching procedure (\ref{proc}) when $\ka = 0$ is given by
\begin{equation}
\mathbb{E}^0\{T\} =
\frac{b+\rho(\Delta) - M/2 - \mathbb{E}\{\min_{t\geq 0} \tilde{S}_t\} + o(1)}{\Delta^2/2}.
\label{DD}
\end{equation}
where
\begin{equation}
\Delta = \|V\transpose \mu\|.  \label{delta_matrix}
\end{equation}
Here $\tilde{S}_t\triangleq \sum_{i=1}^t\delta_i$ is a random walk where the increments $\delta_i$ are independent and identically distributed with mean $\Delta^2/2$ and variance $\Delta^2$, $\mathbb{E}\{\min_{t\geq 0} \tilde{S}_t\} = -\sum_{i=1}^\infty \mathbb{E}\{\tilde{S}_i^-\}$ and $(x)^- = - \min\{x, 0\}$.
\end{thm}


%First, we will provide some intuition why compressive projection may help. The key lies in the fact that: the noise vector is not sparse, and hence projection will reduce energy of the noise vector. However, the signal vector is sparse, with high potability, with proper design of matrix $\A$, the norm of the ${\mu}$ vector will be preserved.


%
The proofs for the above two theorems utilize the fact explained in (\ref{proc2}).
%Recall that $y_t = A x_t$, and $A = U\Sigma V\transpose$. Hence, the log-GLR statistic in (\ref{proc2}) can be written as
%\begin{equation}
%\frac{1}{2}\left\|\Sigma^{-1}U\transpose \left(\frac{\sum_{i=k+1}^t y_i}{\sqrt{t-k}}\right)\right\|^2
%=\frac{1}{2}\left\| \left(\frac{\sum_{i=k+1}^t V\transpose x_i}{\sqrt{t-k}}\right)\right\|^2. \label{white_streams}
%\end{equation}
%
%Let $z_i = V\transpose x_i \in \mathbb{R}^M$, and so that the detection statistic is a function of $z_i$. Under the null hypothesis $\textsf{H}_0$, $x_i \sim \mathcal{N}(0,  I_{N})$; hence, $z_i \triangleq V\transpose x_i \sim \mathcal{N}(0, V\transpose V)$ and $V\transpose V =I_{M}$.
When there is no change-point, the sketching procedure is equivalent to monitoring $M$ data streams consisting of white noise with unit variances, and when there is a change-point, the post-change distribution of the observation vector has mean shifted to $V^{\transpose}\mu$ with the same covariance. This analysis demonstrates that the sketching procedure corresponds to a special case of the so-called ``mixture procedure'' (cf. $T_2$ in  \cite{XieSiegmund2012}) with $p_0 = 1$, $M$ sensors and the post-change mean vector is $V\transpose \mu$. %More detailed analysis are delegate to the appendix.

\begin{remark}
The first order approximation to the EDD is $b/(\Delta^2/2)$, which is the threshold $b$ divided by the Kullback-Leibler (K-L) divergence (see, e.g., \cite{XieSiegmund2012} for the fact that $\Delta^2/2$ is the K-L divergence between $\mathcal{N}(0, I_M)$ and $\mathcal{N}(V\transpose \mu, I_M)$). This makes sense intuitively since the expected increment of the detection statistics (viewed as a random process indexed over time) is roughly the K-L divergence of the test. %And here $\Delta^2/2$ is the Kullback-Leibler divergence between $\mathcal{N}(0, I_M)$ and $\mathcal{N}(V\transpose \mu, I_M)$, which is consistent with our analysis above.
\end{remark}

\begin{remark}
    We may establish the asymptotic optimality of $T'$ in (\ref{proc2}) using the following arguments. This can be based on our earlier work \cite{cao2015multi} that proves the asymptotic optimality for the more complicated non-i.i.d scenarios. Define $C(\gamma) = \{T:  \mathbb{E}^{\infty}\{T\} \geq \gamma \}$ as the set of the stopping times of which the ARL has a prescribed lower bound $\gamma$. We may show that $T'$ minimizes $\mathbb{E}^0\{T\}$ in $C(\gamma)$ as $\gamma \rightarrow \infty$, provided that $b=O(\log \gamma)$ and $\log w = o(\log \gamma)$.
\end{remark}

\begin{figure}
\begin{center}
\includegraphics[width = 0.7\linewidth]{b_versus_M.pdf}
\end{center}
\vspace{-0.2in}
\caption{For a fixed ARL being 5000, the curve for threshold $b$ versus $M$, when $N = 100$ and $w = 200$. Dashed is a tangent line of the curve.}
\label{fig:b_vs_M}
\vspace{-0.1in}
\end{figure}

\begin{remark}\label{remark_b_vs_M}
\yc{For a fixed large ARL, when $M$ increases, the threshold $b$ grows approximately linearly with $M$. Fig. \ref{fig:b_vs_M} is a such example, where $N = 100$, $w = 200$, we fix the ARL to be 5000, and find the corresponding threshold $b$ using Theorem \ref{thmARL} when $M$ increases from 10 to 100. (In Section \ref{sec:numerical}, we verify that Theorem \ref{thmARL} is an accurate approximation.)  Note that the resulting $b$ versus $M$ curve is very close to a straight line (shown as the dashed line). This can also be interpreted using Theorem \ref{thmARL}.
When $b\rightarrow \infty$, in (\ref{ET}), $c(M, b, w)$ and $1/[1-M/(2b)]$ tend to constants, and the remaining terms can be written as
\begin{equation}
\exp\{ b[1 - (1-\log\frac{M}{2b} ) \frac{M}{2b} -\frac{\log M}{2b}] \}.
\label{explanation}
\end{equation}
A common property of sequential change-point detection procedures is that the ARL grows exponentially with the threshold $b$ \cite{changepoint_new_book2014}, since we can achieve a large ARL for a reasonably valued threshold $b$. To achieve this, due to (\ref{explanation}), we need $b$ to grow at least linearly with $M$. }
%and the term that multiplies $b$ in the exponent tends to a positive constant. This property is desirable since we can obtain a large ARL for a reasonably valued threshold $b$.
%Hence, the remaining terms also tend to a constant under such a condition. This leads to a reasonable ARL expression. (\textcolor{red}{also needs further explanation ; the root has reasonable })
\end{remark}

\begin{remark}
\yc{To compare performance of the sketching procedure with that using the original data, we consider the the ratio between EDD using the original data (denoted as EDD($N$)) and EDD using the sketches (denoted as EDD($M$)).} Let %$\Gamma$ to be
\begin{equation}\label{gamma_def}
\Gamma\triangleq \frac{\|V\transpose \mu\|^2}{\|\mu\|^2}.
\end{equation}
From Theorem \ref{thmEDD}, we have that the EDD of the sketching procedure $\mathbb{E}^0\{T\}$ is $2b/\|V\transpose \mu\|^2 + o(1)$. \yc{Remark \ref{remark_b_vs_M} states that to achieve a constant ARL, $b$ should grow linearly with $M$.}  Hence, the ratio EDD($N$)$/$EDD($M$) is proportional to
\begin{equation}
 \frac{N}{M}\cdot
\frac{\|V\transpose\mu\|^2}{\|\mu\|^2} =  \frac{N}{M} \cdot \Gamma.
\label{ratio1}
\end{equation}
%Since the ratio in (\ref{ratio1}) is determined by $\Gamma$, we should study the property of $\Gamma$ for different $A$ to understand the performance of the sketching procedure.
\end{remark}



%
%When $\mu$ is unknown, $\mbox{EDD}(N)\leq \mbox{EDD}(M)$ and we would like their ratio to be as close to 1 as possible.
%
%\textcolor{red}{We will establish the following results..}
In Section \ref{sec:choice}, we will show that when $A$ is a Gaussian random matrix, $\Gamma$ is on the order of $M/N$. Hence, for Gaussian random matrix, this ratio is on the order of $(N/M)\cdot (M/N) = 1$. This result is intriguing, since it implies that there is little loss in performance by Gaussian random projection, as we indeed observe from numerical examples.

Moreover, we can show in Section \ref{sec:expander} that when $A$ is a sparse 0-1 matrix with $d$ non-zero entries on each row (in particular, an expander graph), $\Gamma$ is on the order of ${M(1-\epsilon)/(dN)}$. Hence, the ratio is on the order of $(N/M)\cdot [(M(1-\epsilon))/(dN)]  = (1-\epsilon)/(d)$,
for some small number $\epsilon > 0$. This implies that the sparse matrix may still be used for sketching purposes, but it may perform worse than the Gaussian random matrix by a factor of $1/d$. Such a result is expected as the sparse $A$ maps less information into the sketches in reducing the complexity of the measurements, since only few entries of the signals are required for each measurement. These theoretical predictions are also validated by numerical studies in Section \ref{sec:numerical}.

There is one possible intuitive explanation for the above. Unlike in compressed sensing, where our goal is to recover a sparse signal and hence one needs the projection to preserve norm up to a factor through the restricted isometry property (RIP) \cite{RIPCandes2008}, here our goal is to detect a non-zero vector in Gaussian noise. Hence, even though the projection reduces the norm of the vector, as long as the projection does not diminish the signal normal to be below the noise floor, we are still able to detect it. Since the Gaussian matrix is a dense matrix with probability one, it will achieve this goal; the expander graph, however, since it is sparse, may have some probability of missing the signal through projection.


\subsection{Choice of $A$ and bounding $\Gamma$}
\label{sec:choice}

%\input{expander_notes.tex}

%From the analysis above, note that the key lies in the choice of the sketching matrix $A$, and it affect the performance of sketching procedure through the ratio $\Gamma$ defined in (\ref{gamma_def}).
In the following, we establish bounds for $\Gamma$ when $A$ is a Gaussian matrix and an expander graph, respectively.
%From our expected delay analysis in (\ref{DD}), note that the performance depends on $\Delta$ defined in (\ref{delta_matrix}), which is essentially the Kullback-Leibler divergence.


\subsubsection{Gaussian matrix} \label{sec:Gaussian}

%First we present a simple lower-bound using similar arguments as above. Then we will derive the distribution of the ratio $\|V\transpose \mu\|/\|\mu\|$ and obtain a tighter bound.
%
Consider $A\in \mathbb{R}^{M\times N}$ whose entries are i.i.d. Gaussian
with mean zero and variance $1/M$. %Let $V \in\mathbb{R}^{m\times n}$ corresponds to the  right singular vectors.
%
We show that, provided $M$ and $N$ grow proportionally, $\Gamma$ approaches $\lim_{N\rightarrow\infty} M/N$ at a rate exponential in $N$.
%We first derive the precise distribution of $\Gamma$ under the assumption that $A$ is an i.i.d. Gaussian matrix.
\begin{lemma}[\cite{ruben79}]\label{quotient_dist}
Let $A\in\mathbb{R}^{M\times N}$ have i.i.d. $\mathcal{N}(0, 1)$ entries, and its  SVD be given by $A=U\Sigma V^T$. Then for any fixed vector $\mu$,
\begin{equation}\label{quotient_result}
\Gamma \sim\mbox{Beta}\left(\frac{M}{2},\frac{N-M}{2}\right).
\end{equation}
\end{lemma}
\yc{More related results can be found in \cite{betaDistribution90}.}
%
Since the $\mbox{Beta}(\alpha, \beta)$ distribution has mean $\alpha/(\alpha+\beta)$, we have
\[
\mathbb{E}(\Gamma) = \frac{M/2}{M/2 + (N-M)/2} = \frac{M}{N}.
\]
Next, we may also show that, provided $M$ and $N$ grow proportionally, $\Gamma$ converges to its mean value at a rate exponential in $N$. Define $\delta\in(0,1)$ to be
\begin{equation}\label{delta_def}
\delta\triangleq\displaystyle\lim_{N\rightarrow\infty}\frac{M}{N}.
\end{equation}
%We have the following result.
%
\begin{thm}[Gaussian $A$]\label{mean_concentration}
Let $A\in\mathbb{R}^{M\times N}$ have entries i.i.d. $\mathcal{N}(0, 1)$.  %independent of $\mu$, let $\Gamma$ be defined as in (\ref{gamma_def}), and
Let $N\rightarrow\infty$ such that (\ref{delta_def}) holds. Then  for $0<\epsilon<\min(\delta,1-\delta)$,
\begin{equation}\label{concentration_result}
\mathbb{P}(\delta-\epsilon<\Gamma<\delta+\epsilon)\rightarrow 1,
\end{equation}
at a rate exponential in $N$.
\end{thm}


%\vspace{0.2in}
%{\em Remark:} If ${\mu}$ is Normal and $A$ is any arbitrary matrix, i.e., ${\mu}_i \sim \mathcal{N}(0, 1)$, then $\|{\mu}\|^2 \sim \chi^2_N$, and $(V_0\transpose{\mu})_i \sim \mathcal{N}(0, 1)$ by orthogonal invariance. So $\|V\transpose{\mu}\|^2 \sim \chi_M^2$, $\|W\transpose{\mu}\|^2\sim \chi_{N-M}^2$. But $\|V\transpose{\mu}\|^2 + \|W\transpose{\mu}\|^2 = \|{\mu}\|^2 \sim \chi_N^2$, from which we see that (\ref{quotient_result}), and therefore the concentration result of Theorem~\ref{mean_concentration}, also holds in this case.


\begin{figure}
\begin{center}
\includegraphics[width = .5\linewidth]{expander_graph.pdf}
\end{center}
\caption{Illustration of a bipartite graph with $d=2$ and $c = 3$. }
\label{fig:EG}
\end{figure}

\subsubsection{Sparse sketching matrix $A$}\label{sec:expander}

Interestingly, it can be shown that we may also use a sparse 0-1 matrix $A$ for sketching (even if $\mu$ is not sparse), up to a small performance loss (i.e., longer EDD for fixed ARL), and provided that the post-change mean vector is element-wise positive (corresponds to the so-called ``one-sided'' change scenario typically considered in environmental monitoring \cite{XieSiegmund2012}). These sparse sketching matrices $A$ enable efficient sketching schemes, as each entry in the sketching vector only requires collecting information from few dimensions of the original data vector.

Assume $[\mu]_i\geq 0$ for all $i$. Consider
a matrix $A \in \mathbb{R}^{M\times N}$ consisting of binary entries, which corresponds to a bipartite graph. Following coding theory terminology, we call the left variables nodes (there are $N$ such variables), which correspond to the entries of $x_t$, and we call the right variables parity check nodes (there are $M$ such nodes), which correspond to entries of $y_t$. In a bipartite graph, connections between the variable nodes are not allowed. The adjacency matrix of the bipartite graph corresponds to our $A$ here. We further consider a bipartite graph with regular left degree $c$ (i.e., the number of edges from each variable node is $c$), and regular right degree $d$ (i.e., the number of edges from each parity check node is $d$), as illustrated in Fig. \ref{fig:EG}. Hence, this requires $Nc = Md$. % .


Expander graphs satisfy the above requirements, and they have been used in compressed sensing to sense a sparse vector (e.g. \cite{XuHassibi2007}). In particular, a matrix $A$ corresponds to a $(s,\epsilon)$-expander graph with regular right degree $d$ if and only if each column  of $A$ has exactly $d$ ``1''s, and for any set $S$ of right nodes with $|S|\leq s$, the set of neighbors $\mathcal{N}(S)$ of the left nodes has size $\mathcal{N}(S)\geq (1-\epsilon) d |S|$. If it further holds that each row of $A$ has $c$ ``1''s,  we say $A$ corresponds to a  $(s,\epsilon)$-expander with regular right degree $d$ and regular left degree $c$.
%
The existence of such expander graphs is established in \cite{BM01}:
%
\begin{lemma}[\cite{BM01}]\label{lem:goodexpander}
For any fixed $\epsilon >0$ and $\rho \triangleq M/N <1$, when $N$ is sufficiently large, there always exists an $(\alpha N, \epsilon)$ expander with a regular right degree $d$ and a regular left degree $c$ for some constants $\alpha\in(0,1)$, $d$ and $c$.
\end{lemma}

%, and they have certain good properties (in particular, holds for nonnegative signals as we show below).
\begin{thm}[Expander $A$]\label{thm:expander}
If $A$ corresponds to a $(s, \epsilon)$-expander with regular degree $d$ and regular left degree $c$, for any nonnegative vector $[\mu]_i\geq 0$, we have
\[
\Gamma \geq \frac{M(1-\epsilon)}{dN}.
\]
\end{thm}



\subsection{Approximation for ARL and EDD of time-varying projection}

It is intractable to directly analyze the general procedure $T_2$. However, we may characterize the performance of $T_2$ in (\ref{proc3}) for the ``missing data'' case when $A_i$ consists of $0$-$1$ entries and $M_t=M, t>0$, using the following scheme.
\begin{figure}
\begin{center}
\includegraphics[width = .7\linewidth]{binning.pdf}
\end{center}
\caption{A sampling scheme realization when $A_t$ is a $0$-$1$ matrix. Each dot denotes a sampled observation.}
\label{fig:binning}
\end{figure}
Assume each time we randomly sample $M$ out of $N$ entries to observe. Fig. \ref{fig:binning} explains the idea, which is a realization of the ``sampling'' scheme in procedure $T_2$, when $N=10$ rows, $t-k=17$ columns and we take $M=3$ subsamples each time. Dots denote the observed entries. Since we randomly sample $3$ observations each time, the sum of the number of the dots on each column in Fig. \ref{fig:binning} is $3$.

The idea for approximation is based on relaxing the constraint that each time we observe exactly $M$ out of $N$ entries, to sampling each entry of $x_i$ with probability $M/N$, for $1\leq n\leq N$. This eases finding an approximation to the EDD of $T_2$.
%
Define
$\xi_{ni}$'s as i.i.d. Bernoulli random variables with parameter \[r=\frac{M}{N},\] for $1\leq n\leq N$ and $i\geq 1$ and define
\[
Z_{n,k,t} \triangleq \frac{\sum_{i=k+1}^t [x_i]_n\xi_{ni}}{\sqrt{(t-k)r}}.
\]
Based on this, we introduce an approximate procedure:
\[
\widetilde{T}_2 = \inf \{t\geq 1: \max_{t-w\leq k < t} \frac{1}{2} \sum_{n=1}^N Z_{n,k,t}^2 > b \},
\]
where $b>0$ is the prescribed threshold.

Next, we further approximate $\widetilde{T}_2$ so that we can apply the previous analysis for ARL of (\ref{proc}) in Theorem $1$. Taking the expectation with respect to the Bernoulli random variables, we replace $\sum_{i=k+1}^t [x_i]_n \xi_{nj}$ by its expectation $r\sum_{i=k+1}^t [x_i]_n$ and $\sum_{j=k+1}^t \xi_{nj}$ by its expectation $(t-k)r$.
%The idea for approximation is to assume that for each fixed time window with size $t-k$, the sum of number of solid circles is about $(t-k)M/N$. Even not precise, such crude approximation provides convenient analysis that turns out to be good enough for ARL approximation.
%Based on the above argument, we introduce procedure $T_2'$ that is simpler for the analysis and approximates that performance of $T_2$.
Define
\[
U_{n,k,t} \triangleq \frac{\sum_{i=k+1}^t [x_i]_n}{\sqrt{t-k}},
\]
Then, a further approximation of $\widetilde{T}_2$ is given by
\[
T_2' = \inf \{t\geq 1: \max_{t-w\leq k < t} \frac{r}{2} \sum_{n=1}^N U_{n,k,t}^2 > b \},
\]
where $b>0$ is the prescribed threshold. Following a similar argument as that in Theorem $1$ by replacing $M$ with $N$ and $b$ with $b'=b/r$, we have that
\begin{equation}
\begin{split}
&\mathbb{E}^\infty\{T_2'\} \\
&\approx \frac{2\sqrt{\pi}}{c(N, b', w) } \frac{1}{\sqrt{N}} \frac{1}{1-\frac{N}{2b'}}
\left(\frac{N}{2}\right)^{\frac{N}{2}} b'^{-\frac{N}{2}} e^{b'-\frac{N}{2}},
\end{split} \label{ET2}
\end{equation}
where $c(N, b', w)$ is defined by replacing $b$ with $b'$ in (\ref{c_def}).
Table \ref{table_comp2} in Section \ref{sec:varyingnumerical} verifies that the theoretical threshold $b$ obtained by this approximation is very close to that obtained from simulation.
%The differences between the two approximation of ARLs (\ref{ET}) and (\ref{ET2}) are that we replace $M$ in (\ref{ET}) with $N$ and use $N/M$ as an additional multiplier. %\end{remark}


We can also obtain the first order approximation of EDD for the procedure in (\ref{proc3}) by using the approximation in $\widetilde{T}_2$. Note that
\begin{equation}
\mathbb{E}^0 \{Z_{n,0,t}^2\} = \frac{(rt\mu_n)^2+rt+rt\mu_n^2}{rt} = 1 + (rt+1)\mu_n^2.
\end{equation}
If the window size $w$ is sufficiently large, we have that
\begin{equation}
\mathbb{E}^0 \left\{\max_{k<t} \frac{1}{2}\sum_{n=1}^N Z_{n,k,t}^2\right\} \approx N/2+(rt+1)\sum_{n=1}^N \mu_n^2/2.
\label{lin1}
\end{equation}
Using Wald's identity \cite{Siegmund1985} and ignoring the overshoot of the statistic over the threshold $b$, we obtain a first order approximation by equating the right-hand side of (\ref{lin1}) to $b$, taking expectations on both sides, obtaining the equation
\[
N/2 + \mathbb{E}^0\{ T_2' \} \left(r\sum_{n=1}^N \mu_n^2/2\right) + \sum_{n=1}^N \mu_n^2/2 = b,
\]
from which we can solve for a first-order approximation for the expected detection delay\begin{equation}
\mathbb{E}^0\{T_2'\} \approx \left(\frac{2b-N}{\sum_{n=1}^N \mu_n^2}-1\right)\cdot \frac{N}{M}.
\label{DD2}
\end{equation}
Numerical examples in Section \ref{sec:varyingnumerical} demonstrates that this is a reasonably good approximation.


%Ignoring the overshoot to obtain a simpler result, we have that
%\begin{equation}
%\mathbb{E}^0 \{T_2\} \approx 2N\left(\frac{2b-1}{\sum_{n=1}^N \mu_n^2}-\left(1-\frac{M}{2N}\right)\right)/M,
%\label{DD2}
%\end{equation}
%where $\mu_n$s are the post-change means.





\section{Numerical examples}
\label{sec:numerical}


\yc{In the subsequent examples, we select ARL to be 5000 to represent a low rate of false detection (similar choice is made in other sequential change-point detection work such as \cite{XieSiegmund2012}). In practice, however, the target ARL value depends on sampling rate and how frequent we can tolerate false detection (once a month or once a year). All simulated results are obtained from $10^4$ repeatitions.
% In the following, we use simulation to determine the threshold $b$ corresponds to an ARL 5000. %, we use simulation For each detection procedure, fixing a threshold $b$, we can obtain a simulated ARL by applying the procedure on a standard multivariate normal process. Since the ARL is increasing with the threshold $b$, we apply a bisection search to obtain the threshold $b$ such that the simulated ARL is $5000$.
}

\subsection{Gaussian random matrix}
\label{sec:numerical_gaussian}
Consider $A$ generated as a Gaussian random matrix, with entries i.i.d. $\mathcal{N}(0, 1/N)$.

\subsubsection{Accuracy of theoretical approximations}

We use Theorem \ref{thmARL} to find the threshold $b$ corresponding to an ARL equal to 5000. \yc{Since the ARL is an increasing function of the threshold $b$, we use bisection to find the threshold $b$ that corresponds to a targeted ARL $5000$.} Then we compare it with a threshold $b$ that is found from simulation. As shown in Table \ref{table_comp}, the threshold find using Theorem \ref{thmARL} is very close to that obtained from simulation. \yc{Therefore, even if these theoretical ARL approximation is derived in an asymptotic regime when $N$ tends to infinity, it is applicable when $N$ is large but finite}. This is quite useful in determining a threshold for a targeted ARL, as simulations for large $N$ and $M$ can be quite time-consuming, for a large ARL (e.g., 5000 or 10,000).

Moreover, we also simulate the EDD detecting a signal with a post-change mean vector $\mu$ with entries all equal to $[\mu]_i = 0.5$.
As also shown in Table \ref{table_comp}, the approximations for EDD using Theorem \ref{thmEDD} are also resonably accurate.

\begin{table}[h!]
\caption{$A$ being a fixed Gaussian random matrix. $N$ = 100, $w = 200$, ARL = 5000, for simulated EDD $[\mu]_i = 0.5$.
\yc{Numbers in the parentheses are the standard deviation of the simulated EDD.} }
\begin{center}
\begin{tabular}{c||c|c||c|c}
\hline
$M$ & $b$ (theo)  & $b$ (simu) & EDD (theo) & EDD (simu) \\\hline
100 & 84.65 & 84.44 & 3.4 & 4.3 (0.9) \\
70 & 64.85 & 64.52 & 4.0 & 5.1 (1.2)\\
50 & 51.04 & 50.75 & 4.8 & 5.9 (1.6) \\
30 & 36.36 & 36.43 & 7.7 & 7.6 (2.5)\\
10 & 19.59 & 19.63 & 19.8 & 17.4 (9.8)\\\hline
\end{tabular}
\label{table_comp}
\end{center}
\end{table}

\subsubsection{Comparison of EDD}

Then, we compare EDD of the sketching procedure using the original data (which corresponds to $A = I$) with the sketching procedures using a Gaussian $A$, with $N=500$ and various $M < N$.

First, consider the post-change mean vector $[\mu]_i = \mu_0$. Fig. \ref{fig:EDD_comp}(a) shows EDD versus an increasing signal strength $\mu_0$. Clearly, when $\mu_0$ is sufficiently large, the sketching procedure can approach the performance using the original data as predicted by the theory in Section \ref{sec:choice}. From Fig. \ref{fig:EDD_comp}(a), we can also obtain the minimum $M$ required so that the EDD of the sketching procedure is within one plus the EDD using the original data, which we denote as $M^*$.
The results are shown in Table \ref{tab:minimumM_1}. When $\mu_0$ is sufficiently large, we may use $M$ to be even less than $30$ for a $N = 500$ dimensional signal. Note that we do not require signals to be sparse.

\begin{figure}[h]
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_Gaussian_changeMu.pdf}\\
(a) \\
\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_Gaussian_changeP.pdf}\\
(b)  %\\
%(a) \\
%\includegraphics[width=0.7\linewidth]{detectability}\\
%(b)
\end{tabular}
\end{center}
\caption{$A$ being a fixed Gaussian random matrix.  ${\rm EDD}_{\rm o}$ represents the EDD using the original data vector (can be viewed as setting $A = I$). \yc{The standard deviation of each point is less than half of its value.} (a) EDD versus $\mu_0$, when all $[\mu]_i = \mu_0$; (b) EDD versus $p$ when we randomly select $100p\%$ entries $[\mu]_i$ to be one and set the other entries to be zero; the smallest value of $p$ is 0.05.
%$\mu = $ $\mu\geq 1$, $M \geq 30$; $\mu\geq 0.9$, $M \geq 50$; $\mu\geq 0.7$, $M \geq 100$, $\mu\geq 0.6$, $M \geq 150$, $\mu\geq 0.5$, $M \geq 200$; $\mu\geq 0.3$, $M \geq 300$.
}
\label{fig:EDD_comp}
\end{figure}

\begin{table}[h]
\centering
\caption{$A$ being a fixed Gaussian random matrix. Minimum $M_*$ required for various mean shifts $\mu_0$. $N=500$, $w=200$ and all $[\mu]_i = \mu_0$.
\yc{Numbers in the parentheses are standard deviation of the simulated results.}}
\label{tab:minimumM_1}
\begin{tabular}{c|c|c|c|c|c}
\hline
               $\mu_0$      & $0.3$ & $0.5$ & $0.7$ & $1$ & $1.2$ \\ \hline
$M_*$ & 300        & 150        & 100        & 50       & 30 \\ \hline
$\mbox{EDD}_0$ & {\footnotesize 8.5 (2.0)} & {\footnotesize 4.7 (0.9)} & {\footnotesize 2.7 (0.5)} & {\footnotesize 2.0 (1.1)} & {\footnotesize 2 (0.01)} \\ \hline
\end{tabular}
\end{table}

%\begin{figure}[h]
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_Gaussian_changeP.pdf} %\\
%%(a) \\
%%\includegraphics[width=0.7\linewidth]{detectability}\\
%%(b)
%\end{tabular}
%\end{center}
%\caption{$A$ being a fixed Gaussian random matrix. EDD versus $p$ when we randomly select $100p\%$ entries $[\mu]_i$ to be one and set the other entries to be zero. ${\rm EDD}_o$ represents the EDD using the original data vector (can be viewed as setting $A = I$). \yc{The standard deviation of each point is less than half of its value.}
%\yc{The smallest value of $p$ is 0.05.}}
%\label{fig:EDD_comp_2}
%\end{figure}

\yc{Then, assume that the post-change mean vector has only $100p\%$ entries $\mu_i$ being one and the other entries being zero. Fig. \ref{fig:EDD_comp}(b) shows EDD versus an increasing $p$. As $p$ increases, the signal strength also increases, and the sketching procedure approaches the performance using the original data. Similarly, we can obtain the minimum $M$ required so that the EDD of the sketching procedure is within one plus the EDD using the original data. For example, when $p=0.5$, one can use $M=100$ for a $N=500$ dimensional signal to detect the change in the mean without much performance loss. }
%\textcolor{red}{add a table for Gaussian matrix, and another table for expander graph.}

\subsection{Expander graph}

Now consider $A$ generated as an expander graph. We compare the EDD for the change-point detection procedure using the original data (which corresponds to $A = I$) with the sketching procedures using an expander graph $A$ as above with $N=500$ and various $M < N$. \yc{To demonstrate the performance when $A$ is an expander graph, we run the simulations with the same settings as those in Section \ref{sec:numerical_gaussian}.} We verify the theoretical approximations are accurate for the expander graphs too (details omitted).

\subsubsection{Comparison of EDD}

Again, consider the post-change mean vector $[\mu]_i = \mu_0$. Fig. \ref{fig:EDD_expander}(a) shows EDD with an increasing $\mu_0$. \yc{Note that the simulated EDDs are much smaller than those for the Gaussian random projections in Fig. \ref{fig:EDD_comp}. A possible reason is that the expander graph is better at aggregating the signals when $[\mu]_i$ are all positive. However, when $[\mu]_i$ are can be either positive or negative, the two choices of $A$ have similar performance, as shown in Fig. \ref{fig:EDD_Twotype_M_posi_nega}, where $[\mu]_i$ are drawn i.i.d. uniformly from  $[-3, 3]$.}

Then, assume that the post-change mean vector has only $100p\%$ entries $\mu_i$ being one and the other entries being zero. Fig. \ref{fig:EDD_expander}(b) shows the simulated EDD versus an increasing $p$. As $p$ tends to $1$, the sketching procedure approaches the performance using the original data. From Fig. \ref{fig:EDD_expander}(b) we can obtain the minimum $M$ required so that the EDD of the sketching procedure is within one plus the EDD using the original data, as shown in Table \ref{tab:minimumM_2}. For example, when $p$ is around $0.5$, we may use $M=50$ for a $N=500$ dimensional signal. When $p$ is larger than $0.7$, one may use $M$ less than $30$ for a $N=500$ dimensional signal to detect the change, which implies that we can use a fairly small number of sketches to achieve a quick detection.

\begin{figure}[h]
%\vspace{-0.2in}
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_expander_changeMu.pdf}\\
(a)\\
\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_expander_changeP.pdf}\\
(b)
\end{tabular}
\end{center}
\caption{$A$ being a fixed expander graph. ${\rm EDD}_{\rm o}$ represents the EDD using the original data vector (can be viewed as setting $A = I$). \yc{The standard deviation of each point is less than half of its value.} (a) EDD versus $\mu_0$,  when all $[\mu]_i = \mu_0$; (b) EDD versus $p$ when we randomly select $100p\%$ entries $[\mu]_i$ to be one and set the other entries to be zero; the smallest value of $p$ is 0.05.}
\label{fig:EDD_expander}
\end{figure}

\begin{figure}[h]
%\vspace{-0.2in}
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_Twotype_changeM_posi_nega.pdf}
\end{tabular}
\end{center}
\caption{\yc{Comparison of EDDs for $A$ being a Gaussian random matrix versus an expander graph, when $[\mu]_i$'s are i.i.d. from $[-3,3]$. ${\rm EDD}_{\rm o}$ represents the EDD using the original data vector (can be viewed as setting $A = I$). }}
\label{fig:EDD_Twotype_M_posi_nega}
\end{figure}

%
%\begin{figure}[h]
%%\vspace{-0.2in}
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_expander_changeP.pdf}
%\end{tabular}
%\end{center}
%\caption{$A$ being a fixed expander graph. EDD versus $p$ when we randomly select $100p\%$ entries of mean vector with mean shift $1$. ${\rm EDD}_o$ represents the EDD using the original data vector (can be viewed as setting $A = I$). \yc{The standard deviation of each point is less than half of the value.}
%\yc{The smallest value of $p$ is 0.05.}}
%\label{fig:EDD_expander_p}
%\end{figure}

\begin{table}[h]
\centering
\caption{$A$ being a fixed expander graph. Minimum $M_*$ required for various $p$. $N=500$, $w=200$ and $100p\%$ of entries $[\mu]_i = 1$.
\yc{Number in parentheses are standard deviation of the simulated results.}}
\label{tab:minimumM_2}
\begin{tabular}{c|c|c|c|c|c}
\hline
                  $p$   & $0.1$ & $0.2$ & $0.3$ & $0.5$ & $0.7$ \\ \hline
$M_*$ & 300        & 200        & 100        & 50       & 30 \\ \hline
$\mbox{EDD}_0$ & {\footnotesize 4.4 (0.8)} & {\footnotesize 2.4 (0.5)} & {\footnotesize 2.0 (0.1)} & {\footnotesize 2.0 ($0.02$)} & {\footnotesize 2.0 ($0.01$)} \\ \hline
\end{tabular}
\end{table}


\subsection{Time-varying projection \yc{by 0-1 matrices}}
\label{sec:varyingnumerical}

\begin{table}[h!]
\caption{$A_t'$s being time-varying. $N$ = 100, $w = 200$, ARL = 5000, for simulated EDDs all $[\mu]_i = 0.5$.
\yc{Numbers in the parentheses are standard deviation of the simulated results.}}
\begin{center}
\begin{tabular}{c||c|c||c|c}
\hline
$M$ & $b$ (theo)  & $b$ (simu) & EDD (theo) & EDD (simu) \\\hline
100 & 84.65 & 84.44 & 2.8 & 3.3 (0.8) \\
70 & 83.72 & 83.41 & 3.8 & 4.5 (1.2)\\
50 & 82.84 & 83.02 & 5.3 & 6.1 (1.5) \\
30 & 81.46 & 82.48 & 8.7 & 9.8 (2.4)\\
10 & 78.32 & 79.27 & 23.4 & 26.6 (6.4) \\\hline
\end{tabular}
\label{table_comp2}
\end{center}
\end{table}

\begin{table}[h!]
\caption{$A_t'$s being time-varying. $N$ = 100, $w = 200$, ARL = 5000, for simulated EDDs $25\%$ entries $[\mu]_i = 1$.
\yc{Number in the parentheses are standard deviation of the simulated results.}}
\begin{center}
\begin{tabular}{c||c|c||c|c}
\hline
$M$ & $b$ (theo)  & $b$ (simu) & EDD (theo) & EDD (simu) \\\hline
100 & 84.65 & 84.44 & 2.8 & 3.3 (0.9) \\
70 & 83.72 & 83.41 & 3.8 & 4.5 (1.2)\\
50 & 82.84 & 83.02 & 5.3 & 6.1 (1.6) \\
30 & 81.46 & 82.48 & 8.7 & 9.8 (2.5)\\
10 & 78.32 & 79.27 & 23.4 & 26.7 (6.8) \\\hline
\end{tabular}
\label{table_comp3}
\end{center}
\end{table}


\subsubsection{Accuracy of theoretical approximations}

For time-varying projection, Table \ref{table_comp2} shows the accuracy of the approximation of $b'$s using the ARL in (\ref{ET2}) and the accuracy of the approximation of EDD in (\ref{DD2}) with various $M'$s if $N=100$, $w=200$, when all $[\mu]_i =  0.5$. Note that the approximate is more accurate if $M/N$ is close to $1$. One intuitive explanation is that,  when $M/N$ is close to $0$, replacing the number of observations for each entries by its mean $(t-k)M/N$ is less accurate, and hence, the approximation by using $T_2'$ for $T_2$ is less accurate.

\yc{To understand whether sparsity of the post-change mean vector affects the accuracy of the ARL approximation in (\ref{DD2}), we perform another experiment letting $N=100$, $w=200$ and  $25\%$ of $[\mu]_i = 1$. Note that the total signal energy $\sum_{n=1}^N \mu_n^2$  is the same as above, so the theoretical EDD computed from (\ref{DD2}) should be same. The simulated EDDs are shown in Table \ref{table_comp3}, which shows that the approximation in (\ref{DD2}) is still accurate.}

\subsubsection{Comparison of EDD}

To demonstrate the performance of time-varying projection $A_t$ defined in (\ref{proc3}) with $0$-$1$ entries, we again consider the two cases: the post-change mean vector $[\mu]_i = \mu_0$; the post-change mean vector has $100p\%$ entries $[\mu]_i$ being one and the other entries being zero. The simulated EDDs are shown in Fig. \ref{fig:EDD_randomrow_p}. Note that (\ref{proc3}) is able to detect change quickly with a small subset of observations. Even if EDDs of (\ref{proc3}) are larger than those for the fixed Gaussian matrix and the expander graph in Fig. \ref{fig:EDD_comp} and Fig. \ref{fig:EDD_expander}, this is meaningful, since in practice when performing subsampling for dimensionality reduction we usually scan through the signal by choosing a different subset of signal coordinates to observe.
%From the trend of the curves, we can also say that for each fixed $M$ the EDDs attain its minimum if the magnitude of mean shift is large enough.

\begin{figure}[h]
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_varying_changeMu}\\
(a)\\
\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_varying_changeP.pdf}\\
(b)
\end{tabular}
\end{center}
\caption{$A_t'$s are time-varying projection. ${\rm EDD}_{\rm o}$ represents the EDD using the original data vector (can be viewed as setting $A = I$). \yc{The standard deviation of each point is less than half of its value.} (a) EDD versus $\mu_0$, when all $[\mu]_i = \mu_0$; (b) EDD versus $p$ when we randomly select $100p\%$ entries $[\mu]_i$ to be one and set the other entries to be zero; the smallest value of $p$ is 0.05.}
\label{fig:EDD_randomrow_p}
\end{figure}

\subsection{Comparison with classical Multivariate CUSUM method}

\yc{We compare our sketching method with a benchmark approach based on the classical method. The benchmark method applies the classic multivariate CUSUM procedure \cite{woodall1985multivariate} directly to sketches. By doing so, it does not incoporate the covariance structure of the sketches. Moreover, since the benchmark approach is based on CUSUM, it needs a prescribed post-change mean vector (set to be an all-one vector in our example). Hence, its performance may be affected by parameter misspecification. We compare the performance again in two settings, when all $[\mu]_i$ are equal to a constant, and when $100p\%$ entries of the post-change mean vector are larger. Fig. \ref{fig:EDD_compare_multiCUSUM} demonstrates a significant amount of performance gain of our sketching procedure.}
%First, random projection by $A$ introduces correlations between the rows of the signal, which is not considered by the classical method. Second, the classical multivariate CUSUM procedure needs a prescribed post-change mean vector (which is set to be an all-one vector in our example). These two features are the main reasons that our method performs better than the classical one, as shown in Fig. \ref{fig:EDD_compare_multiCUSUM_Mu} and Fig. \ref{fig:EDD_compare_multiCUSUM_P}. In the examples, we assume  $100p\%$ entries of the post-change mean vector are increased by $\mu_0$. In both figures, for a fixed $M$, our method has smaller EDDs and smaller standard deviations than classical method.

  \begin{figure}[h]
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_Twotype_Classical_Mu.pdf}\\
(a)\\
\includegraphics[width=0.68\linewidth]{New_EDD_mesh_A_Twotype_Classical_P.pdf}\\
(b)
\end{tabular}
\end{center}
\caption{\yc{Comparison of the sketching procedure with a benchmark method based on the classic multivariate CUSUM procedure. (a) EDDs versus various $M$s, when all $[\mu]_i = 0.2$; (b) EDDs versus various $M$s, when we randomly select $10\%$ entries $[\mu]_i$ to be one and set the other entries to be zero.}}
\label{fig:EDD_compare_multiCUSUM}
\end{figure}

\subsection{Solar flare detection }

Consider a video from the Solar Data Observatory (SDO), which demonstrates an abrupt emergence of a solar flare\footnote{The video can be found at http://nislab.ee.duke.edu/MOUSSE/. The Solar Object Locator for the original data is SOL2011-04-30T21-45-49L061C108.}. Each frame is of size $232 \times 292$ pixels, which results in $N=67744$ dimensional streaming data. In this video, the normal states are slowly drifting solar flares, and the anomaly is a much brighter transient solar flare at $t=223$. The true change-point $t=223$ is hand-picked. Fig. \ref{fig:solarflare}(a) provides a snapshot of the original SDO data at $t=150$ when no solar flare exists. Fig. \ref{fig:solarflare}(b) provides a snapshot of the original SDO data at $t=223$ when the solar flare emerges as a brighter curve in the middle of the image, which is not clearly visible. Before running our detection procedures, we preprocess the data by tracking and removing the slowly changing background with the MOUSSE algorithm and detect the change in the residuals (which corresponds to $x_t$ in our formulation).

\yc{We first verify the Gaussianity for the residuals ($x_t$) by the one-sample Kolmogorov-Smirnov (KS) test. The p-value returned by the KS test is $0.47$ for the signal at $t=150$, which indicates that the Gaussian distribution is a reasonable assumption. As an additional proof of Gaussianity, Fig. \ref{fig:qqplot_solar_flare} shows the Quantile-Quantile (Q-Q) plot for the residuals at $t=150$.}


%presented in \cite{xie2013change}, and the results are shown in Fig. \ref{fig:solarflare} (c) and Fig. \ref{fig:solarflare} (d).

We apply the sketching procedure with fixed projection on the MOUSSE residuals, which performs dimensionality reduction. Choosing the sketching matrix $A$ to be an $M$-by-$N$ Gaussian random matrix with entries i.i.d. $\mathcal{N}(0,1/N)$. Note that the signal is deterministic in this case (since we use a deterministic video sequence). To evaluate consistency of the method, we run the procedure $500$ times, each time using a different random Gaussian matrix as the fixed projection $A$. Fig. \ref{fig:EDD_solarflare} shows the error-bars of the EDDs from $500$ runs. As $M$ increases, both the means and standard deviations of the EDDs decrease. When $M$ is larger than $750$, the EDD is less than $3$ with very high probability, which means that our sketching detection procedure can reliably detect the solar flare with only $750$ sketches. This is a significant reduction compared with using the original data with a dimensionality reduction ratio of $750/67744 \approx  0.01$.

\begin{figure}[h!]
\begin{center}
\begin{tabular} {cc}
\includegraphics[width=4cm,height=3cm]{original_noflare.pdf} &
\includegraphics[width=4cm,height=3cm]{original_flare.pdf} \\
(a) & (b) \\
% \includegraphics[width=4cm,height=3cm]{mousse_noflare.pdf} &
 %\includegraphics[width=4cm,height=3cm]{mousse_flare.pdf}\\
%(c) & (d) \\
\end{tabular}
\end{center}
\caption{Snapshot of the original solar flare data. (a) at $t=150$; (b) at $t=223$. The true change-point location is $t=223$.  %(c) MOUSSE tracking residual at t=150; (d) MOUSSE tracking residual at t=223.
}
%\label{fig:solar}
\label{fig:solarflare}
\end{figure}

\begin{figure}[h]
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.7\linewidth]{qqplot_solar_flare.pdf}
\end{tabular}
\end{center}
\caption{Solar flare detection example: Q-Q plot for $x_t$ (track error) at $t=150$ to verify Gaussian distribution.}
\label{fig:qqplot_solar_flare}
\end{figure}


\begin{figure}[h]
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.7\linewidth]{ErrorBar_solarflare.pdf}
\end{tabular}
\end{center}
\caption{Solar flare detection example: EDD versus various $M$ when $A$ is an $M$-by-$N$ Gaussian random matrix. The error-bars are obtained from $10^4$ repetitions with runs with different Gaussian random matrix $A$.}
\label{fig:EDD_solarflare}
\end{figure}


\subsection{Power system}

Finally, we consider a synthetic example based on the real power network topology of the Western States Power Grid of the United States, which consists of $4941$ nodes and $6594$ edges and the minimum degree of a node in the network is $1$, as shown in Fig. \ref{fig:topology}\footnote{The topology of the power network can be downloaded at http://networkdata.ics.uci.edu/data/power/ \cite{watts1998collective}.}. The nodes represent generators, transformers and substations, and edges represent high-voltage transmission lines between them \cite{watts1998collective}.
%
Note that the graph is sparse and that there are many ``communities'' which correspond to densely connected subnetworks. In this example, we simulate a situation for power failure over this large network. \yc{Assume at each time we may observe the real power injection at an edge. When the power system is in steady state, the observation is the true state plus Gaussian observation noise \cite{PowerSystemBook04}. We may estimate the true state (e.g., using techniques in \cite{PowerSystemBook04}), subtract it from the observation vector, and treat the residual vector as our signal $x_i$, which can be assumed to be i.i.d. standard Gaussian distribution.} When failures happen in a power system, there will be a shift in the mean for a small number of affected edges, since in practice, when there is a power failure, usually only a small part of the network is affected simultaneously.

To perform the sketching procedure, we consider $A_t'$s with $N=6594$ and various $M<N$. Each time, we randomly choose $M$ nodes in the network and measure the sum of the quantities over all attached edges as shown in Fig. \ref{fig:powerguide}. Note that in this case we construct a projection matrix that is neither a Gaussian matrix nor an expander graph, but the structure of the projection matrix is constrained by the network topology. Our experiment shows that our detection procedure performs well with such topology constrained $A_t'$s. Although our example is a highly simplified model for power networks, it sheds some light on the potential of our method applied to monitoring real power networks.

%In this example, we simulate a simple power outage scenario, and apply our detection scheme based on our sketching scheme. A recent work done for the power line outage detection is included in \cite{chen2016quickest}. When the power network works well, we can see the signals as drawn from i.i.d standard normal distribution (otherwise we can normalize the signals). To perform the sketching procedure, we consider $A$s with $N=6594$ and various $M<N$. Taking the practical setting into consideration, we can only collect the information on the edges attached to each node. When determining $A$, in other words, we randomly choose $M$ nodes in the network and then set the values of each row to be $0$ or $1$ to compute the sum of the observations of the edges attached to the corresponding node. Therefore, we construct a matrix $A$ that is neither Gaussian matrix nor expander graph, but a topology constrained matrix. Our next experiments show that our detection procedure performs also well with the topology constrained $A$s, which highlights the robustness of our method.

\begin{figure}[h]
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.7\linewidth]{powernetwroktopology.pdf}
\end{tabular}
\end{center}
\caption{Power network topology of the Western States Power Grid of the United States.}
\label{fig:topology}
\end{figure}


In the following experiment, we set the threshold $b$ such that the ARL is $5000$. We assume that the means of $5\%$ of the edges in the network increase by $\mu_0$. Fig. \ref{fig:EDD_topology}(a) shows the simulated EDD versus an increasing signal strength $\mu_0$. Note that the EDD using a small number of sketches is quite small if $\mu_0$ is sufficiently large. For example, when $\mu_0=4$, one may detect the change by observing from only $M=100$ sketches, which is a significant reduction compared with using the original data with a ratio of $100/4941 \approx 0.02$.

Next, we fix the mean shift of affected edges by $1$ and assume that $100p\%$ edges are affected by the change. Fig. \ref{fig:EDD_topology}(b) shows the simulated EDD versus increasing $p$. Similarly, the EDD using a small number of sketches is quite small if $p$ is large enough. If $p$ is around $0.4$, for example, one can detect the change by using only $M=50$ sketches.

\begin{figure}[h]
%\vspace{-0.2in}
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.7\linewidth]{powerguide.pdf}
\end{tabular}
\end{center}
\caption{Illustration of the measurement scheme for a power network. Suppose the physical quantities at edges (e.g, real power flow) at time $i$ form the vector $x_i$, and we can observe the sum of the edge quantities at each node. When there is a power failure, some edges are affected and their means are shifted.}
\label{fig:powerguide}
\end{figure}


\begin{figure}[h]
\begin{center}
\begin{tabular}{c}
\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_topology_changeMu.pdf}\\
(a)\\
\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_topology_changeP.pdf}\\
(b)
\end{tabular}
\end{center}
\caption{Power system example: $A$ being a power network topology constrained sensing matrix. \yc{The standard deviation of each point is less than half of the value.} (a) EDD versus $\mu_0$ when we randomly select $5\%$ edges with mean shift $\mu_0$; (b) EDD versus $p$ when we randomly select $100p\%$ edges with mean shift $1$. The minimum value of $p$ is $0.05$.}
\label{fig:EDD_topology}
\end{figure}

%\begin{figure}[h]
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_topology_changeP.pdf}
%\end{tabular}
%\end{center}
%\caption{$A$ is a power network topology constrained sensing matrix. EDD versus $p$ when we randomly select $100p\%$ edges with mean shift $1$. \yc{The standard deviation of each point is less than half of the value.} Results are from $10000$ Monte Carlo trials. \yc{The minimum value of $p$ is $0.05$.}}
%\label{fig:EDD_topology_p}
%\end{figure}
%\section{Discussions}
%
%If the signal has certain structure we may be able to exploit such structure. For instance, the post-change mean may lie on a low-dimensional subspace: $\mu = V^\star \beta$ for some basis $V^\star \in \mathbb{R}^{M\times N}$ and coefficient vector $\beta$. If we know such a $V^\star$, we could set $A = V^\star$ such that $V = V^\star$, and the ratio becomes
%\[
%\frac{\mbox{EDD}(N)}{\mbox{EDD}(M)} = \mathcal{O}(\frac{N}{M}\cdot
%\frac{\|{V^\star} \transpose V^\star \beta \|^2}{\|V^\star \beta\|^2})
%= \mathcal{O}(\frac{N}{M} \cdot \frac{\|\beta\|^2}{\|\beta\|^2}) = \mathcal{O}(\frac{N}{M}).
%\]
%%
%In general, we have that
%\[
%\frac{\mbox{EDD}(N)}{\mbox{EDD}(M)}
%\approx \mathcal{O}(\frac{N}{M}\cos(\theta_{\min}))
%\]
%where $\theta_{\min} \in [0, \pi/2]$ is the least principle angle between the two subspaces formed by $V$ and $V^\star$. Hence if we know the true subspace to a certain precision such that $\cos(\theta_{\min}) > N/M$ where the signal lies on, the sketching procedure may achieve an even better performance.





%\begin{figure}[h]
%\begin{center}
%\begin{tabular}{c}
%\includegraphics[width=0.7\linewidth]{New_EDD_mesh_A_Twotype_Classical_P.pdf}
%\end{tabular}
%\end{center}
%\caption{\yc{Comparison of our sketching procedure with classic multivariate CUSUM method. EDDs versus various $M$s, where $10\%$ entries of mean vector with mean shift $1$. Results are from $10000$ Monte Carlo trials. }}
%\label{fig:EDD_compare_multiCUSUM_P}
%\end{figure}

\section{Summary and Discussions}
\label{sec:summary}


In this paper, we studied the problem of sequential change-point detection if the observations are linear projections of the high-dimensional signals. We were interested in the problem where the change-point causes an {\it unknown} shift in the mean of the signal, and one would like to detect such a change as quickly as possible. We presented new sketching procedures for fixed and time-varying projections, respectively. The sketching procedures were derived based on the generalized likelihood ratio statistic. We analyzed theoretically the performance of our procedures by deriving theoretical approximations to the average run length (ARL) when there is no change, and the expected detection delay (EDD) when there is a change. Our approximations were shown to be highly accurate numerically. We also explored the choice of the projection matrix theoretically, when it is chosen as the Gaussian random matrix or as a sparse expander graph. We demonstrate the good performance of our procedure using numerical simulations and two real-data examples for solar flare detection and failure detection in power networks.

Thus far, we have assumed that the data streams are independent. In practice, if the data streams are dependent with a {\it known} covariance matrix $\Sigma$, we can whiten the data streams by applying a linear transformation $\Sigma^{-1/2}x_t$. \yc{Otherwise, the covariance matrix $\Sigma$ can also be estimated using a training stage via regularized maximum likelihood methods (see \cite{FanLiaoLiu2016} for an overview). Alternatively, we may estimate the covariance matrix $\Sigma'$ of the sketches $A\Sigma A\transpose$ or $A_t\Sigma A_t\transpose$ directly, which requires fewer samples to estimate due to the lower dimensionality of the covariance matrix. Then we may build statistical change-point detection procedures using $\Sigma'$ (similar to what has been done for the projection Hotelling control chart in \cite{change-point-bookchapt15}), which we leave for future work.}


One additional discussion would be the following. If the signal has certain structure we may exploit such structure. For instance, the post-change mean may lie on a low-dimensional subspace: $\mu = V^\star \beta$ for some basis $V^\star \in \mathbb{R}^{M\times N}$ and coefficient vector $\beta$. If we know such a $V^\star$, we could set $A = V^\star$, and the ratio $\mbox{EDD}(N)/\mbox{EDD}(M)$ is proportional to
\[
\frac{N}{M}\cdot
\frac{\|{V^\star} \transpose V^\star \beta \|^2}{\|V^\star \beta\|^2}
= \frac{N}{M} \cdot \frac{\|\beta\|^2}{\|\beta\|^2} = \frac{N}{M}.
\]
%
In general, we have that the ratio $\mbox{EDD}(N)/\mbox{EDD}(M)$ is proportional to $\mathcal{O}(\frac{N}{M}\cos(\theta_{\min}))$,
where $\theta_{\min} \in [0, \pi/2]$ is the least principle angle between the two subspaces formed by $V$ and $V^\star$. Hence if we know the true subspace in which the signal lies to a certain precision such that $\cos(\theta_{\min}) > N/M$, the sketching procedure may achieve an even better performance. \yc{Exploiting the knowledge of subspace has also been considered in \cite{Bodnar05} in the context of SPC design.}

\section*{Acknowledgement}

%Thank Qingbin Li for running partial results.
This work is partially supported by NSF grants CCF-1442635 and CMMI-1538746.
%\yc{1212}
%\clearpage

\bibliography{RD_proc}

%\clearpage

%\begin{figure}[h!]
%\begin{center}
%\includegraphics[width=.5\textwidth]{dd_Gaussian_m20}
%\includegraphics[width=.5\textwidth]{dd_Gaussian_m50}
%\end{center}
%\caption{Expected detection delay versus the ratio between sparsity and signal dimension $k/p$, using Gaussian random matrix $\A$ of dimension 50 by 200 for projection. Average run length (ARL) is 5000 for both the reduced-dimension and the full-dimension procedures. Typo: in label of x-axis, $k$ should be $s$.}
%\label{fig1}
%\end{figure}

%\newpage
%\begin{table}[h!]
%\caption{Thresholds, $p = 500$, Gaussian random matrix $\A$}
%\begin{center}
%\begin{tabular}{c|c}
%& $b$  (ARL = 5000)\\\hline
%$M = 50$ & 34.4\\
%$M = 100$ & 53.6 \\
%$M = 200$ & 88.1 \\
%$M = 300$ & 120.4 \\
%full (corresponds to $T_2(1)$ with \# sensors = 500) & 84.6 \\
%reduced ($n$ = 50) & 182
%\end{tabular}
%\end{center}
%\end{table}
%
%\begin{figure}[h!]
%\begin{center}
%\includegraphics[width=.4\textwidth]{gaussian_n50}
%\includegraphics[width=.4\textwidth]{gaussian_n100}
%\includegraphics[width=.4\textwidth]{gaussian_n200}
%\includegraphics[width=.4\textwidth]{gaussian_n300}
%\end{center}
%\caption{Expected detection delay versus the ratio between sparsity and signal dimension $s/p$, using Gaussian matrix $\A$ for the reduced dimension procedure. Average run length (ARL) is 5000 for both the reduced-dimension and the full-dimension (which corresponds to $\A = \I$) procedures. The K-L divergence $\Delta^2 = \sum_{i\in \mathcal{S}} \mu_i^2 = 5$ is kept constant with $\mu_i = \sqrt{5/|\mathcal{S}|}$. The subset of affected sensors $\mathcal{S}$ is random.}
%\label{fig1}
%\end{figure}
%
%
%
%
%
%
\clearpage


\appendix

\section{}

We start by deriving the ARL and EDD for the sketching procedure.
\begin{proof}[Proofs for Theorem \ref{thmARL} and Theorem \ref{thmEDD}]
This analysis demonstrates that the sketching procedure corresponds to the so-called mixture procedure (cf. $T_2$ in \cite{XieSiegmund2012}) in a special case of $p_0 = 1$, $M$ sensors, and the post-change mean vector is $V\transpose \mu$. In \cite{XieSiegmund2012}, Theorem 1, it was shown that the ARL of the mixture procedure with parameter $p_0\in[0, 1]$ and $M$ sensors is given by
%
\begin{equation}
\mathbb{E}^\infty\{T\} \sim H(M, \theta_0)/ \underbrace{\int_{[2M\gamma(\theta_0)/m_1]^{1/2}}^{[2M\gamma(\theta_0)/m_0]^{1/2}}  y \nu^2(y) dy}_{c'(M, b, w)}, \label{ET1}
\end{equation}
%
where the detection statistic will search within a time window $m_0\leq t-k\leq m_1$. Let $g(x,p_0) = \log(1-p_0 + p_0e^{x^2/2})$. Then $\psi(\theta) = \log \mathbb{E}\{e^{\theta g(U, p_0)}\}$ is the log moment generating function (MGF) for $g(U, p_0)$, $U\sim \mathcal{N}(0, 1)$,
%
 $\theta_0$ is the solution to $\dot{\psi}(\theta) = b/M$,
\begin{equation} \label{abbrev}
H(M,\theta) = \frac {\theta
[2\pi  \ddot{\psi}(\theta)]^{1/2}}{ \gamma(\theta) M^{1/2}}
\exp\{M[\theta \dot{\psi}(\theta) - \psi(\theta)]\},
\end{equation}
and
\[
\gamma(\theta) = \frac{1}{2}\theta^2
\mathbb{E}\{[\dot{g}(U, p_0)]^2 \exp[\theta g(U, p_0) - \psi(\theta)]\}.
\]%When $p_0=1$, $\phi(\theta)$ is the log moment generating function (MGF) of $U^2/2$, where $U\sim \mathcal{N}(0, 1)$.
Note that $U^2$ is $\chi^2_1$ distributed, whose MGF is given by $\mathbb{E}\{\e^{\theta U^2}\}=1/\sqrt{1-2\theta}$. Hence, when $p_0 = 1$,
\[
\psi(\theta) = \log\mathbb{E}\{e^{\theta U^2/2}\}
%= \log \{(1-2\frac{\theta}{2})^{-1/2}\}
 = -\frac{1}{2}\log(1-\theta).
\]
%
The first-order and second-order derivative of the log MGF are given by, respectively,
\ben
\dot{\psi}(\theta) = \frac{1}{2(1-\theta)}, \quad
\ddot{\psi}(\theta) =  \frac{1}{2(1-\theta)^2}
\een
Set $\dot{\psi}(\theta_0) = b/M$. We obtain the solution that $1-\theta_0 = M/(2b)$, and $\theta_0 = 1-M/(2b)$. Hence, $\ddot{\phi}(\theta_0) = 2b^2/M^2$.
%
%
We have $g(x,1) = x^2/2$, and $\dot{g}(x, 1) = x$.
\begin{align*}
\gamma(\theta) %&= \frac{1}{2}\theta^2 \mathbb{E}\{[\dot{g}(U)]^2 \exp[\theta g(U) - \psi(\theta)]\} \\
& = \frac{\theta^2}{2} \mathbb{E}\{U^2 e^{\frac{\theta U^2}{2}}\} e^{\log\sqrt{1-\theta}} \\
&= \frac{\theta^2}{2} \cdot \frac{1}{(1-\theta)^{3/2}} \cdot \sqrt{1-\theta} = \frac{\theta^2}{2(1-\theta)},
\end{align*}
where
\begin{align*}
\mathbb{E}\{U^2 e^{\frac{\theta U^2}{2}}\}
=& \frac{1}{\sqrt{2\pi}} \int x^2 e^{\frac{\theta x^2}{2}} e^{-\frac{x^2}{2}} dx \\
=& \frac{1}{\sqrt{2\pi}} \int x^2 e^{-\frac{x^2}{2/(1-\theta)}} dx = \frac{1}{(1-\theta)^{3/2}}.
\end{align*}
Combining the above, we have that the ARL of the sketching procedure is given by
\begin{equation}
\begin{split}
\mathbb{E}^\infty \{T\}
%&= \frac{\theta_0 [2\pi \ddot{\psi}(\theta_0)]^{1/2}}{c(N) \gamma(\theta_0)\sqrt{N}} e^{N[\theta_0 \dot{\psi}(\theta_0) - \psi(\theta_0)]} \\
&= \frac{\theta_0 [2\pi\cdot \frac{1}{2(1-\theta_0)^2}]^{1/2}}
{c'(M, b, w)\frac{\theta_0^2}{2(1-\theta_0)}\sqrt{M}} e^{\frac{M\theta_0}{2(1-\theta_0)}}(1-\theta_0)^{M/2} + o(1)\\
&= \frac{\sqrt{\pi}}{c'(M, b, w)}\frac{2}{\theta_0\sqrt{M}} e^{\frac{M\theta_0}{2(1-\theta_0)}}(1-\theta_0)^{M/2} + o(1).
\end{split}
\label{ET_1}
\end{equation}
Next, using the fact that $1/(1-\theta_0) = 2b/M$, we have that the two terms in the above expression can be written as
\[
\frac{M\theta_0}{2(1-\theta_0)} = \frac{M\theta_0}{2} \frac{2b}{M} = \theta_0 b, \quad (1-\theta_0) = \frac{M}{2b},
\]
then (\ref{ET_1}) becomes
\begin{align*}
&\mathbb{E}^\infty\{T\}=\frac{\sqrt{\pi}}{c'(M, b, w)} \frac{2}{\theta_0 \sqrt{M}} e^{\theta_0 b} (\frac{M}{2b}
)^{\frac{N}{2}} + o(1)\\
=& \frac{2\sqrt{\pi}}{c'(M, b, w)} \frac{1}{\sqrt{M}} \frac{1}{1-\frac{M}{2b}} e^{(b-\frac{M}{2})}  (\frac{N}{2b}
)^{\frac{N}{2}} + o(1)
\end{align*}
\begin{align*}
=& \frac{2\sqrt{\pi}}{c'(M, b, w)} \frac{1}{\sqrt{M}} \frac{1}{1-\frac{M}{2b}}
(\frac{M}{2})^{\frac{M}{2}} b^{-\frac{M}{2}} e^{b-\frac{M}{2}} + o(1).
\end{align*}
Finally, note that we can also write
\[\gamma(\theta_0) = \theta_0^2/[2(1-\theta_0)] = (1-M/(2b))^2/(M/b),\] and the constant is
\begin{equation}
\begin{split}
c'(M, b, w) =& \int_{[2M\gamma(\theta_0)/w]^{1/2}}^{[2M\gamma(\theta_0)]^{1/2}}  y \nu^2(y) dy \\
=& \int_{\sqrt{\frac{2b}{w}}(1-\frac{M}{2b})}^{\sqrt{2b}(1-\frac{M}{2b})} y \nu^2(y) dy.
\end{split}
\end{equation}
We are done deriving the ARL. The EDD can be derived by applying Theorem 2 of \cite{XieSiegmund2012} in the case where $\Delta = \|V\transpose\mu\|$, the number of sensors is $M$, and $p_0 = 1$.
\end{proof}

The following proof is for the Gaussian random matrix $A$.

%\begin{proof}[Proof of Lemma \ref{quotient_dist}]
%We first establish the orthogonal invariance of $V V\transpose$. It follows from (\ref{A_SVD}) that
%\begin{equation}\label{alternate}
%A^{\dag}A=V\Sigma^{-1}U\transpose U\Sigma V\transpose=V V\transpose.
%\end{equation}
%Now suppose $Y\in\mathbb{R}^{N\times N}$ is some fixed orthogonal matrix. Then $B=AY\transpose$ also follows a Gaussian distribution since a Gaussian matrix is invariant under multiplication by a fixed orthogonal matrix~\cite[Lemma 4.3]{mythesis}, and we have
%$$YV V\transpose Y\transpose=U A^{\dag}A Y\transpose=(AY\transpose)^{\dag}(AY\transpose)=B^{\dag}B,$$
%which has the same distribution as $V V\transpose$, which proves orthogonal invariance.
%Recalling (\ref{gamma_def}), we next use orthogonal invariance to show that $\Gamma$ shares the same distribution as $[V V\transpose]_{11}$. Choose $Y$ so that its first column is ${\mu}$ normalized, that is
%$$Y=\left[\begin{array}{cc}
%\frac{{\mu}}{\|{\mu}\|}\;\;Y_2\end{array}\right]\in\mathbb{R}^{N\times N}.$$
%Then, by the orthogonal invariance property,
%$$\frac{{\mu}\transpose YVV\transpose Y\transpose{\mu}}{\|{\mu}\|^2}\;\;\;\;\mbox{and}\;\;\;\;\frac{{\mu}\transpose VV\transpose{\mu}}{\|{\mu}\|^2}=\Gamma$$
%share the same distribution. But
%\begin{equation*}
%\begin{split}
%{\mu}\transpose YVV\transpose Y\transpose{\mu}
%=&{\mu}\transpose\left[\begin{array}{cc}
%\frac{{\mu}}{\|{\mu}\|}\;\;Y_2\end{array}\right]VV\transpose\left[\begin{array}{cc}
%\frac{{\mu}}{\|{\mu}\|}\;\;Y_2\end{array}\right]\transpose{\mu}\\
%=& [V V\transpose]_{11}\|{\mu}\|^2,
%\end{split}
%\end{equation*}
%where the last step follows from the orthogonality of the columns of $Y$. Dividing by ${\mu}$, we see that $\Gamma$ does indeed share the same distribution as $[VV\transpose]_{11}$. Finally, we show that $[VV\transpose]_{11}$ has the distribution given in (\ref{quotient_result}). It is well known that the truncated orthogonal matrix $V$ may be extended to a full orthogonal matrix $V_0=[V\;\;W]\in\mathbb{R}^{N\times N}$, such that $V_0$ is distributed according to the uniform Haar measure~\cite{EdelmanRao2005}. It follows that the first row of $V_0$ is uniformly distributed on the sphere in $\mathbb{R}^N$, and we may write
%$$V_0=\left[\begin{array}{c}
%\frac{z\transpose}{\|z\|} \\
%\tilde{V}\transpose
%\end{array}\right],$$
%where $z_i \sim \mathcal{N}(0, 1)$. Let us also write
%$$z=\left[\begin{array}{c}z_1\\z_2\end{array}\right], \quad z_1 \in \mathbb{R}^M, \quad  z_2 \in \mathbb{R}^{N-M}.$$
%Then the first row of $V$ is $z_1\transpose/\|z\|$, and so
%\begin{equation}
%\begin{split}
%[VV\transpose]_{11}&=\frac{z_1\transpose}{\|z\|}\cdot\frac{z_1}{\|z\|} =\frac{\|z_1\|^2}{\|z\|^2}= \frac{\sum_{i=1}^M z_i^2}{\sum_{i=1}^N z_i^2} \\
%&= \frac{\sum_{i=1}^M z_i^2}{\sum_{i=1}^M z_i^2 + \sum_{i =M+1}^N z_i^2}
%= \frac{P}{P+Q},
%\end{split}
%\end{equation}
%where $P\sim \chi_M^2$, $Q\sim \chi_{N-M}^2$. So $[V V\transpose]_{11} \sim \mbox{Beta}\left(\frac{M}{2}, \frac{N-M}{2}\right) $ \cite[26.5.3]{handbook} and the result is proved.
%\end{proof}


\begin{proof}[Proof of Theorem \ref{mean_concentration}]
It follows from (\ref{quotient_result}), and a standard result concerning the distribution function of the beta distribution~\cite[26.5.3]{handbook}, that
\begin{equation}\label{dist_fn}
\mathbb{P}\{\Gamma\le b\}=I_b\left(\frac{M}{2}, \frac{N-M}{2}\right),
\end{equation}
where $I$ is the regularized incomplete beta function (RIBF)~\cite[6.6.2]{handbook}. We first prove the lower bound in (\ref{concentration_result}). Assuming $N\rightarrow\infty$ such that (\ref{delta_def}) holds, we may combine (\ref{dist_fn}) with \cite[Theorem 4.18]{mythesis} to obtain
\begin{equation*}
\begin{split}
&\displaystyle\lim_{(M, N)\rightarrow \infty}\frac{1}{N} \ln\mathbb{P}\{\Gamma\le\delta-\epsilon\}\\
&= - \left[\delta \ln \left(\frac{\delta}{\delta-\epsilon}\right) + (1-\delta)\ln\left(\frac{1-\delta}{1-\delta + \epsilon}\right)\right] = -c < 0,
\end{split}
\end{equation*}
%
from which it follows that there exists $\tilde{N}$ such that, for all $N\geq\tilde{N}$,
$$\frac{1}{N} \ln\mathbb{P}\{\Gamma\le\delta-\epsilon\}<-\frac{c'}{2},$$
which rearranges to give
$$\mathbb{P}\{\Gamma\le\delta-\epsilon\} <e^{\frac{-c'N}{2}},$$
which proves the lower bound in (\ref{concentration_result}). To prove the upper bound, it follows from (\ref{dist_fn}) and a standard property of the RIBF~\cite[6.6.3]{handbook} that
\begin{equation}\label{comp_dist_fn}
\mathbb{P}\{\Gamma\geq b\}=I_{1-b}\left(\frac{N-M}{2}, \frac{M}{2}\right).
\end{equation}
Assuming $N\rightarrow\infty$ such that (\ref{delta_def}) holds, we may combine (\ref{comp_dist_fn}) with \cite[Theorem 4.18]{mythesis} to obtain
\begin{equation*}
\begin{split}
&\displaystyle\lim_{(M, N)\rightarrow \infty}\frac{1}{N} \ln\mathbb{P}\{\Gamma\geq\delta+\epsilon\}\\
&= - \left[(1-\delta)\ln \left(\frac{1-\delta}{1-\delta-\epsilon}\right) + \delta\ln\left(\frac{\delta}{\delta + \epsilon}\right)\right] = -d < 0,
\end{split}
\end{equation*}
and the argument now proceeds analogously to that for the lower bound.
\end{proof}

The following proofs are for the sparse 0-1 matrix $A$. Below, $A_{ij} = [A]_{ij}$.
\begin{lemma}\label{lem:2norm}
If a 0-1 matrix $A$ has constant column sum $d$, %corresponds to a $(s, \epsilon)$ expander with regular right degree $d$,
for every non-negative vector $x$ such that $[x]_i \geq 0$,  we have %(Meng: only hold for nonnegative signals)
%When $A$ is an $(k, \epsilon)$ expander, when $\bm \mu \geq \bf 0$, we have
\begin{equation}\label{eqn:lb}
\|A x\|_2 \geq \sqrt{d}\|x \|_2.
\end{equation}
\end{lemma}
\begin{proof}[Proof of Lemma \ref{lem:2norm}]
\begin{equation*}
\begin{split}
\|A x \|_2^2 & =  \sum_{i=1}^M \left(\sum_{j=1}^N A_{ij}x_j\right)^2 \\  & \geq   \sum_{i=1}^M \sum_{j=1}^N (A_{ij}x_j)^2  = d \|x \|_2^2.
\end{split}
              \end{equation*}
\end{proof}


\begin{lemma}[Bounding $\sigma_{\max}(A)$]\label{lem:sigma}
If $A$ corresponds to a $(s, \epsilon)$-expander with regular degree $d$ and regular left degree $c$, for any nonnegative vector $x$,
\begin{equation}
\frac{\|A x\|_2}{\|x\|_2} \leq d\sqrt{\frac{N}{M}},
\end{equation}
thus,
\begin{equation}\label{eqn:sigmamax}
\sigma_{\max}(A) \leq d \sqrt\frac{N}{M}.
\end{equation}
\end{lemma}



\begin{proof}[Proof of Lemma \ref{lem:sigma}]
For any nonnegative vector $x$,
\begin{eqnarray}
\|Ax\|_2^2 & = & \sum_{i=1}^M (\sum_{j=1}^N A_{ij}x_j)^2 \nonumber \\
              & = & \sum_{i=1}^M \big(\sum_{j=1}^N (A_{ij}x_j)^2 + \sum_{j=1}^N \sum_{l=1, l \leq j}^N (A_{ij}A_{il}x_jx_l) \big)  \nonumber \\
               & \leq & \sum_{i=1}^M \big(\sum_{j=1}^N (A_{ij}x_j)^2 + \sum_{j=1}^N \sum_{l=1, l \leq j}^N \frac{A_{ij}A_{il}}{2}(x_j^2+x_l^2) \big) \nonumber \\
               & = & \sum_{i=1}^M \sum_{j=1}^N \sum_{l=1}^N \frac{A_{ij}A_{il}}{2}(x_j^2+x_l^2)   \nonumber\\
                              & = &  \sum_{j=1}^N \sum_{l=1}^N \sum_{i=1}^M \frac{A_{ij}A_{il}}{2}(x_j^2+x_l^2)   \nonumber \\
                             % & =&\sum_{j=1}^p  \sum_{i=1}^n (A_{ij}x_j)^2 + \sum_{j=1}^p \sum_{l=1, l \leq j}^p  \sum_{i=1}^n \frac{A_{ij}A_{il}}{2}(x_j^2+x_l^2)  \\
                             % & =&\sum_{j=1}^p  d(x_j)^2 + \sum_{j=1}^p \sum_{l=1, l \leq j}^p  \sum_{i=1}^n \frac{A_{ij}A_{il}}{2}(x_j^2+x_l^2)  \\ \label{eqn:right}
                              & \leq&\sum_{j=1}^N  dc(x_j)^2 %+ \sum_{j=1}^p \sum_{l=1, l \leq j}^p  \sum_{i=1}^n \frac{A_{ij}A_{il}}{2}(x_j^2+x_l^2)
                               \label{eqn:left}\\
                               & = & \frac{d^2N}{M}\|x\|_2^2. \label{eqn:2norm}
\end{eqnarray}
%
Above, (\ref{eqn:left}) holds since  %$A$ has exactly $d$ `1's, and all others are `0's for every column $j$. Since $A$ has $c$ `1's in each row,
for a given column $j$, $A_{ij}=1$ holds for exactly $d$ rows. And for each row $i$ of these $d$ rows, $A_{il}=1$ for exactly $c$ columns with $l \in \{1,\ldots,p\}$;  % for a given row $i$ and a given
(\ref{eqn:2norm}) holds since $dN=Mc$. Finally, from the definition of $\sigma_{\max}$, (\ref{eqn:sigmamax}) holds.
\end{proof}




\begin{proof}[Proof for Theorem \ref{thm:expander}]
Note that
\begin{eqnarray} \label{eqn:norm}
\Delta &= &(\mu \transpose VV\transpose \mu)^{1/2}
   =  (\mu\transpose A\transpose U \Sigma^{-2} U\transpose A  \mu)^{1/2}  \nonumber\\
   & \geq & \sigma_{\max}^{-1}(A)\|U\transpose A  \mu\|_2     =  \sigma_{\max}^{-1}(A)\| A  \mu\|_2, \label{eqn:Delta}
\end{eqnarray}
where $\sigma_{\max}=\sigma_{\max}(A)$, and (\ref{eqn:Delta}) holds since $U$ is a unitary matrix.
Thus, in order to bound $\Delta$, we need to characterize $\sigma_{\max}$, as well as $\| A \mu\|_2$  for a $s$ sparse vector $\mu$.
Combining (\ref{eqn:Delta}) with Lemma \ref{lem:2norm} and \ref{lem:sigma}, we have that for every nonnegative vector $\mu$, $[\mu]_i \geq 0$,
\begin{equation}\label{eqn:Deltabound}
\Delta \geq  \frac{1}{d} \sqrt\frac{M}{N}  \sqrt{d(1-\epsilon)}\| \mu \|_2= \sqrt{\frac{M(1-\epsilon)}{dN}}\|\mu \|_2.
\end{equation}
%
Finally, Lemma \ref{lem:goodexpander} characterizes the quantity $[M(1-\epsilon)/(dN)]^{1/2}$ in (\ref{eqn:Deltabound}) and establishes the existence of such an expander graph. When $A$ corresponds to an  $(\alpha N, \epsilon)$ expander described in Lemma \ref{lem:goodexpander}, $\Delta \geq  \|\beta  \mu \|_2$ for all non-negative signals $[\mu]_i \geq 0$ for some constant $\alpha$ and some constant $\beta = (\rho(1-\epsilon)/d)^{1/2}$. Done.
\end{proof}






\end{document}

changepoint-sketching-journal_3rd_arxiv.bbl
\begin{thebibliography}{10}

\bibitem{xie2015sketching}
Yao Xie, Meng Wang, and Andrew Thompson,
\newblock ``Sketching for sequential change-point detection'',
\newblock in {\em Global Conference on Signal and Information Processing
  (GlobalSIP)}, 2015, pp. 78--82.

\bibitem{changepoint_new_book2014}
A.~Tartakovsky, I.~Nikiforov, and M.~Basseville,
\newblock {\em Sequential analysis: Hypothesis Testing and Changepoint
  Detection},
\newblock Chapman and Hall/CRC, 2014.

\bibitem{PoorHadjiliadis2008}
H.~V. Poor and O.~Hadjiliadis,
\newblock {\em Quickest detection},
\newblock Cambridge University Press, Dec. 2008.

\bibitem{sketchingWoodruff2014}
D.~P. Woodruff,
\newblock ``Sketching as a tool for numerical linear algebra'',
\newblock {\em Foundations and Trends in Theoretical Computer Science}, vol.
  10, pp. 1--157, 2014.

\bibitem{XieSiegmund2012}
Yao Xie and David Siegmund,
\newblock ``Sequential multi-sensor change-point detection'',
\newblock {\em The Annals of Statistics}, vol. 41, no. 2, pp. 670--692, 2013.

\bibitem{Runger96}
G.~C. Runger,
\newblock ``Projections and the {U}-squared multivariate control chart'',
\newblock {\em J. of Quality Technology}, vol. 28, no. 3, pp. 313--319, 1996.

\bibitem{Bodnar05}
O.~Bodnar and W.~Schmid,
\newblock ``Multivariate control charts based on a projection approach'',
\newblock {\em Allgemeines Statistisches Archiv}, vol. 89, pp. 75--93, 3005.

\bibitem{Skubalska2013}
E.~Skuba\l{l}ska-Rafajlowicz,
\newblock ``Random projections and {Hotelling's} t2 statistics for change
  detection in high-dimensional data analysis'',
\newblock {\em Int. J. Appl. Math. Comput. Sci.}, vol. 23, no. 2, pp. 447--461,
  2013.

\bibitem{change-point-bookchapt15}
E.~Skuba\l{}ska-Rafajlowicz,
\newblock ``Change-point detection of the mean vector with fewer observations
  than the dimension using instanenous normal random projections'',
\newblock in {\em Stochastic models, statistics and their applications}, vol.
  122. Springer Proc. Math. Stat., 2015.

\bibitem{SQC06}
D.~C. Montgomery,
\newblock {\em Introduction to statistical quality control},
\newblock Wiley, 2008.

\bibitem{compressiveSP2010}
Mark~A Davenport, Petros~T Boufounos, Michael~B Wakin, and Richard~G Baraniuk,
\newblock ``Signal processing with compressive measurements'',
\newblock {\em Selected Topics in Signal Processing, IEEE Journal of}, vol. 4,
  no. 2, pp. 445--460, 2010.

\bibitem{arias2012detecting}
Ery Arias-Castro et~al.,
\newblock ``Detecting a vector based on linear measurements'',
\newblock {\em Electronic Journal of Statistics}, vol. 6, pp. 547--558, 2012.

\bibitem{GengXuLai2013}
J.~Geng, W.~Xu, and L.~Lai,
\newblock ``Quickest search over multiple sequences with mixed observations'',
\newblock {\em arXiv:1302.3826}, Feb. 2013.

\bibitem{XuLai2013}
Weiyu Xu and Lifeng Lai,
\newblock ``Compressed hypothesis testing: To mix or not to mix?'',
\newblock in {\em Communication, Control, and Computing (Allerton), 2013 51st
  Annual Allerton Conference on}. IEEE, 2013, pp. 1443--1449.

\bibitem{compressive_change2014}
George~K Atia,
\newblock ``Change detection with compressive measurements'',
\newblock {\em Signal Processing Letters, IEEE}, vol. 22, no. 2, pp. 182--186,
  2015.

\bibitem{harchaoui2009kernel}
Zaid Harchaoui, Eric Moulines, and Francis~R Bach,
\newblock ``Kernel change-point analysis'',
\newblock in {\em Advances in Neural Information Processing Systems}, 2009, pp.
  609--616.

\bibitem{arlot2012kernel}
Sylvain Arlot, Alain Celisse, and Zaid Harchaoui,
\newblock ``Kernel change-point detection'',
\newblock {\em arXiv preprint arXiv:1202.3878}, 2012.

\bibitem{chen2016quickest}
Yu~Christine Chen, Taposh Banerjee, Alejandro~D Dom{\'\i}nguez-Garc{\'\i}a, and
  Venugopal~V Veeravalli,
\newblock ``Quickest line outage detection and identification'',
\newblock {\em Power Systems, IEEE Transactions on}, vol. 31, no. 1, pp.
  749--758, 2016.

\bibitem{mishin2014real}
Dmitry Mishin, Kieran Brantner-Magee, Ferenc Czako, and Alexander~S Szalay,
\newblock ``Real time change point detection by incremental pca in large scale
  sensor data'',
\newblock in {\em High Performance Extreme Computing Conference (HPEC), 2014
  IEEE}. IEEE, 2014, pp. 1--6.

\bibitem{PCA_changepoint}
F.~Tsung and K.~Wang,
\newblock ``Adaptive charting techniques: Literature review and extensions'',
\newblock in {\em Frontier in Statistical Quality Control}, pp. 19--35.
  Springer-Verlag, 2010.

\bibitem{siegmund2011detecting}
David Siegmund, Benjamin Yakir, and Nancy~R Zhang,
\newblock ``Detecting simultaneous variant intervals in aligned sequences'',
\newblock {\em The Annals of Applied Statistics}, vol. 5, no. 2A, pp. 645--668,
  2011.

\bibitem{basseville1993detection}
Mich{\`e}le Basseville, Igor~V Nikiforov, et~al.,
\newblock {\em Detection of abrupt changes: theory and application}, vol. 104,
\newblock Prentice Hall Englewood Cliffs, 1993.

\bibitem{chen2015information}
Yuxin Chen, Changho Suh, and Andrea~J Goldsmith,
\newblock ``Information recovery from pairwise measurements: A
  shannon-theoretic approach'',
\newblock in {\em Information Theory (ISIT), 2015 IEEE International Symposium
  on}. IEEE, 2015, pp. 2336--2340.

\bibitem{massiminoone}
Andrew~K Massimino and Mark~A Davenport,
\newblock ``One-bit matrix completion for pairwise comparison matrices'',
\newblock in {\em Proc. Workshop on Signal Processing with Adaptive Sparse
  Structured Representations (SPARS)}.

\bibitem{balzano2015local}
Laura Balzano and Stephen~J Wright,
\newblock ``Local convergence of an algorithm for subspace identification from
  partial data'',
\newblock {\em Foundations of Computational Mathematics}, vol. 15, no. 5, pp.
  1279--1314, 2015.

\bibitem{Jackson91}
J.~E. Jackson,
\newblock {\em A user's guide to principle components},
\newblock Wiley, New York, 1991.

\bibitem{logGLR06}
D.~Montiel, H.~Cang, and H.~Yang,
\newblock ``Quantitatve characterization of changes in dynamical behavior for
  single-particle tracking studies'',
\newblock {\em J. of Physical Chemistry B}, vol. 110, no. 40, pp. 19763--19770,
  2006.

\bibitem{SiegmundVenkat1996}
D.~Siegmund and E.~S. Venkatraman,
\newblock ``Using the generalized likelihood ratio statistic for sequential
  detection of a change-point'',
\newblock {\em Annals of Stat.}, vol. 23, no. 1, pp. 255 -- 271, 1995.

\bibitem{Siegmund1985}
D.~Siegmund,
\newblock {\em Sequential Analysis: Test and Confidence Intervals},
\newblock Springer, Aug. 1985.

\bibitem{SiegmundYakir2007}
D.~Siegmund and B.~Yakir,
\newblock {\em The statistics of gene mapping},
\newblock Springer, 2007.

\bibitem{cao2015multi}
Yang Cao, Yao Xie, and Nagi Gebraeel,
\newblock ``Multi-sensor slope change detection'',
\newblock {\em Annals of Operations Research}, pp. 1--27, 2016.

\bibitem{RIPCandes2008}
E.~J. Candes,
\newblock ``The restricted isometry property and its implications for
  compressed sensing'',
\newblock {\em Compte Rendus de l'Academie des Sciences, Paris, Serie I}, vol.
  342, pp. 589--592, 2008.

\bibitem{ruben79}
Harold Ruben,
\newblock ``The volume of an isotropic random parallelotope'',
\newblock {\em Journal of Applied Probability}, vol. 16, no. 1, pp. 84--94,
  1979.

\bibitem{betaDistribution90}
P.~Frankl and H.~Maehara,
\newblock ``Some geometric applications of the beta distribution'',
\newblock {\em Ann. Inst. Statist. Math}, 1990.

\bibitem{XuHassibi2007}
W.~Xu and B.~Hassibi,
\newblock ``Efficient compressive sensing with deterministic guarantees using
  expander graphs'',
\newblock in {\em Info. Theory Workshop}, 2007.

\bibitem{BM01}
D.~Burshtein and G.~Miller,
\newblock ``Expander graph arguments for message-passing algorithms'',
\newblock {\em IEEE Trans. Inf. Theory}, vol. 47, no. 2, pp. 782 --790, 2001.

\bibitem{woodall1985multivariate}
William~H Woodall and Matoteng~M Ncube,
\newblock ``Multivariate cusum quality-control procedures'',
\newblock {\em Technometrics}, vol. 27, no. 3, pp. 285--292, 1985.

\bibitem{watts1998collective}
Duncan~J Watts and Steven~H Strogatz,
\newblock ``Collective dynamics of `small-world' networks'',
\newblock {\em Nature}, vol. 393, no. 6684, pp. 440--442, 1998.

\bibitem{PowerSystemBook04}
A.~Abur and A.~G. Exposito,
\newblock {\em Power system state estimation: Theory and Implementation},
\newblock CRC Press, 2004.

\bibitem{FanLiaoLiu2016}
J.~Fan, Y.~Liao, and H.~Liu,
\newblock ``An overview on the estimation of large covariance and precision
  matrices'',
\newblock {\em Econometrics Journal}, , no. 19, pp. C1--C32, 2016.

\bibitem{handbook}
M.~Abramowitz and I.~Stegun,
\newblock {\em A handbook of mathematical functions, with formulas, graphs and
  mathematical tables},
\newblock Dover, 10th edition, 1964.

\bibitem{mythesis}
A.~Thompson,
\newblock {\em Quantitative analysis of algorithms for compressed signal
  recovery},
\newblock PhD thesis, School of Mathematics, University of Edinburgh, 2012.

\end{thebibliography}
