Separating an Outlier from a Change

TL;DR

This paper introduces a method for change detection that distinguishes true regime shifts from outliers in data with unknown post-change distributions, validated on economic and climate datasets.

Contribution

It characterizes likely outlier events and proposes a distribution testing approach, improving detection accuracy over existing methods like FMA and GLRT.

Findings

Successfully detects regime shifts in market and climate data

Outperforms traditional tests in various performance metrics

Identifies climate change as a significant regime shift

Abstract

We study the change detection problem with an unknown post-change distribution. Under this constraint, the unknown change in the distribution of observations may occur in many ways without much structure on the observations, whereas, before the change point, a false alarm (outlier) is highly structured, following a particular sample path. We first characterize these likely events for the deviation and propose a method to test the empirical distribution, relative to the most likely way for it to occur as an outlier. We benchmark our method with finite moving average (FMA) and generalized likelihood ratio tests (GLRT) under 4 different performance criteria including the run time time complexity. Finally, we apply our method on economic market indicators and climate data. Our method successfully captures the regime shifts during times of historical significance for the markets and identifies the current climate change phenomenon to be a highly likely regime shift rather than a random event.

Tables

Table 1. TABLE I: Time complexity scaling

	IPT	FMA	GLRT
init.	$𝒪\left(n+m\log {ϵ}^{-1}\right)$	$\Theta \left(1\right)$	$\Theta \left(1\right)$
rolling	$\Theta \left(d+m\right)$	$\Theta \left(d\right)$	$𝒪\left(n+m\log m\log {ϵ}^{-1}\right)$
sliding	$\Theta \left(1\right)$	$\Theta \left(1\right)$	$𝒪\left(n+m\log m{\log \left(dϵ\right)}^{-1}\right)$

Equations

$$\underset{\Gamma }{\inf }\underset{f_1}{\sup }$$

$$P_{\infty }\left({\hat{f}}_{X_1^n}\in {\Gamma }_1\right)\le \alpha$$

$$\underset{t_a}{\inf }$$

$$\underset{k}{\sup }{P}_{\infty }\left(k\le t_a<k+n_{\alpha }\right)\le \tilde{\alpha }$$

$$\underset{t_a}{\inf }$$

$$E_{\infty }{t}_a\ge \gamma .$$

$$\underset{n\to \infty }{\lim }-\frac{1}{n}\log P_{\infty }\left(q\left({\hat{f}}_{X_1^n}\right)\ge c\right)$$

$$P_{\infty }\left({\hat{f}}_{X_1^n}\in {\Gamma }_1\right)$$

$$\underset{f_1}{\sup }{P}_{f_1,1}\left({\hat{f}}_{X_1^n}\in {\Gamma }_0\right)$$

$$\underset{k}{\sup }{P}_{\infty }\left(k\le t_I<k+n_{\alpha }\right)$$

$$\underset{f_1,t_1}{\sup }{P}_{f_1,t_1}\left(t_I\ge t_1+n|t_I\ge t_1\right)$$

$$\nu$$

$$\eta$$

$$f_n^{\ast }$$

$$(S_k,i_k)$$

$$D_k$$

$${\tau }_{k+1}$$

$$t_I$$

$$E_{\infty }{t}_I$$

$$t_I$$

$$t_S^{(j)}$$

$$t_I$$

$$E_{f_1,t_1}\left({\left(t_I-t_1+1\right)}^{+}|X_1^{t_1-1}\right)$$

$$E_{f_1,1}{t}_I$$

$$WADD(t_I)$$

$$WADD(t_I)=\frac{\log \gamma }{\underline{q}\left(\frac{v^{\ast }}{2}+\frac{|q_0|}{L^2}\right)}.$$

$$q\left({\hat{f}}_{X_i^k}\right)$$

$$S_k$$

$$=\{\begin{array}{ll}{\left(S_{k-1}+X_k\right)}^{+},& {\tau }_k\ne k\\ {\left(X_k\right)}^{+},& {\tau }_k=k\end{array}$$

$$f_n^{\ast }$$

$$f_n^{\ast }(a)$$

$$D_k$$

$$=I\left(\hat{f}\Vert f_0\right)+\Lambda (r)-r\sum _{a}a\hat{f}(a)$$

$$\Lambda (r)$$

$$E_{\infty }\left(\exp \left(rX-\Lambda \left(r\right)\right)X\right)$$

$$\text{norm}{f-f_1}_1$$

$$q(f)$$

Full text

I Introduction

Not every long-term deviation from the norm should be considered as a result of a change. Such occurrences may also be caused by rare events driven by the system and this distinction is important in many applications from model selection [1, 2] to hardware faults [3, 4], security [5, 6] and even health [7, 8]. To that end, this paper focuses on problem instances of change detection [9] (also see [10] for a wide treatment of change detection and refer to chapter 8 for nonadditive change models and existing approaches) in a time-series with minimax cost with respect to change point and post-change distributions.

The quickest change detection problem with unknown post-change models has been studied widely, both in theory [11, 12, 13, 14, 15] and in application [16, 17] (see [18] for more variants and applications). For some non-Bayesian quickest change detection problems where pre- and post-change distributions have finitely many alternatives it is known that a version of the cumulative sum algorithm (CUSUM) is optimal or asymptotically optimal [19, 20, 21]. For other problems where the pre- or post-change distribution has infinitely many alternatives or is not parametrized, algorithms that form a likelihood ratio (ex. CUSUM) are not applicable. In this case, since the change in the distribution is arbitrary it is hard to perform optimally against all alternatives. A widely used method in the literature is to form a generalized likelihood ratio statistics over the set of alternative distributions and then use it as a drift term [14, 22, 20], which can be computationally demanding or without any optimality guarantees. To address this, we propose a recursive 111if the quasiconcave function $q$ can be estimated recursively for the empirical probability mass function $\hat{f}$ and linear complexity online method, called the information projection test (IPT), that uses an initial convex decision region and modifies it to suppress most likely false alarms. IPT utilizes the separation between pre- and post-change distributions at the first level and compares the relative entropy of the empirical distribution with respect to the distribution of the most likely outlier at the second level. For the quickest change detection problem, we find bounds on the average run length (ARL) and the worst average detection delay (WADD) and prove its asymptotic optimality up to a multiplicative constant and share simulation results which empirically show that IPT performs similar to GLRT.

To motivate our approach, we first define IPT for fixed window size tests and bound the receiver operating characteristic (ROC) of IPT for the composite hypothesis testing problem. Then, we transition to the transient change detection problem to introduce the concept of delay and bound IPT’s performance of detecting changes with bounded delay before studying Lorden’s criteria. Again, we provide empirical comparisons with of FMA [23] and GLRT.

Another important metric considered is the computational complexity of the algorithms. There is already research on low complexity change point estimation techniques, two examples being changes in video content popularity [24] and networks with propagating changes [25]. Even with a similar detection performance, IPT has a few orders of lower complexity ($99.86\%$ reduction in test time for an alphabet with $\sim 6500$ letters and same order of samples) compared to GLRT. Thus, we claim our linear complexity algorithm achieves a better trade-off between computation time and performance compared to other methods in the literature.

Contributions of this paper can be summarized as follows:

• •

We introduce a simple and novel algorithm for detecting changes with known pre-change distribution but unknown post-change distribution.

• •

We prove performance bounds in the non-asymptotic regime and prove asymptotic optimality.

• •

We empirically benchmark our method under (1) composite hypothesis testing problem, (2) worst case post-change distribution transient change detection problem and (3) worst case post-change distribution quickest change detection problem with Lorden’s criteria.

• •

We study special cases to provide insight and application areas for our method.

• •

The time complexity of our approach is compared against existing methods theoretically and empirically.

• •

We draw new insights by applying our idea in analyzing economic and climate time series data to understand whether periods of long deviations from the norm are outliers or change. We identify shifts from the average behavior of market indices and company returns to show that our scheme successfully manages to detect the periods with regime changes like crises. One of the highlights of our analysis is that the climate change phenomenon that is observed in the last 30 years is highly unlikely to be an outlier, giving credence to the hypothesis that it is caused by exogenous factors.

Earlier version of this work has appeared in [26]. In this version we include performance bounds for the transient and quickest change detection problems, include empirical results for the latter, increase the number of applications on real datasets and provide insightful special cases.

II Model and Problem Statement

Let $\mathcal{A}=\left\{a_{1},\dots,a_{m}\right\}$ be a finite alphabet and $\mathcal{P}$ denote the probability simplex of probability mass functions (p.m.f.s) $f$ over $\mathcal{A}$. Consider $X_{1}^{\infty}=\left(X_{1},\dots,X_{k},\dots\right)$, a sequence of observations, where, given change point $t_{1}$, each random variable $X_{1}^{t_{1}-1}$ is independent and identically distributed (i.i.d.) with known pre-change p.m.f. $f_{0}$ and $X_{t_{1}}^{\infty}$ are i.i.d. with the unknown post-change p.m.f. $f_{1}\in\mathcal{P}_{1}$, independent of previous observations $X_{1}^{t_{1}-1}$. Assume there exists $q:\left(\mathcal{P},l_{1}\right)\to\mathbb{R}$ quasiconcave and Lipschitz continuous with Lipschitz constant $L$, that satisfies $q(f_{0})=q_{0}<0<\underline{q}\leq q_{1}=q(f_{1})$. In other words, $\mathcal{P}_{1}$ is a subset of a closed and convex set $q^{-1}\left(\left[\underline{q},\infty\right)\right)$ and $f_{0}\notin q^{-1}\left(\left[\underline{q},\infty\right)\right)$. We use $\hat{f}_{X_{i}^{j}}$ for the empirical p.m.f. of $X_{i}^{j}$, ie. $\hat{f}_{X_{i}^{j}}(x)=\frac{1}{j-i+1}\sum_{k=i}^{j}\mathbb{1}_{x_{i}}(x)$, and $\mathcal{P}^{n}$ for the discrete set of empirical p.m.f.s that can be realized with $n$ samples from $\mathcal{A}$. Finally, $I\left(f\middle\|f^{\prime}\right)=\sum_{a}\log\frac{f(a)}{f^{\prime}(a)}$ denotes the Kullback-Leibler (KL) divergence between distributions $f$ and $f^{\prime}$ in nats and $P_{f_{1},t_{1}}$ denotes the probability law when post-change distribution is $f_{1}$ and change time is $t_{1}$. When change does not occur we denote $P_{f_{1},\infty}$ briefly as $P_{\infty}$.

First, we consider the composite hypothesis testing problem. Given $n$, the null hypothesis, $\mathbf{H}_{0}$, is true if $X_{1}^{n}$ are i.i.d. $f_{0}$ and the alternative hypothesis, $\mathbf{H}_{1}$, is true if i.i.d. $f_{1}$. Our goal is to minimize the maximum probability of misdetection given an upper bound on the probability of false alarm. An example decision rule region for the problem is given in Fig. 1a for $m=3$. Once we define the decision regions $\Gamma_{0}=\left\{\hat{\mathbf{H}}=\mathbf{H}_{0}\right\}$ and $\Gamma_{1}=\mathcal{P}^{n}\setminus\Gamma_{0}=\left\{\hat{\mathbf{H}}=\mathbf{H}_{1}\right\}$ the problem can be stated as follows.

[TABLE]

Second, we consider the transient change detection problem [21] with unknown post-change distribution. The transient change detection problem minimizes the probability of misdetecting a single transient change within a specified window after the change. Without loss of generality (w.l.o.g.), we aasume the change in distribution is permanent but has to be detected in a bounded time $n$. The maximization is subject to an upper bound on the probability of false alarm within any window of size $n_{\alpha}$. If $t_{a}$ denotes an arbitrary alarm (stopping) time of $X_{1}^{\infty}$, the problem can be formulated as

[TABLE]

where we assume the worst post-change distribution. The optimal solution to this problem is unknown even for known post-change distribution except for the special case $n=1$ subject to a minimum ARL [21].

Finally, we consider the quickest change detection problem with unknown post-change distribution. The quickest change detection in Lorden’s sense [27] involves minimizing WADD given a minimum ARL. Assuming the worst post-change distribution, the problem can be formulated as

[TABLE]

Optimal solution to this problem is also unknown. Therefore, we compare our solution of (3), (4) and (5), (6) with frequently used (asymptotically optimal) methods in the literature.

III Motivation

The first problem is an instance of composite hypothesis testing [28], without a parametric model under $\mathbf{H}_{1}$. A typical detector in this case will pick $\mathbf{H}_{1}$ if the empirical distribution, $\hat{f}=\hat{f}_{X_{1}^{n}}$ is closer to a candidate distribution $f_{1}$ than it is to $f_{0}$. Even in that case, however, it is not necessarily true that $X_{1}^{n}$ is drawn from some $f_{1}$. It may still be drawn from $f_{0}$, yet the empirical distribution might look as if it is drawn from some $f_{1}$, leading to a false positive. We call such strings that are drawn from $f_{0}$, but satisfy $\hat{f}\in\mathcal{P}_{1}$ an outlier.

For clarity we provide insight from large deviations theory where the main results are asymptotic, but our contributions apply to the non-asymptotic regime as well. Next, we study the way outliers occur, when they occur using Sanov’s theorem. See [29] for a detailed treatment.

Theorem** (Sanov).**

For any continuous and quasiconcave function $q:\mathcal{P}\to\mathbb{R}$, p.m.f. $f_{0}$ and constant $c\neq\sup_{f}q(f)$,

[TABLE]

The proof of the most general form can be found in [30]. Note that the interesting case is when $q(f_{0})<c<\sup_{f}q(f)$ and that we have the equality in Sanov’s theorem, since $\left\{q(f)\geq c\right\}$ is the closure of its interior as shown in Fig. 1b, where we also illustrate Sanov’s theorem. The solution of the constrained convex optimization $\operatorname*{argmin}_{q(f)\geq c}I\left(f\middle\|f_{0}\right)$ is called the information projection or I-projection of $f_{0}$ onto the superlevel set $\left\{q(f)\geq c\right\}$ and is denoted by $f^{*}$. The theorem gives an asymptotic result in the number of samples and provides the exact characterization of the exponent at which the probability of an outlier decays. This implies that the most likely way for a sustained deviation to occur is the empirical distribution of the associated string to look as if it is drawn from the I-projection $f^{*}$. The probability of this particular way of deviation dominates the probability of all others. An approximate version of this result also holds for a finite string of length $n$ if Stirling’s approximation is satisfied. In the next section, we give our algorithm which is based on using the projected distribution as the most likely outlier distribution. Subsequently, we evaluate the detection performance in the non-asymptotic region to justify its applicability.

IV Information Projection Test

The information projection test (IPT) exploits the observation that outliers occur in a particular way with high probability, while there is no such structure for the deviations. We describe the fixed window size IPT for hypothesis testing and transient change detection here and the variable size IPT for the quickest change detection in Subsection V-C. Given a window size $n$, the fixed window IPT involves the following steps.

Pick $c^{S},c^{D}\geq 0$ and window size $n$. 2. 2.

Find the I-projection $f^{*}=\operatorname*{argmin}_{q(f)\geq c^{S}}I\left(f\middle\|f_{0}\right)$. 3. 3.

At time $k\geq n$, find $S_{k}=q\left(\hat{f}_{X_{k-n+1}^{k}}\right)$. 4. 4.

If $S_{k}\geq c^{S}$, find $D_{k}=I\left(\hat{f}_{X_{k-n+1}^{k}}\middle\|f^{*}\right)$; else, return to Step 3. 5. 5.

If $D_{k}\geq c^{D}$, stop and claim a change has happened; else, return top Step 3.

Note that Step 5 determines if the empirical distribution $\hat{f}$ is close to the most likely fistribution for an outlier, $f^{*}$. This process is illustrated in Fig. 1c. In the rest of the paper, we call the relative entropy $I\left(f\middle\|f^{*}\right)$ also as the relative log likelihood function (RLLF) with respect to $f^{*}$, referring to the comparison of the likelihood of the projection distribution with that of the empirical distribution of the observation. We also use $c$ to denote $\left(c^{S},c^{D}\right)$.

V Performance Guarantees

V-A Composite Hypothesis Testing

The next two theorems bound the probability of false alarm and misdetection as a function of $c$. We show that the false positive probability decays exponentially with $c^{D}+\frac{\left(\left(c^{S}-q_{0}\right)^{+}\right)^{2}}{2L^{2}}$ and that the worst case probability of misdetection decays with rate $\left(\frac{\underline{q}-q_{0}-\left(c^{S}-q_{0}\right)^{+}}{\sqrt{2}L}-\sqrt{c^{D}}\right)^{2}$.

Theorem 1**.**

For any $c$ such that $c^{S}<\sup_{f}q(f)$,

[TABLE]

Theorem 2**.**

For any $c$ such that $\left(c^{S}-q_{0}\right)^{+}+\sqrt{2L^{2}c^{D}}<\underline{q}-q_{0}$,

[TABLE]

V-B Transient Change Detection

In this subsection, we prove corollaries to Theorem 1 and 2 that provide a bound on the transient change detection performance of IPT under the assumption of unknown post-change distribution. The transient change detection is regards the pair of metrics: the worst case probability of false alarm and worst case probability of misdetection within a given length of interval. Since we are given a window of $n$ samples to detect the change and $n\ll n_{\alpha}$, we utilize a fixed window IPT. This problem was first proposed in [31] and is also studied recently in [21, 23].

Corollary 3**.**

For any $c$ such that $\left(c^{S}-q_{0}\right)^{+}+\sqrt{2L^{2}c^{D}}<\underline{q}-q_{0}$ the IPT stopping time $t_{I}$ with fixed window length $\frac{n+1}{2}$ satisfies

[TABLE]

where

[TABLE]

V-C Quickest Change Detection

In this subsection, we apply IPT to the quickest change detection problem with an unknown post-change distribution. To minimize the average detection delay under the worst change time $t_{1}$ and pre-change realizations $X_{1}^{t_{1}-1}$, we modify IPT to have an effective window that only contains the samples with a positive drift from the minimum the random walk has visited, similar to the CUSUM algorithm. We also use a restart mechanism if the empirical distribution is similar to a most likely false alarm, similar to the IPT with fixed window size. Since, at any given time $k$, the effective window size $n_{k}$ is varying, we determine define the most likely deviation as a function of the window size, $f_{n}^{*}$.

[TABLE]

For any $k$, $\tau_{k}$ denotes the last time the IPT has been restarted since the empirical distribution resembles a false alarm. So, $\tau_{1}=1$ and for any $k$,

[TABLE]

We the index $i_{k}$ that maximizes $S_{k}$ to define the effective window size $n_{k}=k-i_{k}+1$. Finally, we compare the proximity of the empirical distribution over the effective window with respect to the most likely false alarm given the window size $n_{k}$. If $D_{k}$ is large, we stop the observation process and raise an alarm, else we conclude a false alarm and reset the algorithm. The threshold for this step is allowed to vary with $n$, so we select a sequence of $c_{n}^{D}\geq 0$ instead of a constant $c^{D}$.

[TABLE]

We also use shorthand notations $c_{n}$ and $Q(i,j)$ to refer to $\left(\frac{c^{S}}{n},c_{n}^{D}\right)$ and $(j-i+1)q\left(\hat{f}_{X_{i}^{j}}\right)$ respectively. We define the auxiliary stopping time $t_{S}=\inf\left\{k\middle|S_{k}\geq c^{S}\right\}$ where only the first statistic crosses the threshold and we define $t_{S}^{k}$ and $t_{I}^{k}$ to be the stopping times starting at $k$, counting only the observations after $X_{k-1}$.

In the following theorem we show that the ARL of our method can be lower bounded in terms of $c$. We prove the theorem where $q$ is the expectation operator but the general result can be obtained by constructing a separating hyperplane between the convex and closed set $\mathcal{P}_{1}$ and the initial distribution $f_{0}\notin\mathcal{P}_{1}$. Then, step 1 of the IPT holds only if the hyperplane is crossed. Thus, a lower bound on $E_{\infty}t_{S}$ can be formulated using this halfspace instead of $\mathcal{P}_{1}$.

Theorem 4**.**

If $c_{n}^{D}=c^{D}\geq\frac{2\left|q_{0}\right|\underline{q}}{(1+\rho)L^{2}}$ for $n>(1+\rho)\frac{c^{S}}{\underline{q}}$ for some $\rho>0$, then, the ARL of IPT is bounded as

[TABLE]

where $\psi(v)=E_{\infty}\exp(vX)$ and $v^{*}>0$ satisfies $\psi(v^{*})=1$.

The following lemma proves a bound on the stopping time for any change time $t_{1}$ and observations $X_{1}^{t_{1}-1}$.

Lemma 5**.**

For any $t_{1},f_{1}$ and $X_{1}^{t_{1}-1}$, the IPT stopping time $t_{I}$ is bounded by

[TABLE]

for any $j=0,\dots$ where $t_{S}^{(j)}$ is the $j$th time $S_{k}$ crosses $c^{S}$, i.e. $t_{S}^{(0)}=0$ and $t_{S}^{(j+1)}=t_{S}^{(j)}+t_{S}^{t_{S}^{(j)}+1}$. Further, for any $j=1,\dots$ and $t_{S}^{(j-1)}<k\leq t_{S}^{(j)}$,

[TABLE]

Thus, for $j^{*}=\min\left\{j\middle|t_{S}^{(j)}\geq t_{1}\right\}$,

[TABLE]

In the next lemma, we prove that the average delay is bounded by a function of the average delay when $t_{1}=1$ using the bound we have found in Lemma 5. This lemma shows that the worst average delay IPT experiences is bounded by the average delay up to a multiplicative constant.

Lemma 6**.**

For any $f_{1},t_{1}$ and $X_{1}^{t_{1}-1}$, the average delay is bounded as follows.

[TABLE]

Before showing a bound on the WADD, we bound the average detection delay when change happens at $t_{1}=1$ asymptotically as $c^{S}\to\infty$.

Lemma 7**.**

For any $f_{1}$, if $c_{n}^{D}=0$ for $n\leq(1+\rho)\frac{c^{S}}{\underline{q}}$ for some $\rho>0$, then, as $c^{S}\to\infty$

[TABLE]

The next theorem is the asymptotic characterization of the WADD of IPT. In the proof we first derive a non-asymptotic bound and then consider the limit as $c^{S}\to\infty$.

Theorem 8**.**

If $c_{n}^{D}=0$ for $n\leq(1+\rho)\frac{c^{S}}{\underline{q}}$ for some $\rho>0$, then, as $c^{S}\to\infty$,

[TABLE]

Lorden showed in his seminal work [27] that the optimal WADD for any change detection algorithm with ARL at least $\gamma$ is $\frac{\log\gamma}{I\left(f_{1}\middle\|f_{0}\right)}$ asymptotically as $\gamma\to\infty$ even if the post-change distribution $f_{1}$ is known. Using Theorem 4 and 8, we prove that IPT is asymptotically optimal in Lorden’s criteria up to a multiplicative constant.

Corollary 9**.**

For $\rho>0,c^{D}\geq\frac{2\left|q_{0}\right|\underline{q}}{(1+\rho)L^{2}}$ and $c_{n}^{D}=\begin{cases}0,&n\leq(1+\rho)\frac{c^{S}}{\underline{q}}\\ c^{D},&\text{otherwise}\end{cases}$, as $\gamma\to\infty$,

[TABLE]

Proof:

Use Theorem 4 and 8. ∎

We end our discussion on the quickest change detection problem by describing two special cases.

V-C1 Special Case: Change in Mean

For the special case where the change is in the mean, IPT reduces to a CUSUM algorithm with restart.

[TABLE]

Further, the most likely deviations $f_{n}^{*}$ become the exponentially tilted distributions with mean $\max\left\{\frac{c^{S}}{n},q_{0}\right\}$.

[TABLE]

where

[TABLE]

V-C2 Special Case: Existence of a Representative $f_{1}\in\mathcal{P}_{1}$ for Log-Likelihood Ratio

For the special case where there exists $f_{1}\in\mathcal{P}_{1}$ and $\rho>0$ such that for all $f\in\mathcal{P}_{1}$,

[TABLE]

the function $q\left(f\right)=I\left(f\middle\|f_{0}\right)-I\left(f\middle\|f_{1}\right)$ is a quasiconcave and Lipschitz continuous function with $q_{0}=q\left(f_{0}\right)=-I\left(f_{0}\middle\|f_{1}\right)<0$ and for any $f\in\mathcal{P}_{1}$,

[TABLE]

Thus, for all $f\in\mathcal{P}_{1}$, $0<\rho I\left(f_{1}\middle\|f_{0}\right)=\underline{q}\leq q(f)$. Further, we can interpret $Q(1,k)=kq\left(\hat{f}_{X_{1}^{k}}\right)$ as a random walk with i.i.d. steps of the log-likelihood ratio $\log\frac{f_{1}\left(X_{i}\right)}{f_{0}\left(X_{i}\right)}$. Finally, $S_{k}=\max_{\tau_{k}\leq i\leq k}\sum_{j=i}^{k}\log\frac{f_{1}\left(X_{i}\right)}{f_{0}\left(X_{i}\right)}$ becomes the CUSUM algorithm with restart.

VI Comparison to Other Detection Schemes

We compare IPT with various change detection algorithms proposed in the literature in terms of ROC and Lorden’s criteria.

VI-A Composite Hypothesis Testing

For the hypothesis testing and transient change detection problems, we describe the simulation procedure and plot the empirical ROC, worst case probability of misdetection versus probability of false alarm, for each of the methods considered.

The ROC associated with our method is illustrated in Fig. 2a and compared with the ROC of an FMA test and GLRT. We use a ternary alphabet $\{-1,0,1\}$, define the pre-change distribution to be uniform and the change in mean to satisfy $\underline{q}=0.25$. The deviations are chosen to have $n=25$. At this stage, we determine the I-projection of $f_{0}$ on $\mathcal{P}_{1}$. To find the optimum performance of our algorithm, we change $c^{D}$ over $[2^{-8},2^{-3}]$. The curves in Fig. 2a show that, our proposed method is able to decrease the false positive rate compared to CUSUM-like FMA filter for the same misdetection rate and performs similar to the GLRT. IPT’s average misdetection performance is within $27\%$ of that of GLRT but performs $34\%$ better than FMA over the measured region of the ROC. The area under the misdetection-false alarm curves are $0.343,0.224$ and $0.184$ for FMA, IPT and GLRT respectively.

VI-B Transient Change Detection

We consider the transient change detection problem under a family of discrete Gaussian distributions. Discrete Gaussian sampling is of particular interest in lattice-based cryptography [32] and thereby in quantum-resilient cryptography [33]. For hardware based solutions proposed to generate fast and true random numbers [34], detecting a change in the distribution of generated discrete Gaussian samples could increase the reliability against failures. Here, we only consider a finite alphabet discrete Gaussian with $\mathcal{A}=\{-5,\dots,5\},m=11,n=80,d_{\alpha}=200$ and $f(k)=N\exp\left(-\frac{k^{2}}{2d^{2}}\right)$. The known pre-change distribution has $\sigma_{0}^{2}=1$ and the post-change distribution satisfies $\sigma_{1}^{2}\geq\underline{\sigma}^{2}=2$. For $q(f)=\sigma_{f}^{2}$, the set of post-change distributions $\mathcal{P}_{1}=q^{-1}([\underline{\sigma}^{2},\infty))$ is a convex subset of $\mathcal{P}$, thus we can utilize the IPT. For each method, we use a rolling window of size $w=20<\frac{n+1}{2}$. The worst change point $t_{1}$ is chosen among $t_{1}=1,\dots,w$ and we randomly specify $100$ different post-change distributions from $\mathcal{P}_{1}$. Each data point in Fig. 2b is the empirical average of $10^{5}$ worst case false alarms and worst case misdetection tests. With $\frac{20}{11}\sim 2$ samples per letter, GLRT overfits the empirical data while estimating $\sigma_{1}^{2}$ and performs worse than IPT and FMA.

VI-C Quickest Change Detection

For the quickest change detection problem, we test the IPT as described in Subsection V-C in terms of WADD versus ARL under the test procedure described. Then, the empirical delay curves of different algorithms are compared.

In Fig. 2c, we consider the quickest change detection performance of IPT, FMA and GLRT. The minimum WADD is plotted against the ARL under the same setting as in Subsection VI-A. IPT provides a tradeoff between the mean test and GLRT with slightly worse detection delay compared to the latter. But in a scenario where computational costs would not allow a user to choose only according to the delay characteristics IPT may outperform both FMA and GLRT. A fast stream of data tracked for change at a central unit is such an example.

VII Complexity

In this section we describe the complexity of the change detection algorithms: IPT, FMA and GLRT. We compare the computational complexity as a function of sample size $n$, alphabet size $m$ or error $\epsilon\sim\norm{\tilde{f^{*}}-f^{*}}_{1}$ involved in I-projection for IPT and maximum likelihood estimation during GLRT. The FMA test using a rolling window has $\Theta(n)$ complexity222$\mathcal{O}(g)$ is the set of functions asymptotically upper bounded by $g$ [35], i.e. $\mathcal{O}(g)=\left\{f\middle|\exists d,u_{0}>0\ \text{such that}\ \forall u\geq u_{0},0\leq f(u)\leq dg(u)\right\}$ where, if $u$ is a vector variable the inequalities are componentwise. Similarly, $\Theta(g)=\left\{f\middle|f\in\mathcal{O}(g)\ \text{and}\ g\in\mathcal{O}(f)\right\}$. to decide on the alarm time during execution. The GLRT forms a likelihood ratio for each new sample and therefore solves the following (nonconvex) optimization problem:

[TABLE]

which (if $\mathcal{P}_{1}$ is convex) is computationally equivalent to binary searching two Lagrange multipliers using an equation of $m$ rational terms, thus it has $\mathcal{O}\left(n+m\log\frac{1}{\epsilon}\log\log\frac{1}{\epsilon}\right)$ complexity. IPT omits this projection step by initially finding the most likely deviations and then comparing any realization to this distribution. Therefore, it has only $\Theta(n+m)$ complexity. Since GLRT does not allow a recursive update, its complexity is the same for a sliding window method. Assuming the function $q\left(\hat{f}\right)$ admits a recursive update in $\Theta(1)$, IPT has overall $\Theta(1)$ complexity for sliding windows. This is because $I\left(\hat{f}\middle\|f^{*}\right)$ has a simple update rule. Similarly FMA with rolling window has $\Theta(1)$ complexity. Thus, whereas IPT and FMA have linear complexity for rolling windows and constant complexity for sliding windows, GLRT has superlinear complexity in both cases.

For a data stream with packet rate $R$ satisfying $\frac{1}{R}\gg n+m$, GLRT can not respond to changes in the time series as fast as IPT or FMA. This is especially true if a worst case detection delay is used as a cost. With a rapidly increasing demand for social media content, video, distributed computation there will be a growing need for fast heuristics than for offline exact computation.

VII-1 Computational Complexity

Using Ohio Supercomputer Center Owens cluster’s single node of 128GB memory and 2.40GHz clock the test times varied as in Fig. 2d [36]. We increased the alphabet and sample sizes proportionally and averaged over the total number of tests. IPT performed $25$ to $670$ times faster than the GLRT where the gain from test time increased with the alphabet size.

VIII Applications to Economics and Climate Data

In this section, we apply our outlier detection scheme, based on checking the RLLF of the most likely outlier distribution to the empirical distribution of the given period. We would like to test whether shifts of historic value are merely outliers or it happens due to an exogenous factor. To obtain the ground truth, we will use the long-term empirical distribution. Thus, the underlying assumption we make in the following data analysis is that, for $r$ shifts of length $n$ the overall cumulative duration of deviations is much smaller compared to the size of the data: $rn\ll T$. To analyze a period of interest, we pick a duration $n$ and a threshold, $c^{S}$, that matches the mean of the underlying data over that period.

VIII-A Daily Returns of Portfolios by Size

We focus on portfolio returns with different market caps [37] to identify segments with a worse performance compared to the average behavior, sustained over a substantial duration of time. The primary question we answer is the following: “is the identified negative sequence a rare event generated by the statistical nature of the market or is it driven by some exogenous factors with the potential to cause a financial crisis?” The monthly averages of the return vary between $0.8-0.9$. We quantized the daily percentage portfolio returns for different market caps over July 1926-March 2019 and selected a threshold $c^{S}$ below the average returns. We then computed the I-projection of the long term empirical distribution over the distributions with mean less than the threshold. The $n=6$ month moving average and the KL divergence of the observation window against the I-projection is given in Fig. 3. With the right choice of $n$ and $c$ we were able to identify historically significant events like the great depression or the 2009 financial crisis. One of the most significant finding is that, even though there are other periods with mean return as low as these periods, the RLLF clearly differentiated these periods from the crises times, showing strength of our scheme.

VIII-B Oil Industry Market Index

Similar to Subsection VIII-A, this time we focus on market index for oil industry. We use monthly return series of Oil Industry Portfolio Index, obtained from Kenneth French Data Library [37] and covers the period from 1926-07-01 to 2019-03-31. The historical monthly return of the index is $0.89$. We have run the experiment with $c^{S}=-0.1\%$ and $n=12$ months. In Fig. 4, we have $7$ first threshold crosses, only $2$ of which classified as change. The greatest divergence is for the $2009$ oil crisis period. Unfortunately, the $1973$ oil embargo is disregarded as a false alarm. This example shows the requirement of tuning for our method.

VIII-C Oxygen Isotope Data

We analyze climate data collected from a polar cap that spans $\sim 1800$ years to gain insight of our method and historic climate data [38]. Fig. 5 uses a moving average window of $n=20$ years of quantized data of oxygen isotope density over the years $226-2009$. We then select thresholds $c^{S}$, half a standard deviation above the mean, and $c^{D}$ and determine the I-projections. Note that the only point IPT classifies as change is the last 20-year period of above norm d18O levels. Although the threshold selections are not unique this intuitive approach with these hyperparameters gives insight that the last 20-year period is less likely to be an outlier rather than an effect of an exogenous change, ie. man made.

IX Conclusion

In a variety of applications, it is important to identify shifts from the typical behavior. However, not every shift from the norm marks a regime change or a deviation; instead it could be an outlier. In this paper, we developed a method that differentiates between outliers and change. To achieve that, we used results from large deviations theory, which gives us statistical characterization for the outliers. Our method tests a given string against the most likely way an outlier occurs and marks deviations using detection-theoretic tools.

Our method uses control parameters $c^{S}$ and $c^{D}$ to shape the decision region to determine the outliers. We have bounded the probability of false positive and misdetection as functions of these control parameters and shown that the detection delay is asymptotically optimal given an ARL.

We also emphasized IPT’s computational simplicity in that it reduces GLRT’s requirement to run non-convex and non-recursive optimization problems to solving convex optimization problems during initialization. This has proven to reduce the average running time per test 25 to 670 fold for varying alphabet sizes.

We also applied our algorithm to a variety of applications and drawed insights from the observed time series data. Our method verified the global warming phenomenon as highly unlikely to be an outlier, giving credence for it to be caused by factors (such as those man-made) leading to a statistical change in the indicator variables. Similarly, we tried to extract historically significant exogeneous events for the market within the data.

We have omitted proves that complicate the exposition whenever necessary but we believe similar results can be extended to the continuous alphabet case via binning.

Future work will be on proposing a metric that naturally and objectively combines sample delay for detection with computational complexity to address true delays as observed in high traffic data centers that detect changes in popularity, inappropriate content, fake accounts or attacks.

Appendix A Inequalities

The inequalities utilized throughout the paper are briefly provided in this subsection for easy reference.

Theorem** (Pythagorean Theorem for Relative Entropy [39]).**

For any closed and convex set $C\subset\mathcal{P}$, let $f^{*}$ denote the I-projection of $f_{0}$ on $C$. Then, for any $f\in C$,

[TABLE]

Theorem** (Pinsker’s Inequality [39]).**

For any pair $f$ and $f^{\prime}$,

[TABLE]

Corollary**.**

For any $f$ and $f^{\prime}$,

[TABLE]

Proof:

The first inequality is the Pinsker’s inequality. The second inequality is due to the Lipschitz property of $q$. ∎

Theorem** (Inequality for the $L_{1}$ Deviation [40]).**

For any $f_{0}$, $n$ and $\delta>0$,

[TABLE]

Appendix B Proofs of Theorems

B-A Proof of Theorem 1

Given $c$ such that $c^{S}<\sup_{f}q(f)$, we have $\Gamma_{1}=\left\{f\in\mathcal{P}^{n}\middle|q\left(f\right)\geq c^{S},I\left(f\middle\|f^{*}\right)\geq c^{D}\right\}$ and

[TABLE]

B-B Proof of Theorem 2

Given $c$ such that $\left(c^{S}-q_{0}\right)^{+}+\sqrt{2L^{2}c^{D}}<\underline{q}-q_{0}$, $\Gamma_{0}=\left\{f\in\mathcal{P}^{n}\middle|q\left(f\right)<c^{S}\ \text{or}\ I\left(f\middle\|f^{*}\right)<c^{D}\right\}$ and for all $f\in\Gamma_{0}$, either $q(f)<c^{S}$ or $q(f)\geq c^{S}$, $I\left(f\middle\|f^{*}\right)<c^{D}$ and

[TABLE]

Thus, for any $f_{1}$,

[TABLE]

and

[TABLE]

B-C Proof of Corollary 3

Consider the probability of false alarm. The alarm time can only occur at $k^{\prime}(l)=\frac{n+1}{2}l$ for positive integers $l=1,\dots$ and each window has $\frac{n+1}{2}$ samples i.i.d. with $f_{0}$. Then,

[TABLE]

where we have used Theorem 1. Thus,

[TABLE]

Consider the probability of misdetection. Since the window size is $\frac{n+1}{2}$ the IPT with fixed size and rolling window will have at least one window where samples are i.i.d. with $f_{1}$. Then, the upper bound is a direct consequence of Theorem 2.

B-D Proof of Theorem 4

Let us express $E_{\infty}t_{I}$ conditional on $t_{S}$.

[TABLE]

Assume $c_{n}^{D}=c^{D}$ for $n>(1+\rho)\frac{c^{S}}{\underline{q}}$, then, we can upper bound the denominator as follows.

[TABLE]

since $n+1\leq\left((1+\rho)\frac{c^{S}}{\underline{q}}+1\right)^{\frac{\underline{q}}{(1+\rho)c^{S}}n}$ for $n>(1+\rho)\frac{c^{S}}{\underline{q}}$. Also, w.l.o.g., we assumed $\left|q_{0}\right|\geq(1+\rho)\underline{q}$. Then,

[TABLE]

for $c^{D}\geq\frac{2\left|q_{0}\right|\underline{q}}{(1+\rho)L^{2}}$. Next, we bound lower bound $E_{\infty}t_{S}$. We assume $q$ is the mean operator and $S_{k}$ is the random walk that is bounded below by zero. Denote the number of zero crossings as $N_{0}$ and the latest upper or lower threshold crossing time as $k^{*}$, then,

[TABLE]

The random walk has i.i.d. and bounded steps. Therefore its steps has moments of all order and we can apply Wald’s identity. For every $v\in\mathbb{R}$,

[TABLE]

and

[TABLE]

Thus,

[TABLE]

Finally we obtain,

[TABLE]

B-E Proof of Lemma 5

Let $j=0,\dots$ be arbitrary. W.l.o.g. assume $t_{I}\geq t_{S}^{(j)}$. If $D_{t_{S}^{(j)}}\geq c^{D}$, we have $t_{I}=t_{S}^{(j)}$. Else, the stopping time is restarted at $t_{S}^{(j)}+1$ and $t_{I}=t_{S}^{(j)}+t_{I}^{t_{S}^{(j)}+1}$. In any case, (20) is satisfied.

Let $j=1,\dots,$ and $t_{S}^{(j-1)}<k\leq t_{S}^{(j)}$ be arbitrary. If $\tau_{k+t_{S}^{k}-1}>k$, then there exists $k^{\prime}$ such that $k<k^{\prime}\leq k+t_{S}^{k}-1$ and $\tau_{k+t_{S}^{k}-1}=\tau_{k^{\prime}}=k^{\prime}$. Thus, $t_{S}^{(j)}\leq k^{\prime}\leq k+t_{S}^{k}-1$. Else, $\tau_{k+t_{S}^{k}-1}\leq k$

[TABLE]

Thus, $t_{S}^{(j)}\leq k+t_{S}^{k}-1$.

Finally, let $j^{*}=\min\left\{j\middle|t_{S}^{(j)}\geq t_{1}\right\}$. Then, using (20) and (21) with $j=j^{*}$ and $k=t_{1}$ we get

[TABLE]

B-F Proof of Lemma 6

For any $f_{1},t_{1}$ and $X_{1}^{t_{1}-1}$,

[TABLE]

Therefore,

[TABLE]

since $t_{S}^{(j^{*})}\geq t_{1}$.

B-G Proof of Lemma 7

First, we express the expected delay conditional to $t_{S}$.

[TABLE]

The nominator is bounded using the probability that $S_{k}$ does not cross the threshold until time $t$.

[TABLE]

Then, for any $\delta>0$,

[TABLE]

The denominator is bounded using the fact that $t_{S}\leq(1+\rho)\frac{c^{S}}{q_{1}}$ with high probability.

[TABLE]

since $c^{D}_{n}=0$ and $t_{I}=t_{S}$ for $n_{t_{S}}\leq(1+\rho)\frac{c^{S}}{\underline{q}}$. So,

[TABLE]

As $c^{S}\to\infty$,

[TABLE]

since $\delta>0$ is arbitrary.

B-H Proof of Theorem 8

Using Lemma 6 and 7, as $c^{S}\to\infty$,

[TABLE]