Robust Change Detection for Large-Scale Data Streams

TL;DR

This paper introduces a robust, scalable change detection method for large-scale data streams that effectively handles outliers and unknown affected streams, with proven robustness and efficiency.

Contribution

It develops a new local robust detection statistic and a global scheme for real-time monitoring, along with novel robustness and efficiency metrics.

Findings

The proposed scheme effectively filters unaffected streams.

It demonstrates robustness against outliers.

Numerical simulations confirm high detection efficiency.

Abstract

Robust change-point detection for large-scale data streams has many real-world applications in industrial quality control, signal detection, biosurveillance. Unfortunately, it is highly non-trivial to develop efficient schemes due to three challenges: (1) the unknown sparse subset of affected data streams, (2) the unexpected outliers, and (3) computational scalability for real-time monitoring and detection. In this article, we develop a family of efficient real-time robust detection schemes for monitoring large-scale independent data streams. For each data stream, we propose to construct a new local robust detection statistic called $L_{\alpha}$-CUSUM statistic that can reduce the effect of outliers by using the Box-Cox transformation of the likelihood function. Then the global scheme will raise an alarm based upon the sum of the shrinkage transformation of these local $L_{\alpha}$-CUSUM statistics so as to filter out unaffected data streams. In addition, we propose a new concept called {\em false alarm breakdown point} to measure the robustness of online monitoring schemes and propose a worst-case detection efficiency score to measure the detection efficiency when the data contain outliers. We then characterize the breakdown point and the efficiency score of our proposed schemes. Asymptotic analysis and numerical simulations are conducted to illustrate the robustness and efficiency of our proposed schemes.

Tables

Table 1. Table 1: A comparison of the detection delays of 8 8 8 schemes with γ = 5000 𝛾 5000 \gamma=5000 under the idealized model. The smallest and largest standard errors of these 8 8 8 schemes are also reported under each post-change hypothesis based on 1000 1000 1000 repetitions in Monte Carlo simulations.

Gross error model with $ϵ=0$
	# affected local data streams
	$1$	$3$	$8$	$10$	$15$	$20$	$50$	$100$
Smallest standard error	$0.29$	$0.12$	$0.05$	$0.04$	$0.03$	$0.03$	$0.01$	$0.00$
Largest standard error	$0.58$	$0.20$	$0.07$	$0.06$	$0.05$	$0.03$	$0.02$	$0.01$
Our proposed robust scheme
$N_{\alpha =0.51}^{\left(soft\right)}\left(b=8.5,d=0.8915\right)$	$41.0$	$18.6$	$10.3$	$9.2$	$7.5$	$6.5$	$4.5$	$3.9$
$N_{\alpha =0.51}^{\left(r=10\right)}\left(b=17.19\right)$	$40.6$	$18.5$	$10.3$	$9.2$	$7.7$	$6.9$	$5.3$	$4.8$
Comparison of other methods
$N_{\alpha =0}^{\left(soft\right)}\left(b=21.52,d=2.3026\right)$	$33.6$	$15.2$	$8.4$	$7.5$	$6.1$	$5.3$	$3.7$	$3.0$
$N_{\alpha =0.51,\max }\left(b=4.3\right)$	$27.7$	$19.6$	$16.2$	$15.6$	$14.8$	$14.2$	$12.7$	$11.9$
$N_{\alpha =0.51,\text{sum}}\left(b=36.85\right)$	$63.7$	$26.9$	$12.5$	$10.5$	$7.8$	$6.4$	$3.3$	$2.0$
$N_{Chan,\alpha =0.51}\left(b=1.04,p_0=0.1\right)$	$31.4$	$17.7$	$10.8$	$9.7$	$7.8$	$6.7$	$4.1$	$3.0$
$N_{Chan,\alpha =0}\left(b=21.6,p_0=0.1\right)$	$32$	$15.2$	$11.2$	$7.5$	$5.3$	$4.2$	$3.3$	$2.3$
$N_{XS}\left(b=19.5,p_0=0.1\right)$	$30.9$	$13.2$	$7.2$	$5.7$	$4.7$	$3.5$	$1.8$	$1.0$

Table 2. Table 2: A comparison of the detection efficiency score of 8 8 8 schemes under the gross error model with ϵ = 0.1 italic-ϵ 0.1 \epsilon=0.1 based on 1000 1000 1000 repetitions in Monte Carlo simulations. The threshold b 𝑏 b is chosen to satisfy γ = 5000 𝛾 5000 \gamma=5000 in the idealized model.

Gross error model with $ϵ=0.1$
	# affected local data streams
	$1$	$3$	$8$	$10$	$15$	$20$	$50$	$100$
Our proposed robust scheme
$N_{\alpha =0.51}^{\left(soft\right)}\left(b=8.5,d=0.8915\right)$	$0.17$	$0.34$	$0.61$	$0.68$	$0.83$	$0.95$	$1.37$	$1.66$
$N_{\alpha =0.51}^{\left(r=10\right)}\left(b=17.09\right)$	$0.17$	$0.35$	$0.62$	$0.68$	$0.82$	$0.92$	$1.2$	$1.35$
Other methods for comparison
$N_{\alpha =0}^{\left(soft\right)}\left(b=21.52,d=2.3026\right)$	$0.27$	$0.32$	$0.4$	$0.43$	$0.48$	$0.53$	$0.7$	$0.8$
$N_{\alpha =0.51,\max }\left(b=4.3\right)$	$0.3$	$0.36$	$0.44$	$0.46$	$0.49$	$0.51$	$0.58$	$0.62$
$N_{\alpha =0.51,\text{sum}}\left(b=36.85\right)$	$0.12$	$0.25$	$0.5$	$0.58$	$0.77$	$0.94$	$1.73$	$2.88$
$N_{Chan,\alpha =0}\left(b=21.6,p_0=0.1\right)$	$0.26$	$0.32$	$0.36$	$0.43$	$0.54$	$0.63$	$0.75$	$1.03$
$N_{Chan,\alpha =0.51}\left(b=1.04,p_0=0.1\right)$	$0.21$	$0.37$	$0.59$	$0.65$	$0.81$	$0.94$	$1.53$	$2.19$
$N_{XS}\left(b=19.5,p_0=0.1\right)$	$0.22$	$0.39$	$0.44$	$0.49$	$0.52$	$0.55$	$0.78$	$0.93$

Equations

$$1:$$

$$2:$$

$$\dots$$

$$K:$$

$$h_0(x)=(1-ϵ)f_0(x)+ϵg_0(x),h_1(x)=(1-ϵ)f_1(x)+ϵg_1(x),$$

$${\mathbf{E}}_{f_0}^{(\infty )}(T)\ge \gamma ,$$

$${ϵ}^{\ast }(T)=\inf \{ϵ\ge 0:\underset{h_0\in {\hslash }_{0,ϵ}}{\inf }\log ({\mathbf{E}}_{h_0}^{(\infty )}T(b_{\gamma }))=o(\log \gamma )\},$$

$${\hslash }_{0,ϵ}=\{h|h=(1-ϵ)f_0+ϵg,g\in \mathcal{G}\},$$

$$\displaystyle\mathbf{D}_{h_{1}}(T)=\quad\sup_{\nu\geq 1}\ \mbox{ess}\sup{\bf E}_{h_{1}}^{(\nu)}\left((T-\nu+1)^{+}\big{|}\mathcal{F}_{\nu-1}\right).$$

$$\text{AE}(T,ϵ;g_0,g_1)=\underset{\gamma \to \infty }{\lim }\frac{\log ({\mathbf{E}}_{h_0}^{(\infty )}(T(b_{\gamma })))}{{\mathbf{D}}_{h_1}(T(b_{\gamma }))},$$

$$\text{WAE}(T,ϵ)=\underset{g_0\in \mathcal{G},g_1\in \mathcal{G}}{\inf }\text{AE}(T,ϵ;g_0,g_1)=\underset{\gamma \to \infty }{\lim }\frac{\underset{g_0\in \mathcal{G}}{\inf }\left[\log ({\mathbf{E}}_{h_0}^{(\infty )}(T(b_{\gamma })))\right]}{\underset{g_1\in \mathcal{G}}{\sup }\left[{\mathbf{D}}_{h_1}(T(b_{\gamma }))\right]}.$$

$$\displaystyle W_{\alpha,k,n}=\max\Big{(}W_{\alpha,k,n-1}+\frac{[f_{1}(X_{k,n})]^{\alpha}-[f_{0}(X_{k,n})]^{\alpha}}{\alpha},0\Big{)},$$

$$\displaystyle W^{*}_{k,n}=\max_{1\leq\nu<\infty}\log\frac{\prod_{i=1}^{\nu-1}f_{0}(X_{k,i})\prod_{i=\nu}^{n}f_{1}(X_{k,i})}{\prod_{i=1}^{n}f_{0}(X_{k,i})}=\max\Big{(}W^{*}_{k,n-1}+\log\frac{f_{1}(X_{k,n})}{f_{0}(X_{k,n})},\ 0\Big{)}.$$

$$N_{\alpha }(b)=\inf \left\{n\ge 1:\sum _{k=1}^{K}h(W_{\alpha ,k,n})\ge b\right\},$$

$$N_{\alpha }^{(soft)}(b,d)$$

$$N_{\alpha }^{(r)}(b)$$

$$N_{Chan,\alpha }(b,p_0)=\inf \left\{n\ge 1:\sum _{k=1}^K\log \left[1-p_0+0.64\ast p_0\exp (W_{\alpha ,k,n}/2)\right]\ge b\right\}.$$

$$N_{\alpha ,\max }(b)$$

$$N_{\alpha ,\text{sum}}(b)$$

$$I_i(ϵ,\alpha ;g_i)$$

$$\displaystyle{\bf E}_{h_{0}}\exp\Big{\{}\lambda(\epsilon,\alpha;g_{0})\frac{[f_{1}(X)]^{\alpha}-[f_{0}(X)]^{\alpha}}{\alpha}\Big{\}}=1.$$

$${\mathbf{E}}_{h_0}^{(\infty )}[N_{\alpha }^{(soft)}(b,d)]\ge \frac{1}{4}\exp \left({\left[\sqrt{\lambda (ϵ,\alpha ;g_0)b}-\sqrt{K\exp \{-\lambda (ϵ,\alpha ;g_0)d\}}\right]}^2\right).$$

$${\mathbf{E}}_{h_0}^{(\infty )}[N_{\alpha }^{(r)}(b)]\ge \frac{1}{4}\exp \left({\left[\sqrt{\lambda (ϵ,\alpha ;g_0)b}-\sqrt{K}\right]}^2\right).$$

$${\mathbf{D}}_{h_1}(N_{\alpha }^{(soft)}(b,d))\le (1+o(1))\frac{1}{I_1(ϵ,\alpha ;g_1)}\left(\frac{b}{m}+d\right),$$

$${\mathbf{D}}_{h_1}(N_{\alpha }^{(r)}(b))\le (1+o(1))\frac{1}{I_1(ϵ,\alpha ;g_1)}\left(\frac{b}{m}\right),$$

$$\text{AE}(N_{\alpha },ϵ;g_0,g_1)\ge m\lambda (ϵ,\alpha ;g_0)I_1(ϵ,\alpha ;g_1).$$

$$\text{AE}(N_{\alpha },ϵ;g_0^{\ast },g_1^{\ast })=0,$$

$$M^{\ast }(\alpha )=\underset{x}{\mathrm{esssup}}|\frac{[f_1(x){]}^{\alpha }-[f_0(x){]}^{\alpha }}{\alpha }|.$$

$$\text{WAE}(N_{\alpha },ϵ)$$

$$d_{opt}=\frac{1}{\lambda ({\alpha }_{opt})}\left(\log \frac{K}{m}+\log \frac{\log \gamma }{m}\right),$$

$$b_{\gamma }=\frac{1}{\lambda ({\alpha }_{opt})}{\left(\sqrt{\log (4\gamma )}+\sqrt{K\exp \{-\lambda ({\alpha }_{opt})d_{opt}\}}\right)}^2,$$

$$M(\alpha )=\underset{x}{\mathrm{esssup}}\frac{[f_1(x){]}^{\alpha }-[f_0(x){]}^{\alpha }}{\alpha },$$

$$d_{\alpha }(f_0,f_1)$$

Full text

Robust Change Detection for Large-Scale Data Streams

**Ruizhi Zhang

**Department of Statistics, University of Nebraska-Lincoln,

Lincoln, Nebraska, USA **Yajun Mei

**H. Milton Stewart School of Industrial and Systems Engineering,

Georgia Institute of Technology, Atlanta, Georgia, USA **Jianjun Shi

**H. Milton Stewart School of Industrial and Systems Engineering,

Georgia Institute of Technology, Atlanta, Georgia, USA

————————————————————————————————-

000Address correspondence to Ruizhi Zhang, Department of Statistics, University of Nebraska-Lincoln, Lincoln, Nebraska, 68583 USA; E-mail: rzhang35@unl.edu

Abstract: Robust change-point detection for large-scale data streams has many real-world applications in industrial quality control, signal detection, biosurveillance. Unfortunately, it is highly non-trivial to develop efficient schemes due to three challenges: (1) the unknown sparse subset of affected data streams, (2) the unexpected outliers, and (3) computational scalability for real-time monitoring and detection. In this article, we develop a family of efficient real-time robust detection schemes for monitoring large-scale independent data streams. For each data stream, we propose to construct a new local robust detection statistic called $L_{\alpha}$-CUSUM statistic that can reduce the effect of outliers by using the Box-Cox transformation of the likelihood function. Then the global scheme will raise an alarm based upon the sum of the shrinkage transformation of these local $L_{\alpha}$-CUSUM statistics so as to filter out unaffected data streams. In addition, we propose a new concept called false alarm breakdown point to measure the robustness of online monitoring schemes and propose a worst-case detection efficiency score to measure the detection efficiency when the data contain outliers. We then characterize the breakdown point and the efficiency score of our proposed schemes. Asymptotic analysis and numerical simulations are conducted to illustrate the robustness and efficiency of our proposed schemes.

Keywords: Robustness; Breakdown Point; Change Detection; Large-scale Data.

Subject Classifications: 62L15; 60G40.

1.  Introduction

Robust statistics have been extensively studied in the offline context when the entire data set is available for decision-making and is contaminated with outliers, e.g., robust estimation [1, 2], robust hypothesis testing [3, 4], and robust regression [5, 6]. Also, see the classical books, [7] or [8], for literature review. In this paper, we propose to develop robust methods in the context of sequential change-point detection when one is interested in detecting sparse, persistent smaller changes in large-scale data streams under the contamination of transient larger outliers. The problem of robust monitoring large-scale data streams in the presence of outliers occurs in many real-world applications such as industrial quality control, biosurveillance, key infrastructure, or internet traffic monitoring, in which sensors are deployed to constantly monitor the changing environment, see [9],[10],[11]. Unfortunately, it is highly non-trivial to develop efficient, robust real-time monitoring schemes or algorithms due to three challenges: (1) the sparsity, where only a few unknown data streams might be affected; (2) the robustness, where we are interested in detecting smaller persistent changes, not the larger transient outliers; and (3) the computational scalability, where the algorithms can be implemented recursively to make real-time decisions.

In the literature of sequential change-point detection for a large number of data streams, to the best of our knowledge, while the sparsity issue has been investigated, no research has been done on the robustness issue. To be more specific, the sparsity has been first addressed by [12] using a semi-Bayesian approach and later by [13] using shrinkage-estimation-based schemes. [14] developed asymptotic optimality theory for large-scale independent Gaussian data streams. Unfortunately, all these methods are sensitive to outliers since they are based on the likelihood function of specific parametric models (e.g. Gaussian) of the observations. Meanwhile, regarding the robustness issue, research is available for monitoring one- or low- dimensional streaming data such as rank-based method in [15, 16], kernel-based method in [17]. However, these nonparametric methodologies generally lose detection efficiency under specific parametric or semi-parametric models. By considering the worst-case of the outlier distribution, [18] formulated the problem of finding the optimal robust change detection procedure by solving a minimax problem. However, the resulting optimal test is based on the least-favorable-pair distributions of two uncertainty sets, which depends on the information of outliers. More importantly, it is unclear how to extend their method from monitoring a single data stream to monitoring multiple data streams when we also need to deal with the sparsity issue in which there is uncertainty on the subset of affected data streams.

In this paper, we develop efficient real-time monitoring schemes that are able to robustly detect smaller persistent changes in the presence of larger transient outliers when online monitoring of large-scale data streams. From the methodology viewpoint, our proposed schemes are semi-parametric and extend two contemporary concepts to the context of online monitoring of data streams: (i) $L_{q}$-likelihood [19, 20] for robustness, and (ii) the sum-shrinkage technique [21, 22] for sparsity. These allow us to develop statistically efficient and computationally simple schemes that can be implemented recursively over time for robust real-time monitoring of a large number of data streams. Moreover, we also extend the concept of breakdown in the offline robust statistics [23] to the sequential change-point detection context and conduct the false alarm breakdown point analysis, which turns out to be useful for the choices of tuning parameters in our proposed schemes.

We should point out that our contribution is not on the optimality theory but on the asymptotic properties of our proposed schemes that include the classical CUSUM-based procedures as a special case. Our research makes four contributions in the statistics field by combining robust statistics with sequential change-point detection for large-scale data streams. First, our proposed method is robust to infrequent outliers as well as the uncertainty of affected data streams. Second, our proposed method can be implemented recursively and distributed via parallel computing and thus is suitable for real-time monitoring over a long time period. Third, inspired by the concept of breakdown point [23] in the offline robust statistics, we propose a novel concept of false alarm breakdown point to quantify the robustness of any online monitoring schemes and show that our proposed schemes indeed have much larger false alarm breakdown point than the classical CUSUM-based schemes. Finally, from the mathematical viewpoint, we use Chebyshev’s inequality to derive non-asymptotic lower bounds on the average run length of false alarm of our proposed methods. The non-asymptotic results hold regardless of dimensionality and allow us to provide a deep insight into the effect of high-dimensionality in change-point detection under the modern asymptotic regime when the dimension or the number of data streams goes to $\infty.$

The remainder of this article is organized as follows. In Section 2, we start with problem formulations and model assumptions. In Section 3, we introduce our proposed family of robust monitoring schemes. In Section 4, the properties of the detection efficiency of our proposed schemes and the guideline to choose tuning parameters in our proposed schemes are provided. Then, we investigate the robustness of our proposed methods by conducting breakdown point analysis in Section 5. Simulation results are presented in Section 6. In Section 7, we conclude our paper with a few remarks The proofs of our main theorems are postponed to Appendix.

2.  Problem Formulation

Suppose we are monitoring $K$ independent data streams in a system.

[TABLE]

Under the classical change-point detection model for monitoring multi-streams (e.g., [12, 14, 24, 21, 22]), one assumes that the data $X_{k,n}$’s are initially independent and identically distributed (i.i.d.) with probability density function (pdf) $f_{0}(x)$. At some unknown time $\nu\geq 1,$ an undesired event occurs, and change the distributions of $m$ out of $K$ data streams, i.e., the affected local streams $X_{k,n}$’s have another distribution $f_{1}(x)$ when $n\geq\nu.$ The objective is to raise an alarm as soon as possible once a change occurs. Here, we refer to this classical model as the idealized model.

In this paper, we investigate the change-point detection problem under Tukey-Huber’s gross error model. As mentioned in the introduction, we want to raise an alarm as quickly as possible if there is a persistent distribution change on the data, but we prefer to take observations without any actions if there are only transient outliers. Mathematically, we assume the distribution of data $X_{k,n}$ might be changed from $h_{0}$ to $h_{1}$ at some change time $\nu,$ the $h_{0}$ and $h_{1}$ are the Tukey-Huber’s gross error model of the mixture densities

[TABLE]

where $\epsilon\in[0,1)$ is referred to as the contamination/outlier ratio, $g_{0}$ and $g_{1}$ are the (unknown) outlier distributions. Denote by ${\bf P}_{h_{0}}^{(\infty)}$ and ${\bf E}_{h_{0}}^{(\infty)}$ the probability measure and expectation when the data $X_{k,n}$’s are i.i.d. with the density $h_{0}$ when no change occurs, and denote by ${\bf P}_{h_{1}}^{(\nu)}$ and ${\bf E}_{h_{1}}^{(\nu)}$ the same when the change occurs at time $\nu$ and $m$ out of $K$ streams $X_{k,n}$’s have the post-change distribution $h_{1}$.

As in the classical sequential change-point problem, a statistical procedure under our setting is defined as a stopping time $T$ that represents the time when we raise an alarm to declare that a change has occurred. Here $T$ is an integer-valued random variable, and the decision $\{T=t\}$ is based only on the observations in the first $t$ time steps. To evaluate the performance of the detection procedure $T$ under Tukey-Huber’s gross error model when the outlier distributions $g_{0}$ and $g_{1}$ in (2.4) are unknown, we first assume the average run length to false alarm of the procedure $T$ is controlled under the idealized model. That is, we assume that the procedure $T$ is designed to satisfy the false alarm constraint

[TABLE]

for some pre-specified value $\gamma>0.$ We then investigate the robustness and the detection efficiency of the monitoring procedure under the gross error model in (2.4).

First, we propose quantifying the robustness of a monitoring procedure $T$ under the gross error model in (2.4) by borrowing the concept of breakdown point analysis from the offline robust statistics literature. To be more specific, we propose to define a new concept called false alarm breakdown point, which characterizes the minimal percentage of outliers that can make the false alarm rate under the gross error model $h_{0}=(1-\epsilon)f_{0}(x)+\epsilon g_{0}(x)$ very different from that under the idealized model $f_{0}.$

The false alarm breakdown point $\epsilon^{*}(T)$ of a family of monitoring schemes $T(b)$’s is defined as

[TABLE]

where ${\bf E}^{(\infty)}_{f_{0}}(T(b_{\gamma}))\sim\gamma$ as $\gamma\to\infty$, and the set $\hbar_{0,\epsilon}$ is the $\epsilon$-contaminated distribution density class of the idealized model $f_{0}(x)$ for given $\epsilon\in[0,1):$

[TABLE]

and $\mathcal{G}$ denotes the class of all probability densities of the data $X_{k,n}$.

Roughly speaking, the false alarm breakdown point characterizes the minimal percentage of outliers that can make the designed average run length to false alarm $\gamma$ unreliable. Thus, a scheme with larger breakdown points is more robust.

Second, we quantify the detection efficiency of the monitoring procedure $T$ under the gross error model in (2.4). For that purpose, recall that under the Lorden’s minimax criteria [25], the worst-case detection delay under $h_{1}$ is defined as

[TABLE]

Here $\mathcal{F}_{\nu-1}=(X_{1,[1,\nu-1]},\ldots,X_{K,[1,\nu-1]})$ denotes past global information at time $\nu.$ $X_{k,[1,\nu-1]}=(X_{k,1},\ldots,X_{k,\nu-1})$ is past local information for the $k$-th data stream. However, since the outlier distribution $g_{1}$ is unknown, we propose to define two quantities on detection efficiency: one is asymptotic efficiency score defined by

[TABLE]

and the other is the worst-case asymptotic efficiency score defined by

[TABLE]

In both definitions, $b_{\gamma}$ is a threshold of $T=T(b_{\gamma})$ so that ${\bf E}^{(\infty)}_{f_{0}}(T(b_{\gamma}))\sim\gamma.$ Clearly, when the data contain outliers, the procedure with a larger asymptotic efficiency score implies more efficiency in detecting the persistent change. Note that the definition of the asymptotic efficiency in (2.9) depends on the outlier distributions $g_{0}$ and $g_{1},$ which are unknown in practice, but the worst-case detection efficiency $\text{WAE}(T,\epsilon)$ in (2.10) measures the worst case among all set $\mathcal{G}$ of outlier distributions $g_{0}$ and $g_{1}.$ Note when we are monitoring a single data stream, i.e., the dimension $K=1,$ the optimal procedure that maximizes $\text{WAE}(T,\epsilon)$ is a CUSUM procedure constructed by a least-favorable-pair $g^{*}_{0},g^{*}_{1},$ as shown in [18]. However, the problem of finding the optimal procedure that minimizes $\text{WAE}(T,\epsilon)$ becomes more complicated when the dimension $K$ is large and the set of affected data streams is unknown.

In this paper, our objective is to develop a family of efficient, robust monitoring schemes that have a large breakdown point $\epsilon^{*}(T)$ in (2.6) and a large worst-case asymptotic efficiency score $\text{WAE}(T,\epsilon)$ in (2.10) subject to the constraints that this family of schemes satisfy the false alarm constraint in (2.5) under the idealized pre-change distribution $f_{0}.$

3.  Our proposed method

In this section, we will present our proposed schemes. At the high-level, our proposed schemes include two components: (i) robust monitoring each local data stream individually in parallel, and then (ii) combining local detection statistics to make an online global-level decision. For the purpose of easy understanding, we split the presentation of our proposed schemes into two subsections, and each subsection focuses on one component of the proposed scheme.

3.1.  Robust local statistics

For the $k^{th}$ data stream, we propose to define a new local $L_{\alpha}$-CUSUM statistic:

[TABLE]

for $n\geq 1,$ and $W_{\alpha,k,0}=0.$ Here $\alpha\geq 0$ is a tuning parameter that can control the tradeoff between statistical efficiency and robustness under the gross error model in (2.4) and its suitable choice will be discussed later.

The motivation of our $L_{\alpha}$-CUSUM statistic in (3.1) is as follows. Recall that when locally monitoring the single $k^{th}$ data stream $X_{k,n}$ with a possible local distribution change from $f_{0}$ to $f_{1},$ the generalized likelihood ratio test becomes the classical CUSUM statistic $W_{k,n}^{*},$ which has a recursive form:

[TABLE]

The CUSUM statistic enjoys nice optimality properties when all models are fully correctly specified [26], but unfortunately it is very sensitive to the outliers as in all other likelihood based methods in offline statistics. One recent idea in offline robust statistics is to replace the log-likelihood statistic $\log f(X)$ by $L_{\alpha}$-likelihood statistic $([f(X)]^{\alpha}-1)/\alpha$ for some $\alpha>0,$ see [19],[20]. At the high-level, $L_{\alpha}$-likelihood statistic $\frac{[f(X)]^{\alpha}-1}{\alpha}$ is always bounded below by $-1/\alpha$ whereas the log-likelihood statistic $\log f(X)$ could go to $-\infty$. Thus, the impact of outliers is bounded for the $L_{\alpha}$-likelihood statistic but unbounded for the log-likelihood statistic. Moreover, as $\alpha\to 0,$ the $L_{\alpha}$-likelihood function converges to the log-likelihood statistic, and thus it keeps statistical efficiencies when $\alpha$ is small. Here we apply this idea to develop our $L_{\alpha}$-CUSUM statistic. More rigorous robust properties will be discussed later in Section 5.

3.2.  Efficient global monitoring statistics

With local $L_{\alpha}$-CUSUM statistics $W_{\alpha,k,n}$ in (3.1) for each local stream, it is important to fuse these local statistics together smartly so as to address the sparsity issue. Here we propose to combine these local statistics together via the sum-shrinkage technique in [21], i.e., we raise a global-level alarm at time

[TABLE]

where $h(\cdot)\geq 0$ are some suitable shrinkage transformation functions, and $b>0$ is a pre-specified constant. Intuitively, the shrinkage functions $h(\cdot)$’s in (3.3) play the role of dimension reduction by automatically filtering out those non-changing local data streams and by keeping only those local streams that might provide information about the changing event. This will allow us to improve the detection power in the sparsity scenario when only a few local features are involved in the change.

For the purpose of illustration, here we focus on two kinds of shrinkage functions: one is the soft-thresholding function $h(x)=\max\{x-d,0\},$ and the other is the order-thresholding function $h(x)=x{\bf 1}\{x\geq w_{(r)}\},$ where $w_{(r)}$ is the $r$-th largest statistic of $w_{1},\cdots,w_{K}.$ Then the corresponding two global monitoring schemes are defined by

[TABLE]

where $W_{\alpha,(1),n}\geq W_{\alpha,(2),n}\geq\ldots\geq W_{\alpha,(K),n}$ are the order statistics of the $K$ local $L_{\alpha}$-CUSUM statistics $W_{\alpha,1,n},\ldots,W_{\alpha,K,n}.$

One can also consider other shrinkage functions such as the detectability score tansformation $h(x)=\log\left[1-p_{0}+0.64p_{0}\exp(x/2)\right]$ proposed in [14]. This yields another global monitoring scheme

[TABLE]

Our extensive numerical simulation experiences illustrate that for a given $\alpha,$ the scheme $N_{Chan,\alpha}(b)$ in (3.2) has the similar statistical/robustness properties to those schemes $N^{(soft)}_{\alpha}(b,d)$ and $N^{(r)}_{\alpha}(b)$ in (3.4) and (3.5) in many interesting sparse post-change scenarios when $p_{0}=r/K.$ This is because all these procedures utilize the same local $L_{\alpha}$-CUSUM statistics $W_{\alpha,k,n}$ in (3.1) and aim to detect the same post-change scenarios (after regularization).

Besides these aforementioned shrinkage transformations, there are other approaches to combine the local detection statistics together to make a global alarm. Two popular approaches in the literature are the “MAX” and the “SUM” schemes, see [27] and [28]:

[TABLE]

On one hand, the “MAX” and the “SUM” schemes could be considered as the special cases of our proposed top-r based scheme $N^{(r)}_{\alpha}(b)$ in (3.5) when $r=1$ and $r=K$ respectively. On the other hand, the “MAX” and “SUM” approaches are generally statistically inefficient unless in extreme cases of very few or many affected local data streams.

Note that there are three tuning parameters in our proposed schemes: $(\alpha,d,b)$ for the schemes $N^{(soft)}_{\alpha}(b,d)$ in (3.4) and $(\alpha,r,b)$ for the scheme $N^{(r)}_{\alpha}(b)$ in (3.5). It is natural to ask what are the “optimal” choices of these tuning parameters. It turns out that the most challenging one is the optimal choice of the common parameter $\alpha$, which is related to the robustness from the gross error models in (2.4), and will be discussed in Section 5. Next, the “optimal” choice of the shrinkage parameter $d$ or $r$ mainly depends on the number of affected local data streams, see our asymptotic properties in the next section. Finally, the choice of the threshold $b$ is straightforward for given two other parameters since it can be chosen to satisfy the false alarm constraint in (2.5).

4.  Worst-case asymptotic efficiency score

In this section, we derive the worst-case asymptotic efficiency score (2.10) of our proposed schemes $N^{(soft)}_{\alpha}(b,d)$ in (3.4) and $N^{(r)}_{\alpha}(b)$ in (3.5). To see that, we first report two standard change-point detection properties of our proposed schemes: the ARL to false alarm and detection delay under the gross error model $h_{i}=(1-\epsilon)f_{i}+\epsilon g_{i},$ where the outlier distributions $g_{i}$ are given and $i=0,1.$ Then, we will look at the worst-case of the outlier distributions to derive the worst-case asymptotic efficiency. It is important to note that our proposed schemes do not involve the contamination ratio $\epsilon$ or the information of outliers $\epsilon,g.$ Finally, based on our detection delay analysis, we provide guidelines on how to choose the tuning parameters in our proposed schemes. The proofs of the theorems are presented in the Appendix.

Let us begin with the definition of the expectation of the $L_{\alpha}$-likelihood ratio statistic $Y=([f_{1}(X)]^{\alpha}-[f_{0}(X)]^{\alpha})/\alpha$ when $X$ is distributed according to $h_{i}=(1-\epsilon)f_{i}+\epsilon g_{i}$ for given $\epsilon,$ $g_{i},$ and $i=0,1$ Note that when $\alpha=0,$ the variable $Y$ should be treated as the log-likelihood ratio $\log(f_{1}(X)/f_{0}(X))$.

Definition 4.1**.**

Given $\epsilon\geq 0$ and $\alpha\geq 0,$ for $i=0,1,$ define

[TABLE]

Note when $\epsilon=0,$ $I_{i}(\epsilon=0,\alpha;g_{i})$ does not depend on $g_{i}.$ So we further denote $I_{i}(\alpha):=I_{i}(\epsilon=0,\alpha;g_{i})$ for simplification. It turns out that the ARL to false alarm and detection delay of our proposed schemes are depending on whether $I_{i}(\epsilon,\alpha;g_{i})<0$ or $>0.$ Next, let us summarize the false alarm properties of our proposed schemes under the gross error model $h_{0}=(1-\epsilon)f_{0}+\epsilon g_{0}.$

Theorem 4.1**.**

Assume $I_{0}(\epsilon,\alpha;g_{0})<0,$ then there exists a unique positive constant $\lambda(\epsilon,\alpha;g_{0})$ depends on $f_{0},f_{1},g_{0},\alpha,\epsilon$ such that

[TABLE]

With the constant $\lambda(\epsilon,\alpha;g_{0})>0$ in (4.2), the ARL to false alarm of our proposed schemes, $N^{(soft)}_{\alpha}(b,d)$ in (3.4) and $N^{(r)}_{\alpha}(b)$ in (3.5), are given as follows under different sufficient conditions:

(a) When $\lambda(\epsilon,\alpha;g_{0})b>K\exp\{-\lambda(\epsilon,\alpha;g_{0})d\},$ we have

[TABLE]

(b) When $\lambda(\epsilon,\alpha;g_{0})b>K,$ we have

[TABLE]

Let us add some comments to better understand the theorem. First, the existence of the unique constant $\lambda(\epsilon,\alpha;g_{0})>0$ in (4.2) is based on the assumption that $I_{0}(\epsilon,\alpha;g_{0})<0$, see Appendix A2 of [29]. Moreover, when $\epsilon=0,$ both $I_{0}(0,\alpha;g_{0})$ and $\lambda(0,\alpha;g_{0})$ only depend on the idealized model $f_{i}(x)$ and $\alpha,$ but do not depend on the information of outliers, i.e., $\epsilon$ and $g_{i}.$ For simplification, we denote $\lambda(\alpha)=\lambda(0,\alpha;g_{0}).$

Second, our rigorous, non-asymptotic results in (4.3) and (4.4) hold no matter how large the number $K$ of data streams is. This allows us to investigate the modern asymptotic regime when the dimension $K$ goes to $\infty.$

Finally, the assumptions of $\lambda(\epsilon,\alpha;g_{0})b>K\exp\{-\lambda(\epsilon,\alpha;g_{0})d\}$ or $\lambda(\epsilon,\alpha;g_{0})b>K$ essentially says that the global threshold $b$ of our proposed schemes should be large enough if one wants to control the global false alarm rate when online monitoring large-scale streams. These results allow us to find a conservative threshold $b$ so as to satisfy the false alarm constraint in (2.5), also see the details of parameter setting below.

Next, the following theorem summarizes the detection delays of our proposed schemes $N^{(soft)}_{\alpha}(b,d)$ in (3.4) and $N^{(r)}_{\alpha}(b)$ in (3.5) when $m$ out of $K$ features are affected by the occurring event for some given $1\leq m\leq K.$ The detailed proof of Theorem 4.2 will be presented in Appendix.

Theorem 4.2**.**

Suppose $I_{1}(\epsilon,\alpha;g_{1})>0$ , and $m$ out of $K$ features are affected.

(a) If $b/m+d$ goes to $\infty$, then the detection delay of $N^{(soft)}_{\alpha}(b,d)$ satisfies

[TABLE]

(b) If $r\geq m$ and $b/m$ goes to $\infty$, then the detection delay of $N^{(r)}_{\alpha}(b)$ satisfies

[TABLE]

where the $o(1)$ term does not depend on the dimension $K$, but might depend on $m$ and $\alpha$ as well as the distributions $h_{1}.$

To simplify the notation, we use $N_{\alpha}$ to denote both the scheme $N^{(soft)}_{\alpha}(b,d)$ and scheme $N^{(r)}_{\alpha}(b).$ By Theorem 4.1 and Theorem 4.2, when $K$ is fixed, if $I_{0}(\epsilon,\alpha;g_{0})<0$ and $I_{1}(\epsilon,\alpha;g_{1})>0,$ we can get a natural lower bound of the asymptotic efficiency score (2.9) of our proposed schemes,

[TABLE]

However, if we can find outlier distributions $g^{*}_{0},g^{*}_{1}$ such that $I_{0}(\epsilon,\alpha;g^{*}_{0})>0$ and $I_{1}(\epsilon,\alpha;g^{*}_{1})<0,$ we will get

[TABLE]

which implies the procedure cannot detect the persistent change from $f_{0}$ to $f_{1}$ at all due to the contamination of outliers.

Now, we are ready to present the worst-case asymptotic efficiency score of our proposed scheme $N_{\alpha}.$ First, assume I_{0}(\alpha)={\bf E}_{f_{0}}\Big{[}\frac{[f_{1}(X)]^{\alpha}-[f_{0}(X)]^{\alpha}}{\alpha}\Big{]}<0 and I_{1}(\alpha)={\bf E}_{f_{1}}\Big{[}\frac{[f_{1}(X)]^{\alpha}-[f_{0}(X)]^{\alpha}}{\alpha}\Big{]}>0, denote

[TABLE]

Then we have the following theorem:

Theorem 4.3**.**

For our proposed scheme $N_{\alpha}(b)$with $\alpha\geq 0,$ suppose $K$ is fixed and $b\to\infty,$

(a) if $\epsilon<-I_{0}(\alpha)/[M^{*}(\alpha)-I_{0}(\alpha)]$ and $\epsilon<I_{1}(\alpha)/[M^{*}(\alpha)+I_{1}(\alpha)],$ we have

[TABLE]

where $\lambda^{*}(\epsilon,\alpha)=\underset{g_{0}\in\mathcal{G}}{\inf}\lambda(\epsilon,\alpha;g_{0})>0.$

(b) Otherwise, $\text{WAE}(N_{\alpha},\epsilon)=0.$

Note if $\log(f_{1}(x)/f_{0}(x))$ is unbounded, we have $M^{*}(0)=+\infty.$ Based on Theorem 4.3, for any $\epsilon>0,\text{WAE}(N_{\alpha=0},\epsilon)=0,$ which implies the CUSUM based method cannot detect the persistent change at all under any percentage of outliers. However, if both $f_{0},f_{1}$ are bounded, for any $\alpha>0,$ we have $M^{*}(\alpha)<+\infty.$ Thus, our proposed schemes $N_{\alpha}$ will always have a positive worst-case asymptotic efficiency score when the contamination ratio $\epsilon$ is small. This implies the detection efficiency of our proposed schemes under the gross error model.

Note that there are three tuning parameters in our proposed schemes: $(\alpha,d,b)$ for the schemes $N^{(soft)}_{\alpha}(b,d)$ in (3.4) and $(\alpha,r,b)$ for the scheme $N^{(r)}_{\alpha}(b)$ in (3.5). It is natural to ask what are the “optimal” choices of these tuning parameters. It turns out that Theorems 4.1 and 4.2 provide the optimal choices of $(d,b)$ or $(r,b)$ that asymptotically minimize the detection delay subject to the false alarm constraint $\gamma$ in (2.5). Below we will report the corresponding results, and detailed proofs and descriptions are postponed to Appendix.

(1) The optimal choice of parameter $\alpha$, which turns out to be the most challenging one, as it is related to the robustness from the gross error models in (2.4). We will discuss in more details in Section 5 through the concept of false alarm breakdown point. The result in Section 5 shows the optimal $\alpha_{opt}$ only depends on the distributions $f_{0},f_{1}$ but independent of other parameters $d,r,b,$ and outliers information $\epsilon,g.$

(2) Given $\alpha_{opt},$ the choice of the shrinkage parameter $d$ or $r$ mainly depends on the number $m$ of affected local feature coefficients. If we want to minimize the detection delay subject to the false alarm constraint $\gamma$ in (2.5), we can set $r=m$ for the scheme $N^{(r)}_{\alpha_{opt}}(b)$ in (3.5). The optimal choice of $d$ for the proposed scheme $N^{(soft)}_{\alpha_{opt}}(b,d)$ in (3.4) is a little complicated, and given by

[TABLE]

where $\lambda(\alpha_{opt})$ is defined in (4.2) and only depends on $f_{0},f_{1}$ and $\alpha_{opt}.$

(3) The choice of the threshold $b$ is straightforward for given two other parameters, since it can be chosen to satisfy the false alarm constraint in (2.5) under the idealized distribution $f_{0}.$ A choice of global detection threshold

[TABLE]

will guarantee that our proposed scheme $N^{(soft)}_{\alpha_{opt}}(b_{\gamma},d_{opt})$ satisfies the global false alarm constraint $\gamma$ in the idealized model as in (2.5).

Note that all these choices of parameters do not depend on the $\epsilon$ or $g$, and only depend on the idealized model $f_{0},f_{1}$ and a prior knowledge on the number $m$ of affected data streams.

5.  Breakdown point analysis

In this section, we will investigate the robustness properties of our proposed schemes, $N^{(soft)}_{\alpha}(b,d)$ in (3.4) and $N^{(r)}_{\alpha}(b)$ in (3.5), through the false alarm breakdown point analysis. This will provide the guideline on how to choose the tuning parameter $\alpha,$ which controls the robustness of our proposed schemes.

In the classical offline robust statistics, the breakdown point is one of the most popular measures of robustness of statistical procedures. At a high-level, in the context of finite samples, the breakdown point is the smallest percentage of contaminations that may cause an estimator or statistical test to be really poor. Since the pioneering work of [23] for the asymptotic definition of breakdown point, much research has been done to investigate the breakdown point for different robust estimators or hypothesis testings in the offline statistics, see [30], [31]. To the best of our knowledge, no research has been done on the breakdown point analysis under the online monitoring or change-point context.

Given the importance of the system-wise false alarm rate for online monitoring large-scale data streams in real-world applications, here we focus on the breakdown point analysis for false alarms. Intuitively, for a family of procedures $T(b)$ that is robust, if it is designed to satisfy the false alarm constraint $\gamma$ in (2.5) under the idealized model $f_{0},$ then its false alarm rate should not be too bad under the gross error model $h_{0}$ with some small amount of outliers. There are two specific technical issues that require further clarification. First, how bad is a “bad” false alarm rate? We propose to follow the sequential change-point detection literature to assess the false alarm rate by $\log{\bf E}^{(\infty)}_{h_{0}}(T(b))$ and deem the false alarm rate unacceptable if $\log{\bf E}^{(\infty)}_{h_{0}}(T(b))$ is much smaller than the designed level of $\log\gamma$, i.e., if $\log{\bf E}^{(\infty)}_{h_{0}}(T(b))=o(\log\gamma).$ Second, what kind of the contamination function $g$ in (2.7) should we consider in the gross error model? Here we propose to follow the offline robust statistics literature to consider the worst-case scenario in the $\epsilon$-contaminated distribution class in [1] that includes any arbitrary contamination functions $g$’s, which leads to the definition of the false alarm breakdown point in (2.6).

Now we are ready to conduct the false alarm breakdown point analysis for our proposed schemes $N^{(soft)}_{\alpha}(b,d)$ and $N^{(r)}_{\alpha}(b)$ with a given tuning parameter $\alpha\geq 0$. To do so, for the densities $f_{0}(x)$ and $f_{1}(x),$ and for any given $\alpha\geq 0,$ we define an intrinsic bound

[TABLE]

and the density power divergence between $f_{0}$ and $f_{1}$:

[TABLE]

Note that $d_{\alpha}(f_{0},f_{1})$ was proposed in [2], which showed that it is always positive when $f_{1}$ and $f_{0}$ are different. Moreover, when $\alpha=0,$ $d_{\alpha=0}(f_{0},f_{1})$ becomes Kullback-Leibler information number $\int f_{0}(x)\log\frac{f_{0}(x)}{f_{1}(x)}dx.$

With these two new notations, the following theorem derives the false alarm breakdown point of our proposed schemes $N^{(soft)}_{\alpha}(b,d)$ and $N^{(r)}_{\alpha}(b)$ as a function of the tuning parameter $\alpha$ for a fixed soft-thresholding parameter $d$ and $r$ when online monitoring a given $K$ number of data streams. Since they have the same breakdown point, to simplify the notation, we use $N_{\alpha}$ to denote both the scheme $N^{(soft)}_{\alpha}(b,d)$ and scheme $N^{(r)}_{\alpha}(b).$

Theorem 5.1**.**

Suppose that $f_{\theta}(x)=f(x-\theta)$ is a location family of density function with continuous probability density function $f(x),$ and assume $f_{\theta_{0}}(x)-f_{\theta_{1}}(x)$ takes both positive and negative values for $x\in(-\infty,+\infty)$. For $\alpha\geq 0,$ and any fixed $d$ and $K,$ the false alarm breakdown point of our proposed schemes $N_{\alpha}$ in (3.4) and (3.5) is the same and given by

[TABLE]

where $M(\alpha)$ and $d_{\alpha}(f_{\theta_{0}},f_{\theta_{1}})$ are defined in (5.1) and (5.2). In particular, $\epsilon^{*}(N_{\alpha})=0$ if $M(\alpha)=\infty$ and $d_{\alpha}(f_{\theta_{0}},f_{\theta_{1}})$ is finite.

The proof of Theorem 5.1 requires the asymptotic properties of our proposed schemes $N^{(soft)}_{\alpha}(b,d)$ in (3.4) and $N^{(r)}_{\alpha}(b)$ in (3.5) under the assumption that $\epsilon$ and $g$ are given, which has been studied in the previous section. The detailed proof of Theorem 5.1 will be presented in the supplementary materials.

Next, let us apply Theorem 5.1 to guide us to choose the optimal robustness parameter $\alpha$. Since the false alarm breakdown point of our proposed schemes do not require any information about the contamination ratio $\epsilon$ and contamination distribution $g,$ one nature idea is to maximize the false alarm breakdown point in (5.3):

[TABLE]

As an illustration, let us see the results of (5.3) and (5.4) for widely used normal distributions, i.e., when $f_{\theta}$ is the pdf of $N(\theta,\sigma^{2}).$ In this case, when $\alpha=0,$ the density power divergence $d_{\alpha=0}(f_{\theta_{0}},f_{\theta_{1}})=\frac{1}{2\sigma^{2}}(\theta_{1}-\theta_{0})^{2}$ is finite, but the bound $M(\alpha=0)$ in (5.1) becomes $+\infty$ since it is the supremum of the log-likelihood ratio $\log f_{\theta_{1}}(x)-\log f_{\theta_{0}}(x)=(\theta_{1}-\theta_{0})x-(\theta_{1}^{2}-\theta_{0}^{2})/2$ over $x\in(-\infty,\infty).$ Hence, $\epsilon^{*}(N_{\alpha=0})=0.$ That is, the false alarm breakdown point of the baseline CUSUM-based scheme $N_{\alpha=0}$ is $0,$ i.e., any amount of outliers will deteriorate the false alarm rate of the classical CUSUM statistics-based schemes. This is consistent with the offline robust statistics literature that the likelihood-function based methods are very sensitive to model assumptions and are generally not robust.

Meanwhile, for any $\alpha>0,$ note that

[TABLE]

and thus it is not difficult to derive from (5.2) that

[TABLE]

Moreover, if we let $M(=1/\sqrt{2\pi\sigma^{2}}),$ then $|f_{\theta}(x)|\leq M$ for all $x.$ By the definition in (5.1), we have $|M(\alpha)|\leq 2M^{\alpha}/\alpha,$ which is finite for any $\alpha>0.$ This implies that for normal distributions, $\epsilon^{*}(N_{\alpha})>0$ for any $\alpha>0.$ Thus our proposed $L_{\alpha}$-CUSUM based scheme with $\alpha>0$ is much more robust than the classical CUSUM scheme.

To see the optimal choice of $\alpha$ based on (5.4), let us consider a concrete numerical example when $f_{\theta_{0}}\sim N(0,1)$ and $f_{\theta_{1}}\sim N(1,1).$ By (5.5), we can compute the value $d_{\alpha}(0,1)$ for any $\alpha\geq 0.$ While we do not have analytic formula for the upper bound $M(\alpha)$ in (5.1), its numerical value can be easily found by brute-force exhaustive search over the real line $x\in(-\infty,\infty).$ The result shows the false alarm breakdown point of our proposed scheme $N_{\alpha}$ will first increase and then decrease as $\alpha$ varies from [math] to $2.$, and yields the optimal choice of $\alpha_{opt}$ as $0.51,$ with corresponding breakdown point as $0.233.$ That means our proposed scheme with the choice of $\alpha=0.51$ could tolerate $23.3\%$ arbitrarily bad observations in terms of keeping the designed false alarm constraint stable.

Finally, we should emphasize that the optimal value $\alpha_{opt}$ in (5.4) and the false alarm breakdown point $\epsilon^{*}(N_{\alpha})$ in (5.3) will generally depend on the change magnitude or signal-to-noise-ratio. To illustrate this, we consider three families: Normal, Laplace, and Logistic distributions with the scale parameter $\sigma=1.$ This yields three families of pdfs, $f_{\theta}(x)=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(x-\theta)^{2}}{2}),$ $\frac{1}{2}\exp(-|x-\theta|)$, or $\frac{\exp(-(x-\theta))}{(1+\exp(-(x-\theta)))^{2}}.$ In each case, we assume that the pre-change parameter $\theta_{0}=0,$ the designed post-change parameter $\theta_{1}$ varies from $1$ to $5.$ In Figures 2 and 2, we plot the optimal value $\alpha_{opt}$ and the corresponding false alarm breakdown point $\epsilon^{*}(N_{\alpha})$ as a function of $\theta_{1}.$ Figure 2 implies that with the increasing of the post-change $\theta_{1}$ or the signal-to-noise-ratio, our proposed robust schemes with optimal $\alpha$ can tolerate more outliers. Also it is interesting to see from Figure 2 that the optimal $\alpha_{opt}$ decreases for normal or logistic distribution as the post-change parameter $\theta_{1}$ increases. A surprising result is that the optimal $\alpha_{opt}=0$ for the Laplace distribution. This implies the classical CUSUM procedure for Laplace distribution is actually optimal in the sense of having the largest breakdown point. One possible explanation is that for the Laplace distribution, the log-likelihood ratio $\log(f_{\theta_{1}}(x)/f_{\theta_{0}}(x))=-|x-\theta_{1}|+|x-\theta_{0}|$ takes values in the interval $[\theta_{0}-\theta_{1},\theta_{1}-\theta_{0}]$ when $\theta_{1}>\theta_{0}.$ Thus the impact of outliers is directly controlled.

6.  Numerical Simulations

In this section we conduct numerical simulation studies to illustrate the robustness and efficiency of our proposed schemes $N^{(soft)}_{\alpha}(b,d)$ and $N^{(r)}_{\alpha}(b)$.

In our simulation studies, we assume there are $K=100$ independent data streams, and at some unknown time, $m=10$ features are affected by the occurring event. Also the change is instantaneous if a stream is affected, and we do not know which subset of streams will be affected. We set $f_{\theta}=$ pdf of $N(\theta,1),$ the pre-change parameter $\theta_{0}=0,$ the post-change parameter $\theta_{1}=1,$ and the contamination densities $g_{0},g_{1}$ are pdfs of $N(0,3^{2}).$ Our proposed schemes $N^{(soft)}_{\alpha}(b,d)$ in (3.4) and $N^{(r)}_{\alpha}(b)$ in (3.5) are constructed by using the density function $f_{\theta_{0}}$ and $f_{\theta_{1}}.$

In the first simulation study, we consider the idealized model when $\epsilon=0.$ In this case, for our proposed robust scheme $N^{(soft)}_{\alpha}(b,d)$ in (3.4), as shown in the previous section, the optimal choices of $\alpha_{opt}=0.51.$ By (4.11), if $\log(\gamma)<<K,$ then the corresponding optimal shrinkage parameters $d\approx\frac{1}{\lambda(\epsilon=0,\alpha=0.51)}\log\frac{K}{m}=0.8915$ for $K=100$ and $m=10,$ since $\lambda(\epsilon=0,\alpha=0.51)=2.5829.$ For our proposed robust scheme $N^{(r)}_{\alpha}(b)$ in (3.5), we choose $\alpha=\alpha_{opt}=0.51$ and $r=10.$ For the baseline CUSUM-based scheme, i.e., $N^{(soft)}_{\alpha=0}(b,d)$ with $\alpha=0$, we choose the shrinkage parameter $d=\frac{1}{\lambda(\epsilon=0,\alpha=0)}\log\frac{K}{m}=2.3026,$ since $\lambda(\epsilon=0,\alpha=0)=1.$

In summary, we will compare the following different schemes.

• •

Our proposed scheme $N^{(soft)}_{\alpha}(b,d)$ in (3.4) with $\alpha_{opt}=0.51$ and $d=0.8915.$

• •

Our proposed scheme $N^{(r)}_{\alpha}(b)$ in (3.5) with $\alpha_{opt}=0.51$ and $r=10.$

• •

The baseline CUSUM-based scheme $N^{(soft)}_{\alpha=0}(b,d)$ with $d=2.3026.$

• •

The MAX scheme $N_{\alpha=0.51,\max}(b)$ in (3.6);

• •

The SUM scheme $N_{\alpha=0.51,\text{sum}}(b)$ in (3.7);

• •

The method $N_{XS}(b,p_{0}=0.1)$ in [12] based on generalized likelihood ratio:

[TABLE]

where for all $1\leq k\leq K,0\leq i<n,$

[TABLE]

• •

The method $N_{Chan,\alpha=0}(b,p_{0}=0.1)$ in [14] under the idealized model that is an extension of the SUM scheme in [28]:

[TABLE]

where $W_{k,n}^{*}$ is the CUSUM statistics in (3.2).

• •

The method $N_{Chan,\alpha=0.51}(b,p_{0}=0.1)$ in (3.2) which is similar to $N_{Chan,\alpha=0}$ but replace the CUSUM statistic by our proposed $L_{\alpha}$-CUSUM statistic.

For each of these schemes $T(b)$, we first find the appropriate values of the threshold $b$ to satisfy the false alarm constraint $\gamma\approx 5000$ under the idealized model with $\epsilon=0$ (within the range of sampling error). Next, using the obtained global threshold value $b,$ we simulate the detection delay when the change-point occurs at time $\nu=1$ under several different post-change scenarios, i.e., different number of affected sensors. All Monte Carlo simulations are based on $1000$ repetitions.

Table 1 summarizes the detection delays of these $8$ schemes under $9$ different post-change hypothesis. Among all schemes, $N_{XS}(b,p_{0})$ generally yields the smallest detection delay. However, we want to emphasize that it is computationally expensive. Specifically, even if we use a time window of size $k$ as in [12] to speed up the implementation of $N_{XS}(b,p_{0})$, at each time $n,$ $O(Kk^{2})$ computations are needed to get the global monitoring statistics, whereas our proposed scheme $N^{(soft)}_{\alpha}(b,d)$ only require $O(K)$ computations to get the global monitoring statistics.

Another interesting observation from Table 1 is that the detection delay of our proposed robust schemes $N^{(soft)}_{\alpha=0.51}(b,d)$ and $N^{(r)}_{\alpha=0.51}(b)$ are not too bad compared with the CUSUM-based scheme $N^{(soft)}_{\alpha=0}(b,d=2.3026)$, and it just takes additional $1.7$ time steps to raise a correct global alarm under the idealized model when $m=10$ data streams are affected.

In the second simulation study, we will examine the detection efficiency of these schemes under the gross error model when $\epsilon=0.1.$ For each of these $8$ schemes, we use the same threshold $b$ obtained from the first simulation to guarantee these schemes satisfy the same false alarm constraint $\gamma=5000$ under the idealized model. Then, we wilsimulate the in-control average run and the detection delay of these schemes when both the pre-change distribution and post-change distribution are the gross error model in (2.4) with $\epsilon=0.1,$ $g_{0},g_{1}$ as pdfs of $N(0,3^{2}).$ We then report the empirical version of the asympototic efficiency score in (2.9) of these schemes under $8$ different post-change hypothesis in Table 2.

First, we can see our proposed scheme $N^{(soft)}_{\alpha=0.51}(b,d=0.8915)$ and $N^{(r=10)}_{\alpha=0.51}(b=18.7)$ have the largest detection efficiency score among all comparison methods when $10$ data streams are affected. Moreover, by using our proposed $L_{\alpha}$-CUSUM statistics with $\alpha=0.51,$ the method $N_{Chan,\alpha=0.51}(b,p_{0}=0.1)$ yields the similar detection efficiency to our proposed schemes. This illustrates that the improvement of $L_{\alpha}$-CUSUM statistics is significant as compared to the baseline CUSUM statistics in the presence of outliers.

It is also interesting to note that the MAX-scheme $N_{\alpha=0.51,\max}(b)$ and the SUM-scheme $N_{\alpha=0.51,\text{sum}}(b)$ are designed for the case when $m=1$ or $m=K$ features are affected, and Table 2 confirmed that their detection efficinecys are indeed the largest in their respective designed scenarios. However, when the number of affected features $m$ is moderate and is arround $10$, our proposed scheme $N^{(soft)}_{\alpha=0.51}(b,d)$ and $N^{(r)}_{\alpha=0.51}(b)$ have larger detection efficiency, which implies our proposed schemes with sum-shrinkage technique could be more robust to the number of affected features.

In the third experiment, we investigate the impact of contamination rate $\epsilon$ on the false alarms of different methods to illustrate the robustness of our proposed $L_{\alpha}$-CUSUM statistics. Since the top-r scheme $N^{(r)}_{\alpha=0.51}(b),$ MAX-scheme $N_{\alpha=0.51,\max}(b),$ the SUM-scheme $N_{\alpha=0.51,\text{sum}}(b)$ and $N_{Chan,\alpha=0.51}(b,p_{0})$ are all based on local $L_{\alpha}$-CUSUM statistics, their robustness properties are similar to our proposed scheme $N^{(soft)}_{\alpha=0.51}(b,d).$ To highlight the robustness of our proposed $L_{\alpha}$-CUSUM statistics, we only compare our proposed scheme $N^{(soft)}_{\alpha=0.51}(b,d)$ with other three schemes: $N^{(soft)}_{\alpha=0}(b,d)$, $N_{Chan,\alpha=0}(b,p_{0}),$ and $N_{XS}(b,p_{0}).$

Figure 3 reports the curve of $\log{\bf E}_{h_{0}}^{(\infty)}(T)$ as the contamination ratio $\epsilon$ varies from $0.02$ to $0.2$ with stepsize $0.02$. Clearly, all curves decrease with the increasing of contaminations, meaning that all schemes will raise false alarm more frequently when there are more outliers. However, the curves for the CUSUM or likelihood-ratio based methods decreased very quickly, whereas our proposed $L_{\alpha}$-CUSUM statistics-based method with $\alpha_{opt}=0.51$ decreases rather slowly. This suggests that our proposed scheme is more robust in the sense of keeping the designed ARL more stable with a small departure from the assumed model.

7.  Conclusion

In this paper, we study the problem of robust monitoring of large-scale data streams when the true observed data follow Huber’s gross error model. We develop a family of efficient and robust detection schemes that can be implemented in real-time. From the worst-case detection efficiency point of view, we show our proposed methods can still have positive detection efficiency under a small proportion of arbitrary outliers. In contrast, the CUSUM-based methods lose all detection efficiency once the data include outliers. From the robustness point of view, we propose a new concept called false alarm breakdown point, which measures the stability of the designed false alarm constraint of any monitoring procedures under the effects of outliers. Our breakdown point analysis implies our proposed methods can have positive breakdown points. We also provide detailed guidelines on the choices of tuning parameters in our detection procedures. However, in this work, we focus on the problem of monitoring homogeneous independent data streams. It is of future interest to extend the problem to nonhomogeneous data streams with some correlation structures.

Appendix

In this online supplementary material, we provide detailed proofs to Theorems 4.1, 4.2, and Theorem 4.3, the optimal parameter choice in Section 4, and the proof of Theorem 5.1.

A. Proof of Theorem 4.1

(a) For any $x\geq 0,$ by Chebyshev’s inequality,

[TABLE]

where $W_{\alpha,k}^{*}=\limsup_{n\rightarrow\infty}W_{\alpha,k,n}.$ We will show that $W_{\alpha,k}^{*}$ exists later, and when it does exist, it is clear that $W_{\alpha,k}^{*}$ are i.i.d. across different $k$ under the pre-change measure ${\bf P}_{h_{0}}^{(\infty)}.$ Now if we define the log-moment generating function of the $W_{\alpha,k}^{*}$’s

[TABLE]

for some $\theta\geq 0,$ then another round application of Chebyshev’s inequality yields

[TABLE]

for $\theta>0.$ Combining (7.1) and (7.6) yields that

[TABLE]

for all $x\geq 0.$ Since $x(1-xu)$ is maximized at $x=1/(2u)$ with the maximum value $1/(4u).$ We conclude from (7.8) that

[TABLE]

for any $\theta>0$ as long as $\psi_{\alpha}(\theta)$ in (7.5) is well-defined.

The remaining proof is to utilize the definition of $\lambda(\epsilon,\alpha;g_{0})>0$ in (4.2) to show that the upper limiting $W_{\alpha,k}^{*}$ of the proposed $L_{\alpha}$-CUSUM statistics is well-defined and derive a careful analysis of $\psi_{\alpha}(\theta)$ in (7.5). When $\alpha=0,$ the $L_{\alpha}$-CUSUM statistics become the classical CUSUM statistics, and the corresponding analysis is well-known, see [21]. Here our main insight is that our proposed $L_{\alpha}$-CUSUM statistics $W_{\alpha,k,n}$ for detecting a change from $h_{0}(x)$ to $h_{1}(x)$ in (2.4) can be thought of as the classical CUSUM statistic for detecting a local change from $h_{0}(x)$ to another new density function $h_{2}(x).$ Hence, under the pre-change hypothesis of $h_{0}(\cdot),$ the false alarm properties of our proposed $L_{\alpha}$-CUSUM statistics can be derived through those of the classical CUSUM statistics.

By the definition of $\lambda(\epsilon,\alpha;g_{0})>0$, if we define a new function

[TABLE]

then $h_{2}(x)$ is a well-defined probability density function. Then in the problem of detection a local change from $h_{0}(x)$ to $h_{2}(x),$ the local CUSUM statistics for the $k$th local data stream is defined recursively by

[TABLE]

Compared with our proposed $L_{\alpha}$-CUSUM statistics $W_{\alpha,k,n}$, it is clear that $W_{k,n}^{\prime}=\lambda(\epsilon,\alpha)W_{\alpha,k,n},$ and thus our proposed $L_{\alpha}$-CUSUM statistics $W_{\alpha,k,n}$’s are equivalent to the standard CUSUM statistics $W_{k,n}^{\prime}$ up to a positive constant $\lambda(\epsilon,\alpha;g_{0}).$ By the classical results on the CUSUM, see Appendix 2 on Page 245 of [32], as $n\rightarrow\infty,$ $W_{k,n}^{\prime}$ converges to a limit and thus $W_{\alpha,k,n}$ also converges to a limit, denoted by $W_{\alpha,k}^{*}.$ Moreover, the tail probability of $W_{\alpha,k}^{*}$ satisfies

[TABLE]

Now we shall use (7.11) to derive information bound of $\psi_{\alpha}(\theta)$ in (7.5). In order to simplify our arguments, we abuse the notation and simply denote $\lambda(\epsilon,\alpha;g_{0})$ by $\lambda$ in the remaining proof of the theorem. By the definition of $\psi_{\alpha,k}(\theta)$ in (7.5) and the tail probability $G(x)$ in (7.11), for $\theta>0,$

[TABLE]

where the second equation is based on the integration by parts. Clearly, relation (A. Proof of Theorem 4.1) holds for any $0<\theta<\lambda=\lambda(\epsilon,\alpha;g_{0}).$

By (7.9) and (A. Proof of Theorem 4.1), we have

[TABLE]

for all $0<\theta<\lambda=\lambda(\epsilon,\alpha;g_{0}).$ When $\lambda b>K\exp\{-\lambda d\},$ relation (4.3) follows at once from (7.13) by letting $\theta=\sqrt{\lambda/b}\left(\sqrt{\lambda b}-\sqrt{K\exp\{-d\lambda\}}\right)\in(0,\lambda).$ This completes the proof of Theorem 4.1 (a).

(b) Note $N^{(soft)}_{\alpha}(b,d=0)\leq N^{(r)}_{\alpha}(b)$ for any $b\geq 0.$ Therefore, (4.4) can be derived directly from (4.3) by letting $d=0$ in (4.3).

B. Proof of Theorem 4.2

First, we will prove the part (a) of Theorem 4.2. To prove the detection delay bound (4.6) in Theorem 4.2, without loss of generality, assume the first $m$ data streams are affected. Consider a new stopping time

[TABLE]

Clearly $N^{(soft)}_{\alpha}(b,d)\leq T^{\prime}(b,d),$ and thus

[TABLE]

Next, by the recursive definition of $W_{\alpha,k,n}$ in (3.1), using the same approach in Theorem 2 of [25] that connects the recursive CUSUM-type scheme to the random walks, we have

[TABLE]

where ${\bf E}_{1}$ denotes the expectation when the change happen at time $\nu=1,$ and $T^{\prime\prime}(b,d)$ is the first passage time when the random walk with i.i.d. increment of mean $mI_{1}(\epsilon,\alpha;g_{1})$ exceeds the bound $b+md,$ and is defined as

[TABLE]

By standard renewal theory, as ($\frac{b}{m}+d$)$\to\infty$, we have

[TABLE]

Relation (4.6) then follows at once from the above relations, which completes the proof of part (a) of Theorem 4.2.

To prove the part (b), we define another stopping time

[TABLE]

Note for the sorted statistics $W_{\alpha,(1),n}\geq W_{\alpha,(2),n}\geq\cdots\geq W_{\alpha,(K),n},$ we have $\sum_{k=1}^{m}W_{\alpha,k,n}\leq\sum_{k=1}^{m}W_{\alpha,(k),n}.$ Thus,when $m\leq r,$ $N^{(r)}_{\alpha}(b)\leq\tau(b).$ By standard renew theory, we have

[TABLE]

which completes the proof of part (b) of Theorem 4.2.

C. Proof of Theorem 4.3

Note if $\epsilon<-I_{0}(\alpha)/[M^{*}(\alpha)-I_{0}(\alpha)],$

[TABLE]

Therefore, by Theorem 4.1, there exists a postive number $\lambda^{*}(\epsilon,\alpha)=\underset{g_{0}\in\mathcal{G}}{\inf}\lambda(\epsilon,\alpha;g_{0})>0$ such that

[TABLE]

Moreover, if $\epsilon<I_{1}(\alpha)/[M^{*}(\alpha)+I_{1}(\alpha)],$

[TABLE]

By Theorem 4.2, we have

[TABLE]

Thus, by the definition of worst-case detection efficiency in (2.10), we have

[TABLE]

D. Parameter Setting in Section 4

The choice of $b=b_{\gamma}$ in (4.12) follows directly from Theorem 4.1 (a). To prove (4.11), we abuse the notation and use $\lambda$ to denote $\lambda(\alpha)$ for simplification. By Theorem 4.2, the optimal $d$ is the non-negative value that minimize the function

[TABLE]

This is an elementary optimization problem, and the optimal $d$ can be found by taking derivative of $\ell(d)$ with respect to $d$, since $\ell(d)$ is a convex function of $d.$ To see this,

[TABLE]

Thus $\ell(d)$ is a convex function on $[0,+\infty),$ and the optimal $d_{opt}$ value can be found by setting $\ell^{\prime}(d)=0:$

[TABLE]

This gives an unique optimal value

[TABLE]

which is equivalent to those in (4.11) under the assumption that $m=m(K)<<\min(\log\gamma,K).$ Plugging $d=d_{opt}$ in (D. Parameter Setting in Section 4) back to (4.12) yields the choice of $b_{\gamma}$.

D. Proof of Theorem 5.1

By Theorems 4.1 and 4.2, the false alarm breakdown point of our proposed method $N_{\alpha}$ can be found by finding the smallest $\epsilon$ value such that $I_{0}(\epsilon,\alpha;g_{0})>0$ for some distribution $g_{0},$ where $I_{0}(\epsilon,\alpha;g_{0})$ is defined in (4.1). That is equivalent to

[TABLE]

The remaining proof is based on a careful analysis of $I_{0}(\epsilon,\alpha;g_{0})$ for any arbitrary outlier density function $g_{0}$. For any $h_{0}(x)=(1-\epsilon)f_{0}(x)+\epsilon g_{0}(x)\in\hbar_{0,\epsilon},$ by (4.1), we have

[TABLE]

where $d_{\alpha}(f_{0},f_{1})$ is defined in (5.2) and is the density power divergence between $f_{0}$ and $f_{1}$ proposed by [2]. Here we use the fact that $\int[f_{1}(x)]^{1+\alpha}dx=\int[f_{0}(x)]^{1+\alpha}dx$ when $f_{0}(x)$ and $f_{1}(x)$ come from the same location family.

By the definition of $M(\alpha)$ in (5.1), it is clear from (7.19) that

[TABLE]

Therefore, by (7.18), if both $d_{\alpha}(f_{0},f_{1})$ and $M(\alpha)$ are finite, the false alarm breakdown point of $N_{\alpha}$ should be

[TABLE]

If $d_{\alpha}(f_{0},f_{1})$ is finite but $M(\alpha)=+\infty$, by (7.18) and (7.20), $\epsilon^{*}(N_{\alpha})=0$. If $d_{\alpha}(f_{0},f_{1})=+\infty$ but $M(\alpha)$ is finite, $\epsilon^{*}(N_{\alpha})=1$. If both $d_{\alpha}(f_{0},f_{1})$ and $M(\alpha)$ are $+\infty$ and $\frac{d_{\alpha}(f_{0},f_{1})}{M(\alpha)}=\rho$, by (7.18) and (7.20), we have $\epsilon^{*}(N_{\alpha})=\frac{\rho}{\rho+(1+\alpha)}$ no matter $\rho$ is finite or not. Therefore, for all cases, the false alarm breakdown point of $N_{\alpha}$ have the same expression in (7.21), which completes the proof of Theorem 5.1.