Change Detection with the Kernel Cumulative Sum Algorithm

TL;DR

This paper introduces a kernel-based CUSUM algorithm for online change detection that operates effectively without strict distributional assumptions, using MMD for non-parametric comparison, and provides performance bounds.

Contribution

It presents the Kernel CUSUM algorithm, a non-parametric change detection method based on MMD, suitable for settings with abundant background data and less restrictive assumptions.

Findings

Derives bounds on expected delay and false alarm time.

Applicable in non-parametric, real-world data scenarios.

Performs well with large background datasets.

Abstract

Online change detection involves monitoring a stream of data for changes in the statistical properties of incoming observations. A good change detector will detect any changes shortly after they occur, while raising few false alarms. Although there are algorithms with confirmed optimality properties for this task, they rely on the exact specifications of the relevant probability distributions and this limits their practicality. In this work we describe a kernel-based variant of the Cumulative Sum (CUSUM) change detection algorithm that can detect changes under less restrictive assumptions. Instead of using the likelihood ratio, which is a parametric quantity, the Kernel CUSUM (KCUSUM) algorithm compares incoming data with samples from a reference distribution using a statistic based on the Maximum Mean Discrepancy (MMD) non-parametric testing framework. The KCUSUM algorithm is applicable in settings where there is a large amount of background data available and it is desirable to detect a change away from this background setting. Exploiting the random-walk structure of the test statistic, we derive bounds on the performance of the algorithm, including the expected delay and the average time to false alarm.

Figures

Equations

$$ARL2FA={\mathbb{E}}_{\infty }[T].$$

$$\begin{array}{rl}& \mathrm{ESADD}=\underset{1\le t<\infty }{\sup }\mathrm{esssup}{\mathbb{E}}_t[(T-t{)}^{+}\mid {\mathcal{F}}_{t-1}].\end{array}$$

$${ARL2FA}_{\mathrm{CUSUM}}\ge \exp (h).$$

$$\begin{array}{rl}& {\mathrm{ESADD}}_{\mathrm{CUSUM}}\le \\ & \frac{h}{d_{\mathrm{KL}}(p_1,p_0)}+\frac{1}{d_{\mathrm{KL}}(p_1,p_0{)}^2}{\mathbb{E}}_1\left[{\left({\left(\log \frac{f_1(x)}{f_0(x)}\right)}^{+}\right)}^2\right].\end{array}$$

$$k(x,y)=\exp (-\Vert x-y{\Vert }^2/2)$$

$$d_k(p_0,p_1)=$$

$$\sqrt{{\mathbb{E}}_{p_0\times p_0}[k(x,x^{\prime })]+{\mathbb{E}}_{p_1\times p_1}[k(y,y^{\prime })]-2{\mathbb{E}}_{p_0\times p_1}[k(x,y)]}.$$

$$\begin{array}{rl}g((x_0,x_1),(y_0,y_1))=& k(x_0,x_1)+k(y_0,y_1)\\ & -k(x_0,y_1)-k(x_1,y_0).\end{array}$$

$$\begin{array}{rl}{\rho }_{\text{L}}(X,Y)& =\frac{1}{\lfloor n/2\rfloor }\sum _{i=1}^{\lfloor n/2\rfloor }g((x_{2i-1},x_{2i}),(y_{2i-1},y_{2i})).\end{array}$$

$$\begin{array}{rl}{g}_{\delta }((x_0,x_1),(y_0,y_1))=& k(x_0,x_1)+k(y_0,y_1)\\ & -k(x_0,y_1)-k(x_1,y_0)-\delta .\end{array}$$

$$z_n=(x_{2n-1},x_{2n},y_{2n-1},y_{2n}).$$

$$c_n=\inf \left\{k\ge n|\sum _{i=n}^k{g}_{\delta }(z_i)>h\right\}.$$

$${ARL2FA}_{\mathrm{KCUSUM}}\ge \frac{2}{{\mathbb{P}}_{\infty }(c_1<\infty )}.$$

$${\mathrm{ESADD}}_{\mathrm{KCUSUM}}\le 2{\mathbb{E}}_1[c_1].$$

$$c=\inf \{n\ge 1\mid b_n>h\}.$$

$$T_{\text{KCUSUM}}=2c.$$

$$c=\underset{n\ge 1}{\inf }{c}_n.$$

$${\mathbb{E}}_{\infty }[c]\ge \frac{1}{{\mathbb{P}}_{\infty }(c_1<\infty )}.$$

$$z_i\sim \{\begin{array}{ll}{p}_0\times p_0\times p_0\times p_0& \text{ for }1\le i<m,\\ p_1\times p_1\times p_0\times p_0& \text{ for }m\le i.\end{array}$$

$$\begin{array}{rl}& {\mathbb{E}}_{2m-1}[(T-(2m-1){)}^{+}\mid {\mathcal{F}}_{2(m-1)}]\\ & \stackrel{\text{A}}{=}{\mathbb{E}}_{2m-1}[(2c-2m+1{)}^{+}\mid {\mathcal{F}}_{2(m-1)}]\\ & \stackrel{\text{B}}{\le }{\mathbb{E}}_{2m-1}[(2c_m-2m+1{)}^{+}\mid {\mathcal{F}}_{2(m-1)}]\\ & \le 2{\mathbb{E}}_{2m-1}[(c_m-m+1{)}^{+}\mid {\mathcal{F}}_{2(m-1)}]\\ & \stackrel{\text{C}}{=}2{\mathbb{E}}_{2m-1}[(c_m-m+1{)}^{+}]\\ & \stackrel{\text{D}}{=}2{\mathbb{E}}_1[(c_1-1+1{)}^{+}]\\ & =2{\mathbb{E}}_1[c_1].\end{array}$$

$$z_i\sim \{\begin{array}{ll}{p}_0\times p_0\times p_0\times p_0& \text{ for }1\le i<m,\\ p_0\times p_1\times p_0\times p_0& \text{ for }i=m,\\ p_1\times p_1\times p_0\times p_0& \text{ for }m<i.\end{array}$$

$$\begin{array}{rl}& {\mathbb{E}}_{2m}[(T-2m{)}^{+}\mid {\mathcal{F}}_{2m-1}]\\ & ={\mathbb{E}}_{2m}[(2c-2m{)}^{+}\mid {\mathcal{F}}_{2m-1}]\\ & \le 2{\mathbb{E}}_{2m}[(c_{m+1}-m{)}^{+}\mid {\mathcal{F}}_{2m-1}]\\ & =2{\mathbb{E}}_{2m}[(c_{m+1}-m{)}^{+}]\\ & =2{\mathbb{E}}_2[(c_2-1{)}^{+}]\\ & =2{\mathbb{E}}_1[(1+c_1-1{)}^{+}]\\ & =2{\mathbb{E}}_1[c_1].\end{array}$$

$${\mathbb{E}}_t[(T-t{)}^{+}\mid {\mathcal{F}}_{t-1}]\le 2{\mathbb{E}}_1[c_1]$$

$${ARL2FA}_{\mathrm{KCUSUM}}\ge 2\exp \left(\frac{h}{4\Vert k{\Vert }_{\infty }}\log \left(1+\frac{\delta }{4\Vert k{\Vert }_{\infty }}\right)\right).$$

$${\mathrm{ESADD}}_{\mathrm{KCUSUM}}\le \frac{2h}{d_k(p_0,p_1{)}^2-\delta }+\frac{8\Vert k{\Vert }_{\infty }^2}{(d_k(p_0,p_1{)}^2-\delta {)}^2}.$$

$${\mathbb{P}}_{\infty }(c_1<\infty )$$

$$={\mathbb{P}}_{\infty }\left(\underset{k\ge 1}{\sup }\sum _{i=1}^k{g}_{\delta }(z_i)>h\right).$$

$${\mathbb{P}}_{\infty }\left(\underset{k\ge 1}{\sup }\sum _{i=1}^k{g}_{\delta }(z_i)>h\right)\le \exp (-rh),$$

$$M(r)=1-r\delta +{\int }_0^r{\int }_0^{\lambda }{\mathbb{E}}_{\infty }[\exp (u{g}_{\delta }(z))g_{\delta }(z{)}^2]\mathrm{d}u\mathrm{d}\lambda .$$

$$M(r)\le 1-r\delta +16\Vert k{\Vert }_{\infty }^2{\int }_0^r{\int }_0^{\lambda }\exp (u4\Vert k{\Vert }_{\infty })\mathrm{d}u\mathrm{d}\lambda$$

$$=1-r\delta +16\Vert k{\Vert }_{\infty }^2{\int }_0^r\frac{1}{4\Vert k{\Vert }_{\infty }}\left(\exp (\lambda 4\Vert k{\Vert }_{\infty })-1\right)\mathrm{d}\lambda$$

$$=1-r\delta +4\Vert k{\Vert }_{\infty }{\int }_0^r\left(\exp (\lambda 4\Vert k{\Vert }_{\infty })-1\right)\mathrm{d}\lambda .$$

Full text

I Introduction

In this work we are interested in the problem of detecting abrupt changes in streams of data. This could mean detecting a change in the average value of the observations, or a change in variance, or, more generally, finding a change in any other distributional property. In particular we are interested in online change detection, where the algorithm should figure out a change has occurred soon after it happens, without waiting for the entire data stream to be observed. Some examples of where this is relevant include intrusion detection, industrial quality control, and others.

In cases where sufficient prior knowledge of the change is available, there are known optimal algorithms for online change detection. If the probability distributions before and after the change are known, then the CUSUM procedure (shown in Algorithm II.1) is known to be optimal for an objective function that takes into account the magnitude of delays and the frequency of false alarms [11, 14]. Aside from optimality, the CUSUM is also simple to program and has an intuitive interpretation in terms of maximum likelihood. However, there are many situations where the relevant probability distributions can not be modeled precisely, and it would be difficult to use the CUSUM in these cases.

Closely related to change detection is statistical hypothesis testing. In particular, in this work we attempt to leverage tools that have been developed for the problem of two-sample testing and adapt them for use in change detection. Two-sample testing is a non-parametric hypothesis testing task where the goal is to determine if two data sets come from the same distribution. For this problem, an approach based on kernel embeddings has been developed, termed Maximum Mean Discrepancy (MMD) [8]. In the MMD approach, two datasets are compared by computing the distance between the corresponding empirical measures, using a distance which is induced by a positive definite kernel function (See Section III for formal definitions.) Compared to other distance measures, MMD distances are appealing because they admit very simple unbiased estimators. Furthermore, being defined through kernels, the methods are not restricted to Euclidean datasets, and are applicable to hypothesis testing problems involving strings, graphs, and other structured data [15].

Motivated by prior work using MMD in hypothesis testing, we introduce the Kernel Cumulative Sum (KCUSUM) algorithm (Algorithm IV.1). Unlike the CUSUM, the KCUSUM does not require exact specifications of the pre- and post-change distributions. Instead, it relies on a database of samples from the pre-change distribution, and continuously compares incoming observations with samples from the database using a kernel function chosen by the user. In this way, the KCUSUM is able to detect a change to any distribution whose distance from the pre-change distribution is above a user-supplied threshold. Our main theoretical results (Theorem IV.2 and Corollary IV.3) concern the delays and false alarms of the KCUSUM. We derive an upper bound on the time to detect a change (that is, the delay) and a lower bound on the time until a false alarm occurs when there is no change. The analysis builds on existing theory for the CUSUM [12, 11].

The rest of this paper is structured as follows. In Section II we review the basic notions of the CUSUM algorithm and in Section III we review the MMD framework. We introduce the KCUSUM algorithm in Section IV, where we also present the analysis. In Section V we present the results of a numerical evaluation.

II Cumulative Sum algorithm

We consider a sequence of random variables $\{x_{n}\}_{n\geq 1}$ and assume that there is an index $t$ such that for all $1\leq i<t$ the variables $x_{i}$ have the distribution $p_{0}$, and for $i\geq t$, the $x_{i}$ have the distribution $p_{1}$. Presently we assume the $x_{i}$ take values in some Euclidean space $\mathbb{R}^{d}$, although the kernel methods that we shall introduce are not restricted to this scenario. The index $t$ is referred to as the change point. An online change detection algorithm tries to identify this change point in real-time, and bases the decision of whether or not a change has occurred by time $n$ on the data available up to time $n$.

The Cumulative Sum algorithm is an online change detection procedure introduced by Page [16]. For the purposes of introducing the CUSUM, assume that the distributions $p_{0}$ and $p_{1}$ have densities $f_{0}$ and $f_{1}$ respectively. The steps of the CUSUM are presented in Algorithm II.1. At step $n$ of the procedure a data point $x_{n}$ is observed, and the log-likelihood ratio $\log\frac{f_{1}(x_{n})}{f_{0}(x_{n})}$ is calculated and the result added onto the statistic $Z_{n}$. If the result would be negative then $Z_{n}$ is set to zero, effectively restarting the algorithm. If $Z_{n}$ crosses a threshold $h$ then a change is declared at time $n$. Formally, the CUSUM stopping time is $T_{\textrm{CUSUM}}=\inf\{n\geq 1\mid Z_{n}\geq h\}$.

For some insight into why the CUSUM works, consider the behavior of the log-likelihood ratio before and after the change. Before the change, it has mean $-d_{KL}(p_{0},p_{1})$, where $d_{\mathrm{KL}}(p_{0},p_{1})=\mathbb{E}_{p_{0}}\left[\log\frac{f_{0}(x)}{f_{1}(x)}\right]$ is the Kullback-Liebler divergence between the distributions $p_{0}$ and $p_{1}$. Since $d_{KL}$ is positive when $p_{0}\neq p_{1}$, the increment term will have a negative mean. This drift combined with the barrier at zero causes the statistic to stay near zero before the change. After the change, the increment has a positive mean equal to $d_{KL}(p_{1},p_{0})$ and $Z_{n}$ begins to increase, eventually crossing any positive threshold $h$ with probability one, which will cause the algorithm to end. Beyond these heuristic arguments, the CUSUM can be shown to be optimal in a certain sense, as we review below.

Example II.1

Denote by $\mathcal{N}(a,b)$ the normal distribution with mean $a$ and variance $b$. Consider detecting a change in variance of normally distributed random variables, where the pre-change distribution is $\mathcal{N}(1,1)$ and the post-change distribution is $\mathcal{N}(1,4)$. A sample sequence $x_{1},\ldots,x_{n}$ of length $n=400$ with a change point at $t=200$ is shown on the left in Fig. 1. The true change time is marked by a dashed line. The log likelihood ratio in this case is $\log\frac{f_{1}(x)}{f_{0}(x)}=\frac{3}{8}x^{2}-\frac{3}{4}x+\log\frac{1}{2}+\frac{3}{8}$. The values of $\log\frac{f_{1}(x_{n})}{f_{0}(x_{n})}$ for $n=1,2,\ldots$ are plotted in the middle of Fig. 1. We can see that the log likelihood ratio has a negative mean before the change and a positive mean after the change. The resulting CUSUM statistic $Z_{n}$ is shown on the right of Fig. 1. Using a threshold of $h=10$ results in detection at time $T_{\textrm{CUSUM}}=212$. $\blacksquare$

Next we review the performance characteristics of the CUSUM. Each possible change time defines a different distribution on the sequences $\{x_{n}\}_{n\geq 1}$. If there is no change, then the variables are independent and identically distributed (i.i.d.) with $x_{i}\sim p_{0}$ for all $i\geq 1$. We denote this distribution on sequences by $\mathbb{P}_{\infty}$, and denote expectations with respect to this distribution by $\mathbb{E}_{\infty}$. In general, a change at time $t\geq 1$ means that for $1\leq i\leq t-1$ the $x_{i}$ are i.i.d. with $x_{i}\sim p_{0}$ and for $i\geq t$ they are i.i.d. with $x_{i}\sim p_{1}$. We let $\mathbb{P}_{t}$ denote the probability distribution on sequences under the assumption of a change at time $t$, and $\mathbb{E}_{t}$ represents the expectation under this distribution. For $n\geq 1$ let $\mathcal{F}_{n}$ be the $\sigma$-algebra $\mathcal{F}_{n}=\sigma(x_{1},x_{2},\ldots,x_{n})$. Intuitively, $\mathcal{F}_{n}$ represents the information contained in the observations up to and including time $n$. Formally, an online change detector can be represented as a stopping time with respect to the filtration $\{\mathcal{F}_{n}\}_{n\geq 1}$, that we denote by $T$, with the interpretation that the value of $T$ is an estimate of the change point.

When running a change detector on a particular sequence, two types of errors may occur. There may be a false alarm, which means the change is detected too early, or there may be a delay, meaning the change was detected late. We formalize the levels of false alarm and delay using the metrics of average run length to false alarm ($\operatorname{ARL2FA}$) and worst case average detection delay ($\operatorname{ESADD}$). These are standard metrics for evaluating change detectors [11, 19, 13].

For a change detector $T$, the average time to false alarm is

[TABLE]

That is, the $\operatorname{ARL2FA}$ is the average amount of time until a change is detected given a sequence of observations with no change.

We measure delay using Lorden’s criterion [11]. If $T$ is a change detector and there is a change at time $t\geq 1$, then the expected delay given the history of the observations up to time $t-1$ is the random variable $\mathbb{E}_{t}[(T-t)^{+}\mid\mathcal{F}_{t-1}].$ 111The notation $(\cdot)^{+}$ refers to the positive part function: $(x)^{+}=\max\{0,x\}$. The worst case delay for a change at time $t$ is obtained by taking an essential supremum over all possible sequences of length $t-1$, denoted by $\operatorname{esssup}\mathbb{E}_{t}[(T-t)^{+}\mid\mathcal{F}_{t-1}]$. Finally, taking the supremum over all change times $t$ we obtain the worst case delay:

[TABLE]

Notably, the CUSUM provides the optimal trade off between the time to false alarm and the worst case delay. This was first proved in an asymptotic form in [11], and the result was later proved in full non-asymptotic form in [14]. A proof of optimality can be found in [14, 19]. Further information on the derivation and properties of the CUSUM may be found in [2], or [19].

The precise relation between the threshold $h$ and the performance levels ARL2FA and ESADD is non-trivial and involves solving numerically intractable integral equations [16, 13]. However, it is possible to derive some upper and lower bounds that may be useful in practice. For the sake of comparison with the analysis of Kernel CUSUM, it will interesting to consider the following quantitative bounds on the performance of the CUSUM.

Proposition II.2

The performance of the $\mathrm{CUSUM}$ (Algorithm II.1) can be bounded as follows. The time to false alarm satisfies

[TABLE]

If it also holds that \mathbb{E}_{1}\left[\big{(}\big{(}\log\frac{f_{1}(x)}{f_{0}(x)}\big{)}^{+}\big{)}^{2}\right]<\infty then

[TABLE]

Proof:

See the appendix. ∎

The intuitive interpretation of these equations is that increasing the threshold $h$ causes an increase the time until false alarm, but it also leads to increased detection delay. From the second equation, we can see that the detection delay increases when the distributions get closer. The term involving the positive part of the log-likelihood ratio is related to the variance of the CUSUM statistic. In Corollary IV.3 below we shall obtain somewhat analogous bounds for the Kernel CUSUM.

III Maximum Mean Discrepancy

Two-sample testing refers to the problem of determining whether two datasets are drawn from the same distribution. One approach to this problem is to consider the empirical measures defined by the datasets, and to compute the distance between the empirical measures using a probability metric. If enough data points are used, then the empirical distance should be close to the true distance. If the distance is large then we can be confident that the datasets are generated by different distributions. This is the idea underlying several classical tests, such as the Kolmogorov-Smirnov test [4], the Cramer-von-Mises test [1] and the Anderson-Darling test [17].

The Maximum Mean Discrepancy (MMD) approach is also based on computing the distance between empirical distributions. In MMD, the datasets are implicitly embedded in a Reproducing Kernel Hilbert Space (RKHS) corresponding to a user-supplied kernel function $k$, and the distance between the embeddings is computed [8]. Compared to classical approaches, there are several features of MMD that make it appealing for non-parametric statistics. First, the MMD distance has a range of simple unbiased estimators (see the definition of $\rho_{L}$ below for one such example.) Second, there is the flexibility offered by choice of kernel, and using kernels means the test can be applied on datasets without a natural Euclidean representation, such as strings, graphs and other structured data [7, 21].

Let $\mathcal{X}$ be a set, and let $k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ be a kernel on this set; this is a symmetric, positive definite function that we regard intuitively as a similarity measure222Symmetric means that $k(x,y)=k(y,x)$ and positive-definite means that for any choice of $n$ elements $x_{1},\ldots,x_{n}$ and $n$ real numbers $a_{1},\ldots,a_{n}$, we have $\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}a_{i}a_{j}k(x_{i},x_{j})\geq 0$. The reader may have in mind the set $\mathcal{X}=\mathbb{R}^{n}$ and the Gaussian kernel

[TABLE]

Some other choices for the kernel function may be found in [15, Table 3.1]. Further suppose that $\mathcal{X}$ has the structure of a measurable space $(\mathcal{X},\Sigma)$ and that $k$ is a measurable function on $\mathcal{X}\times\mathcal{X}$ with the product $\sigma$-algebra. Define $\mathcal{P}(\mathcal{X})$ to be the set of all probability measures on $(\mathcal{X},\Sigma)$, and using the kernel $k$, define the subset $\mathcal{P}_{k}=\{\mu\in\mathcal{P}(\mathcal{X})\mid\mathbb{E}_{x\sim\mu}[\sqrt{k(x,x)}]<\infty\}$. If the kernel has the additional property of being characteristic333We refer the reader to [8] or [15] for a precise definition of characteristic kernel. For instance, if $\mathcal{X}=\mathbb{R}^{n}$ then the Gaussian kernel is characteristic. then we may define the MMD metric on $\mathcal{P}_{k}(\mathcal{X})$, denoted $d_{k}$. This metric is defined as

[TABLE]

See [8] for more details.

One of the unbiased estimators of $d_{k}^{2}$ presented in [8] is the linear statistic $\rho_{\textrm{L}}$ defined below. For convenience, the estimator is expressed using the following function $g$:

[TABLE]

Consider two data sets $X=\{x_{1},x_{2},\ldots,x_{n}\}\subseteq\mathcal{X}$ and $Y=\{y_{1},y_{2},\ldots,y_{n}\}\subseteq\mathcal{X}$. Then $\rho_{\textrm{L}}(X,Y)$ is

[TABLE]

The linear statistic is interesting as it is a sum of i.i.d. terms, which means the central limit theorem may be used to approximate its distributional properties, which can help in tuning the thresholds when MMD is used for hypothesis testing [8]. Furthermore, in the online setting we can study the trajectory of the statistic (that is, as a function of $n$) using the theory of sums of i.i.d. random variables, greatly facilitating analysis. We will exploit this structure in the analysis of KCUSUM below.

III-A *Other related work *

The field of online change detection has its roots in industrial quality control [20], and sequential analysis [22]. In particular, the CUSUM procedure is closely related to the sequential probability ratio test (SPRT), which is a foundational online hypothesis testing algorithm [22].

An alternative to online change detection is offline detection, where the algorithms do not run until the entire sequence is observed, and all data is used to make a decision. Kernel offline change point analysis was explored in [9]. Each hypothetical change time is used to partition the dataset into two groups, consisting of prior observations and later observations, and the two resulting datasets are compared using a kernel based test. This is repeated for each possible change time. If one of the comparisons yields a significant discrepancy between the pre- and post- observation datasets, then a change is declared at the point where the difference was largest. This procedure is of interest as it does not use reference data, instead basing its decisions on comparisons between disjoint sets of observations.

A number of approaches to non-parametric online change detection have been proposed in [3]. As given, they apply only to the case of detecting changes in mean. One straightforward method of applying kernel non-parametric tests to online change detection is a sliding window approach, as pursued in [10]. In this approach, at each time a decision regarding the change is made based on the distance between the most recent fixed-size block of data and a block from the reference distribution, using an MMD distance $d_{k}$. This can be seen as an kernel-based generalization of Shewhart control charts [20].

Besides the CUSUM, there are other algorithms with optimality properties, notably the Shiryaev-Roberts (SR) change detector [18]. While the CUSUM minimizes the worst case delay (ESADD), the SR detector minimizes a form of average delay. Like the CUSUM, the SR test statistic admits a simple recursive form and hence it may be possible to extend our algorithmic construction to SR-type detectors as well, but this is outside the scope of this paper.

IV The Kernel CUSUM Algorithm

The Kernel CUSUM (KCUSUM) algorithm blends features of the CUSUM procedure with the MMD framework. The basic idea is that instead of using each new observation $x_{n}$ to compute the log-likelihood ratio $\log\tfrac{f_{1}(x_{n})}{f_{0}(x_{n})}$, which is an estimator of the KL-divergence $d_{KL}$, we will use the new observation and some other random samples to compute an estimate of an MMD distance $d_{k}$.

The KCUSUM algorithm is defined with the help of a shifted version of the function $g$ used to define the MMD statistics: Given $\delta>0$ (the role of $\delta$ is explained in detail below), define $g_{\delta}$ as

[TABLE]

The details of the Kernel CUSUM are listed in Algorithm IV.1. It is assumed that the change to be detected is from a reference distribution $p_{0}$ to an unknown distribution $p_{1}$. At even numbered iterations, the most recent observation $x_{n}$ is paired with the previous observation $x_{n-1}$ and these two points are compared with two reference points $y_{n},y_{n-1}$ using MMD. We subtract a constant $\delta>0$ from the result to get $v_{n}$. The variable $v_{n}$ is then added onto the statistic $Z_{n-1}$ to get the next value $Z_{n}$. If the new value $Z_{n}$ would be negative, then it resets to [math], effectively restarting the algorithm. At odd numbered iterations, the statistic $Z_{n}$ is unchanged. Note that as a consequence, the algorithm only raises alarms at even numbered iterations.

The reason for subtracting a positive amount $\delta$ at each step of Algorithm IV.1 is to guarantee that the increments $v_{n}$ have negative drift under the pre-change regime and positive drift in the post-change regime. This is a consistency property that enables us to formulate non-trivial bounds on performance. Using the definitions in Algorithm IV.1, it is evident that before a change, $\mathbb{E}[v_{n}]=-\delta<0$ and after a change, $\mathbb{E}[v_{n}]=d_{k}^{2}(p_{0},p_{1})-\delta$. In other words, the algorithm can detect a change to any distribution $p_{1}$ that is at least distance $\sqrt{\delta}$ away from the reference distribution $p_{0}$.

Example IV.1

Consider the problem of detecting a change in variance of normally distributed random variables, as in in Example II.1. The upper left plot of Fig. 2 shows the stream of observations, which are normally distributed with a change in variance at time $t=200$. The KCUSUM is based on the linear statistic $\rho_{L}$, so at each time $n$ the MMD estimate $k(x_{n-1},x_{n})+k(y_{n-1},y_{n})-k(x_{n},y_{n-1})-k(x_{n-1},y_{n})$ is computed, and this quantity is plotted in the middle of Fig. 2. On the right is the resulting $Z_{n}$ sequence. As in the CUSUM, a simple threshold is used to decide that a change has occurred. For this particular realization of the variables, a threshold of $h=5$ results in detection at $T_{\textrm{KCUSUM}}=225$. The value of $\delta$ was $\delta=1/40$. $\blacksquare$

Next we consider the delay and false alarm rate of the KCUSUM. The time to false alarm is defined as in Equation (1). The worst case delay is defined as in Equation (2), using the filtration $\{\mathcal{F}_{n}\}$ where $\mathcal{F}_{n}=\sigma(x_{1},y_{1},\ldots,x_{n},y_{n})$. To prove the bounds we adapt the methods of [11], which allows one to reduce the problem of analyzing the CUSUM to that of analyzing a random walk with i.i.d. terms.

For the analysis it will be convenient to group the variables together as $\{z_{n}\}_{n\geq 1}$ where for $n\geq 1$,

[TABLE]

The grouping reflects how the algorithm process the data in blocks of two. As a consequence, the bounds in this section involve additional factors of two compared the CUSUM (see Proposition II.2).

Associated to the grouped sequence $\{z_{n}\}_{n\geq 1}$, define auxiliary stopping times $c_{1},c_{2},\ldots$ as

[TABLE]

Theorem IV.2

Let $T_{\textrm{KCUSUM}}$ be the change detector corresponding to the Kernel CUSUM (Algorithm IV.1). Let $p_{0}$ be the pre-change, or reference distribution. Then by using a threshold $\delta>0$, the time to false alarm is at least

[TABLE]

If $p_{1}$ is a distribution with $d_{k}(p_{0},p_{1})>\sqrt{\delta}$, then the worst case detection delay is at most

[TABLE]

Proof:

Let $b_{0}=0$ and for $n\geq 1$ let $b_{n}=\max\{0,b_{n-1}+g_{\delta}(z_{n})\}.$ Define the stopping time $c$ as

[TABLE]

The relation between $c$ and the KCUSUM stopping time $T_{\textrm{KCUSUM}}$ is

[TABLE]

Note that, as discussed in [11], the stopping time $c$ can be represented as

[TABLE]

and each stopping time $c_{n}$ uses the same decision rule, the only difference being that they operate on shifted versions of the input sequence $\{x_{i}\}_{i\geq n}$. In this setting Theorem 2 from [11] is applicable, which yields a lower bound on $c$ under $\mathbb{P}_{\infty}$:

[TABLE]

Combining Equations (11) and (13) yields the claim (8).

Now we consider bounding the worst case delay. If the sequence $\{x_{i}\}_{i\geq 1}$ has a change at an odd valued time, say $t=2m-1$ for $m\geq 1$, then the sequence $\{z_{n}\}_{n\geq 1}$ has a change at time $m$. Explicitly,

[TABLE]

From here we reason as in Theorem 2 of [11]:

[TABLE]

Step A follows from Equation (11) and Step B follows since $c$ is the infimum of the $\{c_{n}\}_{n\geq 1}$. Step C follows from the independence of $c_{m}$ from $\mathcal{F}_{2(m-1)}$ and finally Step D follows from the fact that the distribution of $c_{m}-m$ under $t=2m-1$ is the same as the distribution of $c_{1}-1$ under $t=1$.

The situation is slightly more complex if the change occurs a time $t$ that is even, say $t=2m$ for some $m\geq 1$. In this case the grouped sequence $\{z_{n}\}_{n\geq 1}$ does not experience an abrupt change, and instead there are three relevant distributions. Specifically,

[TABLE]

Reasoning as in Equation (14), then,

[TABLE]

Combining (14) and (15), we see that for all $t\geq 1$,

[TABLE]

Combining Equation (16) with the definition of worst case delay (2) yields the inequality (9). ∎

We can combine Theorem IV.2 with certain facts about random walks (Lemma A.2 and Proposition A.1) to get more specific bounds on the delays and false alarms, as shown in the following Corollary. In this corollary, we assume that the kernel $k$ is bounded by a constant $\|k\|_{\infty}$, and also assume that $\delta$ is bounded by $\delta<2\|k\|_{\infty}$. This is a necessary assumption when the kernel is bounded, since if $\delta\geq 2\|k\|_{\infty}$ then $g_{\delta}(z)\leq 0$ for all $z$, and it will not be possible to detect any changes.

Corollary IV.3

Let the assumptions of Theorem IV.2 hold. Further assume that the kernel $k$ is bounded by a constant $\|k\|_{\infty}$ and let $\delta<2\|k\|_{\infty}$. Then the time to false alarm satisfies

[TABLE]

If $p_{1}$ is a distribution with $d_{k}(p_{0},p_{1})>\sqrt{\delta}$, then the worst case detection delay is at most

[TABLE]

Proof:

To start, note that

[TABLE]

To upper-bound the last term in this equation, we will apply Lemma A.2. Note that under our assumption that the kernel is bounded, the moment generating function $M(r)=\mathbb{E}_{\infty}[\exp(rg_{\delta}(z_{1}))]$ is guaranteed to be well defined for all $r\geq 0$. Therefore

[TABLE]

where $r$ is any number satisfying $r>0$ and $M(r)\leq 1$. To identify such an $r$, start with a second order expansion of $M$:

[TABLE]

Under the assumption that $\delta<2\|k\|_{\infty}$ it holds that $|g_{\delta}(z)|\leq 4\|k\|_{\infty}$ and

[TABLE]

Minimizing the right hand side of the final equation above with respect to $r$ yields

[TABLE]

Combining (8), (IV) with (20) and using this definition of $r$ yields the claim (17).

For the delay, note that $\mathbb{E}_{1}[c_{1}]$ is the expected amount of time until a random walk with positive drift crosses an upper boundary. Hence we may apply Proposition A.1. This leads to

[TABLE]

Combining (9), (21), and the bound $g_{\delta}(z)^{+}\leq 2\|k\|_{\infty}$, we obtain the relation (18).

∎

Figure 3 shows the logarithmic relation between false alarm time and delay specified by the theorem. For each level of false alarm $x\in\{1,2,\ldots,10^{4}\}$, we computed the smallest value of $h$ guaranteed to achieve false alarm rate $x$ according to Equation (17), and plug in this threshold to compute delay according to Equation (18). The computations were performed for a hypothetical problem where $d_{k}(p_{0},p_{1})^{2}=1/6,\delta=2^{-5}$ and $\|k\|_{\infty}=0.5$.

V Empirical Results

In this section we evaluate the KCUSUM on several change detection tasks. In each case, the observations consisted of vectors in $\mathbb{R}^{4}$, and the Gaussian kernel (defined in Equation 3) was used with $\sigma^{2}=1$. The pre-change distribution in each task is the normal distribution on $\mathbb{R}^{4}$ with mean zero and a covariance matrix equal to the identity scaled by a factor of $\frac{1}{2}$. The four possible post-change distributions were as follows:

Change in mean: A normal distribution with mean $(1,1,1,1)$ and the pre-change covariance matrix. 2. 2.

Change in variance (all components): A normal distribution with mean [math] and a covariance matrix equal to the identity scaled by a factor of $2$. 3. 3.

Change in variance (random component): The distribution obtained by sampling from the pre-change distribution and scaling a random component by a factor of $2$. 4. 4.

Change to uniform: The distribution on $\mathbb{R}^{4}$ where each component is sampled independently from the uniform distribution on $[-1/(2\sqrt{3}),1/(2\sqrt{3})]$

Note that the interval in Problem 4 was chosen so that the resulting distribution has the same mean and variance as the pre-change distribution.

For each task we used a Monte Carlo approach to estimate the time to false alarm and delay. The time to false alarm was estimated by generating $n=5000$ sequences with no change and running the KCUSUM until a false alarm was detected. We record the time where the false alarm occurs and average theses values to get the estimate for ARL2FA. For the delay, we generated $n=5000$ sequences that had a change at time $t=1$, and recorded the amount of time until the alarm goes off as an estimate of the delay. In the examples, we set $\delta=2^{-7}$ in tasks (1) - (3), and set it $\delta=2^{-9}$ in task (4).

Based on the results in Figure 4 we see that the change detection tasks increase in difficulty as we go from problem 1 to 4. In all of the problems we observed a similar logarithmic growth rate in the the delay as the time to false alarm is allowed to increase. However, the scale of this growth can vary dramatically. For instance, the first two problems are relatively easy, while the fourth problem seems to be quite difficult for the KCUSUM. Note that in each case, the only data used by the algorithm is the incoming observations and samples from the pre-change distribution, and no information or samples about the post-change distribution are used. Overall, the results suggest that the KCUSUM may be a promising approach for change detection problems where less is known about the type of change.

VI Conclusion

This work introduced the Kernel Cumulative Sum Algorithm (KCUSUM), a new approach for online change detection. Unlike the CUSUM algorithm, this approach does not require knowledge of the probability density ratio for its implementation. Instead, it uses incoming observations and samples from the pre-change distribution. The result is that the same algorithm works for detecting many types of changes. Our theoretical analysis establishes the algorithm’s ability to detect changes, and shows a relation between the delay and the MMD distance of the two distributions. These bounds should also be useful in the analysis of other non-parametric change detectors. Finally, we would like to suggest two avenues for future work. First, there are likely variants of KCUSUM that leverage more complex, non i.i.d. statistics, that may lead to improved detection performance. Secondly, the CUSUM has been investigated for detecting changes in scenarios with more complex dependencies among observations [5] and it may also be possible to extend the kernel methods developed in this paper to detect changes in these cases.

Appendix A Appendix

Proposition A.1** (Corollary 1, [12])**

Let $\{a_{i}\}_{i\geq 1}$ be i.i.d. real-valued random variables such that $\mathbb{E}[a_{i}]=\mu>0$ and $\mathbb{E}[(a_{i}^{+})^{2}]<\infty$. Define $S_{n}=\sum_{i=1}^{n}a_{i}$ and for $a\leq 0\leq b$ let $T$ be the stopping time $T=\inf\{n\geq 1\mid S_{n}\notin[a,b]\}$. Then

[TABLE]

where $\alpha=\mathbb{P}(S_{T}<a)$.

Lemma A.2

Let $\{a_{i}\}_{i\geq 1}$ be a sequence of i.i.d. real-valued random variables and for $n\geq 1$ define the partial sums $S_{n}=\sum_{i=1}^{n}a_{i}$. Let $M(r)=\mathbb{E}\left[\exp(ra_{1})\right]$ be the moment generating function of the $a_{i}$ and suppose there is a $q>0$ such that $M(q)\leq 1$. Then for any $h\geq 0$,

[TABLE]

Proof:

For $n\geq 1$ define $Z_{n}=\exp(qS_{n})$. Then for $n\geq 1$,

[TABLE]

Furthermore, for all $n\geq 1$ it holds that $Z_{n}\geq 0$ and $\mathbb{E}[|Z_{n}|]=M(q)^{n}\leq 1$. Therefore $Z_{n}$ is a non-negative supermartingale. Hence

[TABLE]

The first step in the above derivation follows from the monotonicity of the function $x\mapsto\exp(x)$, the second step follows from Theorem 7.8 in [6], and the final step follows from our assumption on $q$. ∎

Proof of Proposition II.2

Define $S_{1},S_{2},\ldots$ as

[TABLE]

and let $\alpha=\mathbb{P}_{\infty}(\sup_{n\geq 1}S_{n}>h)$. It follows from Theorem 2 of [11] that the CUSUM obeys

[TABLE]

Note that $\mathbb{E}_{\infty}[\log\frac{f_{1}(x_{i})}{f_{0}(x_{i})}]=-d_{KL}(p_{0},p_{1})<0$. Hence $S_{n}$ is a random walk with negative drift. Furthermore, the moment generating function of the increments under $\mathbb{P}_{\infty}$ can be expressed as $M(r)=\mathbb{E}_{\infty}\left[\exp\left(r\log\frac{f_{1}(x_{1})}{f_{0}(x_{1})}\right)\right]=\mathbb{E}_{\infty}\left[\left(\frac{f_{1}(x_{1})}{f_{0}(x_{1})}\right)^{r}\right]$ and it is evident that $g(1)=\mathbb{E}_{1}[1]=1$. Then we may apply Lemma A.2 with $q=1$ to conclude that $\alpha\leq\exp(-h)$. This establishes our claim on the time to false alarm.

Let $N$ be the stopping time $N=\inf\{n\geq 1\mid S_{n}>h\}$. Again applying Theorem 2 from [11], it holds that

[TABLE]

Under the assumption of a change point at $t=1$, the variables $\log\frac{f_{1}(x_{i})}{f_{0}(x_{i})};i=1,2,\ldots$ are i.i.d with positive mean, and therefore $\{S_{n}\}_{n\geq 1}$ is a random walk with positive drift $\mu=\mathbb{E}_{1}[\log\frac{f_{1}(x_{1})}{f_{0}(x_{1})}]$. The bound on ESADD then follows from applying Proposition A.1 to the stopping time $N$.

$\blacksquare$