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

Online change-point detection is a classic and fundamental problem in statistics and related areas. The goal is to detect a change in the underlying data distribution as soon as possible once the change has occurred while satisfying the false alarm constraint.
Classic methods are typically based on parametric models which detect changes in distribution parameters, such as Shewhart chart \citep{shewhart1925application},
Cumulative Sum (CUSUM) \citep{page1954continuous} and
Shiryaev-Roberts (SR) \citep{shiryaev1963optimum}.
For comprehensive and systematic surveys on parametric change-point detection, we refer readers to \citet{basseville1993detection,siegmund2013sequential,tartakovsky2014sequential,xie2021sequential}.

In modern applications, especially in high-dimensional data settings, the exact specification of the probability distribution is difficult, where non-parametric and distribution-free methods are needed. Kernel-based methods have attracted much attention due to their flexibility and broad applicability; see, e.g., kernel Maximum Mean Discrepancy (MMD) \citep{gretton2006kernel,smola2007hilbert,gretton2012kernel} and a survey on its application\citep{muandet2017kernel}. As non-parametric methods, the kernel-based methods do not require distributional assumptions and allow flexible kernel choices.

In kernel-based change-point detection, much emphasis has been on offline change-point detection, where the goal is to detect and locate change-points in retrospect. In the offline setting, the detection is typically conducted by a series of two-sample tests for samples divided by candidate change-point locations. For instance, such works include the maximum Kernel Fisher Discriminant Ratio-based procedure \citep{harchaoui2008kernel}, kernel MMD-based scan $B$ procedure \citep{li2019scan}, and a recent kernel-based scan procedure to boost the detection power in high dimensional setting \citep{song2022new}. In contrast, online change-point detection aims to detect the change as soon as possible from sequential/streaming data.
Existing online kernel-based change-point detection has been focusing on the Shewhart chart-type procedure (e.g., \citet{li2019scan}).
The Shewhart chart-type procedure differs from the online CUSUM procedure, which is based on CUSUM recursion and is generally believed to be more sensitive to detecting small changes \citep{xie2021sequential}. The CUSUM-type procedure has yet to be developed for kernel-based methods due to the increased computational and memory cost with time horizon. For existing kernel-based online change-point detection methods, one can see a summary in Table~\ref{table:kernelcpd} and our literature survey in Section~\ref{subsec:related_work}.

\begin{table}[!ht]
\begin{center}
\caption{A summary of existing online kernel change-point detection procedures.}\label{table:kernelcpd}
\vspace{0.1in}
\begin{small}
\resizebox{.85\textwidth}{!}{%
    \begin{tabular}{|c|c|}
\hline
  CUSUM-type procedure & Shewhart chart-type procedure \\
\hline
 \multirow{3}{*}{\citet{flynn2019change}}  &  \citet{huang2014high,li2019scan} \\
   &  \citet{bouchikhi2019kernel,chang2019kernel} \\
    &  \citet{cobb2022sequential} \\
\hline
\end{tabular}
}
\end{small}
\end{center}
\end{table}

In this work, we develop a kernel-based CUSUM procedure for online change-point detection called the {\it Memory Efficient Online Kernel CUSUM}. The proposed CUSUM-type procedure considers a set of possible change-point locations, and for each hypothesized change-point location, a kernel test statistic is formed by comparing the supposed post-change samples (from the change-point location to the current time) with the reference (pre-change) samples. Then the detection statistic takes the maximum over such a set of parallel detection statistics. Compared to the existing Shewhart chart-type of kernel-based online change-point detection procedure such as \citet{li2019scan}, our procedure combines kernel detection statistics with different block sizes, whereas \citet{li2019scan} only considered one fixed block size. We show that this seemingly minor but essential change leads to a significant boost in detection power. Indeed, whether or not searching for the unknown change-point location is the main difference between the Shewhart chart and the CUSUM procedures and leads to performance gain in detecting small changes in the CUSUM-type procedure. While harvesting the performance gain, we also make the proposed CUSUM-type procedure computationally more efficient by presenting a recursive implementation to achieve a constant computational and memory cost at each time step,  which is essential in the online setting.

Our proposed Memory Efficient Online Kernel CUSUM procedure is theoretically grounded. We derive accurate analytic approximations to two fundamental performance metrics: the Average Run Length (ARL), which is useful for calibrating the procedure in practice to control false alarms, and the Expected Detection Delay (EDD), which is related to detection power. Based on this, we further establish the optimal window size to be $\cO(\log {\rm ARL})$ to achieve (asymptotic) no power loss, which is in the same order as the classic result of the window-limited Generalized Likelihood Ratio (GLR) procedure for parametric models \citep{lai1995sequential}. Since window size captures how much data is needed in the memory, this result indicates the minimal computation and memory cost necessary for online kernel change-point detection.

The remainder of this paper is organized as follows. We first briefly review related works.
Section~\ref{sec:formulation} introduces the online change-point detection problem and reviews the literature.
Section~\ref{sec:detect_stat} presents our proposed Memory Efficient Online Kernel CUSUM procedure.
Section~\ref{sec:theory} presents analytic approximations to the ARL and EDD and characterizes the optimal window size.
Section~\ref{sec:numerical} validates our theoretical findings and compares our proposed method with existing methods through extensive numerical experiments.
Section~\ref{sec:conclusion} concludes the paper.

\subsection{Literature}\label{subsec:related_work}

\vspace{.1in}

\noindent \textit{Kernel methods in change-point detection.}
An early work on kernel change-point detection is \citet{harchaoui2007retrospective}, where they incorporated kernels into dynamic programming and developed an offline multiple change-point estimation method. This method was extended by \citet{arlot2019kernel} to handle unknown number of change-points via a model selection penalty; however, its performance heavily depends on the penalty hyperparameter and it is difficult to control the type I error in practice.
Later on, \citet{garreau2018consistent} performed non-asymptotic analysis to show the multiple change-point estimation by model selection \citep{arlot2019kernel} can estimate the correct number of change-points with high probability and locate them at the optimal rate.
The problem of calibrating detection procedures also exists for maximum Kernel Fisher Discriminant Ratio based approach \citep{harchaoui2008kernel}, where they applied the bootstrap resampling method to compute the detection threshold.
In addition, \citet{bouchikhi2019kernel,ferrari2023online} considered a kernel density ratio (KDR) \citep{huang2006correcting,nguyen2010estimating} based Shewhart chart-type method for online change-point detection.

\vspace{.1in}

\noindent \textit{Kernel MMD and its application in change-point detection.}
\citet{gretton2006kernel,smola2007hilbert,gretton2012kernel} proposed to measure the distance/discrepancy between distributions via the Reproducing Kernel Hilbert space embedding of distributions and developed the well-known Maximum Mean Discrepancy (MMD) statistic. Since then, such kernel MMD-based statistical tests have become popular due to their wide range of applicability to arbitrary domains, ranging from
anomaly detection \citep{zou2014nonparametric}
model criticism and goodness-of-fit test \citep{lloyd2015statistical} and the optimality of such kernel-based tests \citep{li2019optimality,balasubramanian2021optimality},
approximate bayesian computation \citep{park2016k2},
transfer learning \citep{van2018transfer,yang2019polynomial},
distributionally robust optimization \citep{staib2019distributionally},
robust hypothesis testing \citep{sun2021data},
optimal sub-sampling \citep{dwivedi2021kernel,dwivedi2021generalized} as well as change-point detection.

In offline kernel MMD-based change-point detection, there are many other notable contributions apart from aforementioned work such as \citet{li2019scan}, including \citet{truong2019greedy}, who developed a greedy single change-point estimation procedure, \citet{jones2020kernel}, who developed a multiple change-point estimation method by incorporating MMD in the least-squares framework developed by \citet{harchaoui2007retrospective} and leveraging Nyström kernel approximation \citep{williams2001using} to achieve linear complexity, and so on.

Kernel MMD-based online change-point detection procedures typically test for changes between a window of reference data and a window of streaming data.
Generally speaking, there are three types of procedures in online change-point detection, i.e., Shewhart chart, CUSUM, and SR \citep{xie2021sequential}.
Recently, there has been much effort focusing on kernel MMD-based online Shewhart chart-type procedures, which consider fixed-window scan statistics.
In addition to the Scan $B$ procedure \citep{li2019scan} introduced above, notable contributions include \citet{huang2014high}, who first introduced MMD to online change-point detection and applied bootstrap to calibrate the detection;
\citet{chang2019kernel}, who leveraged neural network-based data-driven kernel selection \citep{gretton2012optimal} to boost detection power; \citet{cobb2022sequential}, who considered a memoryless geometric distribution of the stopping time to control the false alarm rate and calibrated the detection via bootstrap.

However, the fixed-length scanning window of Shewhart chart-type procedures fails to account for the change-point location, which may decrease the detection power since the window of streaming data could contain samples from pre-change (in-control) distribution. \citet{cobb2022sequential} indeed tackled this problem by considering a time-varying threshold, but such approach was a heuristic and did not have strong theoretical guarantee.
A much more sensible approach is CUSUM-type procedure, but such procedures are rarely considered in kernel MMD online change-point detection, potentially due to its higher computational and memory cost compared to Shewhart chart-type procedures.
One notable contribution on this aspect is a recent work by \citet{flynn2019change}, who proposed a Kernel CUSUM procedure by replacing the likelihood ratio in classic CUSUM with the linear-time MMD statistic and derived bounds for both ARL and EDD. Unfortunately, such ARL bound fails to help calibrate the detection in practice. In this work, we fill this void by presenting a kernel MMD-based CUSUM-type procedure, explicitly characterizing its important statistical properties, i.e., ARL and EDD, and proving its optimal window size to give a memory efficient procedure yet maintain the detection power.

\vspace{.1in}

\noindent \textit{Other non-parametric change-point detection methods.}
In addition to the aforementioned kernel-based approaches, there is also rich literature on other non-parametric change-point detection procedures.
Notable contributions include
Mann–Whitney-type statistic \citep{hawkins2010nonparametric},
total variation regularization to estimate a piecewise constant structure, where the jumps correspond to multiple change-points \citep{harchaoui2010multiple},
direct density-ratio estimation \citep{kawahara2012sequential},
K Nearest Neighbor based CUSUM-type statistic \citep{chen2019sequential},
Kolmogorov–Smirnov distance \citep{padilla2021optimal} and so on.
Besides, in addition to classic independent and identically distributed (i.i.d.) setting, it is also of interest to handle data with serial correlation in change-point detection.
For example, there is a line of research considering the total variation regularization to estimate the changes in high-dimensional regression parameters \citep{rinaldo2021localizing}, autoregressive parameters in time series \citep{wang2019localizing}, and piecewise constant dynamic background in autoregressive time series \citep{wei2021inferring}.
Another example is \citet{liu2013change}, which extended the direct density-ratio estimation approach to the time series data.
However, how to handle data with serial correlation, e.g., autoregressive time series, is outside the scope of this work and we leave this topic for further study.

\section{Set-up and Background}\label{sec:formulation}

Consider sufficient reference samples from domain $\cX$ (typically $\cX = \RR^d$) with an unknown pre-change distribution $\PP$:
$$X_1,\dots,X_M \overset{i.i.d.}{\sim} \PP.$$
At time $t$, we receive i.i.d. streaming data $Y_1,\dots,Y_t$. The online change-point detection problem tests the null hypothesis
$$H_0: Y_1,\dots,Y_t \overset{i.i.d.}{\sim} \PP,$$ against the alternative hypothesis
$$H_1: \exists \  r<t, \  Y_1,\dots,Y_r \overset{i.i.d.}{\sim} \PP, \  Y_{r+1},\dots,Y_t \overset{i.i.d.}{\sim} \QQ,$$
for some unknown pre- and post-change distributions $\QQ \not= \PP$.

In online change-point detection, the goal is to detect the change as soon as possible subject to the Average Run Length constraint when there is no change. Typically, we control the ARL by properly choosing the detection threshold $b$, and we study the Expected Detection Delay under such ARL constraint. For the same ARL, the smaller the EDD is, the more sensitive/better the detection procedure is.

\vspace{0.1in}

\noindent
\textit{Notations.}
Our notations are conventional. Under the null hypothesis where there is no change, we use $\mathbb{P}_{\infty}$ and $\mathbb{E}_{\infty}$ to denote the probability and expectation in this case. When there is a change from $\PP$ to $\QQ$ at time $t = r$, all the reference samples and the streaming data before the change-point, i.e., $t < r$, follow the pre-change distribution $\PP$ whereas the rest of the streaming data will follow the post-change distribution $\QQ$. We use $\mathbb{P}_{r}$ and $\mathbb{E}_{r}$ to denote the probability and expectation in this case.
In particular, we use $\mathbb{P}_{0}$ and $\mathbb{E}_{0}$ to denote the probability and expectation when all streaming samples follow post-change distribution $\QQ$.
We denote $\mathbf{1}_d = (1,\dots,1)^{ \mathrm{\scriptscriptstyle T} } \in \RR^d$, $\mathbf{0}_d = (0,\dots,0)^{ \mathrm{\scriptscriptstyle T} } \in \RR^d$ and $I_d \in \RR^{d \times d}$ as the identity matrix, where superscript $^{ \mathrm{\scriptscriptstyle T} }$ standards for the vector/matrix transpose.
For integers $0 < m \leq n$, we denote $[m:n] = \{m,\dots,n\}$; in a special case where $m=1$, we denote $[n] = \{1,\dots,n\}$.
In addition, we denote $a \wedge b = \min\{a,b\}$ and $a \vee b = \max\{a,b\}$. In particular, we denote $(a)^+ = \max\{a,0\}$.
Lastly, for asymptotic notations, we use the standard definitions as follows:
\begin{itemize}
    \item $f(m) = o(g(m))$ means for all $c > 0$ there exists $k > 0$ such that $0 \leq f(m) < cg(m)$ for all $m \geq k$;
    \item $f(m) =  \cO(g(m))$ means there exist positive constants $c$ and $k$, such that $0 \leq f(m) \leq cg(m)$ for all $m \geq k$;
    \item $f(m) \sim g(m)$ means there exist positive constants $c$ and $k$, such that $0 \leq cg(m) \leq f(m)$ for all $m \geq k$ and at the same time $f(m) =  \cO(g(m))$.
\end{itemize}

\subsection{Classic parametric procedures}
Let us first introduce three popular types of parametric procedures for online change-point detection: Shewhart chart, Cumulative Sum, and Generalized Likelihood Ratio.

\vspace{.1in}

\noindent \textit{Shewhart Chart.} One of the earliest online change-point detection procedures is the Shewhart chart \citep{shewhart1925application}, where the pre-change and post-change distributions need to be completely specified. Suppose the pre-change distribution $\PP$ has density function $p(\cdot)$ and the post-change distribution $\QQ$ has density function $q(\cdot)$. At time step $t$, the Shewhart chart considers the log-likelihood ratio on the most recent steaming observation $Y_t$ as the detection statistic. Formally, the Shewhart chart detection procedure is defined as
\begin{equation}\label{eq:T_Sh}
    T_{\rm Sh} = \inf \left\{t>0: \log \frac{q(Y_t)}{p(Y_t)} \geq b\right\},
\end{equation}
where $b>0$ is a pre-selected constant detection threshold to meet the ARL constraint. However, the Shewhart chart ignores the past streaming observations, and this ``information loss'' could lead to power loss \citep{xie2021sequential}.

\vspace{.1in}

\noindent \textit{Cumulative Sum.}
CUSUM procedure was first introduced by \citet{page1954continuous}, where again the densities of pre-change distribution $\PP$ and the post-change distribution $\QQ$, i.e., $p(\cdot)$ and $q(\cdot)$, are assumed completely known. At time step $t$, observing steaming data $Y_1,\dots,Y_t$, CUSUM considers the cumulative log-likelihood ratio as the detection statistic, i.e.,
\begin{equation}\label{eq:cusum_original_form}
\begin{split}
    S_t &= \sum_{s=1}^t \log \frac{q(Y_s)}{p(Y_s)} - \min_{0 \leq r \leq t} \sum_{s=1}^r \log \frac{q(Y_s)}{p(Y_s)} \\
    &= \max_{1 \leq r \leq t+1} \sum_{s=r}^{t} \log \frac{q(Y_s)}{p(Y_s)}.
\end{split}
\end{equation}
By considering the complete information, the CUSUM procedure tackles the ``information loss'' issue of the Shewhart chart.

Since the maximum in \eqref{eq:cusum_original_form} is attained at the maximum likelihood estimator of the unknown change-point location, the CUSUM statistic has the following recursive computation rule:
\begin{equation}\label{eq:cusum}
    S_t = \left(S_{t-1} + \log \frac{q(Y_t)}{p(Y_t)}\right)^+ , \quad S_0 = 0.
\end{equation}
The CUSUM procedure for threshold $b$ is then defined by the following stopping time:
\begin{equation}\label{eq:T_cusum}
    T_{\rm CUSUM} = \inf \{t>0: S_t \geq b\}.
\end{equation}
Although CUSUM was first introduced as a heuristic, \citet{moustakides1986optimal} later showed that CUSUM minimizes the worst-case expected detection delay under the ARL constraint.
As the ARL goes to infinity, the CUSUM procedure \eqref{eq:T_cusum} has
\begin{equation}\label{eq:lorden_lower_bound}
    \mathbb{E}_{0}[T_{\rm CUSUM}] \sim \frac{\log \text{\rm ARL}}{\text{\rm KL}(\mathbb{P}\|\mathbb{Q})},
\end{equation}
where $\text{\rm KL}(\mathbb{P}\|\mathbb{Q}) = \mathbb{E}_{0}[\log ({q(Y_1)}/{p(Y_1)})]$ is the Kullback-Leibler (KL) divergence. This is known as the Lorden's Lower Bound \citep{lorden1971procedures}.

\vspace{.1in}
\noindent \textit{Generalized Likelihood Ratio.} GLR procedure \citep{lorden1971procedures} is another popular parametric approach. Similarly, the pre-change density $p(\cdot)$ is assumed known, whereas the post-change density $q(\cdot;\theta)$ has a known form but unknown parameter $\theta$, which takes values from $\Theta$.
The GLR detection procedure is defined as follows:
\begin{equation}\label{eq:GLR}
    T_{\rm GLR} = \inf \left\{t>0: \max_{1\leq r \leq t-1} \sup_{\theta \in \Theta}\sum_{s= r }^{t} \log \frac{q(Y_s;\theta)}{p(Y_s)} \geq b\right\}.
\end{equation}

The detection statistic used by the GLR procedure involves maximum likelihood estimation of both change-point location $r$, which is similar to CUSUM, and density parameter $\theta$, which accounts for the unknown post-change distribution parameter. Unfortunately, the GLR statistic \eqref{eq:GLR} cannot be updated recursively as CUSUM \eqref{eq:cusum} does. To reduce the computational cost, the following {\it window-limited} GLR (W-GLR)
is adopted:
\begin{equation}\label{eq:GLR_window_limited}
    T_{\rm W-GLR} = \inf \left\{t>0: \max_{t-w\leq r \leq t-1} \sup_{\theta \in \Theta}\sum_{s= r }^{t} \log \frac{q(Y_s;\theta)}{p(Y_s)} \geq b\right\},
\end{equation}
which only considers the potential change-point in a window with size $w$.

The choice of window size $w$ is the key to the success of the window-limited GLR procedure. \citet{lai1995sequential} proved that the optimal window size $w$ under the Gaussian mean shift assumption is
\begin{equation}\label{eq:lai_opt_window_size}
    w \sim \log \text{\rm ARL}.
\end{equation}
According to \citet{lai2001sequential}, by choosing such window size, ``there is little loss of efficiency in reducing the computational complexity by the window-limited modification''.

\subsection{Scan $B$ procedure}
In this subsection, we will introduce the Scan $B$ procedure, which is a kernel change-point detection procedure built on the empirical estimator of MMD, which we call {\it MMD statistic}.
We will omit the formal definition of Reproducing Kernel Hilbert space as well as the MMD, which we call {\it population MMD}, and refer readers to \citet{gretton2006kernel,smola2007hilbert,gretton2012kernel} for their formal definitions. Now, let us begin with introducing the MMD statistic.

\vspace{.1in}
\noindent \textit{Maximum Mean Discrepancy statistic.} For two probability distributions $\PP$ and $\QQ$ supported on domain $\cX$, the population MMD, which we denote as $\text{\rm MMD}^2(\PP,\QQ)$, is a non-parametric metric to measure their distance/discrepancy. Given two sets of samples $\mathbf X$ and $\mathbf Y$ with equal sample size $n$, where $\mathbf X = (X_1,\dots,X_n)$ are i.i.d. samples from $\PP$ and
$\mathbf Y = (Y_1,\dots,Y_n)$ are i.i.d. samples from $\QQ$, an unbiased estimator of MMD, i.e., MMD statistic, is given by using U-statistic as follows \citep{gretton2012kernel}:
\begin{equation}\label{eq:MMD_estimator}
    \text{\rm MMD}_{\rm u}^2(\mathbf X, \mathbf Y)= \frac{1}{n(n-1)} \sum_{i=1}^{n} \sum_{j \neq i}^{n} h\left(X_{i}, X_{j}, Y_{i}, Y_{j}\right),
\end{equation}
where
\if\mycolumn1
 \begin{equation}\label{eq:kernel_h}
    h(x_{1}, x_{2}, y_{1}, y_{2})=k\left(x_{1}, x_{2}\right)+k\left(y_{1}, y_{2}\right)-k\left(x_{1}, y_{2}\right)-k\left(x_{2}, y_{1}\right),
 \end{equation}
\else
 \begin{equation}\label{eq:kernel_h}
 \begin{split}
    h&(x_{i}, x_{j}, y_{i}, y_{j})\\
    &=k\left(x_{i}, x_{j}\right)+k\left(y_{i}, y_{j}\right)-k\left(x_{i}, y_{j}\right)-k\left(x_{j}, y_{i}\right),
 \end{split}
 \end{equation}
\fi
and $k(\cdot,\cdot)$ is the user-specified kernel function:
\begin{equation}\label{eq:kernel}
    k(\cdot,\cdot): \cX \times \cX \rightarrow \RR.
\end{equation}
Commonly used kernel functions include Laplace kernel $k(x,y) = \exp \{-{\norm{x-y}}/{\gamma}\},$ and Gaussian radial basis function (RBF) kernel $k(x,y) = \exp \{-{\norm{x-y}^2}/{\gamma^2}\},$ where $\gamma$ is the bandwidth parameter. As one will see later in the next section, our theoretical analysis is general and allows flexible kernel choices.

\vspace{.1in}
\noindent \textit{Scan $B$ procedure.}
One drawback of the above unbiased MMD estimator \eqref{eq:MMD_estimator} comes from its quadratic computational cost. This is undesirable in the online setting, where we need to evaluate the detection statistic in real time as we receive the streaming data. To fix this issue, \citet{li2019scan} leveraged the $B$-test statistic \citep{zaremba2013b} and developed the Scan $B$ procedure.

In Scan $B$ procedure, they randomly sample $N B$ samples from reference data $X_1,\dots,X_M$ without replacement, in order to build $N$ {\it pre-change blocks} with equal block size $B$, i.e., $\mathbf{X}_B^{(i)}, \ i \in [N]$. In both Scan $B$ procedure and our proposed procedure below, by saying we have sufficient amount of reference samples, we assume reference sample size $M$ is large enough such that $M > N B$. At time $t$, the {\it post-change block} consists of the most recent $B$ streaming data, i.e., $\mathbf{Y}_B(t) = (Y_{t-B+1},\dots, Y_t)$.
They compute unbiased MMD statistics between the post-change block and each pre-change block and take their average as the $B$-test statistic.
By reusing the post-change block $\mathbf{Y}_B(t)$, Scan $B$ procedure achieves overall linear computational complexity with respect to both time $t$ and sample size $N$. Formally, the detection statistic, which they call Scan $B$-statistic, is defined as follows:
\begin{equation}\label{eq:scan-b}
    {Z}_B(t) = \frac{{Z}^\prime_B(t)}{\sqrt{\operatorname{Var}_{\infty}({Z}^\prime_B(t))}},
\end{equation}
where $\operatorname{Var}_{\infty} (\cdot)$ denotes the variance under $H_0$ and
\begin{equation}\label{eq:scan-b_unmoralized}
    {Z}^\prime_B(t) = \frac{1}{N}\sum_{i=1}^N {\text{\rm MMD}_{\rm u}^2(\mathbf{X}_B^{(i)},\mathbf{Y}_B(t))}.
\end{equation}

For detection threshold $b$, the Scan $B$ procedure is defined by the following stopping time:
\begin{equation}\label{eq:stopping_time_scanb}
    T_{{\rm Sc}B} = \inf \{t : {Z}_B(t) > b\}.
\end{equation}

The above normalizing constant, i.e., the standard deviation of ${Z}^\prime_B(t)$ under $H_0$, remains unchanged with increasing $t$. Moreover, it has an an analytic expression (see Lemma~\ref{lma:var_H0} in Appendix~\ref{appendix:moment}) and can be computed using the reference samples before the detection procedure.

\begin{remark}[Effect of the normalizing constant]
Given that $\operatorname{Var}_{\infty}({Z}^\prime_B(t))$ is held constant in the online Scan $B$ procedure and does NOT influence the online detection, one may ask why we include this normalizing constant here. In fact, \citet{li2019scan} introduced this normalizing constant under the offline setting, where they proposed to use window-limited maximization with respect to (w.r.t.) the block size the estimate the change-point location.
To be precise, consider a fixed number of observations $Y_1,\dots,Y_t$, for a hypothesized change-point at $r < t$, the un-normalized statistics ${Z}^\prime_{2}(t),\dots,{Z}^\prime_{t-r+1}(t)$ \eqref{eq:scan-b_unmoralized} all have the same expectation, i.e., the population MMD, making this un-normalized version undesirable in change-point localization. By scaling those statistics with a factor that decreases w.r.t. the block size $B$ (e.g., $\operatorname{Var}_{\infty}^{-1/2}({Z}^\prime_B(t))$) one can expect that the normalized statistic takes maximum at the change-point location. In the next section, we will see the same effect of normalization when we perform window-limited maximization to locate the change-point in online setting. For brevity, we will use the term Scan $B$-statistic (or procedure) for online Scan $B$-statistic (or procedure) for simplicity, as we do not consider offline setting in this work.
\end{remark}

\section{Proposed Detection Procedure} \label{sec:detect_stat}

In this section, we will formally present our proposed Memory Efficient Online Kernel CUSUM procedure for online change-point detection, which includes a recursive computation of the detection statistic. Compared with Scan $B$ procedure, our proposed procedure significantly achieves quicker change detection while maintains the same computational and memory complexity. Now, we begin with presenting our proposed detection procedure.

\subsection{Memory Efficient Online Kernel CUSUM}

Inspired by the original form of CUSUM \eqref{eq:cusum_original_form}, we propose to consider a parallel set of Scan $B$-statistics ${Z}_B(t)$ \eqref{eq:scan-b} with block size $B$ taking values from $[2:w]$ and take their maximum.
For a pre-selected threshold $b$, the {\it Memory Efficient Online Kernel CUSUM} procedure is defined via the following stopping time:
\begin{equation}\label{eq:stopping_time}
    T_b^w = \inf \left\{t : \max_{B \in [2:w]} {Z}_B(t) > b\right\}.
\end{equation}
where $ w $ is the {\it window size} parameter.

\begin{remark}[Comparison with Scan $B$ procedure]
The window-limited maximization coupled with the normalization can remove pre-change samples from the post-change block, which gives our procedure an edge over Scan $B$ procedure both empirically and theoretically --- it can not only result in quicker change detection (see an numerical evidence in Figure~\ref{fig:EDDARL_compare_w_scanB}), but also help us derive an analytic EDD approximation (see Section~\ref{subsec:arledd} later for details).
In contrast, Scan $B$ procedure include samples from $\PP$ in its post-change block, making it harder to be differentiated from the pre-change blocks. Moreover, the portion of pre-change samples in the post-change block varies at each time step, making it difficult to approximate the EDD in closed-form.
Moreover, as we will see in the next subsection, even though the window-limited maximization introduces more computations, the computational and memory complexity remains the same with Scan $B$ procedure and does not grow with time horizon, which is essential in online setting.

\begin{figure}[!htp]
\centerline{
\includegraphics[width = .95\textwidth]{fig/new exp/fig1_illus_power.pdf}
}
\caption{Comparison of Expected Detection Delay for given Average Run Length between our proposed (red) and the Scan $B$ (orange) detection procedures. ARL is in base $10$ logarithm. The pre-change distribution is $\PP = N(\mathbf{0}_{20}, I_{20})$ and the post-change distribution $\QQ$ is specified on the top of each figure (with $\sigma = 2$). We can observe that our proposed procedure achieves quicker change detection under all three settings.}
\label{fig:EDDARL_compare_w_scanB}
\end{figure}
\end{remark}

\subsection{Recursive computation of detection statistics}

\begin{algorithm}[!htp]
 \caption{Memory Efficient Online Kernel CUSUM}\label{algo:kernel_CUSUM}
 \begin{algorithmic}[1]
 \renewcommand{\algorithmicrequire}{\textbf{Input:} }

 \REQUIRE Pre-change blocks' number $N$, window size $w$, ARL, reference samples $X_1,\dots,X_M$.

  \STATE Initialization: use the most recent $w$ observations to form the post-change block $\mathbf{Y}_t \in \cX^w$; randomly sample $N \times w$ reference samples to form $N$ pre-change blocks $\mathbf{X}^{(i)} \in \cX^w, \ i = 1,\dots,N$; select the threshold $b$ to meet the ARL constraint; compute $C$ \eqref{eq:cpq} using reference samples; calculate the Gram matrices $G_{XX}^{(i)} = k(\mathbf{X}^{(i)},\mathbf{X}^{(i)})$, $G_{XY}^{(i)} = k(\mathbf{X}^{(i)},\mathbf{Y}_t), \ i=1,\dots,N,$ and $G_{YY} = k(\mathbf{Y}_t,\mathbf{Y}_t)$
 \\
  \WHILE {True}
  \STATE receive sample $Y_t$
  \STATE update post-change block:
      $$\mathbf{Y}_t (1:w-1) \leftarrow \mathbf{Y}_t (2:w), \quad
      \mathbf{Y}_t (w) \leftarrow Y_t$$
      \vspace{-0.2in}
  \STATE update Gram matrix $G_{YY}$:
  \vspace{-0.1in}
  \begin{align*}
      &G_{YY} (1:w-1, 1:w-1)  \leftarrow G_{YY} (2:w, 2:w), \\
      &G_{YY} (:, w)  \leftarrow k (\mathbf{Y}_t, Y_t), \quad G_{YY} (w, :)  \leftarrow G_{YY}^{ \mathrm{\scriptscriptstyle T} } (:, w)
  \end{align*}
  \vspace{-0.2in}
  \FOR {$i = 1,\dots,N$}
  \STATE update Gram matrix $G_{XY}^{(i)}$:
  \vspace{-0.1in}
  \begin{align*}
      G_{XY}^{(i)} (:, 1:w-1) & \leftarrow G_{XY}^{(i)} (:, 2:w), \\
      G_{XY}^{(i)} (:, w) & \leftarrow k (\mathbf{X}^{(i)}, Y_t)
  \end{align*}
  \ENDFOR
  \STATE calculate $z = \sum_{i=1}^N \sum_{j=1}^2 \sum_{k \not= j}^2 G_{XX}^{(i)} (j, k) + G_{YY} (j, k) - G_{XY}^{(i)} (j, k) - G_{XY}^{(i)} (k, j)$ and $Z_2(t) = C z / (\sqrt{2} N)$; let $Z_t = Z_2(t)$
  \FOR {$B = 3, \dots, w$}
  \STATE update
  \vspace{-0.1in}
  \if\mycolumn1
\begin{align*}
      z \leftarrow z + \sum_{i=1}^N  \sum_{j = 1}^B & G_{XX}^{(i)} (B, k) + G_{YY} (B, k) -  G_{XY}^{(i)} (B, k) - G_{XY}^{(i)} (k, B)
  \end{align*}
\else
\begin{align*}
      z \leftarrow z + \sum_{i=1}^N  \sum_{j = 1}^B & G_{XX}^{(i)} (B, k) + G_{YY} (B, k) \\
       &-  G_{XY}^{(i)} (B, k) - G_{XY}^{(i)} (k, B)
  \end{align*}
\fi
  \vspace{-0.1in}
  \STATE calculate $Z_B(t) = {C z}/(\sqrt{B(B-1)}N)$
  \STATE update $Z_t \leftarrow \max\{Z_t, Z_B(t)\}$
  \ENDFOR
  \IF {$Z_t > b$}
  \STATE raise an alarm and set $T_b^w = t$
  \STATE \textbf{break}
  \ENDIF
  \ENDWHILE

 \end{algorithmic}
 \end{algorithm}

Before presenting the algorithm, we introduce some notations which only appear in this subsection: for a matrix $G \in \RR^{m \times n}$, $G(i_1:i_2,j_1:j_2)$ represents a slice of matrix $G$ containing its $i_1$-th to $i_2$-th rows and $j_1$-th to $j_2$-th columns; in particular, when we want to include all rows (or columns) in this slice, we use ``$:$'' instead of ``$1:m$'' (or ``$1:n$'').
Similarly, for a vector $g \in \RR^{m}$, $g(i:j)$ represents a slice of $g$, containing its $i$-th to $j$-th elements. Now, we give step-by-step instructions on how to preform our proposed detection procedure in Algorithm~\ref{algo:kernel_CUSUM}.

In Algorithm~\ref{algo:kernel_CUSUM}, the recursive computation of detection statistics comes from two parts: recursive update of Gram matrices (lines $5 - 8$; see a graphical illustration in Figure~\ref{fig:online_update_illus}) and a recursive formulation to find the window-limited maximum (lines $9 - 14$), which are direct benefits of CUSUM and closely resemble that of the original CUSUM \eqref{eq:cusum}.

\begin{figure}[!htp]
\centerline{
\includegraphics[width = 0.9\textwidth]{fig/illus/fig0_illus_onlineupdate.pdf}
}
\caption{Graphical illustration of the recursive update of the Gram matrix in lines $5 - 8$ in Algorithm~\ref{algo:kernel_CUSUM}. Such recursive update was first considered in \citet{li2019scan}. The dark gray part corresponds to the Gram matrix used at the corresponding time step. The highlighted Gram matrix is readily obtained from last time step, which results in $ \cO(N  w )$ complexity for the update.}
\label{fig:online_update_illus}
\end{figure}

Even though the online update of Gram matrices only takes $ \cO(Nw)$ computational complexity, the whole procedure still takes $ \cO(N w^2)$ computations to calculate the detection statistic at each time step (lines $9 - 14$), which is the same with Scan $B$ procedure.
In addition, at each time step, our detection procedure only requires additional $ \cO(w)$ computations in lines $12 - 13$ and does not require additional memory to store those intermediate calculations.
We summarize the comparison among Shewhart chart CUSUM, W-GLR, Scan $B$, our proposed Memory Efficient Online Kernel CUSUM and its oracle variant \eqref{eq:stopping_time_noWL} in Table~\ref{table:complexity}. In particular, we report (i) whether or not the method is robust to the post-change distribution specification, and the (ii) computation and (iii) memory complexity at each time step $t$ for all six methods.
As we can see, the additional window-limit maximization in our proposed Memory Efficient Online Kernel CUSUM does NOT increase the computational and memory complexity compared with the Scan $B$ procedure. We need to clarify that Algorithm~\ref{algo:kernel_CUSUM} actually requires $\cO(Nw^2)$ memory to store the Gram matrix. However, in theory, we can only use $\cO(Nw)$ memory to store raw observations and re-compute the whole Gran matrix with $\cO(Nw^2)$ computations at each time step.

\begin{table}[!ht]
\begin{center}
\caption{Comparison of online change-point detection methods. Here, $w$ is the window size parameter and $N$ is the number of pre-change blocks.}\label{table:complexity}
\vspace{0.1in}
\begin{small}
\resizebox{\textwidth}{!}{%
    \begin{tabular}{|l|c|c|c|}
\hline
 ~ & Robustness & Computation & Memory  \\ \hline
 Shewhart chart \eqref{eq:T_Sh} & \multirow{2}{*}{$\PP$ and $\QQ$ are both known} & $ \cO(1)$ & $-$  \\ \cline{1-1} \cline{3-4}
 CUSUM \eqref{eq:T_cusum} &  & $ \cO(1)$ & $ \cO(1)$  \\ \hline
 W-GLR \eqref{eq:GLR_window_limited} & $\PP$ is known and $\QQ$ has known density form & $ \cO(w^2)$ & $ \cO(w)$  \\ \hline
 Scan $B$ \eqref{eq:stopping_time_scanb} & \multirow{3}{*}{\shortstack{$\PP$ and $\QQ$ are both unknown, but sufficient \\ amount of reference samples is available} } &  $ \cO(Nw^2)$ & $ \cO(Nw)$  \\ \cline{1-1} \cline{3-4}
 Proposed \eqref{eq:stopping_time} &   & $ \cO(Nw^2)$ & $ \cO(Nw)$  \\ \cline{1-1} \cline{3-4}
 Oracle \eqref{eq:stopping_time_noWL} &  & $ \cO(Nt^2)$ & $ \cO(Nt)$  \\ \hline
\end{tabular}
}
\end{small}
\end{center}
\end{table}

\section{Theoretical Analysis} \label{sec:theory}

In this section, we will derive analytic approximations to two fundamental performance metrics in online change-point detection:
(i) the Average Run Length $\mathbb{E}_{\infty}[T_b^w]$, which is the expected stopping time when there is no change, and
(ii) the Expected Detection Delay $\mathbb{E}_{0}[T_b^w]$, which is the expected stopping time when the change occurs immediately at $r = 0$.
We will numerically show that our ARL approximation is very accurate, which helps us to calibrate the detection procedure without time-consuming Monte Carlo simulation.
Furthermore, these results help us to establish the asymptotic optimal window size $ w = w^\star $ in exactly the same sense with \citet{lai1995sequential}.
Now, let us start with the necessary assumptions for our theoretical analysis.

\subsection{Assumptions}
We consider the probability measure $\PP$ (i.e., the pre-change distribution) on a measurable space $(\cX,\cB)$. For symmetric, positive (semi-)definite, and square integrable kernel function $k(\cdot,\cdot)$ \eqref{eq:kernel}. For such kernels, Mercer's theorem gives the following decomposition:
\begin{equation*}
    k(x,x^\prime) = \sum_{j = 1}^\infty \lambda_j \varphi_j(x) \varphi_j(x^\prime),
\end{equation*}
where the limit is in $L^2(\PP)$, $\lambda_1 \geq \lambda_2 \geq \dots \geq 0$ are eigenvalues of
the integral operator induced by kernel $k$, and $\{\varphi_j(\cdot): j \geq 1\}$ are the corresponding orthonormal eigenfunctions such that
\begin{equation*}
    \int_{\cX} \varphi_j(x) \varphi_{j^\prime} (x) d \PP(x) = \left\{\begin{array}{ll}
1 & \text{ if } j = j^\prime, \\
0 & \text{ otherwise}.
\end{array}\right.
\end{equation*}
For concreteness, we further assume kernel $k$ is universal in that $\{\varphi_j(\cdot): j \geq 1\}$ forms an orthonormal basis of $L^2(\PP)$.

Consider the probability measure $\QQ$ (i.e., the post-change distribution) dominated by $\PP$ in that the Radon-Nikodym derivative $d \QQ/ d \PP$ is well-defined. The (squared) population MMD, denoted by $\text{\rm MMD}^2(\PP,\QQ)$, can be expressed as follows:
\begin{align}
    \text{\rm MMD}^2(\PP,\QQ) &= \int_{\cX} \int_{\cX} k(x,x^\prime) d (\PP - \QQ) (x) d (\PP - \QQ) (x^\prime) \nonumber \\
    & = \sum_{j = 1}^\infty \lambda_j \left(\int_{\cX} \varphi_j(x) (p(x) - q(x)) dx \right)^2, \label{eq:eigen_decomp}
\end{align}
where $p$ and $q$ are the densities of $\PP$ and $\QQ$, respectively. To ensure the departure $p - q$ is detectable by a given kernel $k$, we assume there exists $\lambda_j > 0$ such that
\begin{equation}\label{assumption:detection}\tag*{(A1)}
    a_j = \int_{\cX} \varphi_j(x) (p(x) - q(x)) dx \not= 0.
\end{equation}
Here, $a_j$ represents the ``projection'' of the departure $p - q$ along the ``direction'' $\varphi_j$, and Assumption~\ref{assumption:detection} guarantees that
$$\text{\rm MMD}^2(\PP,\QQ) = \sum_{j = 1}^\infty \lambda_j a_j^2 > 0.$$

In addition to the above {\it detectability} assumption, we further impose a technical assumption that all eigenfunctions are uniformly bounded, i.e.,
$$\sup_{j \geq 1} \norm{\varphi_j}_{\infty} < \infty,$$
where $\norm{\varphi_j}_{\infty} = \sup_{x \in \cX} |\varphi_j(x)|$.

This also ensures that the kernel $k(\cdot,\cdot)$ is uniformly bounded on domain $\cX \times \cX$, i.e., there exists a constant $K > 0$ such that
\begin{equation}
    | k(x,y) | \leq K, \quad \forall x, y \in \cX. \label{assumption:bounded_kernel}\tag*{(A2)}
\end{equation}
As one can see, Assumption~\ref{assumption:bounded_kernel} guarantees the existence of the moments of $k(X,Y)$, where $X$ and $Y$ can follow either $\PP$ or $\QQ$.

\begin{remark}
The biased MMD statistic, different from the unbiased estimator $\text{\rm MMD}_{\rm u}^2(\mathbf X, \mathbf Y)$ \eqref{eq:MMD_estimator}, can be expressed by replacing $\PP$ and $\QQ$ with their empirical versions $\hat \PP_n$ and $\hat \QQ_n$ in \eqref{eq:eigen_decomp}, respectively, i.e.,
\begin{align}
    \text{\rm MMD}_{\rm b}^2(\mathbf X, \mathbf Y)& = \int \int k(x,x^\prime) d (\hat \PP_n - \hat \QQ_n) (x) d (\hat \PP_n - \hat \QQ_n) (x^\prime) \nonumber \\
    & = \sum_{j = 1}^\infty \lambda_j \left(\frac{1}{n} \sum_{i = 1}^n \varphi_j(X_i) - \varphi_j(Y_i) \right)^2. \nonumber
\end{align}
Thus, in practice the detectability depends on the empirical venison of the ``projection'', i.e.,
$$\hat a_j = \frac{1}{n} \sum_{i = 1}^n \varphi_j(X_i) - \varphi_j(Y_i).$$
In the large sample limit $n \rightarrow \infty$, both MMD statistics will converge to the population MMD \eqref{eq:eigen_decomp}; see, e.g., Lemma~\ref{lma:scanb_concentration} in Appendix~\ref{appendix:thm2}. Thus, here we assume the population level detectability for simplicity.
\end{remark}

\subsection{ARL and EDD approximations}\label{subsec:arledd}

\subsubsection{Average Run Length}\label{sec:arl_approx}

Approximating ARL $\mathbb{E}_{\infty}[T_b^w]$ in closed-form requires studying the behaviors of extremes of random fields.
Typically, solving this problem involves three steps \citep{yakir2013extremes}: (i) change-of-measure via exponential tilting, (ii) applying likelihood ratio identity to change the probability of interest into expectation, and (iii) invoking localization theorem. Examples using this technique include \citet{siegmund1995using,siegmund2008detecting,siegmund2013sequential,li2019scan}, and so on.
We defer the technical details to Appendix~\ref{appendix:thm1} and give the ARL approximation as follows:

\begin{lemma}[ARL approximation]\label{lma:arl}
As $b \rightarrow \infty$, an approximation to ARL for our proposed Memory Efficient Online Kernel CUSUM procedure $T_b^w$ \eqref{eq:stopping_time} is given by:
\if\mycolumn1
\begin{align}
    \mathbb{E}_{\infty}[T_b^w]=\frac{\sqrt{2 \pi}}{b} \cdot\Bigg\{  \sum_{B = 2}^w e^{\psi_B(\theta_B) - \theta_B b} \frac{\left(2 B-1\right)}{ B\left(B-1\right)} \label{eq:arl_approx}  \cdot \nu\left(\theta_B \sqrt{\frac{2\left(2 B-1\right)}{B\left(B-1\right)}}\right)\Bigg\}^{-1} \cdot[1+o(1)],
\end{align}
\else
\begin{align}
    \mathbb{E}_{\infty}[T_b^w]=\frac{\sqrt{2 \pi}}{b} \cdot\Bigg\{ & \sum_{B = 2}^w e^{\psi_B(\theta_B) - \theta_B b} \frac{\left(2 B-1\right)}{ B\left(B-1\right)} \label{eq:arl_approx} \\
    & \cdot \nu\left(\theta_B \sqrt{\frac{2\left(2 B-1\right)}{B\left(B-1\right)}}\right)\Bigg\}^{-1} \cdot[1+o(1)], \nonumber
\end{align}
\fi
where $\psi_B(\cdot)$ is defined as:
\begin{equation*}
    \psi_B(\theta) = \sum_{n=1}^\infty \frac{\mathbb{E}_{\infty} [Z_B^n(t)]}{n!} \theta^n,
\end{equation*}
and $\theta_B$ is obtained via solving $\dot{\psi_B}(\theta_B) = b$. Moreover, the function $\nu(\cdot)$ can be approximated as:
\begin{equation*}
    \nu(\mu) \approx \frac{(2 / \mu)(\Phi(\mu / 2)-0.5)}{(\mu / 2) \Phi(\mu / 2)+\phi(\mu / 2)},
\end{equation*}where $\phi(\cdot)$ and $\Phi(\cdot)$ are the probability density function and the cumulative distribution function of the standard
normal distribution, respectively.
\end{lemma}

\begin{remark}[Comparison with \citet{li2019scan}]
The difference between our ARL approximation \eqref{eq:arl_approx} and that of Scan $B$ procedure (cf. Theorem 4.2 \citep{li2019scan}) is the additional summation over $B \in [2:w]$. This is straightforward --- In the exponential approximation \eqref{eq:exp_approx}, there is an additional maximization term over $B \in [2:w]$. When applying the change-of-measure by likelihood ratio identity in step (ii), we will see this maximization changes into summation, which leads to an additional summation over $B \in [2:w]$ in \eqref{eq:arl_approx}.
\end{remark}

In practice, we typically incorporate the first two order moments' information to approximate ARL, i.e.,
\begin{equation}\label{eq:theta_B_Gaussian_assumption}
    \psi_B(\theta_B) = \theta_B^2/2, \quad \theta_B = b.
\end{equation}
This gives us
\begin{equation}\label{eq:ARL_Gaussian_assumption}
    \mathbb{E}_{\infty}[T_b^w]=\sqrt{2 \pi} \frac{b  e^{b^2/2}}{ w } [1+o(1)].
\end{equation}
Detailed derivation of the approximation in \eqref{eq:ARL_Gaussian_assumption} can be found in Appendix~\ref{appendix:thm1}.

To obtain a more accurate ARL approximation, we could consider incorporating higher order moments when solving for $\theta_B$ via equation $\dot{\psi_B}(\theta_B) = b$. However, pre-computing/estimating higher order moments accurately from the reference sample pool could be a challenging practical problem. Nevertheless, we introduce {\it skewness correction}, which additionally considers the third order moment $\mathbb{E}_{\infty} [Z_B^3(t)]$, as follows:
\begin{equation}\label{eq:theta_B_skewness_correction}
\begin{split}
    \psi_B(\theta_B) &= \theta_B^2/2 + \mathbb{E}_{\infty} [Z_B^3(t)] \theta_B^3/6, \\
    \theta_B &= \frac{-1 + \sqrt{1+2b\mathbb{E}_{\infty} [Z_B^3(t)]}}{\mathbb{E}_{\infty} [Z_B^3(t)]}.
\end{split}
\end{equation}
Similar to the second order moment, the third order moment $\mathbb{E}_{\infty} [Z_B^3(t)]$ can be pre-computed using the reference data; please see Lemma~\ref{lma:third_moment_H0} in Appendix~\ref{appendix:moment} for its analytic expression.

\begin{figure}[!htp]
\centerline{
\includegraphics[width = 0.7\textwidth]{fig/exp/fig1_ARL_approx.pdf} }
\caption{Comparison of the detection threshold $b$ obtained via (i) Monte Carlo simulation and (ii) theoretical approximation. By this comparison, we can observe that our theoretical ARL approximation, especially after skewness correction, is very close to the Monte Carlo simulation result.}
\label{fig:arl_approx}
\end{figure}

To demonstrate our ARL approximation's practical usefulness, we compare it with Monte Carlo simulation results for three types of pre-change distributions (Gaussian, Exponential, and Laplace distributions) in Figure~\ref{fig:arl_approx}. From this figure, we can observe that, when the ARL is not very large, the approximations are fairly close to the simulated ones for all three cases.
Another key finding is that skewness correction can help improve the ARL approximation accuracy for all three cases, highlighting its practical usefulness in calibrating the detection.
In particular, for Gaussian pre-change distribution, the ARL approximation coupled with skewness correction remains to be very accurate even when the ARL becomes very large.
Later in Appendix~\ref{appendix:opT_b^wmax}, we will theoretically show that considering only finite order moments can lead to underestimating threshold $b$ asymptotically; this underestimation of $b$ phenomena is again validated in our numerical example above.
Lastly, we want to clarify that Gaussian pre-change distribution $\PP$ does not imply Gaussian distribution of the Scan $B$-statistic $Z_B(t)$, meaning that considering its first two orders' moments will not completely characterize its distribution. As we can see in Figure~\ref{fig:arl_approx}, skewness correct does help improve the ARL approximation accuracy for Gaussian pre-change distribution case.

Lastly, we want to emphasize that our proposed skewness correction results in a more accurate ARL approximation than that of the Scan $B$ procedure (cf. Section 6 \citep{li2019scan}), even though our derivations are very similar.
To be precise, $\theta_B$ is different before \eqref{eq:theta_B_Gaussian_assumption} and after \eqref{eq:theta_B_skewness_correction} the skewness correction. To obtain an accurate ARL approximation, we should update $\theta_B$ with its skewness correction version in \eqref{eq:arl_approx}. However, \citet{li2019scan} did NOT consider such update of $\theta_B$ in $\nu\left(\theta_B \sqrt{\frac{2\left(2 B-1\right)}{B\left(B-1\right)}}\right)$, which is as important as updating $\theta_B$ in $e^{\psi_B(\theta_B) - \theta_B b}$, since both terms involve exponentially large parts w.r.t. $\theta_B$.
As a result, our proposed skewness correction gives a more accurate ARL approximation. We numerically validate this finding in Figure~\ref{fig:arl_approx_sw} for the Scan $B$ procedure.

\begin{figure*}[!htp]
\centerline{
\includegraphics[width = \textwidth]{fig/exp/fig1_better_skewness.pdf}
}
\caption{Comparison of the detection threshold $b$ for the Scan $B$ procedure between our proposed skewness correction \eqref{eq:theta_B_skewness_correction} and the original one in Section 6 \citep{li2019scan}. The data is generated from pre-change distribution $\PP$, which is specified on top of each panel. We can observe that ARL approximation with our proposed skewness correction is more accurate under different settings.}
\label{fig:arl_approx_sw}
\end{figure*}

\subsubsection{Expected Detection Delay}

The difficulty in approximating the EDD $\mathbb{E}_{0}[T_b^w]$ comes from the form of the detection statistic, which is the maximum of several correlated r.v.s. and has increment NOT independent from its own past. Therefore, we cannot naively apply Wald's identity to approximate the EDD.
Nevertheless, we will leverage the Azuma's concentration inequality for martingales with bounded difference to study the concentration of the Scan $B$-statistic; see Lemma~\ref{lma:scanb_concentration} in Appendix~\ref{appendix:thm2}. We denote
\begin{equation}\label{eq:cpq}
    C = \sqrt{2}\left(\frac{\mathbb{E}\left[h^{2}\left(X, X^{\prime}, Y, Y^{\prime}\right)\right]}{N} + \frac{N-1}{N}\operatorname{Cov}\left[h(X,X',Y,Y'), h\left(X^{\prime \prime}, X^{\prime \prime \prime}, Y, Y^{\prime}\right)\right]\right)^{-1/2},
\end{equation}
where $X, X^{\prime}, X^{\prime \prime}, X^{\prime \prime \prime}, Y, Y^{\prime} $ are i.i.d. r.v.s with $\PP$. The EDD approximation for sufficiently large $b$ is then given as follows:

\begin{lemma}[EDD approximation]\label{lma:EDD}
Under Assumptions~\ref{assumption:detection} and~\ref{assumption:bounded_kernel}, for window size $w$ chosen as:
\begin{equation}\label{eq:Bmax_choice}
     w  \geq \frac{7 b}{C \cdot \text{\rm MMD}^2(\PP,\QQ)},
\end{equation}
where constant $C$ is defined in \eqref{eq:cpq},
for any $N > 0$, as $b \rightarrow \infty$, the first-order approximation to EDD for our proposed Memory Efficient Online Kernel CUSUM $T_b^w$ \eqref{eq:stopping_time} is given by:
\begin{equation}\label{eq:first_edd_approx}
    \mathbb{E}_{0} [T_b^w] = \frac{b}{C \cdot \text{\rm MMD}^2(\PP,\QQ)  } \ (1+o(1)).
\end{equation}
\end{lemma}

The above first-order EDD approximation agrees with the intuition --- when there is a change-point, the average increment of the detection statistic at each time step will be the increment of including one more sample from $\QQ$ to the post-change block, which is $C \cdot \text{\rm MMD}^2(\PP,\QQ)$ (as evidenced by Lemma~\ref{lma:mean_under_H1} in Appendix~\ref{appendix:moment}). Therefore, the EDD will be approximately the detection threshold $b$ divided by the average one-step increment $C \cdot \text{\rm MMD}^2(\PP,\QQ)$.
We defer the proof of Lemma~\ref{lma:EDD} to Appendix~\ref{appendix:thm2}. In that proof, one pitfall is that we cannot precisely characterize the $o(1)$ overshoot term, and thus we cannot characterize the $o(1)$ term in the EDD approximation \eqref{eq:first_edd_approx} even though the rest of our analysis is non-asymptotic.

\begin{remark}[Effect of pre-change block number $N$]\label{rmk:Neffect}
On one hand, $N$ controls the remainder term in \eqref{eq:EZTb}, and only when $N$ is large enough can we ignore this remainder term. On the other hand, it directly controls the EDD \eqref{eq:first_edd_approx} by controlling the normalizing constant $C$ \eqref{eq:cpq}. For example, consider $\PP = N(\mu \mathbf{1}_{d}, \sigma^2 I_{d})$ and Gaussian RBF kernel with bandwidth parameter $\gamma$ selected by median heuristic. The Taylor expansion (see detailed derivations in Appendix~\ref{appendix:moment}) will give us:
$$0 < \mathbb{E}\left[h^{2}\left(X, X^{\prime}, Y, Y^{\prime}\right)\right] - \operatorname{Cov}\left[h(X,X',Y,Y'), h\left(X^{\prime \prime}, X^{\prime \prime \prime}, Y, Y^{\prime}\right)\right] = \cO(1/d).$$
This tells us that $C$ {\it monotonically increases} to ${\sqrt{2}}/{\sqrt{\operatorname{Cov}\left[h(X,X',Y,Y'), h\left(X^{\prime \prime}, X^{\prime \prime \prime}, Y, Y^{\prime}\right)\right]}}$ as $N \rightarrow \infty$, i.e., larger $N$ leads to smaller EDD \eqref{eq:first_edd_approx}.
In other words, when $N$ goes from one to infinity, the variance of ${Z}^\prime_B(t)$ \eqref{eq:scan-b_unmoralized} shrinks from
${B \choose 2}^{-1} \mathbb{E}\left[h^{2}\left(X, X^{\prime}, Y, Y^{\prime}\right)\right]$
to ${B \choose 2}^{-1}\operatorname{Cov}[h(X,X',Y,Y'), h(X^{\prime \prime},X^{\prime \prime \prime}, Y, Y^{\prime})]$ (as evidenced by Lemma~\ref{lma:var_H0} in Appendix~\ref{appendix:moment}), which enables quicker change detection.
Moreover, the fact that $C$ cannot grow unboundedly with increasing $N$
agrees with the intuition --- in two-sample test, if we only increase the size of one set of samples (i.e., more pre-change blocks) but hold another set of samples (i.e., post-change block) fixed, we cannot expect that the variance of the test statistic shrinks to zero.
\end{remark}

\subsection{Optimal window size}\label{subsec:optimal_Bmax}
In this subsection, we will establish the optimal choice of window size parameter $w$, where the optimality is in the sense of \citet{lai1995sequential} by comparing our procedure with an oracle procedure. To begin with, let us introduce the oracle detection procedure.

\subsubsection{Oracle detection procedure}
The oracle detection procedure, for a threshold $b$, is defined by the following stopping time:
\begin{equation}\label{eq:stopping_time_noWL}
    T_b^{\infty} = \inf \left\{t : \max_{B \in [2:t]} {Z}_B(t) > b \right\},
\end{equation}
where we assume infinite memory capacity and memorize all observations up to time $t$. In contrast, our proposed window-limited version \eqref{eq:stopping_time} only memorizes the most recent $w$ streaming samples.

\subsubsection{Formal definition of the optimality}

There is a trade-off between computation and memory efficiency and the detection power. Here, the memory efficiency is characterized by window size $w$ and we study the ``detection power'' via the EDD for a given ARL. On one hand, as we will see later in Figure~\ref{fig:EDDARL_compare_bminbmax},
larger $w$ leads to smaller EDD and ``better power'', but such ``power growth'' w.r.t. increasing $w$ has a limit, which is characterized by the above oracle procedure.
On the other hand, the above oracle procedure involves in total $ \cO(t^2)$ computation and $ \cO(t)$ memory up to time $t$, which is NOT favorable in online setting, whereas our proposed memory efficient version \eqref{eq:stopping_time} only requires linear computation complexity and constant memory complexity.
Therefore, it is natural to ask:
{\it when holding constant ARL, what is the optimal window size $w = w^\star$ in order to achieve a target performance loss compared with the oracle procedure at the minimal computation and memory cost?}

Inspired by \citet{lai1995sequential}, we define the performance loss of our proposed procedure $T_b^w$ \eqref{eq:stopping_time} compared with the oracle procedure $T_b^{\infty}$ \eqref{eq:stopping_time_noWL} as their EDD difference, i.e., $\mathbb{E}_0[T_b^w] - \mathbb{E}_0[T_b^{\infty}]$.
Here, the detection threshold $b$ is typically chosen such that the corresponding procedure lies in the following constant ARL procedures set:
\begin{equation}\label{eq:constant_arl_procedure}
    \cC_{\rm ARL} = \{T : \mathbb{E}_{\infty} [T] \geq \text{\rm ARL}\}.
\end{equation}
Two fundamental questions are:
\begin{itemize}
    \item [(i)] For a given ARL, how do we choose $b$ to ensure $T_b^w, T_b^{\infty} \in \cC_{\rm ARL}$?
    \item [(ii)] For any fixed tolerance $\varepsilon > 0$ and $b$ choice from question (i) above, what is the optimal window size $w^\star$ to achieve the most memory efficient procedure subject to $\varepsilon$-performance loss? To be precise, $w^\star$ is defined as follows:
    \begin{align}
w^\star = \ & \min w,  \nonumber \\
   & \ {\rm s.t. \ } 0 \leq \mathbb{E}_0[T_b^w] - \mathbb{E}_0[T_b^{\infty}] \leq \varepsilon. \label{eq:opt_w_def}
\end{align}
\end{itemize}

\subsubsection{Main results}

To approach question (i), we leverage the analytic ARL approximation of our proposed procedure in Lemma~\ref{lma:arl} to choose $b$ such that $T_b^w \in \cC_{\rm ARL}$. Unfortunately, it is difficult to obtain such ARL approximation for the oracle procedure. Instead, we apply the Azuma's inequality to show the $b$ choice for our proposed procedure also ensures $T_b^{\infty} \in \cC_{\rm ARL}$. Answering question (ii) requires deriving analytic expression of the loss, i.e., $\mathbb{E}_0[T_b^w] - \mathbb{E}_0[T_b^{\infty}]$, which is a very difficult problem. Alternatively, we give a tight upper bound on the loss using the Azuma's inequality, which in turn gives us a tight upper bound on the optimal window size. Here, we say the upper bound on some quantity is tight in the sense that as the ARL and $b$ approach infinity, the upper bound has the same order with the original quantity up to a constant factor.

Now, we answer both questions formally as follows:

\begin{theorem}[Optimal window size]\label{thm:opT_b^wmax}
Under Assumptions~\ref{assumption:detection} and~\ref{assumption:bounded_kernel}, for any $N > 0$, we will have:
(i) As the ARL goes to infinity, if we choose
\begin{equation}\label{eq:opt_choice}
    b \sim {\log \text{\rm ARL}},
\end{equation}
then the ARLs of our propose Memory Efficient Online Kernel CUSUM $T_b^w$ \eqref{eq:stopping_time} and its oracle procedure $T_b^{\infty}$ \eqref{eq:stopping_time_noWL} will both satisfy the ARL constraint, i.e.,
\begin{equation*
    T_b^w, \  T_b^{\infty} \in \cC_{\rm ARL}.
\end{equation*}

\noindent
(ii) For detection threshold $b$ chosen as in \eqref{eq:opt_choice}, for any fixed tolerance $\varepsilon > 0$,
the optimal window size $w^\star$ \eqref{eq:opt_w_def} has an upper bound $w^\star \leq \bar w$ up to a constant factor, i.e.,
\begin{equation}\label{eq:w_star}
    w^\star \sim \bar w = \frac{6b}{C \cdot \text{\rm MMD}^2(\PP,\QQ) } + \frac{512K^2 \log {3}/{\varepsilon}}{b^2 \left(\frac{N}{4} \wedge \frac{b}{C \cdot \text{\rm MMD}^2(\PP,\QQ)} \right)},
\end{equation}
where constant $C$ is defined in \eqref{eq:cpq}.
\end{theorem}

Here, the notation $\sim$ in \eqref{eq:opt_choice} and \eqref{eq:w_star} means same order up to a constant factor,
which is formally defined at the beginning of Section~\ref{sec:formulation}.
Unfortunately, we cannot give a precise characterization of the constants since (i) our ARL approximation \eqref{eq:arl_approx} has a very complex form and (ii) it is difficult to evaluate the performance loss $\mathbb{E}_0[T_b^w] - \mathbb{E}_0[T_b^{\infty}]$ in closed-form.
Nevertheless, the $b$ and $w$ choices in \eqref{eq:opt_choice} and \eqref{eq:w_star} are informative enough for us to make connections with classic parametric results.
In addition, our non-asymptotic analysis does enable us to derive an
upper bound on the performance loss as follows:
\if\mycolumn1
\begin{align*}
\begin{split}
        &\mathbb{E}_0[T_b^w]  - \mathbb{E}_0[T_b^{\infty}]  \leq 3\exp\left\{-\frac{b^2}{512K^2} \left(\frac{N}{4} \wedge \frac{b}{C \cdot \text{\rm MMD}^2(\PP,\QQ) }\right)\left(w- \frac{6b}{C \cdot \text{\rm MMD}^2(\PP,\QQ)} \right)\right\}.
\end{split}
\end{align*}
\else
\begin{align*}
\begin{split}
        &\mathbb{E}_0[T_b^w]  - \mathbb{E}_0[T_b^{\infty}]  \\
        &\leq 3\exp\left\{-\frac{b^2}{512K^2} \left(\frac{N}{4} \wedge \frac{b}{C \cdot \text{\rm MMD}^2(\PP,\QQ) }\right)(w-6b/(C \cdot \text{\rm MMD}^2(\PP,\QQ)) )\right\}.
\end{split}
\end{align*}
\fi
This upper bound tells us that choosing $w \sim b$ can achieve zero power loss asymptotically. Furthermore, it enables us to derive a non-asymptotic upper bound on $w^\star$ to guarantee $\varepsilon$-performance loss \eqref{eq:opt_w_def} for any fixed tolerance $\varepsilon > 0$ as follows:
\begin{equation*}
    w \geq \bar w,
\end{equation*}
where $\bar w$ is given in \eqref{eq:w_star}. Here, the upper bound $\bar w$ on the optimal window size $w^\star$,
is tight up to a constant factor, i.e. $w^\star \sim \bar w$, as we can show that any $w = o(b)$ choice cannot meet the $\varepsilon$-performance loss constraint \eqref{eq:opt_w_def}.
The proof of Theorem~\ref{thm:opT_b^wmax}
is deferred to Appendix~\ref{appendix:opT_b^wmax}.

\subsubsection{Interpretations of main results}\label{rmk:connection_to_classic}
The above analysis tells us that large $w$ does shrink the performance loss compared with the oracle procedure. However, for large enough ARL, the performance improvement by increasing $w$ will no longer be significant when $w \sim b$; this is also validated by our EDD approximation in Lemma~\ref{lma:EDD} and later in our numerical simulation.

The main message from Theorem~\ref{thm:opT_b^wmax} is that, for any fixed $\PP$ and $\QQ$, the optimal window size \eqref{eq:opt_w_def} is
\begin{equation*}
    w^\star  \sim {\log \text{\rm ARL}},
\end{equation*}
which is the same with the classic result \eqref{eq:lai_opt_window_size} for window-limited GLR under the Gaussian mean shift assumption \citep{lai1995sequential}. Moreover, under the above optimal choice, the corresponding EDD of our Memory Efficient Online Kernel CUSUM will be
\begin{equation}\label{eq:opt_EDD}
    \mathbb{E}_{0} [T_b^{w^\star}] \sim \frac{{\log \text{\rm ARL}}}{\text{\rm MMD}^2(\mathbb{P},\mathbb{Q})},
\end{equation}
which closely resembles the EDD of the classic CUSUM procedure \eqref{eq:cusum}.

Before ending our theoretical analysis, we want to give two remarks highlighting potential future improvements of the current work.

\begin{remark}[Sub-optimality]
Since the KL divergence $\text{\rm KL}(\mathbb{P}\|\mathbb{Q})$ appearing in the denominator of the Lorden's Lower Bound \eqref{eq:lorden_lower_bound} is always larger than $\text{\rm MMD}^2(\mathbb{P},\mathbb{Q})$ \citep{sriperumbudur2012empirical}, \eqref{eq:opt_EDD} indicates that our proposed procedure is not minimax optimal and the amount of sub-optimality is quantified by $\text{\rm MMD}^2(\mathbb{P},\mathbb{Q})$.
However, \citet{balasubramanian2021optimality} showed that such gap can be largely bridged by leveraging a modified kernel embedding \citep{balasubramanian2021optimality} and proved asymptotic minimax optimality for the kernel Goodness-of-Fit test.
Furthermore, it is very important to notice that the kernel embedding modification is sample-size dependent \citep{balasubramanian2021optimality}, and such sample-size dependent modifications to prove minimax optimality are commonly seen in literature, such as
\citet{li2019optimality}, who considered sample-size dependent bandwidth parameter for Gaussaian RBF kernel,
and \citet{zhu2019asymptotically}, who considered a threshold decaying at the rate of square root of the sample size and proved minimax optimality via Sanov's theorem.
As one can see, in online change-point detection, the ``ARL goes to infinity'' limiting process will lead to the ``sample size goes to infinity'' limiting process under our increasing window size scheme $w \sim b$.
Existing literature suggests that sample-size dependent remedy on kernel choice or calibration method might close the gap theoretically.
However, such sample-size dependent scheme is less conventional in online change-point detection and we leave this topic as a future development of the current work.
\end{remark}

\begin{remark}[Curse of dimensionality]
It is well known that MMD, i.e., $\text{\rm MMD}^2(\mathbb{P},\mathbb{Q})$, suffers from the curse of dimensionality as it drops to zero exponentially fast with increasing dimension \citep{ramdas2015decreasing,reddi2015high}. Clearly, based on \eqref{eq:opt_EDD}, EDD will grow exponentially fast w.r.t. dimensionality and our proposed method cannot alleviate this issue. A good news is, even under high dimensional setting, if the difference between $\PP$ and $\QQ$ lies on a low-dimensional manifold, the testing power could still be maintained in the two-sample test problem \citep{cheng2021kernel}. How to explore the effect of such low-dimensional structure in online manifold change-point detection problem remains an open problem and we leave it to future discussion.
\end{remark}

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

In this section, we conduct comprehensive numerical experiments to demonstrate the good performance of our proposed method by comparing it with various benchmark procedures under various settings. Moreover, we will study the influence of hyperparameters, such as $N$ and $ w $, to guide their choices in practice and validate our theoretical findings.

\subsection{Comparison with benchmarks}
To begin with, we will demonstrate our method is capable of achieving quicker change detection. In addition to the (a) Scan $B$ procedure considered in Figure~\ref{fig:EDDARL_compare_w_scanB}, we consider (b) Kernel CUSUM (KCUSUM) and (c) Hotelling's $T^2$ procedures as the benchmark procedures here. We will numerically compare their EDDs for given ARLs under various settings. Before getting into details, let us briefly introduce the detection statistics of the aforementioned benchmark procedures as follows:

\vspace{.1in}
\noindent
{\it (a) Scan $B$-statistic:} For a fixed block size $B_0 \geq 2$, the Scan $B$-statistic $Z_{B_0}(t)$ is defined in \eqref{eq:scan-b}. Here, we choose $B_0 = w$, where $ w $ is the window size of our proposed method.

\vspace{.1in}
\noindent
{\it (b) KCUSUM statistic:} \citet{flynn2019change} replaced the GLR statistic in the classic CUSUM recursive update \eqref{eq:cusum} with the linear-time MMD statistic \citep{gretton2012kernel} and proposed the KCUSUM procedure.
\if\mycolumn1
The KCUSUM detection statistic has the following recursive update rule:
\begin{align*}
    S^{\rm K}_0 = 0, \quad S^{\rm K}_t=\left\{\begin{array}{ll}
\left(S^{\rm K}_{t-1} + h\left(x_{1,t}, x_{2,t}, y_{t-1}, y_{t}\right) - \delta \right)^+ & \text{ if $t$ is even}, \\
S^{\rm K}_{t-1} & \text{ if $t$ is odd},
\end{array}\right.
\end{align*}
\else
The KCUSUM detection statistic has the following recursive update rule: if $t$ is even, then
\begin{align*}
    & S^{\rm K}_t \\
    = & 0 \vee (S^{\rm K}_{t-1} + h\left(x_{1,t}, x_{2,t}, y_{t-1}, y_{t}\right) - \delta),
\end{align*}
if $t$ is odd, then
\begin{align*}
    S^{\rm K}_t=  S^{\rm K}_{t-1},
\end{align*}
and $S^{\rm KCUSUM}(0) = 0$.
\fi
where $h(\cdot,\cdot,\cdot,\cdot)$ is defined in \eqref{eq:kernel_h} and can be viewed as a linear-time MMD statistic with sample size $n = 2$. The hyperparameter $\delta > 0$ makes sure a negative drift and is chosen to be $\delta = 1/50$ as suggested by the original paper.

\vspace{.1in}
\noindent
{\it (c) Hotelling's $T^2$ statistic:}
At time step $t$, for hypothetical change-point $r < t$, we split the sample points into two parts: $U = (X_1,X_2,\dots,X_M,Y_1,\dots,Y_{r-1}), V = (Y_{r},\dots,Y_t).$ For notational simplicity, we denote the elements in $U$ and $V$ by $U_1,\dots,U_{M + r - 1}$ and $V_1,\dots,V_{t - r + 1}$, respectively.
We define
$$
T^{2}_t(r)=\frac{(M + r - 1)(t - r + 1)}{M + t}\left(\bar{U}-\bar{V}\right)^{{ \mathrm{\scriptscriptstyle T} }} \widehat{\Sigma}^{-1}\left(\bar{U}-\bar{V}\right),
$$where $\bar{U} = \sum_{i=1}^{M+r-1} U_i/(M + r - 1)$, $\bar{V} =  \sum_{i=1}^{t-r+1} V_i/(t - r + 1)$ and $\widehat{\Sigma}$ is the pooled covariance matrix:
\if\mycolumn1
\begin{align*}
    \widehat{\Sigma}=(M + t-2)^{-1}\Bigg(&\sum_{i=1}^{M+r-1}\left(U_{i}-\bar{U}\right)\left(U_{i}-\bar{U}\right)^{{ \mathrm{\scriptscriptstyle T} }} +\sum_{i=1}^{t-r+1}\left(V_{i}-\bar{V}\right)\left(V_{i}-\bar{V}\right)^{{ \mathrm{\scriptscriptstyle T} }}\Bigg).
\end{align*}
\else
\begin{align*}
    \widehat{\Sigma}=(M + t-2)^{-1}\Bigg(&\sum_{i=1}^{M+r-1}\left(U_{i}-\bar{U}\right)\left(U_{i}-\bar{U}\right)^{{ \mathrm{\scriptscriptstyle T} }} \\
    &+\sum_{i=1}^{t-r+1}\left(V_{i}-\bar{V}\right)\left(V_{i}-\bar{V}\right)^{{ \mathrm{\scriptscriptstyle T} }}\Bigg).
\end{align*}
\fi
The Hotelling's $T^2$ detection statistic is defined as follows:
\begin{equation*
    S^{{\rm HT}^2}_t = \max_{1\leq r \leq t-1} T^{2}_t(r).
\end{equation*}

\vspace{.1in}
\noindent \textit{Behavior of our proposed statistic.} We first compare the trajectory of our proposed detection statistic with that of the Scan $B$-statistic.
For our proposed detection procedure, we will choose $N=15$, $w = 50$ in the following experiments unless explicitly specified, as suggested by the experiments in \citet{li2019scan} and Figure~\ref{fig:N_effect} below. For MMD-based kernel methods, we choose the Gaussian RBF kernel and use the median heuristic to select the bandwidth parameter in the following experiments.
We plot the mean trajectories over $100$ independent trials and highlight the mean $\pm$ standard deviation area in Figure~\ref{fig:detect_stat_illus_more}.

\begin{figure}[!htp]
\centerline{
\includegraphics[width = 0.7\textwidth]{fig/new exp/fig1_stat_traj_more.pdf}
}
\caption{The trajectories of the mean (solid line) and standard deviation (shaded area) of our proposed and Scan $B$ detection statistics over 100 independent trials. Change from $\mathbb{P}=N(\mathbf{0}_{20}, I_{20})$ to Gaussian mixture with $N(\mathbf{0}_{20}, I_{20}) \ \text{\rm w.p. } 0.3$ and $N(\mathbf{0}_{20},4 \ I_{20}) \ \text{\rm w.p. } 0.7$ occurs at $t=51$.
}
\label{fig:detect_stat_illus_more}
\end{figure}

From Figure~\ref{fig:detect_stat_illus_more}, we can observe that our detection statistic has a ``larger slope'' after the change-point, which is a result of the window-limited maximization in the CUSUM-type procedure. This illustrates that our proposed procedure achieves quicker change detection by removing pre-change samples from the post-change block.
Moreover, we can observe that after $ w  = 50$ steps from the change-point, both statistics stop increasing and oscillate around a ``platform'', and they have almost the same value with each other. This validates that choosing such large $ w $ \eqref{eq:Bmax_choice} in our EDD analysis is necessary, otherwise it will take much longer to stop the procedure and leads to overly large EDD.

We need to additionally mention that larger detection statistic does not necessarily induce quicker detection. This is because the window-limited maximization will also enlarge the statistic under $H_0$ and the threshold $b$ for a given ARL in our proposed procedure. To numerically validate this practical merit, i.e., quicker change detection, of our proposed method, we conduct more experiments by comparing their EDDs directly.

\vspace{.1in}
\noindent \textit{EDD comparison.}
As we mentioned before, for a given ARL, the detection threshold $b$ should be chosen accordingly and thus it can be different for different detection procedures. This means Figure~\ref{fig:detect_stat_illus_more} cannot serve as evidence that our proposed procedure outperforms the Scan $B$ procedure. Here, we choose threshold $b$ to meet the ARL constraint independently for each detection procedure via $1000$ Monte Carlo trials. Next, for each procedure, the EDDs are obtained as the average detection delays over another $1000$ Monte Carlo trials with $10000$ reference samples following pre-change distribution $\PP$ and $50$ streaming samples following post-change distribution $\QQ$. We plot the EDD against ARL for those four procedures in Figure~\ref{fig:EDDARL_compare_comprenhensive}.

\begin{figure*}[!htp]
\centerline{
\includegraphics[width = \textwidth]{fig/new exp/fig1_power.pdf}
}
\caption{Comparison of EDD for our proposed (red), Scan $B$ (orange), KCUSUM (blue) and Hotelling's $T^2$ (brown) procedures. The setting is $\PP = N(\mathbf{0}_{20}, I_{20})$ and  $\QQ$ is Gaussian mixture:  $N(\mathbf{0}_{20}, I_{20}) \ \text{\rm w.p. } 0.3; N(\mu \mathbf{1}_{20}, \sigma^2 \ I_{20}) \ \text{\rm w.p. } 0.7$.
}
\label{fig:EDDARL_compare_comprenhensive}
\end{figure*}

In Figure~\ref{fig:EDDARL_compare_comprenhensive}, the absence of the dot indicates the corresponding procedure fails to detect the change before time $t = 50$ under the corresponding setting. We have the following observations:
(i) First and foremost, our proposed method is not only the most powerful, i.e., it achieves the quickest detection, but also the most robust, i.e., it can detect very small changes which KCUSUM and parametric Hotelling's $T^2$ procedures fail to detect, in the experimental settings considered here.
(ii) Secondly, even though KCUSUM achieves quicker detection than the other two benchmark procedures, it easily fails when the change is subtle or the ARL is large, meaning it is not robust enough.
(iii) Lastly, since the mean shift is obvious under the settings considered in Figure~\ref{fig:EDDARL_compare_comprenhensive}, Hotelling's $T^2$ prceodure performs pretty well compared with Scan $B$ procedure.

\begin{table*}
\caption{EDDs for given ARLs. The change is from standard Gaussian distribution  $\mathbb{P}=N(\mathbf{0}_{20}, I_{20})$ to Gaussian mixture distribution.
Non-value ``$-$'' indicates the corresponding detection procedure fails to detect the change within the first 50 time steps under the corresponding setting. Smallest EDDs are highlighted and we can observe that our proposed method outperforms all existing approaches under the following settings.}\label{table:EDD_Gmix}

\vspace{-0.05in}

\begin{center}
\begin{small}
\begin{sc}
\resizebox{.9\textwidth}{!}{%
\begin{tabular}{lccccccccccccccccc}
\multicolumn{16}{c}{\large{Setting I: $N(\mathbf{0}_{20}, I_{20})$ v.s.
Gaussian mixture with $N(\mathbf{0}_{20}, I_{20}) \ \text{\rm w.p. } 0.3; N(\mu \mathbf{1}_{20},\sigma^2  I_{20}) \ \text{\rm w.p. } 0.7$.}} \\
\toprule[1pt]
& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 0.1, \ \sigma^2 = 0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ \sigma^2 = 0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ \sigma^2 = 1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ \sigma^2 = 4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ \sigma^2 = 9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{19.2} & \textbf{19.55} & \textbf{21.57} & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & \textbf{6.99} & \textbf{7.08} & \textbf{7.66} & \textbf{3.47} & \textbf{3.49} & \textbf{3.6} \\
Scan $B$ & 28.1 & 28.7 & 30.83 & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 15.48 & 15.91 & 17.36 & 9.63 & 9.94 & 10.98 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 6.03 & 6.89 & $-$ \\
Hotelling $T^2$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 10.04 & 10.46 & 12.17 \\
\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 0.3, \ \sigma^2 = 0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ \sigma^2 = 0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ \sigma^2 = 1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ \sigma^2 = 4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ \sigma^2 = 9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{12.55} & \textbf{12.77} & \textbf{14.19} & \textbf{17.26} & \textbf{17.53} & \textbf{19.46} & $-$ & $-$ & $-$ & \textbf{6.69} & \textbf{6.77} & \textbf{7.33} & \textbf{3.46} & \textbf{3.47} & \textbf{3.61} \\
Scan $B$ & 22.36 & 22.93 & 24.84 & 24.87 & 25.5 & 27.61 & $-$ & $-$ & $-$ & 14.85 & 15.28 & 16.69 & 9.53 & 9.85 & 10.91 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 6.04 & 6.71 & 8.18 \\
Hotelling $T^2$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 17.34 & 17.88 & $-$ & 8.88 & 9.28 & 10.68 \\
\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 1, \ \sigma^2 = 0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1, \ \sigma^2 = 0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1, \ \sigma^2 = 1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1, \ \sigma^2 = 4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1, \ \sigma^2 = 9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{3.33} & \textbf{3.34} & \textbf{3.45} & \textbf{3.63} & \textbf{3.67} & \textbf{4.05} & \textbf{4.79} & \textbf{4.85} & \textbf{5.26} & \textbf{4.65} & \textbf{4.7} & \textbf{5.15} & \textbf{3.36} & \textbf{3.37} & \textbf{3.46} \\
Scan $B$ & 9.64 & 10.01 & 11.11 & 10.11 & 10.46 & 11.69 & 11.39 & 11.81 & 13.23 & 11.16 & 11.56 & 12.94 & 9.07 & 9.36 & 10.41 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 5.16 & 5.85 & 6.73 \\
Hotelling $T^2$ & 9.25 & 9.56 & 10.53 & 9.12 & 9.36 & 10.39 & 8.72 & 8.96 & 9.91 & 7.13 & 7.42 & 8.34 & 5.52 & 5.71 & 6.43 \\
\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 1.5, \ \sigma^2 = 0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ \sigma^2 = 0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ \sigma^2 = 1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ \sigma^2 = 4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ \sigma^2 = 9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{3.12} & \textbf{3.13} & \textbf{3.18} & \textbf{3.16} & \textbf{3.17} & \textbf{3.21} & \textbf{3.3} & \textbf{3.3} & \textbf{3.39} & \textbf{3.62} & \textbf{3.64} & \textbf{3.86} & \textbf{3.32} & \textbf{3.32} & \textbf{3.38} \\
Scan $B$ & 7.18 & 7.31 & 8.24 & 7.33 & 7.48 & 8.5 & 8.06 & 8.39 & 9.35 & 9.27 & 9.62 & 10.72 & 8.75 & 8.98 & 10.01 \\
KCUSUM & 3.99 & 4.05 & 4.08 & 4.1 & 4.18 & 4.28 & 4.95 & 5.3 & 5.86 & 7.05 & 8.37 & $-$ & 4.78 & 5.07 & 5.62 \\
Hotelling $T^2$ & 5.74 & 5.86 & 6.59 & 5.72 & 5.87 & 6.51 & 5.63 & 5.8 & 6.34 & 5.11 & 5.24 & 5.77 & 4.37 & 4.49 & 4.98 \\
\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 2, \ \sigma^2 = 0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 2, \ \sigma^2 = 0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 2, \ \sigma^2 = 1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 2, \ \sigma^2 = 4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 2, \ \sigma^2 = 9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{2.89} & \textbf{2.89} & \textbf{2.97} & \textbf{2.95} & \textbf{2.96} & \textbf{3.03} & \textbf{3.11} & \textbf{3.12} & \textbf{3.16} & \textbf{3.31} & \textbf{3.33} & \textbf{3.38} & \textbf{3.28} & \textbf{3.29} & \textbf{3.33} \\
Scan $B$ & 5.94 & 6.44 & 7.11 & 6.39 & 6.74 & 7.19 & 7.02 & 7.17 & 7.87 & 8.32 & 8.62 & 9.49 & 8.53 & 8.71 & 9.68 \\
KCUSUM & 3.98 & 3.98 & 4.01 & 3.98 & 4 & 4.01 & 4.03 & 4.03 & 4.08 & 4.77 & 5.07 & 5.56 & 4.4 & 4.59 & 4.89 \\
Hotelling $T^2$ & 4.31 & 4.38 & 5.17 & 4.34 & 4.44 & 5.08 & 4.37 & 4.45 & 4.92 & 4.14 & 4.24 & 4.6 & 3.72 & 3.82 & 4.16 \\
\bottomrule[1pt]
\end{tabular}
}
\end{sc}
\end{small}
\end{center}

\end{table*}

\begin{table*}
\caption{EDDs for given ARLs under more settings. The change is from standard Gaussian distribution  $\mathbb{P}=N(\mathbf{0}_{20}, I_{20})$ to Laplace and Uniform distributions.
}\label{table:EDD_more}

\vspace{-0.05in}

\begin{center}
\begin{small}
\begin{sc}
\resizebox{.9\textwidth}{!}{%
\begin{tabular}{lccccccccccccccccc}
\multicolumn{16}{c}{\large{Setting II: $N(\mathbf{0}_{20}, I_{20})$ v.s.
${\rm Laplace}(\mu \mathbf{1}_{20}, b \mathbf{1}_{20})$.}} \\
\toprule[1pt]
& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 0.1, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{10.16} & \textbf{10.33} & \textbf{11.33} & \textbf{15.91} & \textbf{16.17} & \textbf{17.73} & $-$ & $-$ & $-$ & \textbf{4.64} & \textbf{4.69} & \textbf{5.01} & \textbf{2.37} & \textbf{2.38} & \textbf{2.54} \\
Scan $B$ & 21.75 & 22.1 & 23.64 & 25.73 & 26.44 & 28.69 & $-$ & $-$ & $-$ & 11.25 & 11.6 & 12.99 & 6.82 & 7.03 & 7.7 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 3.67 & 4.19 & 5.06 \\
Hotelling $T^2$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 7.03 & 7.35 & 8.47 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 0.3, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{6.43} & \textbf{6.51} & \textbf{7.23} & \textbf{8.92} & \textbf{9.04} & \textbf{9.97} & \textbf{17.9} & \textbf{18.26} & \textbf{20.22} & \textbf{4.4} & \textbf{4.46} & \textbf{4.76} & \textbf{2.34} & \textbf{2.36} & \textbf{2.5} \\
Scan $B$ & 15.88 & 16.21 & 17.92 & 18.01 & 18.56 & 20.48 & 23.82 & 24.42 & 26.67 & 10.53 & 10.9 & 12.22 & 6.72 & 6.96 & 7.59 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 3.62 & 4.12 & 4.83 \\
Hotelling $T^2$ & 24.95 & 25.4 & 27.12 & 24.11 & 24.56 & 26.1 & 20.89 & 21.38 & 23.08 & 11.83 & 12.31 & 13.99 & 6.12 & 6.43 & 7.31 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 0.5, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.5, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.5, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.5, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.5, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{4.15} & \textbf{4.18} & \textbf{4.6} & \textbf{5.08} & \textbf{5.17} & \textbf{5.65} & \textbf{7.51} & \textbf{7.61} & \textbf{8.32} & \textbf{4.05} & \textbf{4.07} & \textbf{4.35} & \textbf{2.3} & \textbf{2.31} & \textbf{2.45} \\
Scan $B$ & 12.07 & 12.56 & 13.49 & 13.06 & 13.41 & 14.49 & 15.1 & 15.51 & 16.97 & 9.58 & 9.98 & 11.19 & 6.62 & 6.86 & 7.48 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 3.45 & 3.97 & 4.65 \\
Hotelling $T^2$ & 14.66 & 14.83 & 15.45 & 14.32 & 14.52 & 15.28 & 13.04 & 13.32 & 14.35 & 8.42 & 8.75 & 9.94 & 5.28 & 5.52 & 6.27 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 0.7, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.7, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.7, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.7, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.7, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{3.81} & \textbf{3.86} & \textbf{3.99} & \textbf{3.9} & \textbf{3.93} & \textbf{4.03} & \textbf{4.58} & \textbf{4.63} & \textbf{4.96} & \textbf{3.61} & \textbf{3.64} & \textbf{3.89} & \textbf{2.25} & \textbf{2.25} & \textbf{2.38} \\
Scan $B$ & 8.94 & 9.23 & 10.45 & 9.5 & 9.87 & 11.02 & 10.97 & 11.36 & 12.71 & 8.6 & 8.94 & 10.11 & 6.53 & 6.72 & 7.35 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 3.28 & 3.73 & 4.39 \\
Hotelling $T^2$ & 9.69 & 9.98 & 10.94 & 9.33 & 9.66 & 10.78 & 8.51 & 8.82 & 9.93 & 6.41 & 6.62 & 7.49 & 4.6 & 4.79 & 5.4 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 1.5, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2.03} & \textbf{2.03} & \textbf{2.06} & \textbf{2.25} & \textbf{2.27} & \textbf{2.43} & \textbf{2.03} & \textbf{2.04} & \textbf{2.07} \\
Scan $B$ & 5 & 5 & 5.77 & 5 & 5.03 & 5.95 & 5.47 & 5.68 & 6.28 & 6.26 & 6.44 & 7.1 & 6.05 & 6.15 & 7 \\
KCUSUM & 2.01 & 2.03 & 2.05 & 2.07 & 2.11 & 2.18 & 2.54 & 2.72 & 2.94 & 3.52 & 3.99 & 4.81 & 2.48 & 2.75 & 3.06 \\
Hotelling $T^2$ & 4.01 & 4.05 & 4.93 & 4.06 & 4.12 & 4.76 & 4.06 & 4.14 & 4.51 & 3.66 & 3.74 & 4.11 & 3.12 & 3.21 & 3.57 \\

\bottomrule[1pt]
\end{tabular}
}
\end{sc}
\end{small}
\end{center}

\begin{center}
\begin{small}
\begin{sc}
\resizebox{.9\textwidth}{!}{%
\begin{tabular}{lccccccccccccccccc}
\multicolumn{16}{c}{\large{Setting III: $N(\mathbf{0}_{20}, I_{20})$ v.s.
$U[(a-b) \mathbf{1}_{20},  (a+b) \mathbf{1}_{20}]$.}} \\
\toprule[1pt]
& \multicolumn{3}{c}{{$\boldsymbol{(a =  0.1, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.1, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.1, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.1, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.1, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{6} & \textbf{6} & \textbf{6.93} & \textbf{5.26} & \textbf{5.39} & \textbf{5.99} & \textbf{4} & \textbf{4} & \textbf{4} & \textbf{2.06} & \textbf{2.07} & \textbf{2.18} & \textbf{2} & \textbf{2} & \textbf{2} \\
Scan $B$ & 16 & 16.02 & 18 & 13.96 & 14.22 & 15.23 & 10.36 & 10.74 & 11.66 & 6.17 & 6.41 & 7.08 & 5.07 & 5.21 & 5.98 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 3.67 & 4.2 & 5.08 & 2.18 & 2.23 & 2.33 \\
Hotelling $T^2$ & 29.1 & 29.46 & 30.77 & 20.28 & 20.61 & 22.21 & 12.3 & 12.65 & 13.41 & 5.19 & 5.31 & 5.95 & 3.95 & 4 & 4.26 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(a =  0.3, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.3, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.3, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.3, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.3, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{4} & \textbf{4} & \textbf{4.51} & \textbf{4} & \textbf{4} & \textbf{4} & \textbf{2.92} & \textbf{2.99} & \textbf{3.52} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} \\
Scan $B$ & 12.63 & 13 & 13.97 & 10.66 & 10.98 & 11.91 & 7.54 & 7.98 & 9.06 & 5.48 & 5.79 & 6.25 & 5 & 5.01 & 5.61 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 2.4 & 2.55 & 2.83 & 2.05 & 2.07 & 2.1 \\
Hotelling $T^2$ & 15.77 & 15.99 & 16.89 & 13.03 & 13.17 & 14.03 & 7.55 & 7.92 & 9.17 & 4.8 & 4.89 & 5.05 & 3.5 & 3.64 & 3.98 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(a =  0.7, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.7, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.7, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.7, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.7, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{2.01} & \textbf{2.02} & \textbf{2.95} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} \\
Scan $B$ & 7 & 7 & 8.1 & 6.28 & 6.96 & 7.03 & 5.67 & 5.98 & 6.27 & 5 & 5 & 5.21 & 4.97 & 5 & 5 \\
KCUSUM & 6.71 & $-$ & $-$ & 3.88 & 4.77 & 6.2 & 2.37 & 2.55 & 2.84 & 2.01 & 2.02 & 2.03 & \textbf{2} & \textbf{2} & \textbf{2} \\
Hotelling $T^2$ & 7 & 7 & 8.07 & 6 & 6.06 & 7 & 5 & 5 & 5.4 & 3.92 & 3.97 & 4.04 & 3 & 3.01 & 3.16 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(a =  1.5, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(a =  1.5, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  1.5, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  1.5, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  1.5, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} \\
Scan $B$ & 5 & 5 & 5 & 5 & 5 & 5 & 4.12 & 4.99 & 5 & 4 & 4.09 & 5 & 4.05 & 4.69 & 5 \\
KCUSUM & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} \\
Hotelling $T^2$ & 4 & 4 & 4 & 4 & 4 & 4 & 3 & 3.02 & 3.97 & 2.96 & 2.99 & 3 & 2.05 & 2.1 & 2.6 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(a =  2, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(a =  2, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  2, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  2, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  2, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} \\
Scan $B$ & 4 & 4 & 5 & 4 & 4 & 5 & 4 & 4 & 5 & 4 & 4 & 5 & 4.01 & 4.45 & 5 \\
KCUSUM & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} \\
Hotelling $T^2$ & 3 & 3 & 3 & 3 & 3 & 3 & 3 & 3 & 3 & 2.01 & 2.04 & 2.74 & \textbf{2} & \textbf{2} & 2.01 \\

\bottomrule[1pt]
\end{tabular}
}
\end{sc}
\end{small}
\end{center}

\end{table*}

To further demonstrate the good performance of our proposed method, we compare it with the above three benchmark methods under more scenarios.
To be precise, we additionaly consider more settings for change from Gaussian to Gaussian mixture as well as change from Gaussian distribution to Laplace and uniform distributions.
The results are reported in Tables~\ref{table:EDD_Gmix} and \ref{table:EDD_more}.
From these two tables, we can observe that our proposed Memory Efficient Online Kernel CUSUM achieves the quickest detection among all benchmarks, even under the settings where the changes are too small for parametric procedures to detect.
In addition to aforementioned observations, we can see (iv) Hotelling's $T^2$ procedure performs pretty well when the mean shift is significant, but fails easily otherwise.
Under those small mean shift circumstances, kernel methods can easily achieve better performance than Hotelling's $T^2$ procedure; see the first two rows in the table for evidence.
This observation agrees with our intuition since Hotelling's $T^2$ procedure is a parametric method which is designed to detect the mean shift.

\subsection{Study on the hyperparameter choices}
In this subsection, we will study how the choices of the hyperparameters of our detection procedure influence the EDD. As mentioned above, the choice of window size $w$ is important and we will begin our numerical study with its effect on the EDD.

\vspace{.1in}
\noindent \textit{Effect of window size.} Recall that we use $[2:w]$ as the search region to locate the change-point. Here, in addition to right end-point $w$, i.e., window size, we also study the effect of left end-point. Here, we use $[B_{\text{\rm min}}: B_{\text{\rm max}}]$ as the region in which we optimize block size parameter $B$; we use $w$ and $B_{\text{\rm max}}$ interchangeably to denote the window size parameter. Intuitively, the wider this search region is, the quicker the detection should be, but this merit comes with higher computation and memory cost. We validate this by comparing EDDs for a given ARL under different $B_{\text{\rm min}}$ and $B_{\text{\rm max}}$ choices and plot the results in Figure~\ref{fig:EDDARL_compare_bminbmax}.

\begin{figure*}[!htp]
\centerline{
\includegraphics[width = \textwidth]{fig/bminbmax/fig3_search_region.pdf}
}
\caption{Comparison of EDD for different search regions $[B_{\text{\rm min}}: B_{\text{\rm max}}]$.
The setting is $\mathbb{P}=N(\mathbf{0}_{20}, I_{20})$, and $\QQ$ is, from left to right, (a) $N(\mathbf{0}_{20},0.01 \ I_{20})$; (b) $N(\mathbf{0}_{20},0.2 \ I_{20})$; (c) $N(\mathbf{0}_{20},9 \ I_{20})$.
}
\label{fig:EDDARL_compare_bminbmax}
\end{figure*}

From Figure~\ref{fig:EDDARL_compare_bminbmax}, we observe that,
on one hand, $B_{\text{\rm max}}$ {\it controls the success of detection} --- the detection procedure typically fails to raise an alarm within a reasonable time if EDD of the oracle procedure \eqref{eq:detect_stat_noWL} exceeds $B_{\text{\rm max}}$. However, this does not tell users to choose an overly large window size, since (i) there is a trade-off between the EDD and the computation and memory cost and (ii) EDD improvement by choosing a larger window size $w$ for a given ARL is no longer significant after $w \sim \log$ ARL, as shown by our non-asymptotic study of the performance loss in Section~\ref{subsec:optimal_Bmax}.
Most importantly, (ii) can be numerically verified by the leftmost panel in Figure~\ref{fig:EDDARL_compare_bminbmax} --- As we can see, for $B_{\text{\rm min}} = 2$ case, the improvement of $B_{\text{\rm max}}$ increasing from 40 to 50 is smaller than that of $B_{\text{\rm max}}$ increasing from 30 to 40. In addition, in the rightmost panel, we can see the improvement of EDD by increasing $B_{\text{\rm max}}$ is marginal.
On the other hand, $B_{\text{\rm min}}$ {\it is important to guarantee quick detection} --- ``power loss'' will occur if $B_{\text{\rm min}}$ is greater than EDD of the optimal procedure \eqref{eq:detect_stat_noWL}; see $B_{\text{\rm min}} = 10$ case in the rightmost panel for evidence.
Otherwise, the effect of $B_{\text{\rm min}}$ is negligible when the change is subtle. Here, a simple and safe choice would be $B_{\text{\rm min}} = 2$.

\begin{figure}[!htp]
\centerline{
\includegraphics[width = 0.7\textwidth]{fig/bandwidth/fig3_bandwidth.pdf}
}
\caption{Comparison of EDD for different kernel bandwidths of Gaussian RBF kernel function. We consider different bandwidth choices $\gamma \cdot {\rm med}_{\PP}$ where $\gamma \in \{0.5,1,10,100\}$, where ${\rm med}_{\PP}$ is chosen via median heuristic. The change is from $\mathbb{P}=N(\mathbf{0}_{20}, I_{20})$ to $\QQ$: (a) top left: $N(\mathbf{0}_{20},9 \ I_{20})$; (b) top right: $N(\mathbf{1}_{20},9 \ I_{20})$;
(c) bottom left: Gaussian mixture: $N(\mathbf{0}_{20}, I_{20}) \ \text{\rm w.p. } 0.3; N(\mathbf{0}_{20},9 \ I_{20}) \ \text{\rm w.p. } 0.7$;
(d) bottom right: Gaussian mixture: $N(\mathbf{0}_{20}, I_{20}) \ \text{\rm w.p. } 0.3; N(\mathbf{1}_{20},9 \ I_{20}) \ \text{\rm w.p. } 0.7$.
}
\label{fig:EDDARL_compare_bandwidth}
\end{figure}

\vspace{.1in}
\noindent \textit{Effect of kernel bandwidth.}
Now, we study the effect of kernel bandwidth on the EDD.
As mentioned in \citet{ramdas2015decreasing}, the median heuristic gives the largest possible MMD (see observation 2 therein) and over- as well as under-estimated kernel bandwidth will result in much faster decay in MMD (see observations 1 and 3 therein). As mentioned in Theorem~\ref{thm:opT_b^wmax} above, the optimal EDD is proportional to the inverse of MMD and therefore median heuristic should result in the quickest detection. We demonstrate this empirically by plotting EDD v.s. $\log$ ARL, where threshold $b$ is selected via Monte Carlo simulation to satisfy the corresponding ARL in Figure~\ref{fig:EDDARL_compare_bandwidth}.

Here, we slightly abuse the notation $\gamma$. The kernel bandwidth is chosen to be $\gamma \cdot {\rm med}_{\PP}$, where ${\rm med}_{\PP}$ denotes the median of the $\ell_2$ distance matrix for pre-change samples: $\gamma < 1$ corresponds to under-estimation of kernel bandwidth whereas $\gamma > 1$ represents over-estimation. Figure~\ref{fig:EDDARL_compare_bandwidth} shows under- and over-estimation of kernel bandwidth would both lead to larger EDD and therefore median heuristic is indeed the best. Therefore, in practice, it is safe to apply the median heuristic to select kernel bandwidth.

\vspace{.1in}
\noindent \textit{Effect of number of pre-change blocks $N$.} Lastly, we study the behavior of detection statistics when choosing different $N$'s.
Intuitively, by the Law of Large Numbers, the larger the $N$ is, the better the performance/convergence of block-version MMD statistic, i.e., $B$-test statistic \citep{zaremba2013b}, will be. This is theoretically demonstrated in Remark~\ref{rmk:Neffect}. However, the question is, empirically speaking, how large should $N$ be to achieve a satisfying performance? To answer this question, we plot the trajectory of the averaged detection statistic (over 50 independent runs) for our procedure in Figure~\ref{fig:N_effect}.

\begin{figure}[!htp]
\centerline{
\includegraphics[width = 0.6\textwidth]{fig/new exp/fig2_effect_N.pdf}
}
\caption{Comparison of mean trajectories over 50 independent trials for different pre-change blocks' numbers $N$. The change is from $\mathbb{P}=N(\mathbf{0}_{20}, I_{20})$ to Gaussian mixture with $N(\mathbf{0}_{20}, I_{20}) \ \text{\rm w.p. } 0.3$ and $N(\mathbf{0}_{20},4 \ I_{20}) \ \text{\rm w.p. } 0.7$ occurs at $t=16$.
}
\label{fig:N_effect}
\end{figure}

As we can observe in Figure~\ref{fig:N_effect}, the slopes for $N \geq 15$ are large enough and almost the same after the change-point. We can focus on the first few steps after the change-point, since (i) we are interested in the quick detection of the change-point and (ii) the convergence of the statistic to its mean is controlled by both $N$ and $B$. By taking the trade-off between the computation time and memory cost and the performance into consideration, we believe $N=15$ suffices for our numerical simulation. However, we need to remark that, in practice, this $N$ choice is typically determined by how many reference samples we have.

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

In this work, we suggest a simple remedy by considering the CUSUM-type online kernel change-point detection procedure to tackle the potential power loss of existing Shewhart chart-type procedures, such as Scan $B$ procedure.
We show that this seemingly simple modification is necessary both theoretically and empirically --- on one hand, it enables analytic approximations to both ARL and EDD, and helps us establish the optimal window size which coincides with classic parametric results; on the other hand, our extensive numerical experiments demonstrate the competitive performance of our proposed procedure by comparing it with various existing procedures.

\newpage
\bibliographystyle{plainnat}

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

Online change-point detection is a classic and fundamental problem in statistics and related areas. The goal is to detect a change in the underlying data distribution as soon as possible once the change has occurred while satisfying the false alarm constraint.
Classic methods are typically based on parametric models which detect changes in distribution parameters, such as Shewhart chart \citep{shewhart1925application},
Cumulative Sum (CUSUM) \citep{page1954continuous} and
Shiryaev-Roberts (SR) \citep{shiryaev1963optimum}.
For comprehensive and systematic surveys on parametric change-point detection, we refer readers to \citet{basseville1993detection,siegmund2013sequential,tartakovsky2014sequential,xie2021sequential}.

In modern applications, especially in high-dimensional data settings, the exact specification of the probability distribution is difficult, where non-parametric and distribution-free methods are needed. Kernel-based methods have attracted much attention due to their flexibility and broad applicability; see, e.g., kernel Maximum Mean Discrepancy (MMD) \citep{gretton2006kernel,smola2007hilbert,gretton2012kernel} and a survey on its application\citep{muandet2017kernel}. As non-parametric methods, the kernel-based methods do not require distributional assumptions and allow flexible kernel choices.

In kernel-based change-point detection, much emphasis has been on offline change-point detection, where the goal is to detect and locate change-points in retrospect. In the offline setting, the detection is typically conducted by a series of two-sample tests for samples divided by candidate change-point locations. For instance, such works include the maximum Kernel Fisher Discriminant Ratio-based procedure \citep{harchaoui2008kernel}, kernel MMD-based scan $B$ procedure \citep{li2019scan}, and a recent kernel-based scan procedure to boost the detection power in high dimensional setting \citep{song2022new}. In contrast, online change-point detection aims to detect the change as soon as possible from sequential/streaming data.
Existing online kernel-based change-point detection has been focusing on the Shewhart chart-type procedure (e.g., \citet{li2019scan}).
The Shewhart chart-type procedure differs from the online CUSUM procedure, which is based on CUSUM recursion and is generally believed to be more sensitive to detecting small changes \citep{xie2021sequential}. The CUSUM-type procedure has yet to be developed for kernel-based methods due to the increased computational and memory cost with time horizon. For existing kernel-based online change-point detection methods, one can see a summary in Table~\ref{table:kernelcpd} and our literature survey in Section~\ref{subsec:related_work}.

\begin{table}[!ht]
\begin{center}
\caption{A summary of existing online kernel change-point detection procedures.}\label{table:kernelcpd}
\vspace{0.1in}
\begin{small}
\resizebox{.85\textwidth}{!}{%
    \begin{tabular}{|c|c|}
\hline
  CUSUM-type procedure & Shewhart chart-type procedure \\
\hline
 \multirow{3}{*}{\citet{flynn2019change}}  &  \citet{huang2014high,li2019scan} \\
   &  \citet{bouchikhi2019kernel,chang2019kernel} \\
    &  \citet{cobb2022sequential} \\
\hline
\end{tabular}
}
\end{small}
\end{center}
\end{table}

In this work, we develop a kernel-based CUSUM procedure for online change-point detection called the {\it Memory Efficient Online Kernel CUSUM}. The proposed CUSUM-type procedure considers a set of possible change-point locations, and for each hypothesized change-point location, a kernel test statistic is formed by comparing the supposed post-change samples (from the change-point location to the current time) with the reference (pre-change) samples. Then the detection statistic takes the maximum over such a set of parallel detection statistics. Compared to the existing Shewhart chart-type of kernel-based online change-point detection procedure such as \citet{li2019scan}, our procedure combines kernel detection statistics with different block sizes, whereas \citet{li2019scan} only considered one fixed block size. We show that this seemingly minor but essential change leads to a significant boost in detection power. Indeed, whether or not searching for the unknown change-point location is the main difference between the Shewhart chart and the CUSUM procedures and leads to performance gain in detecting small changes in the CUSUM-type procedure. While harvesting the performance gain, we also make the proposed CUSUM-type procedure computationally more efficient by presenting a recursive implementation to achieve a constant computational and memory cost at each time step,  which is essential in the online setting.

Our proposed Memory Efficient Online Kernel CUSUM procedure is theoretically grounded. We derive accurate analytic approximations to two fundamental performance metrics: the Average Run Length (ARL), which is useful for calibrating the procedure in practice to control false alarms, and the Expected Detection Delay (EDD), which is related to detection power. Based on this, we further establish the optimal window size to be $\cO(\log {\rm ARL})$ to achieve (asymptotic) no power loss, which is in the same order as the classic result of the window-limited Generalized Likelihood Ratio (GLR) procedure for parametric models \citep{lai1995sequential}. Since window size captures how much data is needed in the memory, this result indicates the minimal computation and memory cost necessary for online kernel change-point detection.

The remainder of this paper is organized as follows. We first briefly review related works.
Section~\ref{sec:formulation} introduces the online change-point detection problem and reviews the literature.
Section~\ref{sec:detect_stat} presents our proposed Memory Efficient Online Kernel CUSUM procedure.
Section~\ref{sec:theory} presents analytic approximations to the ARL and EDD and characterizes the optimal window size.
Section~\ref{sec:numerical} validates our theoretical findings and compares our proposed method with existing methods through extensive numerical experiments.
Section~\ref{sec:conclusion} concludes the paper.

\subsection{Literature}\label{subsec:related_work}

\vspace{.1in}

\noindent \textit{Kernel methods in change-point detection.}
An early work on kernel change-point detection is \citet{harchaoui2007retrospective}, where they incorporated kernels into dynamic programming and developed an offline multiple change-point estimation method. This method was extended by \citet{arlot2019kernel} to handle unknown number of change-points via a model selection penalty; however, its performance heavily depends on the penalty hyperparameter and it is difficult to control the type I error in practice.
Later on, \citet{garreau2018consistent} performed non-asymptotic analysis to show the multiple change-point estimation by model selection \citep{arlot2019kernel} can estimate the correct number of change-points with high probability and locate them at the optimal rate.
The problem of calibrating detection procedures also exists for maximum Kernel Fisher Discriminant Ratio based approach \citep{harchaoui2008kernel}, where they applied the bootstrap resampling method to compute the detection threshold.
In addition, \citet{bouchikhi2019kernel,ferrari2023online} considered a kernel density ratio (KDR) \citep{huang2006correcting,nguyen2010estimating} based Shewhart chart-type method for online change-point detection.

\vspace{.1in}

\noindent \textit{Kernel MMD and its application in change-point detection.}
\citet{gretton2006kernel,smola2007hilbert,gretton2012kernel} proposed to measure the distance/discrepancy between distributions via the Reproducing Kernel Hilbert space embedding of distributions and developed the well-known Maximum Mean Discrepancy (MMD) statistic. Since then, such kernel MMD-based statistical tests have become popular due to their wide range of applicability to arbitrary domains, ranging from
anomaly detection \citep{zou2014nonparametric}
model criticism and goodness-of-fit test \citep{lloyd2015statistical} and the optimality of such kernel-based tests \citep{li2019optimality,balasubramanian2021optimality},
approximate bayesian computation \citep{park2016k2},
transfer learning \citep{van2018transfer,yang2019polynomial},
distributionally robust optimization \citep{staib2019distributionally},
robust hypothesis testing \citep{sun2021data},
optimal sub-sampling \citep{dwivedi2021kernel,dwivedi2021generalized} as well as change-point detection.

In offline kernel MMD-based change-point detection, there are many other notable contributions apart from aforementioned work such as \citet{li2019scan}, including \citet{truong2019greedy}, who developed a greedy single change-point estimation procedure, \citet{jones2020kernel}, who developed a multiple change-point estimation method by incorporating MMD in the least-squares framework developed by \citet{harchaoui2007retrospective} and leveraging Nyström kernel approximation \citep{williams2001using} to achieve linear complexity, and so on.

Kernel MMD-based online change-point detection procedures typically test for changes between a window of reference data and a window of streaming data.
Generally speaking, there are three types of procedures in online change-point detection, i.e., Shewhart chart, CUSUM, and SR \citep{xie2021sequential}.
Recently, there has been much effort focusing on kernel MMD-based online Shewhart chart-type procedures, which consider fixed-window scan statistics.
In addition to the Scan $B$ procedure \citep{li2019scan} introduced above, notable contributions include \citet{huang2014high}, who first introduced MMD to online change-point detection and applied bootstrap to calibrate the detection;
\citet{chang2019kernel}, who leveraged neural network-based data-driven kernel selection \citep{gretton2012optimal} to boost detection power; \citet{cobb2022sequential}, who considered a memoryless geometric distribution of the stopping time to control the false alarm rate and calibrated the detection via bootstrap.

However, the fixed-length scanning window of Shewhart chart-type procedures fails to account for the change-point location, which may decrease the detection power since the window of streaming data could contain samples from pre-change (in-control) distribution. \citet{cobb2022sequential} indeed tackled this problem by considering a time-varying threshold, but such approach was a heuristic and did not have strong theoretical guarantee.
A much more sensible approach is CUSUM-type procedure, but such procedures are rarely considered in kernel MMD online change-point detection, potentially due to its higher computational and memory cost compared to Shewhart chart-type procedures.
One notable contribution on this aspect is a recent work by \citet{flynn2019change}, who proposed a Kernel CUSUM procedure by replacing the likelihood ratio in classic CUSUM with the linear-time MMD statistic and derived bounds for both ARL and EDD. Unfortunately, such ARL bound fails to help calibrate the detection in practice. In this work, we fill this void by presenting a kernel MMD-based CUSUM-type procedure, explicitly characterizing its important statistical properties, i.e., ARL and EDD, and proving its optimal window size to give a memory efficient procedure yet maintain the detection power.

\vspace{.1in}

\noindent \textit{Other non-parametric change-point detection methods.}
In addition to the aforementioned kernel-based approaches, there is also rich literature on other non-parametric change-point detection procedures.
Notable contributions include
Mann–Whitney-type statistic \citep{hawkins2010nonparametric},
total variation regularization to estimate a piecewise constant structure, where the jumps correspond to multiple change-points \citep{harchaoui2010multiple},
direct density-ratio estimation \citep{kawahara2012sequential},
K Nearest Neighbor based CUSUM-type statistic \citep{chen2019sequential},
Kolmogorov–Smirnov distance \citep{padilla2021optimal} and so on.
Besides, in addition to classic independent and identically distributed (i.i.d.) setting, it is also of interest to handle data with serial correlation in change-point detection.
For example, there is a line of research considering the total variation regularization to estimate the changes in high-dimensional regression parameters \citep{rinaldo2021localizing}, autoregressive parameters in time series \citep{wang2019localizing}, and piecewise constant dynamic background in autoregressive time series \citep{wei2021inferring}.
Another example is \citet{liu2013change}, which extended the direct density-ratio estimation approach to the time series data.
However, how to handle data with serial correlation, e.g., autoregressive time series, is outside the scope of this work and we leave this topic for further study.

\section{Set-up and Background}\label{sec:formulation}

Consider sufficient reference samples from domain $\cX$ (typically $\cX = \RR^d$) with an unknown pre-change distribution $\PP$:
$$X_1,\dots,X_M \overset{i.i.d.}{\sim} \PP.$$
At time $t$, we receive i.i.d. streaming data $Y_1,\dots,Y_t$. The online change-point detection problem tests the null hypothesis
$$H_0: Y_1,\dots,Y_t \overset{i.i.d.}{\sim} \PP,$$ against the alternative hypothesis
$$H_1: \exists \  r<t, \  Y_1,\dots,Y_r \overset{i.i.d.}{\sim} \PP, \  Y_{r+1},\dots,Y_t \overset{i.i.d.}{\sim} \QQ,$$
for some unknown pre- and post-change distributions $\QQ \not= \PP$.

In online change-point detection, the goal is to detect the change as soon as possible subject to the Average Run Length constraint when there is no change. Typically, we control the ARL by properly choosing the detection threshold $b$, and we study the Expected Detection Delay under such ARL constraint. For the same ARL, the smaller the EDD is, the more sensitive/better the detection procedure is.

\vspace{0.1in}

\noindent
\textit{Notations.}
Our notations are conventional. Under the null hypothesis where there is no change, we use $\mathbb{P}_{\infty}$ and $\mathbb{E}_{\infty}$ to denote the probability and expectation in this case. When there is a change from $\PP$ to $\QQ$ at time $t = r$, all the reference samples and the streaming data before the change-point, i.e., $t < r$, follow the pre-change distribution $\PP$ whereas the rest of the streaming data will follow the post-change distribution $\QQ$. We use $\mathbb{P}_{r}$ and $\mathbb{E}_{r}$ to denote the probability and expectation in this case.
In particular, we use $\mathbb{P}_{0}$ and $\mathbb{E}_{0}$ to denote the probability and expectation when all streaming samples follow post-change distribution $\QQ$.
We denote $\mathbf{1}_d = (1,\dots,1)^{ \mathrm{\scriptscriptstyle T} } \in \RR^d$, $\mathbf{0}_d = (0,\dots,0)^{ \mathrm{\scriptscriptstyle T} } \in \RR^d$ and $I_d \in \RR^{d \times d}$ as the identity matrix, where superscript $^{ \mathrm{\scriptscriptstyle T} }$ standards for the vector/matrix transpose.
For integers $0 < m \leq n$, we denote $[m:n] = \{m,\dots,n\}$; in a special case where $m=1$, we denote $[n] = \{1,\dots,n\}$.
In addition, we denote $a \wedge b = \min\{a,b\}$ and $a \vee b = \max\{a,b\}$. In particular, we denote $(a)^+ = \max\{a,0\}$.
Lastly, for asymptotic notations, we use the standard definitions as follows:
\begin{itemize}
    \item $f(m) = o(g(m))$ means for all $c > 0$ there exists $k > 0$ such that $0 \leq f(m) < cg(m)$ for all $m \geq k$;
    \item $f(m) =  \cO(g(m))$ means there exist positive constants $c$ and $k$, such that $0 \leq f(m) \leq cg(m)$ for all $m \geq k$;
    \item $f(m) \sim g(m)$ means there exist positive constants $c$ and $k$, such that $0 \leq cg(m) \leq f(m)$ for all $m \geq k$ and at the same time $f(m) =  \cO(g(m))$.
\end{itemize}

\subsection{Classic parametric procedures}
Let us first introduce three popular types of parametric procedures for online change-point detection: Shewhart chart, Cumulative Sum, and Generalized Likelihood Ratio.

\vspace{.1in}

\noindent \textit{Shewhart Chart.} One of the earliest online change-point detection procedures is the Shewhart chart \citep{shewhart1925application}, where the pre-change and post-change distributions need to be completely specified. Suppose the pre-change distribution $\PP$ has density function $p(\cdot)$ and the post-change distribution $\QQ$ has density function $q(\cdot)$. At time step $t$, the Shewhart chart considers the log-likelihood ratio on the most recent steaming observation $Y_t$ as the detection statistic. Formally, the Shewhart chart detection procedure is defined as
\begin{equation}\label{eq:T_Sh}
    T_{\rm Sh} = \inf \left\{t>0: \log \frac{q(Y_t)}{p(Y_t)} \geq b\right\},
\end{equation}
where $b>0$ is a pre-selected constant detection threshold to meet the ARL constraint. However, the Shewhart chart ignores the past streaming observations, and this ``information loss'' could lead to power loss \citep{xie2021sequential}.

\vspace{.1in}

\noindent \textit{Cumulative Sum.}
CUSUM procedure was first introduced by \citet{page1954continuous}, where again the densities of pre-change distribution $\PP$ and the post-change distribution $\QQ$, i.e., $p(\cdot)$ and $q(\cdot)$, are assumed completely known. At time step $t$, observing steaming data $Y_1,\dots,Y_t$, CUSUM considers the cumulative log-likelihood ratio as the detection statistic, i.e.,
\begin{equation}\label{eq:cusum_original_form}
\begin{split}
    S_t &= \sum_{s=1}^t \log \frac{q(Y_s)}{p(Y_s)} - \min_{0 \leq r \leq t} \sum_{s=1}^r \log \frac{q(Y_s)}{p(Y_s)} \\
    &= \max_{1 \leq r \leq t+1} \sum_{s=r}^{t} \log \frac{q(Y_s)}{p(Y_s)}.
\end{split}
\end{equation}
By considering the complete information, the CUSUM procedure tackles the ``information loss'' issue of the Shewhart chart.

Since the maximum in \eqref{eq:cusum_original_form} is attained at the maximum likelihood estimator of the unknown change-point location, the CUSUM statistic has the following recursive computation rule:
\begin{equation}\label{eq:cusum}
    S_t = \left(S_{t-1} + \log \frac{q(Y_t)}{p(Y_t)}\right)^+ , \quad S_0 = 0.
\end{equation}
The CUSUM procedure for threshold $b$ is then defined by the following stopping time:
\begin{equation}\label{eq:T_cusum}
    T_{\rm CUSUM} = \inf \{t>0: S_t \geq b\}.
\end{equation}
Although CUSUM was first introduced as a heuristic, \citet{moustakides1986optimal} later showed that CUSUM minimizes the worst-case expected detection delay under the ARL constraint.
As the ARL goes to infinity, the CUSUM procedure \eqref{eq:T_cusum} has
\begin{equation}\label{eq:lorden_lower_bound}
    \mathbb{E}_{0}[T_{\rm CUSUM}] \sim \frac{\log \text{\rm ARL}}{\text{\rm KL}(\mathbb{P}\|\mathbb{Q})},
\end{equation}
where $\text{\rm KL}(\mathbb{P}\|\mathbb{Q}) = \mathbb{E}_{0}[\log ({q(Y_1)}/{p(Y_1)})]$ is the Kullback-Leibler (KL) divergence. This is known as the Lorden's Lower Bound \citep{lorden1971procedures}.

\vspace{.1in}
\noindent \textit{Generalized Likelihood Ratio.} GLR procedure \citep{lorden1971procedures} is another popular parametric approach. Similarly, the pre-change density $p(\cdot)$ is assumed known, whereas the post-change density $q(\cdot;\theta)$ has a known form but unknown parameter $\theta$, which takes values from $\Theta$.
The GLR detection procedure is defined as follows:
\begin{equation}\label{eq:GLR}
    T_{\rm GLR} = \inf \left\{t>0: \max_{1\leq r \leq t-1} \sup_{\theta \in \Theta}\sum_{s= r }^{t} \log \frac{q(Y_s;\theta)}{p(Y_s)} \geq b\right\}.
\end{equation}

The detection statistic used by the GLR procedure involves maximum likelihood estimation of both change-point location $r$, which is similar to CUSUM, and density parameter $\theta$, which accounts for the unknown post-change distribution parameter. Unfortunately, the GLR statistic \eqref{eq:GLR} cannot be updated recursively as CUSUM \eqref{eq:cusum} does. To reduce the computational cost, the following {\it window-limited} GLR (W-GLR)
is adopted:
\begin{equation}\label{eq:GLR_window_limited}
    T_{\rm W-GLR} = \inf \left\{t>0: \max_{t-w\leq r \leq t-1} \sup_{\theta \in \Theta}\sum_{s= r }^{t} \log \frac{q(Y_s;\theta)}{p(Y_s)} \geq b\right\},
\end{equation}
which only considers the potential change-point in a window with size $w$.

The choice of window size $w$ is the key to the success of the window-limited GLR procedure. \citet{lai1995sequential} proved that the optimal window size $w$ under the Gaussian mean shift assumption is
\begin{equation}\label{eq:lai_opt_window_size}
    w \sim \log \text{\rm ARL}.
\end{equation}
According to \citet{lai2001sequential}, by choosing such window size, ``there is little loss of efficiency in reducing the computational complexity by the window-limited modification''.

\subsection{Scan $B$ procedure}
In this subsection, we will introduce the Scan $B$ procedure, which is a kernel change-point detection procedure built on the empirical estimator of MMD, which we call {\it MMD statistic}.
We will omit the formal definition of Reproducing Kernel Hilbert space as well as the MMD, which we call {\it population MMD}, and refer readers to \citet{gretton2006kernel,smola2007hilbert,gretton2012kernel} for their formal definitions. Now, let us begin with introducing the MMD statistic.

\vspace{.1in}
\noindent \textit{Maximum Mean Discrepancy statistic.} For two probability distributions $\PP$ and $\QQ$ supported on domain $\cX$, the population MMD, which we denote as $\text{\rm MMD}^2(\PP,\QQ)$, is a non-parametric metric to measure their distance/discrepancy. Given two sets of samples $\mathbf X$ and $\mathbf Y$ with equal sample size $n$, where $\mathbf X = (X_1,\dots,X_n)$ are i.i.d. samples from $\PP$ and
$\mathbf Y = (Y_1,\dots,Y_n)$ are i.i.d. samples from $\QQ$, an unbiased estimator of MMD, i.e., MMD statistic, is given by using U-statistic as follows \citep{gretton2012kernel}:
\begin{equation}\label{eq:MMD_estimator}
    \text{\rm MMD}_{\rm u}^2(\mathbf X, \mathbf Y)= \frac{1}{n(n-1)} \sum_{i=1}^{n} \sum_{j \neq i}^{n} h\left(X_{i}, X_{j}, Y_{i}, Y_{j}\right),
\end{equation}
where
\if\mycolumn1
 \begin{equation}\label{eq:kernel_h}
    h(x_{1}, x_{2}, y_{1}, y_{2})=k\left(x_{1}, x_{2}\right)+k\left(y_{1}, y_{2}\right)-k\left(x_{1}, y_{2}\right)-k\left(x_{2}, y_{1}\right),
 \end{equation}
\else
 \begin{equation}\label{eq:kernel_h}
 \begin{split}
    h&(x_{i}, x_{j}, y_{i}, y_{j})\\
    &=k\left(x_{i}, x_{j}\right)+k\left(y_{i}, y_{j}\right)-k\left(x_{i}, y_{j}\right)-k\left(x_{j}, y_{i}\right),
 \end{split}
 \end{equation}
\fi
and $k(\cdot,\cdot)$ is the user-specified kernel function:
\begin{equation}\label{eq:kernel}
    k(\cdot,\cdot): \cX \times \cX \rightarrow \RR.
\end{equation}
Commonly used kernel functions include Laplace kernel $k(x,y) = \exp \{-{\norm{x-y}}/{\gamma}\},$ and Gaussian radial basis function (RBF) kernel $k(x,y) = \exp \{-{\norm{x-y}^2}/{\gamma^2}\},$ where $\gamma$ is the bandwidth parameter. As one will see later in the next section, our theoretical analysis is general and allows flexible kernel choices.

\vspace{.1in}
\noindent \textit{Scan $B$ procedure.}
One drawback of the above unbiased MMD estimator \eqref{eq:MMD_estimator} comes from its quadratic computational cost. This is undesirable in the online setting, where we need to evaluate the detection statistic in real time as we receive the streaming data. To fix this issue, \citet{li2019scan} leveraged the $B$-test statistic \citep{zaremba2013b} and developed the Scan $B$ procedure.

In Scan $B$ procedure, they randomly sample $N B$ samples from reference data $X_1,\dots,X_M$ without replacement, in order to build $N$ {\it pre-change blocks} with equal block size $B$, i.e., $\mathbf{X}_B^{(i)}, \ i \in [N]$. In both Scan $B$ procedure and our proposed procedure below, by saying we have sufficient amount of reference samples, we assume reference sample size $M$ is large enough such that $M > N B$. At time $t$, the {\it post-change block} consists of the most recent $B$ streaming data, i.e., $\mathbf{Y}_B(t) = (Y_{t-B+1},\dots, Y_t)$.
They compute unbiased MMD statistics between the post-change block and each pre-change block and take their average as the $B$-test statistic.
By reusing the post-change block $\mathbf{Y}_B(t)$, Scan $B$ procedure achieves overall linear computational complexity with respect to both time $t$ and sample size $N$. Formally, the detection statistic, which they call Scan $B$-statistic, is defined as follows:
\begin{equation}\label{eq:scan-b}
    {Z}_B(t) = \frac{{Z}^\prime_B(t)}{\sqrt{\operatorname{Var}_{\infty}({Z}^\prime_B(t))}},
\end{equation}
where $\operatorname{Var}_{\infty} (\cdot)$ denotes the variance under $H_0$ and
\begin{equation}\label{eq:scan-b_unmoralized}
    {Z}^\prime_B(t) = \frac{1}{N}\sum_{i=1}^N {\text{\rm MMD}_{\rm u}^2(\mathbf{X}_B^{(i)},\mathbf{Y}_B(t))}.
\end{equation}

For detection threshold $b$, the Scan $B$ procedure is defined by the following stopping time:
\begin{equation}\label{eq:stopping_time_scanb}
    T_{{\rm Sc}B} = \inf \{t : {Z}_B(t) > b\}.
\end{equation}

The above normalizing constant, i.e., the standard deviation of ${Z}^\prime_B(t)$ under $H_0$, remains unchanged with increasing $t$. Moreover, it has an an analytic expression (see Lemma~\ref{lma:var_H0} in Appendix~\ref{appendix:moment}) and can be computed using the reference samples before the detection procedure.

\begin{remark}[Effect of the normalizing constant]
Given that $\operatorname{Var}_{\infty}({Z}^\prime_B(t))$ is held constant in the online Scan $B$ procedure and does NOT influence the online detection, one may ask why we include this normalizing constant here. In fact, \citet{li2019scan} introduced this normalizing constant under the offline setting, where they proposed to use window-limited maximization with respect to (w.r.t.) the block size the estimate the change-point location.
To be precise, consider a fixed number of observations $Y_1,\dots,Y_t$, for a hypothesized change-point at $r < t$, the un-normalized statistics ${Z}^\prime_{2}(t),\dots,{Z}^\prime_{t-r+1}(t)$ \eqref{eq:scan-b_unmoralized} all have the same expectation, i.e., the population MMD, making this un-normalized version undesirable in change-point localization. By scaling those statistics with a factor that decreases w.r.t. the block size $B$ (e.g., $\operatorname{Var}_{\infty}^{-1/2}({Z}^\prime_B(t))$) one can expect that the normalized statistic takes maximum at the change-point location. In the next section, we will see the same effect of normalization when we perform window-limited maximization to locate the change-point in online setting. For brevity, we will use the term Scan $B$-statistic (or procedure) for online Scan $B$-statistic (or procedure) for simplicity, as we do not consider offline setting in this work.
\end{remark}

\section{Proposed Detection Procedure} \label{sec:detect_stat}

In this section, we will formally present our proposed Memory Efficient Online Kernel CUSUM procedure for online change-point detection, which includes a recursive computation of the detection statistic. Compared with Scan $B$ procedure, our proposed procedure significantly achieves quicker change detection while maintains the same computational and memory complexity. Now, we begin with presenting our proposed detection procedure.

\subsection{Memory Efficient Online Kernel CUSUM}

Inspired by the original form of CUSUM \eqref{eq:cusum_original_form}, we propose to consider a parallel set of Scan $B$-statistics ${Z}_B(t)$ \eqref{eq:scan-b} with block size $B$ taking values from $[2:w]$ and take their maximum.
For a pre-selected threshold $b$, the {\it Memory Efficient Online Kernel CUSUM} procedure is defined via the following stopping time:
\begin{equation}\label{eq:stopping_time}
    T_b^w = \inf \left\{t : \max_{B \in [2:w]} {Z}_B(t) > b\right\}.
\end{equation}
where $ w $ is the {\it window size} parameter.

\begin{remark}[Comparison with Scan $B$ procedure]
The window-limited maximization coupled with the normalization can remove pre-change samples from the post-change block, which gives our procedure an edge over Scan $B$ procedure both empirically and theoretically --- it can not only result in quicker change detection (see an numerical evidence in Figure~\ref{fig:EDDARL_compare_w_scanB}), but also help us derive an analytic EDD approximation (see Section~\ref{subsec:arledd} later for details).
In contrast, Scan $B$ procedure include samples from $\PP$ in its post-change block, making it harder to be differentiated from the pre-change blocks. Moreover, the portion of pre-change samples in the post-change block varies at each time step, making it difficult to approximate the EDD in closed-form.
Moreover, as we will see in the next subsection, even though the window-limited maximization introduces more computations, the computational and memory complexity remains the same with Scan $B$ procedure and does not grow with time horizon, which is essential in online setting.

\begin{figure}[!htp]
\centerline{
\includegraphics[width = .95\textwidth]{fig/new exp/fig1_illus_power.pdf}
}
\caption{Comparison of Expected Detection Delay for given Average Run Length between our proposed (red) and the Scan $B$ (orange) detection procedures. ARL is in base $10$ logarithm. The pre-change distribution is $\PP = N(\mathbf{0}_{20}, I_{20})$ and the post-change distribution $\QQ$ is specified on the top of each figure (with $\sigma = 2$). We can observe that our proposed procedure achieves quicker change detection under all three settings.}
\label{fig:EDDARL_compare_w_scanB}
\end{figure}
\end{remark}

\subsection{Recursive computation of detection statistics}

\begin{algorithm}[!htp]
 \caption{Memory Efficient Online Kernel CUSUM}\label{algo:kernel_CUSUM}
 \begin{algorithmic}[1]
 \renewcommand{\algorithmicrequire}{\textbf{Input:} }

 \REQUIRE Pre-change blocks' number $N$, window size $w$, ARL, reference samples $X_1,\dots,X_M$.

  \STATE Initialization: use the most recent $w$ observations to form the post-change block $\mathbf{Y}_t \in \cX^w$; randomly sample $N \times w$ reference samples to form $N$ pre-change blocks $\mathbf{X}^{(i)} \in \cX^w, \ i = 1,\dots,N$; select the threshold $b$ to meet the ARL constraint; compute $C$ \eqref{eq:cpq} using reference samples; calculate the Gram matrices $G_{XX}^{(i)} = k(\mathbf{X}^{(i)},\mathbf{X}^{(i)})$, $G_{XY}^{(i)} = k(\mathbf{X}^{(i)},\mathbf{Y}_t), \ i=1,\dots,N,$ and $G_{YY} = k(\mathbf{Y}_t,\mathbf{Y}_t)$
 \\
  \WHILE {True}
  \STATE receive sample $Y_t$
  \STATE update post-change block:
      $$\mathbf{Y}_t (1:w-1) \leftarrow \mathbf{Y}_t (2:w), \quad
      \mathbf{Y}_t (w) \leftarrow Y_t$$
      \vspace{-0.2in}
  \STATE update Gram matrix $G_{YY}$:
  \vspace{-0.1in}
  \begin{align*}
      &G_{YY} (1:w-1, 1:w-1)  \leftarrow G_{YY} (2:w, 2:w), \\
      &G_{YY} (:, w)  \leftarrow k (\mathbf{Y}_t, Y_t), \quad G_{YY} (w, :)  \leftarrow G_{YY}^{ \mathrm{\scriptscriptstyle T} } (:, w)
  \end{align*}
  \vspace{-0.2in}
  \FOR {$i = 1,\dots,N$}
  \STATE update Gram matrix $G_{XY}^{(i)}$:
  \vspace{-0.1in}
  \begin{align*}
      G_{XY}^{(i)} (:, 1:w-1) & \leftarrow G_{XY}^{(i)} (:, 2:w), \\
      G_{XY}^{(i)} (:, w) & \leftarrow k (\mathbf{X}^{(i)}, Y_t)
  \end{align*}
  \ENDFOR
  \STATE calculate $z = \sum_{i=1}^N \sum_{j=1}^2 \sum_{k \not= j}^2 G_{XX}^{(i)} (j, k) + G_{YY} (j, k) - G_{XY}^{(i)} (j, k) - G_{XY}^{(i)} (k, j)$ and $Z_2(t) = C z / (\sqrt{2} N)$; let $Z_t = Z_2(t)$
  \FOR {$B = 3, \dots, w$}
  \STATE update
  \vspace{-0.1in}
  \if\mycolumn1
\begin{align*}
      z \leftarrow z + \sum_{i=1}^N  \sum_{j = 1}^B & G_{XX}^{(i)} (B, k) + G_{YY} (B, k) -  G_{XY}^{(i)} (B, k) - G_{XY}^{(i)} (k, B)
  \end{align*}
\else
\begin{align*}
      z \leftarrow z + \sum_{i=1}^N  \sum_{j = 1}^B & G_{XX}^{(i)} (B, k) + G_{YY} (B, k) \\
       &-  G_{XY}^{(i)} (B, k) - G_{XY}^{(i)} (k, B)
  \end{align*}
\fi
  \vspace{-0.1in}
  \STATE calculate $Z_B(t) = {C z}/(\sqrt{B(B-1)}N)$
  \STATE update $Z_t \leftarrow \max\{Z_t, Z_B(t)\}$
  \ENDFOR
  \IF {$Z_t > b$}
  \STATE raise an alarm and set $T_b^w = t$
  \STATE \textbf{break}
  \ENDIF
  \ENDWHILE

 \end{algorithmic}
 \end{algorithm}

Before presenting the algorithm, we introduce some notations which only appear in this subsection: for a matrix $G \in \RR^{m \times n}$, $G(i_1:i_2,j_1:j_2)$ represents a slice of matrix $G$ containing its $i_1$-th to $i_2$-th rows and $j_1$-th to $j_2$-th columns; in particular, when we want to include all rows (or columns) in this slice, we use ``$:$'' instead of ``$1:m$'' (or ``$1:n$'').
Similarly, for a vector $g \in \RR^{m}$, $g(i:j)$ represents a slice of $g$, containing its $i$-th to $j$-th elements. Now, we give step-by-step instructions on how to preform our proposed detection procedure in Algorithm~\ref{algo:kernel_CUSUM}.

In Algorithm~\ref{algo:kernel_CUSUM}, the recursive computation of detection statistics comes from two parts: recursive update of Gram matrices (lines $5 - 8$; see a graphical illustration in Figure~\ref{fig:online_update_illus}) and a recursive formulation to find the window-limited maximum (lines $9 - 14$), which are direct benefits of CUSUM and closely resemble that of the original CUSUM \eqref{eq:cusum}.

\begin{figure}[!htp]
\centerline{
\includegraphics[width = 0.9\textwidth]{fig/illus/fig0_illus_onlineupdate.pdf}
}
\caption{Graphical illustration of the recursive update of the Gram matrix in lines $5 - 8$ in Algorithm~\ref{algo:kernel_CUSUM}. Such recursive update was first considered in \citet{li2019scan}. The dark gray part corresponds to the Gram matrix used at the corresponding time step. The highlighted Gram matrix is readily obtained from last time step, which results in $ \cO(N  w )$ complexity for the update.}
\label{fig:online_update_illus}
\end{figure}

Even though the online update of Gram matrices only takes $ \cO(Nw)$ computational complexity, the whole procedure still takes $ \cO(N w^2)$ computations to calculate the detection statistic at each time step (lines $9 - 14$), which is the same with Scan $B$ procedure.
In addition, at each time step, our detection procedure only requires additional $ \cO(w)$ computations in lines $12 - 13$ and does not require additional memory to store those intermediate calculations.
We summarize the comparison among Shewhart chart CUSUM, W-GLR, Scan $B$, our proposed Memory Efficient Online Kernel CUSUM and its oracle variant \eqref{eq:stopping_time_noWL} in Table~\ref{table:complexity}. In particular, we report (i) whether or not the method is robust to the post-change distribution specification, and the (ii) computation and (iii) memory complexity at each time step $t$ for all six methods.
As we can see, the additional window-limit maximization in our proposed Memory Efficient Online Kernel CUSUM does NOT increase the computational and memory complexity compared with the Scan $B$ procedure. We need to clarify that Algorithm~\ref{algo:kernel_CUSUM} actually requires $\cO(Nw^2)$ memory to store the Gram matrix. However, in theory, we can only use $\cO(Nw)$ memory to store raw observations and re-compute the whole Gran matrix with $\cO(Nw^2)$ computations at each time step.

\begin{table}[!ht]
\begin{center}
\caption{Comparison of online change-point detection methods. Here, $w$ is the window size parameter and $N$ is the number of pre-change blocks.}\label{table:complexity}
\vspace{0.1in}
\begin{small}
\resizebox{\textwidth}{!}{%
    \begin{tabular}{|l|c|c|c|}
\hline
 ~ & Robustness & Computation & Memory  \\ \hline
 Shewhart chart \eqref{eq:T_Sh} & \multirow{2}{*}{$\PP$ and $\QQ$ are both known} & $ \cO(1)$ & $-$  \\ \cline{1-1} \cline{3-4}
 CUSUM \eqref{eq:T_cusum} &  & $ \cO(1)$ & $ \cO(1)$  \\ \hline
 W-GLR \eqref{eq:GLR_window_limited} & $\PP$ is known and $\QQ$ has known density form & $ \cO(w^2)$ & $ \cO(w)$  \\ \hline
 Scan $B$ \eqref{eq:stopping_time_scanb} & \multirow{3}{*}{\shortstack{$\PP$ and $\QQ$ are both unknown, but sufficient \\ amount of reference samples is available} } &  $ \cO(Nw^2)$ & $ \cO(Nw)$  \\ \cline{1-1} \cline{3-4}
 Proposed \eqref{eq:stopping_time} &   & $ \cO(Nw^2)$ & $ \cO(Nw)$  \\ \cline{1-1} \cline{3-4}
 Oracle \eqref{eq:stopping_time_noWL} &  & $ \cO(Nt^2)$ & $ \cO(Nt)$  \\ \hline
\end{tabular}
}
\end{small}
\end{center}
\end{table}

\section{Theoretical Analysis} \label{sec:theory}

In this section, we will derive analytic approximations to two fundamental performance metrics in online change-point detection:
(i) the Average Run Length $\mathbb{E}_{\infty}[T_b^w]$, which is the expected stopping time when there is no change, and
(ii) the Expected Detection Delay $\mathbb{E}_{0}[T_b^w]$, which is the expected stopping time when the change occurs immediately at $r = 0$.
We will numerically show that our ARL approximation is very accurate, which helps us to calibrate the detection procedure without time-consuming Monte Carlo simulation.
Furthermore, these results help us to establish the asymptotic optimal window size $ w = w^\star $ in exactly the same sense with \citet{lai1995sequential}.
Now, let us start with the necessary assumptions for our theoretical analysis.

\subsection{Assumptions}
We consider the probability measure $\PP$ (i.e., the pre-change distribution) on a measurable space $(\cX,\cB)$. For symmetric, positive (semi-)definite, and square integrable kernel function $k(\cdot,\cdot)$ \eqref{eq:kernel}. For such kernels, Mercer's theorem gives the following decomposition:
\begin{equation*}
    k(x,x^\prime) = \sum_{j = 1}^\infty \lambda_j \varphi_j(x) \varphi_j(x^\prime),
\end{equation*}
where the limit is in $L^2(\PP)$, $\lambda_1 \geq \lambda_2 \geq \dots \geq 0$ are eigenvalues of
the integral operator induced by kernel $k$, and $\{\varphi_j(\cdot): j \geq 1\}$ are the corresponding orthonormal eigenfunctions such that
\begin{equation*}
    \int_{\cX} \varphi_j(x) \varphi_{j^\prime} (x) d \PP(x) = \left\{\begin{array}{ll}
1 & \text{ if } j = j^\prime, \\
0 & \text{ otherwise}.
\end{array}\right.
\end{equation*}
For concreteness, we further assume kernel $k$ is universal in that $\{\varphi_j(\cdot): j \geq 1\}$ forms an orthonormal basis of $L^2(\PP)$.

Consider the probability measure $\QQ$ (i.e., the post-change distribution) dominated by $\PP$ in that the Radon-Nikodym derivative $d \QQ/ d \PP$ is well-defined. The (squared) population MMD, denoted by $\text{\rm MMD}^2(\PP,\QQ)$, can be expressed as follows:
\begin{align}
    \text{\rm MMD}^2(\PP,\QQ) &= \int_{\cX} \int_{\cX} k(x,x^\prime) d (\PP - \QQ) (x) d (\PP - \QQ) (x^\prime) \nonumber \\
    & = \sum_{j = 1}^\infty \lambda_j \left(\int_{\cX} \varphi_j(x) (p(x) - q(x)) dx \right)^2, \label{eq:eigen_decomp}
\end{align}
where $p$ and $q$ are the densities of $\PP$ and $\QQ$, respectively. To ensure the departure $p - q$ is detectable by a given kernel $k$, we assume there exists $\lambda_j > 0$ such that
\begin{equation}\label{assumption:detection}\tag*{(A1)}
    a_j = \int_{\cX} \varphi_j(x) (p(x) - q(x)) dx \not= 0.
\end{equation}
Here, $a_j$ represents the ``projection'' of the departure $p - q$ along the ``direction'' $\varphi_j$, and Assumption~\ref{assumption:detection} guarantees that
$$\text{\rm MMD}^2(\PP,\QQ) = \sum_{j = 1}^\infty \lambda_j a_j^2 > 0.$$

In addition to the above {\it detectability} assumption, we further impose a technical assumption that all eigenfunctions are uniformly bounded, i.e.,
$$\sup_{j \geq 1} \norm{\varphi_j}_{\infty} < \infty,$$
where $\norm{\varphi_j}_{\infty} = \sup_{x \in \cX} |\varphi_j(x)|$.

This also ensures that the kernel $k(\cdot,\cdot)$ is uniformly bounded on domain $\cX \times \cX$, i.e., there exists a constant $K > 0$ such that
\begin{equation}
    | k(x,y) | \leq K, \quad \forall x, y \in \cX. \label{assumption:bounded_kernel}\tag*{(A2)}
\end{equation}
As one can see, Assumption~\ref{assumption:bounded_kernel} guarantees the existence of the moments of $k(X,Y)$, where $X$ and $Y$ can follow either $\PP$ or $\QQ$.

\begin{remark}
The biased MMD statistic, different from the unbiased estimator $\text{\rm MMD}_{\rm u}^2(\mathbf X, \mathbf Y)$ \eqref{eq:MMD_estimator}, can be expressed by replacing $\PP$ and $\QQ$ with their empirical versions $\hat \PP_n$ and $\hat \QQ_n$ in \eqref{eq:eigen_decomp}, respectively, i.e.,
\begin{align}
    \text{\rm MMD}_{\rm b}^2(\mathbf X, \mathbf Y)& = \int \int k(x,x^\prime) d (\hat \PP_n - \hat \QQ_n) (x) d (\hat \PP_n - \hat \QQ_n) (x^\prime) \nonumber \\
    & = \sum_{j = 1}^\infty \lambda_j \left(\frac{1}{n} \sum_{i = 1}^n \varphi_j(X_i) - \varphi_j(Y_i) \right)^2. \nonumber
\end{align}
Thus, in practice the detectability depends on the empirical venison of the ``projection'', i.e.,
$$\hat a_j = \frac{1}{n} \sum_{i = 1}^n \varphi_j(X_i) - \varphi_j(Y_i).$$
In the large sample limit $n \rightarrow \infty$, both MMD statistics will converge to the population MMD \eqref{eq:eigen_decomp}; see, e.g., Lemma~\ref{lma:scanb_concentration} in Appendix~\ref{appendix:thm2}. Thus, here we assume the population level detectability for simplicity.
\end{remark}

\subsection{ARL and EDD approximations}\label{subsec:arledd}

\subsubsection{Average Run Length}\label{sec:arl_approx}

Approximating ARL $\mathbb{E}_{\infty}[T_b^w]$ in closed-form requires studying the behaviors of extremes of random fields.
Typically, solving this problem involves three steps \citep{yakir2013extremes}: (i) change-of-measure via exponential tilting, (ii) applying likelihood ratio identity to change the probability of interest into expectation, and (iii) invoking localization theorem. Examples using this technique include \citet{siegmund1995using,siegmund2008detecting,siegmund2013sequential,li2019scan}, and so on.
We defer the technical details to Appendix~\ref{appendix:thm1} and give the ARL approximation as follows:

\begin{lemma}[ARL approximation]\label{lma:arl}
As $b \rightarrow \infty$, an approximation to ARL for our proposed Memory Efficient Online Kernel CUSUM procedure $T_b^w$ \eqref{eq:stopping_time} is given by:
\if\mycolumn1
\begin{align}
    \mathbb{E}_{\infty}[T_b^w]=\frac{\sqrt{2 \pi}}{b} \cdot\Bigg\{  \sum_{B = 2}^w e^{\psi_B(\theta_B) - \theta_B b} \frac{\left(2 B-1\right)}{ B\left(B-1\right)} \label{eq:arl_approx}  \cdot \nu\left(\theta_B \sqrt{\frac{2\left(2 B-1\right)}{B\left(B-1\right)}}\right)\Bigg\}^{-1} \cdot[1+o(1)],
\end{align}
\else
\begin{align}
    \mathbb{E}_{\infty}[T_b^w]=\frac{\sqrt{2 \pi}}{b} \cdot\Bigg\{ & \sum_{B = 2}^w e^{\psi_B(\theta_B) - \theta_B b} \frac{\left(2 B-1\right)}{ B\left(B-1\right)} \label{eq:arl_approx} \\
    & \cdot \nu\left(\theta_B \sqrt{\frac{2\left(2 B-1\right)}{B\left(B-1\right)}}\right)\Bigg\}^{-1} \cdot[1+o(1)], \nonumber
\end{align}
\fi
where $\psi_B(\cdot)$ is defined as:
\begin{equation*}
    \psi_B(\theta) = \sum_{n=1}^\infty \frac{\mathbb{E}_{\infty} [Z_B^n(t)]}{n!} \theta^n,
\end{equation*}
and $\theta_B$ is obtained via solving $\dot{\psi_B}(\theta_B) = b$. Moreover, the function $\nu(\cdot)$ can be approximated as:
\begin{equation*}
    \nu(\mu) \approx \frac{(2 / \mu)(\Phi(\mu / 2)-0.5)}{(\mu / 2) \Phi(\mu / 2)+\phi(\mu / 2)},
\end{equation*}where $\phi(\cdot)$ and $\Phi(\cdot)$ are the probability density function and the cumulative distribution function of the standard
normal distribution, respectively.
\end{lemma}

\begin{remark}[Comparison with \citet{li2019scan}]
The difference between our ARL approximation \eqref{eq:arl_approx} and that of Scan $B$ procedure (cf. Theorem 4.2 \citep{li2019scan}) is the additional summation over $B \in [2:w]$. This is straightforward --- In the exponential approximation \eqref{eq:exp_approx}, there is an additional maximization term over $B \in [2:w]$. When applying the change-of-measure by likelihood ratio identity in step (ii), we will see this maximization changes into summation, which leads to an additional summation over $B \in [2:w]$ in \eqref{eq:arl_approx}.
\end{remark}

In practice, we typically incorporate the first two order moments' information to approximate ARL, i.e.,
\begin{equation}\label{eq:theta_B_Gaussian_assumption}
    \psi_B(\theta_B) = \theta_B^2/2, \quad \theta_B = b.
\end{equation}
This gives us
\begin{equation}\label{eq:ARL_Gaussian_assumption}
    \mathbb{E}_{\infty}[T_b^w]=\sqrt{2 \pi} \frac{b  e^{b^2/2}}{ w } [1+o(1)].
\end{equation}
Detailed derivation of the approximation in \eqref{eq:ARL_Gaussian_assumption} can be found in Appendix~\ref{appendix:thm1}.

To obtain a more accurate ARL approximation, we could consider incorporating higher order moments when solving for $\theta_B$ via equation $\dot{\psi_B}(\theta_B) = b$. However, pre-computing/estimating higher order moments accurately from the reference sample pool could be a challenging practical problem. Nevertheless, we introduce {\it skewness correction}, which additionally considers the third order moment $\mathbb{E}_{\infty} [Z_B^3(t)]$, as follows:
\begin{equation}\label{eq:theta_B_skewness_correction}
\begin{split}
    \psi_B(\theta_B) &= \theta_B^2/2 + \mathbb{E}_{\infty} [Z_B^3(t)] \theta_B^3/6, \\
    \theta_B &= \frac{-1 + \sqrt{1+2b\mathbb{E}_{\infty} [Z_B^3(t)]}}{\mathbb{E}_{\infty} [Z_B^3(t)]}.
\end{split}
\end{equation}
Similar to the second order moment, the third order moment $\mathbb{E}_{\infty} [Z_B^3(t)]$ can be pre-computed using the reference data; please see Lemma~\ref{lma:third_moment_H0} in Appendix~\ref{appendix:moment} for its analytic expression.

\begin{figure}[!htp]
\centerline{
\includegraphics[width = 0.7\textwidth]{fig/exp/fig1_ARL_approx.pdf} }
\caption{Comparison of the detection threshold $b$ obtained via (i) Monte Carlo simulation and (ii) theoretical approximation. By this comparison, we can observe that our theoretical ARL approximation, especially after skewness correction, is very close to the Monte Carlo simulation result.}
\label{fig:arl_approx}
\end{figure}

To demonstrate our ARL approximation's practical usefulness, we compare it with Monte Carlo simulation results for three types of pre-change distributions (Gaussian, Exponential, and Laplace distributions) in Figure~\ref{fig:arl_approx}. From this figure, we can observe that, when the ARL is not very large, the approximations are fairly close to the simulated ones for all three cases.
Another key finding is that skewness correction can help improve the ARL approximation accuracy for all three cases, highlighting its practical usefulness in calibrating the detection.
In particular, for Gaussian pre-change distribution, the ARL approximation coupled with skewness correction remains to be very accurate even when the ARL becomes very large.
Later in Appendix~\ref{appendix:opT_b^wmax}, we will theoretically show that considering only finite order moments can lead to underestimating threshold $b$ asymptotically; this underestimation of $b$ phenomena is again validated in our numerical example above.
Lastly, we want to clarify that Gaussian pre-change distribution $\PP$ does not imply Gaussian distribution of the Scan $B$-statistic $Z_B(t)$, meaning that considering its first two orders' moments will not completely characterize its distribution. As we can see in Figure~\ref{fig:arl_approx}, skewness correct does help improve the ARL approximation accuracy for Gaussian pre-change distribution case.

Lastly, we want to emphasize that our proposed skewness correction results in a more accurate ARL approximation than that of the Scan $B$ procedure (cf. Section 6 \citep{li2019scan}), even though our derivations are very similar.
To be precise, $\theta_B$ is different before \eqref{eq:theta_B_Gaussian_assumption} and after \eqref{eq:theta_B_skewness_correction} the skewness correction. To obtain an accurate ARL approximation, we should update $\theta_B$ with its skewness correction version in \eqref{eq:arl_approx}. However, \citet{li2019scan} did NOT consider such update of $\theta_B$ in $\nu\left(\theta_B \sqrt{\frac{2\left(2 B-1\right)}{B\left(B-1\right)}}\right)$, which is as important as updating $\theta_B$ in $e^{\psi_B(\theta_B) - \theta_B b}$, since both terms involve exponentially large parts w.r.t. $\theta_B$.
As a result, our proposed skewness correction gives a more accurate ARL approximation. We numerically validate this finding in Figure~\ref{fig:arl_approx_sw} for the Scan $B$ procedure.

\begin{figure*}[!htp]
\centerline{
\includegraphics[width = \textwidth]{fig/exp/fig1_better_skewness.pdf}
}
\caption{Comparison of the detection threshold $b$ for the Scan $B$ procedure between our proposed skewness correction \eqref{eq:theta_B_skewness_correction} and the original one in Section 6 \citep{li2019scan}. The data is generated from pre-change distribution $\PP$, which is specified on top of each panel. We can observe that ARL approximation with our proposed skewness correction is more accurate under different settings.}
\label{fig:arl_approx_sw}
\end{figure*}

\subsubsection{Expected Detection Delay}

The difficulty in approximating the EDD $\mathbb{E}_{0}[T_b^w]$ comes from the form of the detection statistic, which is the maximum of several correlated r.v.s. and has increment NOT independent from its own past. Therefore, we cannot naively apply Wald's identity to approximate the EDD.
Nevertheless, we will leverage the Azuma's concentration inequality for martingales with bounded difference to study the concentration of the Scan $B$-statistic; see Lemma~\ref{lma:scanb_concentration} in Appendix~\ref{appendix:thm2}. We denote
\begin{equation}\label{eq:cpq}
    C = \sqrt{2}\left(\frac{\mathbb{E}\left[h^{2}\left(X, X^{\prime}, Y, Y^{\prime}\right)\right]}{N} + \frac{N-1}{N}\operatorname{Cov}\left[h(X,X',Y,Y'), h\left(X^{\prime \prime}, X^{\prime \prime \prime}, Y, Y^{\prime}\right)\right]\right)^{-1/2},
\end{equation}
where $X, X^{\prime}, X^{\prime \prime}, X^{\prime \prime \prime}, Y, Y^{\prime} $ are i.i.d. r.v.s with $\PP$. The EDD approximation for sufficiently large $b$ is then given as follows:

\begin{lemma}[EDD approximation]\label{lma:EDD}
Under Assumptions~\ref{assumption:detection} and~\ref{assumption:bounded_kernel}, for window size $w$ chosen as:
\begin{equation}\label{eq:Bmax_choice}
     w  \geq \frac{7 b}{C \cdot \text{\rm MMD}^2(\PP,\QQ)},
\end{equation}
where constant $C$ is defined in \eqref{eq:cpq},
for any $N > 0$, as $b \rightarrow \infty$, the first-order approximation to EDD for our proposed Memory Efficient Online Kernel CUSUM $T_b^w$ \eqref{eq:stopping_time} is given by:
\begin{equation}\label{eq:first_edd_approx}
    \mathbb{E}_{0} [T_b^w] = \frac{b}{C \cdot \text{\rm MMD}^2(\PP,\QQ)  } \ (1+o(1)).
\end{equation}
\end{lemma}

The above first-order EDD approximation agrees with the intuition --- when there is a change-point, the average increment of the detection statistic at each time step will be the increment of including one more sample from $\QQ$ to the post-change block, which is $C \cdot \text{\rm MMD}^2(\PP,\QQ)$ (as evidenced by Lemma~\ref{lma:mean_under_H1} in Appendix~\ref{appendix:moment}). Therefore, the EDD will be approximately the detection threshold $b$ divided by the average one-step increment $C \cdot \text{\rm MMD}^2(\PP,\QQ)$.
We defer the proof of Lemma~\ref{lma:EDD} to Appendix~\ref{appendix:thm2}. In that proof, one pitfall is that we cannot precisely characterize the $o(1)$ overshoot term, and thus we cannot characterize the $o(1)$ term in the EDD approximation \eqref{eq:first_edd_approx} even though the rest of our analysis is non-asymptotic.

\begin{remark}[Effect of pre-change block number $N$]\label{rmk:Neffect}
On one hand, $N$ controls the remainder term in \eqref{eq:EZTb}, and only when $N$ is large enough can we ignore this remainder term. On the other hand, it directly controls the EDD \eqref{eq:first_edd_approx} by controlling the normalizing constant $C$ \eqref{eq:cpq}. For example, consider $\PP = N(\mu \mathbf{1}_{d}, \sigma^2 I_{d})$ and Gaussian RBF kernel with bandwidth parameter $\gamma$ selected by median heuristic. The Taylor expansion (see detailed derivations in Appendix~\ref{appendix:moment}) will give us:
$$0 < \mathbb{E}\left[h^{2}\left(X, X^{\prime}, Y, Y^{\prime}\right)\right] - \operatorname{Cov}\left[h(X,X',Y,Y'), h\left(X^{\prime \prime}, X^{\prime \prime \prime}, Y, Y^{\prime}\right)\right] = \cO(1/d).$$
This tells us that $C$ {\it monotonically increases} to ${\sqrt{2}}/{\sqrt{\operatorname{Cov}\left[h(X,X',Y,Y'), h\left(X^{\prime \prime}, X^{\prime \prime \prime}, Y, Y^{\prime}\right)\right]}}$ as $N \rightarrow \infty$, i.e., larger $N$ leads to smaller EDD \eqref{eq:first_edd_approx}.
In other words, when $N$ goes from one to infinity, the variance of ${Z}^\prime_B(t)$ \eqref{eq:scan-b_unmoralized} shrinks from
${B \choose 2}^{-1} \mathbb{E}\left[h^{2}\left(X, X^{\prime}, Y, Y^{\prime}\right)\right]$
to ${B \choose 2}^{-1}\operatorname{Cov}[h(X,X',Y,Y'), h(X^{\prime \prime},X^{\prime \prime \prime}, Y, Y^{\prime})]$ (as evidenced by Lemma~\ref{lma:var_H0} in Appendix~\ref{appendix:moment}), which enables quicker change detection.
Moreover, the fact that $C$ cannot grow unboundedly with increasing $N$
agrees with the intuition --- in two-sample test, if we only increase the size of one set of samples (i.e., more pre-change blocks) but hold another set of samples (i.e., post-change block) fixed, we cannot expect that the variance of the test statistic shrinks to zero.
\end{remark}

\subsection{Optimal window size}\label{subsec:optimal_Bmax}
In this subsection, we will establish the optimal choice of window size parameter $w$, where the optimality is in the sense of \citet{lai1995sequential} by comparing our procedure with an oracle procedure. To begin with, let us introduce the oracle detection procedure.

\subsubsection{Oracle detection procedure}
The oracle detection procedure, for a threshold $b$, is defined by the following stopping time:
\begin{equation}\label{eq:stopping_time_noWL}
    T_b^{\infty} = \inf \left\{t : \max_{B \in [2:t]} {Z}_B(t) > b \right\},
\end{equation}
where we assume infinite memory capacity and memorize all observations up to time $t$. In contrast, our proposed window-limited version \eqref{eq:stopping_time} only memorizes the most recent $w$ streaming samples.

\subsubsection{Formal definition of the optimality}

There is a trade-off between computation and memory efficiency and the detection power. Here, the memory efficiency is characterized by window size $w$ and we study the ``detection power'' via the EDD for a given ARL. On one hand, as we will see later in Figure~\ref{fig:EDDARL_compare_bminbmax},
larger $w$ leads to smaller EDD and ``better power'', but such ``power growth'' w.r.t. increasing $w$ has a limit, which is characterized by the above oracle procedure.
On the other hand, the above oracle procedure involves in total $ \cO(t^2)$ computation and $ \cO(t)$ memory up to time $t$, which is NOT favorable in online setting, whereas our proposed memory efficient version \eqref{eq:stopping_time} only requires linear computation complexity and constant memory complexity.
Therefore, it is natural to ask:
{\it when holding constant ARL, what is the optimal window size $w = w^\star$ in order to achieve a target performance loss compared with the oracle procedure at the minimal computation and memory cost?}

Inspired by \citet{lai1995sequential}, we define the performance loss of our proposed procedure $T_b^w$ \eqref{eq:stopping_time} compared with the oracle procedure $T_b^{\infty}$ \eqref{eq:stopping_time_noWL} as their EDD difference, i.e., $\mathbb{E}_0[T_b^w] - \mathbb{E}_0[T_b^{\infty}]$.
Here, the detection threshold $b$ is typically chosen such that the corresponding procedure lies in the following constant ARL procedures set:
\begin{equation}\label{eq:constant_arl_procedure}
    \cC_{\rm ARL} = \{T : \mathbb{E}_{\infty} [T] \geq \text{\rm ARL}\}.
\end{equation}
Two fundamental questions are:
\begin{itemize}
    \item [(i)] For a given ARL, how do we choose $b$ to ensure $T_b^w, T_b^{\infty} \in \cC_{\rm ARL}$?
    \item [(ii)] For any fixed tolerance $\varepsilon > 0$ and $b$ choice from question (i) above, what is the optimal window size $w^\star$ to achieve the most memory efficient procedure subject to $\varepsilon$-performance loss? To be precise, $w^\star$ is defined as follows:
    \begin{align}
w^\star = \ & \min w,  \nonumber \\
   & \ {\rm s.t. \ } 0 \leq \mathbb{E}_0[T_b^w] - \mathbb{E}_0[T_b^{\infty}] \leq \varepsilon. \label{eq:opt_w_def}
\end{align}
\end{itemize}

\subsubsection{Main results}

To approach question (i), we leverage the analytic ARL approximation of our proposed procedure in Lemma~\ref{lma:arl} to choose $b$ such that $T_b^w \in \cC_{\rm ARL}$. Unfortunately, it is difficult to obtain such ARL approximation for the oracle procedure. Instead, we apply the Azuma's inequality to show the $b$ choice for our proposed procedure also ensures $T_b^{\infty} \in \cC_{\rm ARL}$. Answering question (ii) requires deriving analytic expression of the loss, i.e., $\mathbb{E}_0[T_b^w] - \mathbb{E}_0[T_b^{\infty}]$, which is a very difficult problem. Alternatively, we give a tight upper bound on the loss using the Azuma's inequality, which in turn gives us a tight upper bound on the optimal window size. Here, we say the upper bound on some quantity is tight in the sense that as the ARL and $b$ approach infinity, the upper bound has the same order with the original quantity up to a constant factor.

Now, we answer both questions formally as follows:

\begin{theorem}[Optimal window size]\label{thm:opT_b^wmax}
Under Assumptions~\ref{assumption:detection} and~\ref{assumption:bounded_kernel}, for any $N > 0$, we will have:
(i) As the ARL goes to infinity, if we choose
\begin{equation}\label{eq:opt_choice}
    b \sim {\log \text{\rm ARL}},
\end{equation}
then the ARLs of our propose Memory Efficient Online Kernel CUSUM $T_b^w$ \eqref{eq:stopping_time} and its oracle procedure $T_b^{\infty}$ \eqref{eq:stopping_time_noWL} will both satisfy the ARL constraint, i.e.,
\begin{equation*
    T_b^w, \  T_b^{\infty} \in \cC_{\rm ARL}.
\end{equation*}

\noindent
(ii) For detection threshold $b$ chosen as in \eqref{eq:opt_choice}, for any fixed tolerance $\varepsilon > 0$,
the optimal window size $w^\star$ \eqref{eq:opt_w_def} has an upper bound $w^\star \leq \bar w$ up to a constant factor, i.e.,
\begin{equation}\label{eq:w_star}
    w^\star \sim \bar w = \frac{6b}{C \cdot \text{\rm MMD}^2(\PP,\QQ) } + \frac{512K^2 \log {3}/{\varepsilon}}{b^2 \left(\frac{N}{4} \wedge \frac{b}{C \cdot \text{\rm MMD}^2(\PP,\QQ)} \right)},
\end{equation}
where constant $C$ is defined in \eqref{eq:cpq}.
\end{theorem}

Here, the notation $\sim$ in \eqref{eq:opt_choice} and \eqref{eq:w_star} means same order up to a constant factor,
which is formally defined at the beginning of Section~\ref{sec:formulation}.
Unfortunately, we cannot give a precise characterization of the constants since (i) our ARL approximation \eqref{eq:arl_approx} has a very complex form and (ii) it is difficult to evaluate the performance loss $\mathbb{E}_0[T_b^w] - \mathbb{E}_0[T_b^{\infty}]$ in closed-form.
Nevertheless, the $b$ and $w$ choices in \eqref{eq:opt_choice} and \eqref{eq:w_star} are informative enough for us to make connections with classic parametric results.
In addition, our non-asymptotic analysis does enable us to derive an
upper bound on the performance loss as follows:
\if\mycolumn1
\begin{align*}
\begin{split}
        &\mathbb{E}_0[T_b^w]  - \mathbb{E}_0[T_b^{\infty}]  \leq 3\exp\left\{-\frac{b^2}{512K^2} \left(\frac{N}{4} \wedge \frac{b}{C \cdot \text{\rm MMD}^2(\PP,\QQ) }\right)\left(w- \frac{6b}{C \cdot \text{\rm MMD}^2(\PP,\QQ)} \right)\right\}.
\end{split}
\end{align*}
\else
\begin{align*}
\begin{split}
        &\mathbb{E}_0[T_b^w]  - \mathbb{E}_0[T_b^{\infty}]  \\
        &\leq 3\exp\left\{-\frac{b^2}{512K^2} \left(\frac{N}{4} \wedge \frac{b}{C \cdot \text{\rm MMD}^2(\PP,\QQ) }\right)(w-6b/(C \cdot \text{\rm MMD}^2(\PP,\QQ)) )\right\}.
\end{split}
\end{align*}
\fi
This upper bound tells us that choosing $w \sim b$ can achieve zero power loss asymptotically. Furthermore, it enables us to derive a non-asymptotic upper bound on $w^\star$ to guarantee $\varepsilon$-performance loss \eqref{eq:opt_w_def} for any fixed tolerance $\varepsilon > 0$ as follows:
\begin{equation*}
    w \geq \bar w,
\end{equation*}
where $\bar w$ is given in \eqref{eq:w_star}. Here, the upper bound $\bar w$ on the optimal window size $w^\star$,
is tight up to a constant factor, i.e. $w^\star \sim \bar w$, as we can show that any $w = o(b)$ choice cannot meet the $\varepsilon$-performance loss constraint \eqref{eq:opt_w_def}.
The proof of Theorem~\ref{thm:opT_b^wmax}
is deferred to Appendix~\ref{appendix:opT_b^wmax}.

\subsubsection{Interpretations of main results}\label{rmk:connection_to_classic}
The above analysis tells us that large $w$ does shrink the performance loss compared with the oracle procedure. However, for large enough ARL, the performance improvement by increasing $w$ will no longer be significant when $w \sim b$; this is also validated by our EDD approximation in Lemma~\ref{lma:EDD} and later in our numerical simulation.

The main message from Theorem~\ref{thm:opT_b^wmax} is that, for any fixed $\PP$ and $\QQ$, the optimal window size \eqref{eq:opt_w_def} is
\begin{equation*}
    w^\star  \sim {\log \text{\rm ARL}},
\end{equation*}
which is the same with the classic result \eqref{eq:lai_opt_window_size} for window-limited GLR under the Gaussian mean shift assumption \citep{lai1995sequential}. Moreover, under the above optimal choice, the corresponding EDD of our Memory Efficient Online Kernel CUSUM will be
\begin{equation}\label{eq:opt_EDD}
    \mathbb{E}_{0} [T_b^{w^\star}] \sim \frac{{\log \text{\rm ARL}}}{\text{\rm MMD}^2(\mathbb{P},\mathbb{Q})},
\end{equation}
which closely resembles the EDD of the classic CUSUM procedure \eqref{eq:cusum}.

Before ending our theoretical analysis, we want to give two remarks highlighting potential future improvements of the current work.

\begin{remark}[Sub-optimality]
Since the KL divergence $\text{\rm KL}(\mathbb{P}\|\mathbb{Q})$ appearing in the denominator of the Lorden's Lower Bound \eqref{eq:lorden_lower_bound} is always larger than $\text{\rm MMD}^2(\mathbb{P},\mathbb{Q})$ \citep{sriperumbudur2012empirical}, \eqref{eq:opt_EDD} indicates that our proposed procedure is not minimax optimal and the amount of sub-optimality is quantified by $\text{\rm MMD}^2(\mathbb{P},\mathbb{Q})$.
However, \citet{balasubramanian2021optimality} showed that such gap can be largely bridged by leveraging a modified kernel embedding \citep{balasubramanian2021optimality} and proved asymptotic minimax optimality for the kernel Goodness-of-Fit test.
Furthermore, it is very important to notice that the kernel embedding modification is sample-size dependent \citep{balasubramanian2021optimality}, and such sample-size dependent modifications to prove minimax optimality are commonly seen in literature, such as
\citet{li2019optimality}, who considered sample-size dependent bandwidth parameter for Gaussaian RBF kernel,
and \citet{zhu2019asymptotically}, who considered a threshold decaying at the rate of square root of the sample size and proved minimax optimality via Sanov's theorem.
As one can see, in online change-point detection, the ``ARL goes to infinity'' limiting process will lead to the ``sample size goes to infinity'' limiting process under our increasing window size scheme $w \sim b$.
Existing literature suggests that sample-size dependent remedy on kernel choice or calibration method might close the gap theoretically.
However, such sample-size dependent scheme is less conventional in online change-point detection and we leave this topic as a future development of the current work.
\end{remark}

\begin{remark}[Curse of dimensionality]
It is well known that MMD, i.e., $\text{\rm MMD}^2(\mathbb{P},\mathbb{Q})$, suffers from the curse of dimensionality as it drops to zero exponentially fast with increasing dimension \citep{ramdas2015decreasing,reddi2015high}. Clearly, based on \eqref{eq:opt_EDD}, EDD will grow exponentially fast w.r.t. dimensionality and our proposed method cannot alleviate this issue. A good news is, even under high dimensional setting, if the difference between $\PP$ and $\QQ$ lies on a low-dimensional manifold, the testing power could still be maintained in the two-sample test problem \citep{cheng2021kernel}. How to explore the effect of such low-dimensional structure in online manifold change-point detection problem remains an open problem and we leave it to future discussion.
\end{remark}

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

In this section, we conduct comprehensive numerical experiments to demonstrate the good performance of our proposed method by comparing it with various benchmark procedures under various settings. Moreover, we will study the influence of hyperparameters, such as $N$ and $ w $, to guide their choices in practice and validate our theoretical findings.

\subsection{Comparison with benchmarks}
To begin with, we will demonstrate our method is capable of achieving quicker change detection. In addition to the (a) Scan $B$ procedure considered in Figure~\ref{fig:EDDARL_compare_w_scanB}, we consider (b) Kernel CUSUM (KCUSUM) and (c) Hotelling's $T^2$ procedures as the benchmark procedures here. We will numerically compare their EDDs for given ARLs under various settings. Before getting into details, let us briefly introduce the detection statistics of the aforementioned benchmark procedures as follows:

\vspace{.1in}
\noindent
{\it (a) Scan $B$-statistic:} For a fixed block size $B_0 \geq 2$, the Scan $B$-statistic $Z_{B_0}(t)$ is defined in \eqref{eq:scan-b}. Here, we choose $B_0 = w$, where $ w $ is the window size of our proposed method.

\vspace{.1in}
\noindent
{\it (b) KCUSUM statistic:} \citet{flynn2019change} replaced the GLR statistic in the classic CUSUM recursive update \eqref{eq:cusum} with the linear-time MMD statistic \citep{gretton2012kernel} and proposed the KCUSUM procedure.
\if\mycolumn1
The KCUSUM detection statistic has the following recursive update rule:
\begin{align*}
    S^{\rm K}_0 = 0, \quad S^{\rm K}_t=\left\{\begin{array}{ll}
\left(S^{\rm K}_{t-1} + h\left(x_{1,t}, x_{2,t}, y_{t-1}, y_{t}\right) - \delta \right)^+ & \text{ if $t$ is even}, \\
S^{\rm K}_{t-1} & \text{ if $t$ is odd},
\end{array}\right.
\end{align*}
\else
The KCUSUM detection statistic has the following recursive update rule: if $t$ is even, then
\begin{align*}
    & S^{\rm K}_t \\
    = & 0 \vee (S^{\rm K}_{t-1} + h\left(x_{1,t}, x_{2,t}, y_{t-1}, y_{t}\right) - \delta),
\end{align*}
if $t$ is odd, then
\begin{align*}
    S^{\rm K}_t=  S^{\rm K}_{t-1},
\end{align*}
and $S^{\rm KCUSUM}(0) = 0$.
\fi
where $h(\cdot,\cdot,\cdot,\cdot)$ is defined in \eqref{eq:kernel_h} and can be viewed as a linear-time MMD statistic with sample size $n = 2$. The hyperparameter $\delta > 0$ makes sure a negative drift and is chosen to be $\delta = 1/50$ as suggested by the original paper.

\vspace{.1in}
\noindent
{\it (c) Hotelling's $T^2$ statistic:}
At time step $t$, for hypothetical change-point $r < t$, we split the sample points into two parts: $U = (X_1,X_2,\dots,X_M,Y_1,\dots,Y_{r-1}), V = (Y_{r},\dots,Y_t).$ For notational simplicity, we denote the elements in $U$ and $V$ by $U_1,\dots,U_{M + r - 1}$ and $V_1,\dots,V_{t - r + 1}$, respectively.
We define
$$
T^{2}_t(r)=\frac{(M + r - 1)(t - r + 1)}{M + t}\left(\bar{U}-\bar{V}\right)^{{ \mathrm{\scriptscriptstyle T} }} \widehat{\Sigma}^{-1}\left(\bar{U}-\bar{V}\right),
$$where $\bar{U} = \sum_{i=1}^{M+r-1} U_i/(M + r - 1)$, $\bar{V} =  \sum_{i=1}^{t-r+1} V_i/(t - r + 1)$ and $\widehat{\Sigma}$ is the pooled covariance matrix:
\if\mycolumn1
\begin{align*}
    \widehat{\Sigma}=(M + t-2)^{-1}\Bigg(&\sum_{i=1}^{M+r-1}\left(U_{i}-\bar{U}\right)\left(U_{i}-\bar{U}\right)^{{ \mathrm{\scriptscriptstyle T} }} +\sum_{i=1}^{t-r+1}\left(V_{i}-\bar{V}\right)\left(V_{i}-\bar{V}\right)^{{ \mathrm{\scriptscriptstyle T} }}\Bigg).
\end{align*}
\else
\begin{align*}
    \widehat{\Sigma}=(M + t-2)^{-1}\Bigg(&\sum_{i=1}^{M+r-1}\left(U_{i}-\bar{U}\right)\left(U_{i}-\bar{U}\right)^{{ \mathrm{\scriptscriptstyle T} }} \\
    &+\sum_{i=1}^{t-r+1}\left(V_{i}-\bar{V}\right)\left(V_{i}-\bar{V}\right)^{{ \mathrm{\scriptscriptstyle T} }}\Bigg).
\end{align*}
\fi
The Hotelling's $T^2$ detection statistic is defined as follows:
\begin{equation*
    S^{{\rm HT}^2}_t = \max_{1\leq r \leq t-1} T^{2}_t(r).
\end{equation*}

\vspace{.1in}
\noindent \textit{Behavior of our proposed statistic.} We first compare the trajectory of our proposed detection statistic with that of the Scan $B$-statistic.
For our proposed detection procedure, we will choose $N=15$, $w = 50$ in the following experiments unless explicitly specified, as suggested by the experiments in \citet{li2019scan} and Figure~\ref{fig:N_effect} below. For MMD-based kernel methods, we choose the Gaussian RBF kernel and use the median heuristic to select the bandwidth parameter in the following experiments.
We plot the mean trajectories over $100$ independent trials and highlight the mean $\pm$ standard deviation area in Figure~\ref{fig:detect_stat_illus_more}.

\begin{figure}[!htp]
\centerline{
\includegraphics[width = 0.7\textwidth]{fig/new exp/fig1_stat_traj_more.pdf}
}
\caption{The trajectories of the mean (solid line) and standard deviation (shaded area) of our proposed and Scan $B$ detection statistics over 100 independent trials. Change from $\mathbb{P}=N(\mathbf{0}_{20}, I_{20})$ to Gaussian mixture with $N(\mathbf{0}_{20}, I_{20}) \ \text{\rm w.p. } 0.3$ and $N(\mathbf{0}_{20},4 \ I_{20}) \ \text{\rm w.p. } 0.7$ occurs at $t=51$.
}
\label{fig:detect_stat_illus_more}
\end{figure}

From Figure~\ref{fig:detect_stat_illus_more}, we can observe that our detection statistic has a ``larger slope'' after the change-point, which is a result of the window-limited maximization in the CUSUM-type procedure. This illustrates that our proposed procedure achieves quicker change detection by removing pre-change samples from the post-change block.
Moreover, we can observe that after $ w  = 50$ steps from the change-point, both statistics stop increasing and oscillate around a ``platform'', and they have almost the same value with each other. This validates that choosing such large $ w $ \eqref{eq:Bmax_choice} in our EDD analysis is necessary, otherwise it will take much longer to stop the procedure and leads to overly large EDD.

We need to additionally mention that larger detection statistic does not necessarily induce quicker detection. This is because the window-limited maximization will also enlarge the statistic under $H_0$ and the threshold $b$ for a given ARL in our proposed procedure. To numerically validate this practical merit, i.e., quicker change detection, of our proposed method, we conduct more experiments by comparing their EDDs directly.

\vspace{.1in}
\noindent \textit{EDD comparison.}
As we mentioned before, for a given ARL, the detection threshold $b$ should be chosen accordingly and thus it can be different for different detection procedures. This means Figure~\ref{fig:detect_stat_illus_more} cannot serve as evidence that our proposed procedure outperforms the Scan $B$ procedure. Here, we choose threshold $b$ to meet the ARL constraint independently for each detection procedure via $1000$ Monte Carlo trials. Next, for each procedure, the EDDs are obtained as the average detection delays over another $1000$ Monte Carlo trials with $10000$ reference samples following pre-change distribution $\PP$ and $50$ streaming samples following post-change distribution $\QQ$. We plot the EDD against ARL for those four procedures in Figure~\ref{fig:EDDARL_compare_comprenhensive}.

\begin{figure*}[!htp]
\centerline{
\includegraphics[width = \textwidth]{fig/new exp/fig1_power.pdf}
}
\caption{Comparison of EDD for our proposed (red), Scan $B$ (orange), KCUSUM (blue) and Hotelling's $T^2$ (brown) procedures. The setting is $\PP = N(\mathbf{0}_{20}, I_{20})$ and  $\QQ$ is Gaussian mixture:  $N(\mathbf{0}_{20}, I_{20}) \ \text{\rm w.p. } 0.3; N(\mu \mathbf{1}_{20}, \sigma^2 \ I_{20}) \ \text{\rm w.p. } 0.7$.
}
\label{fig:EDDARL_compare_comprenhensive}
\end{figure*}

In Figure~\ref{fig:EDDARL_compare_comprenhensive}, the absence of the dot indicates the corresponding procedure fails to detect the change before time $t = 50$ under the corresponding setting. We have the following observations:
(i) First and foremost, our proposed method is not only the most powerful, i.e., it achieves the quickest detection, but also the most robust, i.e., it can detect very small changes which KCUSUM and parametric Hotelling's $T^2$ procedures fail to detect, in the experimental settings considered here.
(ii) Secondly, even though KCUSUM achieves quicker detection than the other two benchmark procedures, it easily fails when the change is subtle or the ARL is large, meaning it is not robust enough.
(iii) Lastly, since the mean shift is obvious under the settings considered in Figure~\ref{fig:EDDARL_compare_comprenhensive}, Hotelling's $T^2$ prceodure performs pretty well compared with Scan $B$ procedure.

\begin{table*}
\caption{EDDs for given ARLs. The change is from standard Gaussian distribution  $\mathbb{P}=N(\mathbf{0}_{20}, I_{20})$ to Gaussian mixture distribution.
Non-value ``$-$'' indicates the corresponding detection procedure fails to detect the change within the first 50 time steps under the corresponding setting. Smallest EDDs are highlighted and we can observe that our proposed method outperforms all existing approaches under the following settings.}\label{table:EDD_Gmix}

\vspace{-0.05in}

\begin{center}
\begin{small}
\begin{sc}
\resizebox{.9\textwidth}{!}{%
\begin{tabular}{lccccccccccccccccc}
\multicolumn{16}{c}{\large{Setting I: $N(\mathbf{0}_{20}, I_{20})$ v.s.
Gaussian mixture with $N(\mathbf{0}_{20}, I_{20}) \ \text{\rm w.p. } 0.3; N(\mu \mathbf{1}_{20},\sigma^2  I_{20}) \ \text{\rm w.p. } 0.7$.}} \\
\toprule[1pt]
& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 0.1, \ \sigma^2 = 0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ \sigma^2 = 0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ \sigma^2 = 1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ \sigma^2 = 4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ \sigma^2 = 9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{19.2} & \textbf{19.55} & \textbf{21.57} & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & \textbf{6.99} & \textbf{7.08} & \textbf{7.66} & \textbf{3.47} & \textbf{3.49} & \textbf{3.6} \\
Scan $B$ & 28.1 & 28.7 & 30.83 & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 15.48 & 15.91 & 17.36 & 9.63 & 9.94 & 10.98 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 6.03 & 6.89 & $-$ \\
Hotelling $T^2$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 10.04 & 10.46 & 12.17 \\
\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 0.3, \ \sigma^2 = 0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ \sigma^2 = 0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ \sigma^2 = 1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ \sigma^2 = 4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ \sigma^2 = 9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{12.55} & \textbf{12.77} & \textbf{14.19} & \textbf{17.26} & \textbf{17.53} & \textbf{19.46} & $-$ & $-$ & $-$ & \textbf{6.69} & \textbf{6.77} & \textbf{7.33} & \textbf{3.46} & \textbf{3.47} & \textbf{3.61} \\
Scan $B$ & 22.36 & 22.93 & 24.84 & 24.87 & 25.5 & 27.61 & $-$ & $-$ & $-$ & 14.85 & 15.28 & 16.69 & 9.53 & 9.85 & 10.91 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 6.04 & 6.71 & 8.18 \\
Hotelling $T^2$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 17.34 & 17.88 & $-$ & 8.88 & 9.28 & 10.68 \\
\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 1, \ \sigma^2 = 0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1, \ \sigma^2 = 0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1, \ \sigma^2 = 1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1, \ \sigma^2 = 4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1, \ \sigma^2 = 9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{3.33} & \textbf{3.34} & \textbf{3.45} & \textbf{3.63} & \textbf{3.67} & \textbf{4.05} & \textbf{4.79} & \textbf{4.85} & \textbf{5.26} & \textbf{4.65} & \textbf{4.7} & \textbf{5.15} & \textbf{3.36} & \textbf{3.37} & \textbf{3.46} \\
Scan $B$ & 9.64 & 10.01 & 11.11 & 10.11 & 10.46 & 11.69 & 11.39 & 11.81 & 13.23 & 11.16 & 11.56 & 12.94 & 9.07 & 9.36 & 10.41 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 5.16 & 5.85 & 6.73 \\
Hotelling $T^2$ & 9.25 & 9.56 & 10.53 & 9.12 & 9.36 & 10.39 & 8.72 & 8.96 & 9.91 & 7.13 & 7.42 & 8.34 & 5.52 & 5.71 & 6.43 \\
\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 1.5, \ \sigma^2 = 0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ \sigma^2 = 0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ \sigma^2 = 1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ \sigma^2 = 4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ \sigma^2 = 9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{3.12} & \textbf{3.13} & \textbf{3.18} & \textbf{3.16} & \textbf{3.17} & \textbf{3.21} & \textbf{3.3} & \textbf{3.3} & \textbf{3.39} & \textbf{3.62} & \textbf{3.64} & \textbf{3.86} & \textbf{3.32} & \textbf{3.32} & \textbf{3.38} \\
Scan $B$ & 7.18 & 7.31 & 8.24 & 7.33 & 7.48 & 8.5 & 8.06 & 8.39 & 9.35 & 9.27 & 9.62 & 10.72 & 8.75 & 8.98 & 10.01 \\
KCUSUM & 3.99 & 4.05 & 4.08 & 4.1 & 4.18 & 4.28 & 4.95 & 5.3 & 5.86 & 7.05 & 8.37 & $-$ & 4.78 & 5.07 & 5.62 \\
Hotelling $T^2$ & 5.74 & 5.86 & 6.59 & 5.72 & 5.87 & 6.51 & 5.63 & 5.8 & 6.34 & 5.11 & 5.24 & 5.77 & 4.37 & 4.49 & 4.98 \\
\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 2, \ \sigma^2 = 0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 2, \ \sigma^2 = 0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 2, \ \sigma^2 = 1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 2, \ \sigma^2 = 4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 2, \ \sigma^2 = 9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{2.89} & \textbf{2.89} & \textbf{2.97} & \textbf{2.95} & \textbf{2.96} & \textbf{3.03} & \textbf{3.11} & \textbf{3.12} & \textbf{3.16} & \textbf{3.31} & \textbf{3.33} & \textbf{3.38} & \textbf{3.28} & \textbf{3.29} & \textbf{3.33} \\
Scan $B$ & 5.94 & 6.44 & 7.11 & 6.39 & 6.74 & 7.19 & 7.02 & 7.17 & 7.87 & 8.32 & 8.62 & 9.49 & 8.53 & 8.71 & 9.68 \\
KCUSUM & 3.98 & 3.98 & 4.01 & 3.98 & 4 & 4.01 & 4.03 & 4.03 & 4.08 & 4.77 & 5.07 & 5.56 & 4.4 & 4.59 & 4.89 \\
Hotelling $T^2$ & 4.31 & 4.38 & 5.17 & 4.34 & 4.44 & 5.08 & 4.37 & 4.45 & 4.92 & 4.14 & 4.24 & 4.6 & 3.72 & 3.82 & 4.16 \\
\bottomrule[1pt]
\end{tabular}
}
\end{sc}
\end{small}
\end{center}

\end{table*}

\begin{table*}
\caption{EDDs for given ARLs under more settings. The change is from standard Gaussian distribution  $\mathbb{P}=N(\mathbf{0}_{20}, I_{20})$ to Laplace and Uniform distributions.
}\label{table:EDD_more}

\vspace{-0.05in}

\begin{center}
\begin{small}
\begin{sc}
\resizebox{.9\textwidth}{!}{%
\begin{tabular}{lccccccccccccccccc}
\multicolumn{16}{c}{\large{Setting II: $N(\mathbf{0}_{20}, I_{20})$ v.s.
${\rm Laplace}(\mu \mathbf{1}_{20}, b \mathbf{1}_{20})$.}} \\
\toprule[1pt]
& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 0.1, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.1, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{10.16} & \textbf{10.33} & \textbf{11.33} & \textbf{15.91} & \textbf{16.17} & \textbf{17.73} & $-$ & $-$ & $-$ & \textbf{4.64} & \textbf{4.69} & \textbf{5.01} & \textbf{2.37} & \textbf{2.38} & \textbf{2.54} \\
Scan $B$ & 21.75 & 22.1 & 23.64 & 25.73 & 26.44 & 28.69 & $-$ & $-$ & $-$ & 11.25 & 11.6 & 12.99 & 6.82 & 7.03 & 7.7 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 3.67 & 4.19 & 5.06 \\
Hotelling $T^2$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 7.03 & 7.35 & 8.47 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 0.3, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.3, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{6.43} & \textbf{6.51} & \textbf{7.23} & \textbf{8.92} & \textbf{9.04} & \textbf{9.97} & \textbf{17.9} & \textbf{18.26} & \textbf{20.22} & \textbf{4.4} & \textbf{4.46} & \textbf{4.76} & \textbf{2.34} & \textbf{2.36} & \textbf{2.5} \\
Scan $B$ & 15.88 & 16.21 & 17.92 & 18.01 & 18.56 & 20.48 & 23.82 & 24.42 & 26.67 & 10.53 & 10.9 & 12.22 & 6.72 & 6.96 & 7.59 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 3.62 & 4.12 & 4.83 \\
Hotelling $T^2$ & 24.95 & 25.4 & 27.12 & 24.11 & 24.56 & 26.1 & 20.89 & 21.38 & 23.08 & 11.83 & 12.31 & 13.99 & 6.12 & 6.43 & 7.31 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 0.5, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.5, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.5, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.5, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.5, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{4.15} & \textbf{4.18} & \textbf{4.6} & \textbf{5.08} & \textbf{5.17} & \textbf{5.65} & \textbf{7.51} & \textbf{7.61} & \textbf{8.32} & \textbf{4.05} & \textbf{4.07} & \textbf{4.35} & \textbf{2.3} & \textbf{2.31} & \textbf{2.45} \\
Scan $B$ & 12.07 & 12.56 & 13.49 & 13.06 & 13.41 & 14.49 & 15.1 & 15.51 & 16.97 & 9.58 & 9.98 & 11.19 & 6.62 & 6.86 & 7.48 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 3.45 & 3.97 & 4.65 \\
Hotelling $T^2$ & 14.66 & 14.83 & 15.45 & 14.32 & 14.52 & 15.28 & 13.04 & 13.32 & 14.35 & 8.42 & 8.75 & 9.94 & 5.28 & 5.52 & 6.27 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 0.7, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.7, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.7, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.7, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 0.7, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{3.81} & \textbf{3.86} & \textbf{3.99} & \textbf{3.9} & \textbf{3.93} & \textbf{4.03} & \textbf{4.58} & \textbf{4.63} & \textbf{4.96} & \textbf{3.61} & \textbf{3.64} & \textbf{3.89} & \textbf{2.25} & \textbf{2.25} & \textbf{2.38} \\
Scan $B$ & 8.94 & 9.23 & 10.45 & 9.5 & 9.87 & 11.02 & 10.97 & 11.36 & 12.71 & 8.6 & 8.94 & 10.11 & 6.53 & 6.72 & 7.35 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 3.28 & 3.73 & 4.39 \\
Hotelling $T^2$ & 9.69 & 9.98 & 10.94 & 9.33 & 9.66 & 10.78 & 8.51 & 8.82 & 9.93 & 6.41 & 6.62 & 7.49 & 4.6 & 4.79 & 5.4 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(\mu = 1.5, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(\mu = 1.5, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2.03} & \textbf{2.03} & \textbf{2.06} & \textbf{2.25} & \textbf{2.27} & \textbf{2.43} & \textbf{2.03} & \textbf{2.04} & \textbf{2.07} \\
Scan $B$ & 5 & 5 & 5.77 & 5 & 5.03 & 5.95 & 5.47 & 5.68 & 6.28 & 6.26 & 6.44 & 7.1 & 6.05 & 6.15 & 7 \\
KCUSUM & 2.01 & 2.03 & 2.05 & 2.07 & 2.11 & 2.18 & 2.54 & 2.72 & 2.94 & 3.52 & 3.99 & 4.81 & 2.48 & 2.75 & 3.06 \\
Hotelling $T^2$ & 4.01 & 4.05 & 4.93 & 4.06 & 4.12 & 4.76 & 4.06 & 4.14 & 4.51 & 3.66 & 3.74 & 4.11 & 3.12 & 3.21 & 3.57 \\

\bottomrule[1pt]
\end{tabular}
}
\end{sc}
\end{small}
\end{center}

\begin{center}
\begin{small}
\begin{sc}
\resizebox{.9\textwidth}{!}{%
\begin{tabular}{lccccccccccccccccc}
\multicolumn{16}{c}{\large{Setting III: $N(\mathbf{0}_{20}, I_{20})$ v.s.
$U[(a-b) \mathbf{1}_{20},  (a+b) \mathbf{1}_{20}]$.}} \\
\toprule[1pt]
& \multicolumn{3}{c}{{$\boldsymbol{(a =  0.1, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.1, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.1, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.1, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.1, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{6} & \textbf{6} & \textbf{6.93} & \textbf{5.26} & \textbf{5.39} & \textbf{5.99} & \textbf{4} & \textbf{4} & \textbf{4} & \textbf{2.06} & \textbf{2.07} & \textbf{2.18} & \textbf{2} & \textbf{2} & \textbf{2} \\
Scan $B$ & 16 & 16.02 & 18 & 13.96 & 14.22 & 15.23 & 10.36 & 10.74 & 11.66 & 6.17 & 6.41 & 7.08 & 5.07 & 5.21 & 5.98 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 3.67 & 4.2 & 5.08 & 2.18 & 2.23 & 2.33 \\
Hotelling $T^2$ & 29.1 & 29.46 & 30.77 & 20.28 & 20.61 & 22.21 & 12.3 & 12.65 & 13.41 & 5.19 & 5.31 & 5.95 & 3.95 & 4 & 4.26 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(a =  0.3, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.3, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.3, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.3, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.3, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{4} & \textbf{4} & \textbf{4.51} & \textbf{4} & \textbf{4} & \textbf{4} & \textbf{2.92} & \textbf{2.99} & \textbf{3.52} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} \\
Scan $B$ & 12.63 & 13 & 13.97 & 10.66 & 10.98 & 11.91 & 7.54 & 7.98 & 9.06 & 5.48 & 5.79 & 6.25 & 5 & 5.01 & 5.61 \\
KCUSUM & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & $-$ & 2.4 & 2.55 & 2.83 & 2.05 & 2.07 & 2.1 \\
Hotelling $T^2$ & 15.77 & 15.99 & 16.89 & 13.03 & 13.17 & 14.03 & 7.55 & 7.92 & 9.17 & 4.8 & 4.89 & 5.05 & 3.5 & 3.64 & 3.98 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(a =  0.7, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.7, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.7, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.7, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  0.7, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{2.01} & \textbf{2.02} & \textbf{2.95} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} \\
Scan $B$ & 7 & 7 & 8.1 & 6.28 & 6.96 & 7.03 & 5.67 & 5.98 & 6.27 & 5 & 5 & 5.21 & 4.97 & 5 & 5 \\
KCUSUM & 6.71 & $-$ & $-$ & 3.88 & 4.77 & 6.2 & 2.37 & 2.55 & 2.84 & 2.01 & 2.02 & 2.03 & \textbf{2} & \textbf{2} & \textbf{2} \\
Hotelling $T^2$ & 7 & 7 & 8.07 & 6 & 6.06 & 7 & 5 & 5 & 5.4 & 3.92 & 3.97 & 4.04 & 3 & 3.01 & 3.16 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(a =  1.5, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(a =  1.5, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  1.5, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  1.5, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  1.5, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} \\
Scan $B$ & 5 & 5 & 5 & 5 & 5 & 5 & 4.12 & 4.99 & 5 & 4 & 4.09 & 5 & 4.05 & 4.69 & 5 \\
KCUSUM & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} \\
Hotelling $T^2$ & 4 & 4 & 4 & 4 & 4 & 4 & 3 & 3.02 & 3.97 & 2.96 & 2.99 & 3 & 2.05 & 2.1 & 2.6 \\

\cmidrule(l){1-16}

& \multicolumn{3}{c}{{$\boldsymbol{(a =  2, \ b^2 =  0.1)}$}} & \multicolumn{3}{c}{$\boldsymbol{(a =  2, \ b^2 =  0.3)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  2, \ b^2 =  1)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  2, \ b^2 =  4)}$} & \multicolumn{3}{c}{$\boldsymbol{(a =  2, \ b^2 =  9)}$}\\
ARL & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 & 500 & 1000 & 2000 \\% Column names row
\cmidrule(l){2-4} \cmidrule(l){5-7} \cmidrule(l){8-10} \cmidrule(l){11-13} \cmidrule(l){14-16}
Proposed & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} \\
Scan $B$ & 4 & 4 & 5 & 4 & 4 & 5 & 4 & 4 & 5 & 4 & 4 & 5 & 4.01 & 4.45 & 5 \\
KCUSUM & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} & \textbf{2} \\
Hotelling $T^2$ & 3 & 3 & 3 & 3 & 3 & 3 & 3 & 3 & 3 & 2.01 & 2.04 & 2.74 & \textbf{2} & \textbf{2} & 2.01 \\

\bottomrule[1pt]
\end{tabular}
}
\end{sc}
\end{small}
\end{center}

\end{table*}

To further demonstrate the good performance of our proposed method, we compare it with the above three benchmark methods under more scenarios.
To be precise, we additionaly consider more settings for change from Gaussian to Gaussian mixture as well as change from Gaussian distribution to Laplace and uniform distributions.
The results are reported in Tables~\ref{table:EDD_Gmix} and \ref{table:EDD_more}.
From these two tables, we can observe that our proposed Memory Efficient Online Kernel CUSUM achieves the quickest detection among all benchmarks, even under the settings where the changes are too small for parametric procedures to detect.
In addition to aforementioned observations, we can see (iv) Hotelling's $T^2$ procedure performs pretty well when the mean shift is significant, but fails easily otherwise.
Under those small mean shift circumstances, kernel methods can easily achieve better performance than Hotelling's $T^2$ procedure; see the first two rows in the table for evidence.
This observation agrees with our intuition since Hotelling's $T^2$ procedure is a parametric method which is designed to detect the mean shift.

\subsection{Study on the hyperparameter choices}
In this subsection, we will study how the choices of the hyperparameters of our detection procedure influence the EDD. As mentioned above, the choice of window size $w$ is important and we will begin our numerical study with its effect on the EDD.

\vspace{.1in}
\noindent \textit{Effect of window size.} Recall that we use $[2:w]$ as the search region to locate the change-point. Here, in addition to right end-point $w$, i.e., window size, we also study the effect of left end-point. Here, we use $[B_{\text{\rm min}}: B_{\text{\rm max}}]$ as the region in which we optimize block size parameter $B$; we use $w$ and $B_{\text{\rm max}}$ interchangeably to denote the window size parameter. Intuitively, the wider this search region is, the quicker the detection should be, but this merit comes with higher computation and memory cost. We validate this by comparing EDDs for a given ARL under different $B_{\text{\rm min}}$ and $B_{\text{\rm max}}$ choices and plot the results in Figure~\ref{fig:EDDARL_compare_bminbmax}.

\begin{figure*}[!htp]
\centerline{
\includegraphics[width = \textwidth]{fig/bminbmax/fig3_search_region.pdf}
}
\caption{Comparison of EDD for different search regions $[B_{\text{\rm min}}: B_{\text{\rm max}}]$.
The setting is $\mathbb{P}=N(\mathbf{0}_{20}, I_{20})$, and $\QQ$ is, from left to right, (a) $N(\mathbf{0}_{20},0.01 \ I_{20})$; (b) $N(\mathbf{0}_{20},0.2 \ I_{20})$; (c) $N(\mathbf{0}_{20},9 \ I_{20})$.
}
\label{fig:EDDARL_compare_bminbmax}
\end{figure*}

From Figure~\ref{fig:EDDARL_compare_bminbmax}, we observe that,
on one hand, $B_{\text{\rm max}}$ {\it controls the success of detection} --- the detection procedure typically fails to raise an alarm within a reasonable time if EDD of the oracle procedure \eqref{eq:detect_stat_noWL} exceeds $B_{\text{\rm max}}$. However, this does not tell users to choose an overly large window size, since (i) there is a trade-off between the EDD and the computation and memory cost and (ii) EDD improvement by choosing a larger window size $w$ for a given ARL is no longer significant after $w \sim \log$ ARL, as shown by our non-asymptotic study of the performance loss in Section~\ref{subsec:optimal_Bmax}.
Most importantly, (ii) can be numerically verified by the leftmost panel in Figure~\ref{fig:EDDARL_compare_bminbmax} --- As we can see, for $B_{\text{\rm min}} = 2$ case, the improvement of $B_{\text{\rm max}}$ increasing from 40 to 50 is smaller than that of $B_{\text{\rm max}}$ increasing from 30 to 40. In addition, in the rightmost panel, we can see the improvement of EDD by increasing $B_{\text{\rm max}}$ is marginal.
On the other hand, $B_{\text{\rm min}}$ {\it is important to guarantee quick detection} --- ``power loss'' will occur if $B_{\text{\rm min}}$ is greater than EDD of the optimal procedure \eqref{eq:detect_stat_noWL}; see $B_{\text{\rm min}} = 10$ case in the rightmost panel for evidence.
Otherwise, the effect of $B_{\text{\rm min}}$ is negligible when the change is subtle. Here, a simple and safe choice would be $B_{\text{\rm min}} = 2$.

\begin{figure}[!htp]
\centerline{
\includegraphics[width = 0.7\textwidth]{fig/bandwidth/fig3_bandwidth.pdf}
}
\caption{Comparison of EDD for different kernel bandwidths of Gaussian RBF kernel function. We consider different bandwidth choices $\gamma \cdot {\rm med}_{\PP}$ where $\gamma \in \{0.5,1,10,100\}$, where ${\rm med}_{\PP}$ is chosen via median heuristic. The change is from $\mathbb{P}=N(\mathbf{0}_{20}, I_{20})$ to $\QQ$: (a) top left: $N(\mathbf{0}_{20},9 \ I_{20})$; (b) top right: $N(\mathbf{1}_{20},9 \ I_{20})$;
(c) bottom left: Gaussian mixture: $N(\mathbf{0}_{20}, I_{20}) \ \text{\rm w.p. } 0.3; N(\mathbf{0}_{20},9 \ I_{20}) \ \text{\rm w.p. } 0.7$;
(d) bottom right: Gaussian mixture: $N(\mathbf{0}_{20}, I_{20}) \ \text{\rm w.p. } 0.3; N(\mathbf{1}_{20},9 \ I_{20}) \ \text{\rm w.p. } 0.7$.
}
\label{fig:EDDARL_compare_bandwidth}
\end{figure}

\vspace{.1in}
\noindent \textit{Effect of kernel bandwidth.}
Now, we study the effect of kernel bandwidth on the EDD.
As mentioned in \citet{ramdas2015decreasing}, the median heuristic gives the largest possible MMD (see observation 2 therein) and over- as well as under-estimated kernel bandwidth will result in much faster decay in MMD (see observations 1 and 3 therein). As mentioned in Theorem~\ref{thm:opT_b^wmax} above, the optimal EDD is proportional to the inverse of MMD and therefore median heuristic should result in the quickest detection. We demonstrate this empirically by plotting EDD v.s. $\log$ ARL, where threshold $b$ is selected via Monte Carlo simulation to satisfy the corresponding ARL in Figure~\ref{fig:EDDARL_compare_bandwidth}.

Here, we slightly abuse the notation $\gamma$. The kernel bandwidth is chosen to be $\gamma \cdot {\rm med}_{\PP}$, where ${\rm med}_{\PP}$ denotes the median of the $\ell_2$ distance matrix for pre-change samples: $\gamma < 1$ corresponds to under-estimation of kernel bandwidth whereas $\gamma > 1$ represents over-estimation. Figure~\ref{fig:EDDARL_compare_bandwidth} shows under- and over-estimation of kernel bandwidth would both lead to larger EDD and therefore median heuristic is indeed the best. Therefore, in practice, it is safe to apply the median heuristic to select kernel bandwidth.

\vspace{.1in}
\noindent \textit{Effect of number of pre-change blocks $N$.} Lastly, we study the behavior of detection statistics when choosing different $N$'s.
Intuitively, by the Law of Large Numbers, the larger the $N$ is, the better the performance/convergence of block-version MMD statistic, i.e., $B$-test statistic \citep{zaremba2013b}, will be. This is theoretically demonstrated in Remark~\ref{rmk:Neffect}. However, the question is, empirically speaking, how large should $N$ be to achieve a satisfying performance? To answer this question, we plot the trajectory of the averaged detection statistic (over 50 independent runs) for our procedure in Figure~\ref{fig:N_effect}.

\begin{figure}[!htp]
\centerline{
\includegraphics[width = 0.6\textwidth]{fig/new exp/fig2_effect_N.pdf}
}
\caption{Comparison of mean trajectories over 50 independent trials for different pre-change blocks' numbers $N$. The change is from $\mathbb{P}=N(\mathbf{0}_{20}, I_{20})$ to Gaussian mixture with $N(\mathbf{0}_{20}, I_{20}) \ \text{\rm w.p. } 0.3$ and $N(\mathbf{0}_{20},4 \ I_{20}) \ \text{\rm w.p. } 0.7$ occurs at $t=16$.
}
\label{fig:N_effect}
\end{figure}

As we can observe in Figure~\ref{fig:N_effect}, the slopes for $N \geq 15$ are large enough and almost the same after the change-point. We can focus on the first few steps after the change-point, since (i) we are interested in the quick detection of the change-point and (ii) the convergence of the statistic to its mean is controlled by both $N$ and $B$. By taking the trade-off between the computation time and memory cost and the performance into consideration, we believe $N=15$ suffices for our numerical simulation. However, we need to remark that, in practice, this $N$ choice is typically determined by how many reference samples we have.

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

In this work, we suggest a simple remedy by considering the CUSUM-type online kernel change-point detection procedure to tackle the potential power loss of existing Shewhart chart-type procedures, such as Scan $B$ procedure.
We show that this seemingly simple modification is necessary both theoretically and empirically --- on one hand, it enables analytic approximations to both ARL and EDD, and helps us establish the optimal window size which coincides with classic parametric results; on the other hand, our extensive numerical experiments demonstrate the competitive performance of our proposed procedure by comparing it with various existing procedures.

\newpage
\bibliographystyle{plainnat}

redpajama_set_name: RedPajamaArXiv