Sequential Change-point Detection for High-dimensional and non-Euclidean Data

In many applications, it is often of practical and scientific interest to detect anomaly events in a streaming sequence of high-dimensional or non-Euclidean observations. We study a non-parametric framework that utilizes nearest neighbor information among the observations to detect changes in an online setting. It can be applied to data in arbitrary dimension and non-Euclidean data as long as a similarity measure on the sample space can be defined. We consider new test statistics under this framework that can detect anomaly events more effectively than the existing test while keeping the false discovery rate controlled at a fixed level. Analytic formulas approximating the average run lengths of the new approaches are derived to make them fast applicable to modern datasets. Simulation studies are provided to support theoretical results. The proposed approach is illustrated with an analysis of the NYC taxi dataset.


Sequential Change-point Detection for High-dimensional and non-Euclidean Data Lynna Chu and Hao Chen
Abstract-In many applications, it is often of practical and scientific interest to detect anomaly events in a streaming sequence of high-dimensional or non-Euclidean observations. We study a non-parametric framework that utilizes nearest neighbor information among the observations to detect changes in an online setting. It can be applied to data in arbitrary dimension and non-Euclidean data as long as a similarity measure on the sample space can be defined. We consider new test statistics under this framework that can detect anomaly events more effectively than the existing test while keeping the false discovery rate controlled at a fixed level. Analytic formulas approximating the average run lengths of the new approaches are derived to make them fast applicable to modern datasets. Simulation studies are provided to support theoretical results. The proposed approach is illustrated with an analysis of the NYC taxi dataset.
Index Terms-Anomaly detection; Online change-point; Streaming data; Graph-based tests; Non-parametric.

I. INTRODUCTION
S EQUENTIAL change-point detection aims to detect abrupt anomalies, observations that deviate from regular behavior, in a streaming sequence of observations as quickly as possible, while controlling the number of false alarms. In many modern applications, the sequence of observations may consist of high-dimensional observations (where the dimensional of each observation is larger than the sample size) or non-Euclidean objects (for example, sequences of networks or images). Examples include fraud detection involving large amounts of customer transactions [1]; disease surveillance or medical monitoring using sequences of images or multiple diagnostic measures [2], [3]; and consumer-based data streams such as wearable and smart device monitoring, internet of things sensors, network security systems, cybersecurity, and other web applications [4]- [6]. In all these examples, the goal is to detect a change (if a change is present) as quickly as possible once it occurs, while also limiting the risk of false discovery.
The sequential change-point setting can be formulated as follows: let the observation at t be denoted as Y t , t = 1, 2, . . . , n, . . .. Here, t could be the time index or some other meaningful indices, Y t could be a vector, image, or network, and n is the index for the observation currently being observed. When there is no change-point, Y t 's are identically distributed The work was supported in part by the NSF awards DMS-1513653 and DMS-1848579. Lynna  from an unknown distribution, denoted as F 0 . If there is a change-point at τ , the observation after τ are from a different (unknown) distribution: where F 0 and F 1 are two different probability measures.
This formulation is very general. F 0 and F 1 are unknown and not specified. We do not impose any constraints on how the change happens here. For example, if the observation is a high-dimensional vector, the change may occur in a subset of (unknown) data streams and the subset may be of size one.
Our aim is to construct a stopping rule (denoted as T ) that will alert us to a change as quickly as possible, once it occurs, but keep the number of false discoveries at a fixed level. Moreover, in keeping with the spirit of the problem formulation, the stopping rule should be able to handle modern data types and should be relatively fast and easy to implement. To be precise, this means we would like to define T such that the detection delay, E τ (T − τ |T > τ ), is small, subject to a fixed average run length (formally E ∞ (T ) ≥ c, where c is a pre-specified large value) without making assumptions on the underlying sequence of observations. Here E τ denotes the expectation under the hypothesis that the true changepoint happens at τ and E ∞ denotes the expectation under the hypothesis of no change.

A. Related Works
When the data is univariate (or scalar), the sequential change-point detection has been studied extensively (see [7] and [8] for a review). For low-dimensional data, likelihood based methods have been explored which require knowledge or parameter estimates of the probability density functions (see for example: [9]- [12]). For high-dimensional data, the available methods are somewhat limited and may impose strict assumptions on the data. As an example, many existing methods have the assumption that the different data streams are independent [13]- [17]. Other works allow for more flexible application; for example kernel-based methods [18] and a modified sliding window algorithm [19]. Computationally efficient methods have also been proposed which combine summary statistics based on geometric entropy minimization (GEM) with the cumulative sum (CUSUM) algorithm [20], [21]. For network data, [22] proposed a sequential approach that acts by embedding each graph into a vector domain, where a conventional multivariate change-point detection procedure can be then applied. In general, for non-parametric methods applicable to high-dimensional and non-Euclidean data, theoretical analysis establishing false discovery control is very difficult to carry out.
Recently, [23] proposed a new non-parametric framework that utilizes nearest neighbor information to detect changes in an online setting. They also provided a general, analytical formula for false discovery control. This method can be applied to data in arbitrary dimensions (with no assumption that the different data streams are independent) and to non-Euclidean data. The author proposed to use the following stopping rule: where n 0 , n 1 , and L are pre-specified values, N 0 is the number of historical observations with no change-point, n > N 0 , and Z L|y (t, n) is a two-sample test statistic that tests whether {Y n−L+1 , . . . , Y t } and {Y t+1 , . . . , Y n } are from the same distribution. We refer to Z L|y (t, n) as the edge-count twosample test based on k-nearest neighbor (k-NN) in the following for simplicity. For more details of this test, please see Section II. The author also provided an analytic formula to compute b Z such that the average run length E ∞ (T Z (b Z )) is controlled at a pre-determined value. Simulation studies show that this method beats likelihood-based methods when the dimension is high.
Despite these nice properties, we find that the edge-count two-sample test on k-NN can have low power for some common types of changes when dimension is moderate to high, causing the stopping rule (1) to behave unexpectedly. To illustrate, consider a simple scenario where data are from a d-dimensional Gaussian distribution and there is a change at τ = 201: with Σ(i, j) = 0.3 |i−j| . We consider two types of changes: • Scenario 1 (only mean differs): ||µ|| 2 = ∆. • Scenario 2 (both mean and variance differ): ||µ|| 2 = ∆ and σ. Tables I presents the performance of T Z (b Z ) for both scenarios based on 1,000 simulation runs. Here k = 5, L = 200, n 0 = 25, and n 1 = 175. The table reports the fraction of trials (out of 1,000) to successfully detect the change within 30 (or 50) observations after the change occurs. The average detection delay (EDD) is estimated as the average elapsed time between when the change occurs (τ = 201) and when T Z (b Z ) detects a change. False alarms are not counted here. In each scenario, the threshold b Z is computed by formulas given in [23] such that E ∞ (T Z (b Z )) = 2000.
In Table I, we see that the method performs worse in Scenario 2 than Scenario 1 under two different comparison criteria: the fraction of trials that can be detected given a fixed time is smaller under Scenario 2 and the average detection delay (EDD) is larger under Scenario 2. However, common sense tells us that the additional change in variance should make the two distributions more different and the change in Scenario 2 easier to detect. This phenomenon is due to the curse-of-dimensionality (a more detailed explanation can be found in Section II-B) and results in the existing method having diminished power to detect general changes.

B. Our contributions
To address the problem of the stopping rule T Z (b Z ), we propose three new stopping rules: The definitions of S L|y (t, n), W L|y (t, n), and M L|y (t, n) are provided in Sections II-C, II-D and II-E, respectively. Under the same setup detailed in Table I, Table II shows that these new stopping rules are more successful in detecting the change quickly after it has occurred. They also have shorter detection delays than T Z (b Z ) under both above scenarios and all have shorter detection delays in Scenario 2 than that in Scenario 1. Further comparisons between the stopping rules can be found in Section IV; these demonstrate that the new stopping rule have improved power and detection delay compared to existing methods over a range of general scenarios. To construct these new stopping rules, we propose new twosample tests on k-NN. Specifically, we extend the generalized (S) /weighted (W ) /max-type (M ) edge-count test defined on an undirected similarity graph [24], [25] to the directed k-NN graph. The generalized edge-count and max-type tests on k-NN are well defined except for a particular construction of a k-NN graph (see Theorem 1). The detailed definitions of these stopping rules are given in Section II.
To make the new stopping rules useful for real-data applications, we provide analytical formulas to compute the stopping thresholds, b S , b W , and b M , such that the average run length (ARL) for each new stopping rule is controlled at a pre-determined value. This involves studying how the k-NN graph updates and obtaining expressions that explicitly characterize these dynamics, which lead to the development of new theoretical treatments. Specifically, for all the stopping rules, more accurate expressions of the graph updates are derived and the techniques used are an improvement over the approach utilized in [23]. We demonstrate that the analytical formulas for the stopping thresholds are reasonable to use and we further improve upon their accuracy for finite sample sizes by implementing a skewness correction technique on the thresholds.
In general, each test statistic has its own niche where it dominates. When interested in general change (for example, both mean and variance change), the stopping rules T S (b S ) and T M (b M ) are recommended. The stopping rule T M (b M ) has an advantage over T S (b S ) in that we can obtain more accurate analytical expression for the ARL for false discovery control. If the change of interest is in mean only, the stopping rule based on T W (b W ) is recommended. See Section IV for a comparison of their performance. In this paper, the types of changes we explore are confined to mean and/or variance change. However, the approach can be used to detect other changes in distribution, such as changes in covariance. For illustration, Table III shows the performance of the stopping rules under covariance change only. The data are again generated from d-dimensional Gaussian with d = 100, k = 5, L = 200, n 0 = 25, n 1 = 175, and there is change at τ = 201: with Σ 0 (i, j) = 0.3 |i−j| and Σ 1 (i, j) = 0.68 |i−j| . We see here again that the new stopping rules based on W , S, and M perform better than the stopping rule based on Z.
These new approaches are implemented in an R package gStream.
The organization of the rest of the paper is as follows. Section II discusses the new stopping rules in details. Section III studies the asymptotic properties of the proposed stopping rules and analytic ways to determine the thresholds. The performance of the new methods are further explored in Section IV and the new methods are illustrated on a real data set in Section V.

II. NEW TESTS
The test statistics in the stopping rules (2) -(4) trace from the offline version of the problem studied in [25] where the statistics were defined on an undirected similarity graph. For online detection, observations keep arriving and the similarity graph updates as a new observation arrives. Therefore one needs to understand the dynamics of the series of similarity graphs. [23] studied the directed nearest neighbor graphs as the dynamics of nearest neighbor graphs can be well understood. In this work, we continue to use the directed nearest neighbor graphs to construct improved statistics. The extension to other types of graphs is saved for future work.

A. Notation
First we define a random variable which indicates whether or not an observation Y i is among the rth nearest neighbor to another observation Y j among the observations in n L . Specifically, for any n > N 0 and i, j ∈ n L ∆ = {n − L + 1, . . . , n}, we let A (r) is the indicator function that takes value 1 if the event is true and 0 otherwise. In terms of graph construction, each observation points to its k nearest neighbors. For example, if A (r) to be the indicator function that Y j is one of the first k NNs of Y i among the observations in n L . We use y i 's to denote the realizations of Y i 's and let a + n L ,ij = k r=1 a (r) n L ,ij with a (r) n L ,ij = I(y j is the rth NN of y i among y n−L+1 , . . . , y n ). For any n, each t ∈ {n − L + 1, . . . , n} divides the data sequence into two groups: one group being the observations before t : {Y n−L+1 , . . . , Y t } (Group 1) and the other group being the observations after t : {Y t+1 , . . . , Y n } (Group 2). Define, b 0,ij (t, n L ) =I(n − L + 1 ≤ i ≤ t and t < j ≤ n or t < i ≤ n and n − L + 1 ≤ j ≤ t), b 1,ij (t, n L ) =I(n − L + 1 ≤ i ≤ t and n − L + 1 ≤ j ≤ t), b 2,ij (t, n L ) =I(t < i ≤ n and t < j ≤ n).
Then b 0,ij is the indicator function that Y i and Y j belong to different groups, b 1,ij is the indicator function that Y i and Y j both belong to Group 1, and b 2,ij is the indicator function that Y i and Y j both belong to Group 2.
We define our test statistics as follows:  B 0,ij (t, n L ), B 1,ij (t, n L ), and B 2,ij (t, n L ) are the random variable versions of b 0,ij (t, n L ), b 1,ij (t, n L ), and b 2,ij (t, n L ) such that the distribution of these random variables is defined to be the permutation distribution. The permutation distribution is the distribution induced by all L! possible permutations of observation indices among the observations in n L . The null hypothesis and independence assumption imply that the observation indices are exchangeable, and therefore under the permutation distribution, the true (unknown) distribution of the observations remains unchanged when the null hypothesis is true.
It is clear that R 0,L (t, n) is twice the number of edges in the k-NN graph connecting observations before t and after t, R 1,L (t, n) is twice the number of edges connecting observations prior to t, and R 2,L (t, n) is twice the number of edges that connect observations after t. The notation of the graph-based test quantities emphasizes their dependency on the graph which is constructed on the L most recent observations, with the most recent observation indexed at n. Figure 1 illustrates these test statistics constructed for different times t and n. In the top row of Figure 1, n = 30 and we construct the graph on the L = 20 most recent observations: y 11 , y 12 , . . . , y 30 . The edge-counts are calculated for different values of t ∈ {11, . . . , 30}. R 0,L (t, n) is twice the number of directed black edges, R 1,L (t, n) is twice the number of directed red edges, and R 2,L (t, n) is twice the number of directed blue edges. It is clear for a fixed n, the graph does not change: each graph in the first row of Figure 1 is the same. But as new observations continually arrive, the graph itself will update. For example when n = 40, we now construct the graph on observations y 21 , y 22 , . . . , y 40 (see the second row of Figure 1). For n = 40, the edge-counts are then calculated for different values of t ∈ {21, . . . , 40}.

B. Limitations of the method based on the edge-count test (Z)
The stopping rule T Z (b Z ) (1) is based on the edge-count two-sample test statistic Z L|y (t, n) [23]. To obtain Z L|y (t, n), first a directed k-NN graph is constructed based on a similarity measure (for example, Euclidean distance). Then the number of edges in the k-NN graph that connect observations before t (Y n−L+1 , . . . , Y t ) and after t (Y t+1 , . . . , Y n ) is counted (we refer to this as the between-sample edge-count). Please see Figure 1 for an illustration of how R 0, . Analytical expressions of E(R 0,L (t, n)) and Var(R 0,L (t, n)) can be obtained under the permutation distribution, which is defined as all L! possible rearrangements of the observation indices among n L , and are omitted here for brevity. A relatively low between-sample edge-count (or a large Z L|y (t, n)) indicates the observations before and after t are less mixed and this is evidence against the null hypothesis of no change. The intuition is that observations before t tend to find their nearest neighbor among other observations before t (and similarly for observations after t), which implies a distributional difference between the two groups of observations.
The rationale of a relatively small between sample edgecount (R 0,L (t, n))) holds particularly well under the scenario of mean change only and/or low-dimensional data. However, when the dimension of the data is moderate/high and the change in distribution is not only in mean, for example an additional variance change is present, this rationale breaks down. This is because in the high-dimensional setting observations that are similar (i.e. from the same distribution) are not necessarily close in distance. Consider two distributions that differ in both mean and variance and the dimension of the data is d = 100. The observations would be separated into two layers: observations with the smaller variance in the inner layer and observations with the larger variance in the outer layer. Since the volume of a d-dimensional space increases exponentially in d, for practical sample sizes it is not uncommon for observations in the outer layer to find themselves to be closer to observations in the inner layer compared to other observations in the outer layer. In this scenario the between-sample edge-count would be relatively large under the alternative, rendering Z L|y (t, n) ineffective. The new two-sample test statistics S L|y (t, n) and M L|y (t, n) address this curse-of-dimensionality problem.
Moreover, even for mean change only the stopping rule based on Z L|y (t, n) can still suffer from increased detection delay. This limitation of Z L|y (t, n) is due to a variance boosting problem under unequal sample sizes in the twosample test setting. For more details on the variance boosting issue in the two-sample test setting see [26]. In the sequential setting, this means that using Z L|y (t, n) can lead to increased detection delay since it may only be able to detect the change when it is near the middle of the sequence. To resolve this, we propose the weighted edge-count two-sample test W L|y (t, n).

C. The method based on the generalized edge-count test (S)
The generalized edge-count two-sample test at t ∈ {n−L+ 1, . . . , n} under k-NN can be defined as and Σ y (t, n) = Var((R 1,L (t, n), R 2,L (t, n)) T |y) such that E(·) and Var(·) denote the expectation and variance taken under the permutation distribution.
If a change-point τ > N 0 occurs in the sequence, we would expect S L|y (t, n) to be large when n > τ and t close to τ . The test statistic is defined in this way so that either direction of deviations in the number of within-sample edges from its null expectation would contribute to the test statistic. For example, under the location alternatives, we would expect both R 1,L (t, n) and R 2,L (t, n) to be larger than their null expectations, which would lead to a large S L|y (t, n). Under the scale alternatives, the group with the smaller variance would have a within-edge count larger than its null expectation and the group with a larger variance would have a within-edge count smaller than its null expectation (due to the curse-ofdimensionality), which would also lead to a large S L|y (t, n). Therefore, this test is powerful for both location and scale alternatives.
The following theorem ensures that the statistic is well defined except in very rare scenarios. These scenarios can be checked by calculating the node in-degree of each observation in the k-NN graph. A node's in-degree is the number of other observations that find that node to be its nearest neighbors. Theorem 1. For L ≥ 5, the generalized edge-count twosample test statistic under k-NN is well defined except for when all nodes have an in-degree of exactly k, i.e. d i = k ∀i, where d i = j∈n L a + ji . Based on Theorem 1, as long as each node in the graph does not have an in-degree of k and the graph is constructed on at least 5 observations, then Σ y (t, n) is invertible and S L|y (t, n) is well-defined. A proof of the above theorem is in the supplement.
This leads us to the stopping rule based on the generalized edge-count test under k-NN (2): D. The method based on the weighted edge-count test (W ) Following the same notations in Section II-C, for each t ∈ {n − L + 1, . . . , n}, the weighted edge-count two-sample test statistic under k-NN can be defined as Since it is more difficult for the sample with a smaller sample size to form an edge within the same sample, R 1,L (t, n) and R 2,L (t, n) are weighted by the inverse of their corresponding sample sizes. The test statistic defined in this way resolves the variance boosting problem described in Section II-B.
Relatively large values of R w,L (t, n) are evidence against the null hypothesis of no change. Let Var(R w,L (t, n)|y) .
Under the permutation distribution, analytical formulas for E(R w,L (t, n)) and var(R w,L (t, n)|y) can be calculated based on E(R 1,L (t, n)), E(R 2,L (t, n)), and Σ y (t, n) provided in Section II-C: The stopping rule based on the weighted edge-count test under k-NN is (3): E. The method based on the max-type edge-count test (M ) We can define the max-type test statistic under k-NN based on the following lemma: Lemma 1. The generalized edge-count two-sample test under k-NN can be expressed as The proof of this lemma is in the supplement. The analytical expressions for the expectation and variance of R diff,L (t, n) under the permutation null are: The variance of R diff,L (t, n)|y is well-defined as long as d i = k ∀ i, in other words as long as each node does not have an in-degree of k.
From the above lemma, S L|y (t, n) is the sum of squares of two uncorrelated quantities (these two quantities are further asymptotically independent; details given in Section III). Here, W L|y (t, n) is sensitive to location changes: when the change is in mean, W L|y (t, n) tends to be large. On the other hand, D L|y (t, n) is more sensitive to scale changes: when the change is in variance, |D L|y (t, n)| tends to be large. The sign of D L|y (t, n) depends on whether the distribution after the change has a larger spread or not. This leads to the following max-type edge-count two-sample test statistic under k-NN: M L|y (t, n) = max(|D L|y (t, n)|, W L|y (t, n)).
When there is a change in location and/or scale, depending on the signal of interest, it is useful to consider an extended version of the max-type edge-count two-sample test: where ξ ≥ 0. Different choices of ξ lead to different focuses of the alternatives. For example, if we are more interested in locational changes, we could choose a large ξ. On the other hand, setting ξ to be small would favor detecting scale changes. When ξ = 1, the test reduces to the plain max-type edge-count test. For more detailed discussion on how to select ξ, see Supplement H in [25] (under the offline change-point detection setting, but similar arguments apply to the online setting). The stopping rule based on the max-type edge-count test under k-NN is as follows : where n 1 = n − n 1 and n 0 = n − n 0 . This reduces to (4) when ξ = 1.

III. AVERAGE RUN LENGTH
Given the new stopping rules presented in Section II, we would like to determine the thresholds b S , b W , and b M ξ in an analytic way such that the false discovery rate is controlled at a pre-specified value. A common way to measure the false discovery rate under the online change-point detection is the average run length, i.e., the expected time to stop when there is no change-point, which we denote as In the comparisons in Section I-B (Table II), the thresholds were chosen such that the average run lengths are 2,000 based on simulation runs. This is doable when the underlying distribution of the sequence is known. However, in many applications, the distribution of the sequence is unknown. Furthermore, since new observations keep arriving, resampling based methods, such as permutation and bootstrap, are not appropriate here and even if they were, directly resampling could be very time consuming. Therefore, to make the method fast applicable, we seek to derive analytical expressions for the average run lengths. Given the non-parametric nature of the proposed method, we would not be able to get exact analytic formulas for the average run lengths under finite L. In the following, we first approach the problem asymptotically (Section III-A), and then make adjustments for finite samples (Section III-B).

A. Asymptotic results
To derive analytical expressions for the average run length, we must study the asymptotic distribution of the stopping rules. Since these stopping rules are composed of the random fields {S L|y (t, n L )}, {W L|y (t, n)}, and {M ξ,L|y (t, n)}, we study their asymptotic properties. To obtain the limiting distribution of these random fields, we only need to focus on {D L|y (t, n)} and {W L|y (t, n)}. We show that the limiting distribution of the random fields converge to independent twodimensional Gaussian random fields (Theorem 2). This proof utilizes Stein's method [27]. To fully specify the Gaussian processes, we must derive the covariance functions of the new processes (Theorem 3). Since the k-NN graph updates each time new observations arrive, we must study the dynamics of the k-NN series for the new test statistics. While [23] laid a framework for graph-based sequential detection, the techniques developed in [23] were not directly applicable to the stopping rules proposed in this manuscript and needed to be adapted. Integration of techniques from [25] were useful in this development. However, since the new stopping rules can cover more types of change, we found that a direct extension of these previous works were not sufficient in developing accurate analytical formulas approximating ARL. To push forward the theory, we carefully studied the dynamics of directed k-NN graphs and considered additional ways the graph can update. A more in-depth understanding of the graph updates led us to derive analytical expressions for these updates and incorporate them into our analytical formulas for the thresholds. A comparison of the improvement over a direct extension of previous works can be found in the supplement. Putting together Theorems 2 and 3 with results from [28], allow us to obtain analytical expressions for the average run lengths of the new stopping rules.
In the following, we provide a sketch of the key steps to obtain each result. We defer readers to the supplement for a more technical treatment. The subsequent results are derived under the following condition: Condition 1. There is a positive constant C, 1 ≤ C < ∞, depending only on k, such that In k-NN, each observation points to its first k NNs, so the out-degree of each observation is k. However the in-degree of each observation can vary. This condition states that the indegree of each observation is bounded. It is satisfied almost surely for multivariate data [29], [30]. Let Var(R diff,L (t, n)) , Var(R w,L (t, n)) , which replaces the conditional variances of the test statistics with the unconditional variances. See the supplement for more details.
Our proof that the limiting distributions converge to independent two-dimensional random fields depends on utilizing Stein's method [27]. The general idea of Stein's method is to show for a random variable W and a standard normal random variable Z, that for some family of functions h ∈ Lip(1), the following bound holds: where δ depends on the structure of W . The specific form of Stein's method used here requires A + n L ,ij to be locally dependent. However, even for different i, j, l, r, A + n L ,ij and A + n L ,lr are dependent due to the constraint that j∈n L A + n L ,ij = k for all i ∈ n L . Following [23], we relax these dependencies by considering a similar set of Bernoulli random variables {Ã + n L ,ij } i,j∈n L . We keep the following probabilities unchanged: P(Ã + n L ,ij = 1) = P(A + n L ,ij = 1), P(Ã + n L ,ij = 1,Ã + n L ,ji = 1) = P(A + n L ,ij = 1, A + n L ,ji = 1), P(Ã + n L ,ji = 1,Ã + n L ,li = 1) = P(A + n L ,ji = 1, A + n L ,li = 1), but relax the other dependencies such thatÃ + n L ,ij is independent of {Ã + n L ,il ,Ã + n L ,li } l =j , andÃ + n L ,ij andÃ + n L ,lr are independent when i, j, l, r are all different. ThenÃ + n L ,ij are only locally dependent and can be analyzed through the Stein's method [27].
We are now ready to present the main results.
converge to independent two-dimensional Gaussian random fields, which we denote as Here [x] denotes the largest integer smaller than or equal to x for any real number x.
The detailed proof for Theorem 2 is in the supplement. Based on Theorem 2, we can approximate E ∞ (T S (b S )), E ∞ (T W (b W )), and E ∞ (T M ξ (b M ξ )) by examining the asymptotic behavior of our stopping rules: where n 1 = n − n 1 and n 0 = n − n 0 . Our approximations involve the function ν(x) defined as This function is closely related to the Laplace transform of the overshoot over the boundary of a random walk. A simple approximation given in [31] is sufficient for numerical purposes: where Φ(·) is the cumulative distribution function of the standard normal distribution and φ(·) the density function of the standard normal distribution. According to [28], We now have approximations of the average run lengths for the new stopping rules. The only remaining unspecified quantities are the directional partial derivatives of the covariance functions of the Gaussian random fields. Their analytical expressions are derived in the following theorem.
Here, p k,∞ is the limiting expected number of mutual NNs a node has in k-NN, q k,∞ is the limiting expected number of nodes that share a NN with another node in k-NN, p ∞ (r, s) is the limiting expected number of mutual NNs shared between the rth and sth NNs, and similarly q ∞ (r, s) is the limiting expected number of nodes shared between the rth and sth NNs. Explicitly, p k,∞ = k r=1 k s=1 p ∞ (r, s), and q k, To derive these partial derivatives, we studied the dynamics of the k-NN series as new observations are added. It turns out that a few key quantities are enough to characterize the dynamics in the asymptotic domain. The proof of this theorem is in the supplement.

B. Finite L
We now consider the practical scenario where L is finite. Based on results in Section III-A, E ∞ (T S (b S )), E ∞ (T W (b W )), and E ∞ (T M ξ (b M )) can be approximated by with H L () and G L () are finite sample versions of H() and G(), respectively. In practice, when L is finite we use g W,2 (L, x) and g D,2 (L, x) in place of g W,2 (x) and g D,2 (x) in the above formulas, respectively, where Here, p k,L , p k+1,∞ and can be estimated in a data-driven way.

C. Skewness correction
Analytical approximations provided in Section III-B become less precise for finite L when n 0 is relatively small. This is mainly because the convergence of W L (t, n) and D L (t, n) to normal is slow if (n − t)/L is close to 0 or 1. This problem becomes more severe when dimension is high. To improve upon the analytic approximations for finite sample sizes, we perform skewness correction. We adopt a skewness correction approach discussed in [32] that does the correction up to different extents based on the amount of skewness at each value of t. In particular, we provide better approximations to the marginal probabilities P(W (u − x, w) ∈ b + du) and P(D (u − x, w) ∈ b + du). Following the method based on cumulant-generating functions and change of measure (details refer to [32]), we can approximate the marginal probability by where θ b is chosen such thatψ(θ b ) = b. By a third Taylor approximation, we get where γ L (t, n) := E P (Z L (t, n) 3 ) and explicit expressions are derived using combinatorial analysis. Skewness corrected thresholds are only obtained for W L|y (t, n) and D L|y (t, n), but not S L|y (t, n). This is because for S L|y (t, n) the integrand can easily be non-finite and the approach depends heavily on extrapolation. Therefore, the stopping rule based on M L|y (t, n) is often recommended because it can detect general changes but we can obtain more accurate stopping thresholds.

D. Checking accuracy of analytic formulas for the average run lengths
Here, we check the accuracy of the analytic formulas for the average run lengths. For all three new tests, we have analytic formulas based on asymptotic results (19), (20) and (21), and for the tests based on the weighted/max-type edge-count tests, we have analytic formulas after skewness correction, provided in the supplement. We compare the empirical ARL obtained from these analytic formulas (with skewness correction with applicable) to those obtained from 1, 000 Monte Carlo simulations. The analytical thresholds are obtained so that the average run length is 2, 000. We generated data from three different settings: multivariate normal with d = 10 (denoted by C1), multivariate t 5 with d = 100 (denoted by C2), and multivariate log-normal with d = 1000 (denoted by C3).
Results for different choices of n 0 are shown in Tables IV -VIII. We set n 1 = L − n 0 and k = 5. The asymptotic analytic results are denoted by 'A1' and the skewness corrected approximations are denoted by 'A2'. We see in general that the empirical ARLs obtained from the asymptotic approximations ('A1') are not very close, illustrating the need for skewness correction here. After skewness correction, the empirical ARLs are much closer to 2, 000. It is clear that the accuracy of the skewness corrected approximations depends on n 0 : in general, when n 0 = 40, the skewness corrected approximations do well across most dimensions.
We also investigate how the analytical threshold approximations perform compare to the Monte Carlo thresholds in the presence of change. In Tables V, VII, and IX we report the power, defined as the fraction of trials able to detect the change within 50 observations after the change occurs, and the average detection delay (reported in parenthesis) obtained using analytical formulas ('A1'), analytical formulas with skewness correction ('A2') and Monte Carlo thresholds ('MC') . The thresholds are obtained so that the ARL is set to be 2, 000. The amount of signal used correspond to the power setup described for Tables X, XII, and XI in Section IV, respectively. It is clear that in the presence of change, the analytical formulas (with skewness correction) and Monte Carlo thresholds all lead to similar results in power and detection delay.

IV. POWER ASSESSMENT A. Multivariate data
To examine the performance of the three new test statistics, we compare them to the existing approach in [23] (max t Z L|y (t, n)) and two parametric likelihood-based approaches: Hotelling's T 2 test when there is change in mean and the generalized likelihood ratio test when there is variance change (both these two-sample tests are adapted to the scan statistic setting). The simulation setup is as follows: there are N 0 = 200 historical observations and a change occurs at t = 400 (200 new observations after the start of the test). The observations are independent and follow a d-dimensional distribution. When there is a change in mean, the observations are shifted from 0 by amount ∆ in Euclidean distance. When the covariance matrix changes, to make the change less significant, only the first d/5 of the diagonal elements change with a multiple of σ, and the rest are unchanged. The amount of change is chosen so that the tests have moderate power to be comparable. For fair comparison, we use Monte Carlo simulations to determine the threshold for each of the test so that their average run lengths are all 2,000. Power is reported as the fraction of trials for which the change-point is detected within 50 observations after the change occurred. In the following, we use 'HT' to refer to the scan statistic over the Hotelling's T 2 statistic and use 'GLR' to refer to the scan statistic over the generalized likelihood ratio statistic. For d > 100, in order for HT and GLR to be applicable in higher dimensions, we treat each data stream as if it is independent so that the covariance matrix's inverse and determinant are well-defined.
When there is both mean and variance change, Table X shows the results under the Gaussian setting. When d ≤ 500, the GLR dominates in power. However, when the dimension increases, GLR is no longer able to retain competitive power compared to S and M . On the other hand, this setting is not well-suited for W which is meant to capture mean change only and its performance is the worst here.
To consider other distributions, we also compared the tests for multivariate log-normal data and multivariate t 5 data. The results for the log-normal data are shown in Table XI. Here there is a change in mean parameter only and ∆ is chosen such that the location change dominates. In this setting, W 's performance dominates. We see that all the new tests outperform Z and the parametric tests for d > 10.
The result for the multivariate t 5 data are shown in Table  XII. When there is a change in both mean and variance, Z is unable to outperform the new test statistics. Among the new test statistics, their performance depends on whether the mean or variance signal dominates. When the mean change is stronger (for example, when d = 100), W performs comparably well. However, when the variance change is stronger (for example, when d ≥ 1000), M and S dominate.
Based on the results of these tables, we see that the new graph-based methods perform well under various scenarios and have improved detection delay over the existing method in [23]. In general, if one is certain that the change is locational, the test based on W is recommended; while for more general changes, the tests based on S and M are recommended.

B. Network data
Non-Euclidean data, also referred to as object data, is simply data that does not lie in the Euclidean space. Examples include networks, images, shapes, and trees. Fundamental statistical tools involving vector space analysis are no longer applicable for object data, making it challenging to directly analyze these data types. To demonstrate the new test statistics' power on non-Euclidean data, we evaluate the proposed methods on a sequence of networks. We generate random networks using the configuration model, which is a specific method that allows us to generate random networks with a specified degree sequence (please see [33], [34] for an overview). For every node with degree k i , we create k half-edges (referred to as 'stubs'). The network is created by iteratively selecting two stubs uniformly at random and connecting them to form an edge. This is done until no stubs remain and will result in a network with a predefined degree sequence. For the graph-based results, we show the results for graphs constructed using two different similarity measures. Specifically, for a network at time t, we encode the network using an adjacency matrix A t with 1 for element (i, j) if node i and j are connected, and 0 otherwise. The similarity measures are: 1) Similarity 1: ||A t − A s || F , 2) Similarity 2: Under this setting, we compare the graph-based approach to another method designed specifically to detect changes in a stream of networks/graphs [22]. The approach proposed in [22] is quite general: it does not make explicit assumptions on the network/graph and can be applied to a stream of graphs of varying sizes. Their approach consists of two steps: (1) each graph G t is mapped to a vector y t through a prototype-based embedding, and (2) a change is then detected in the stream of vectors y t using any conventional multivariate changedetection procedure. Specifically, embedding is carried out by assessing the dissimilarity between graphs and a selection of prototypes. In the paper, they rely on the graph edit distance (GED), which count and weights the edit operations that are needed in order to make two input graphs equal and is applicable to graphs where the nodes are unidentified. To carry out their approach, we follow the implementation proposed in [22]; specific details can be found in Section IV of [22]. Since computing the GED is quite computationally expensive, we allow each network to have only 6 nodes. We set a change at τ = 101. Before and after the change, the total node degree remains the same, but the sequence of node degree changes. Before the change, each network is constructed such that nodes 1 and 4 have node degree of 1 and the remaining nodes have node degree 2. After the change, nodes 3 and 6 have node degree of 1 and the remaining nodes have node degree 2. Observe that in this setting Similarity 1 and 2 are equivalent and so we only report results for Similarity 1. The first 50 observations are treated as training observations for [22]. We set L = 50 for the graph-based approach. For fair comparison with the method in [22], we set the window size to be 50 and we implement their method using both GED and the Forbenius norm (Similarity 1) as the similarity measure between networks. The ARL is set to be 2,000. Power is defined as the number of trials (out of 100) where the change is detected within 50 observations. Among those trials where the change is detected, we also report the expected detection delay (EDD) and its standard deviation when applicable. We note that the approach in [22] depends on the choice of embedding and its hyper-parameters, which should be chosen carefully. A thorough discussion is provided in [22]. We carry out simulations under a variety of settings for the number of dimensions to embed (d) and the number of prototypes to select (nproto). Table XIII show the power performance of the graph-based methods compared to the approach in [22]. We can see that all of the graph-based statistics do better in terms of power and detection delay compared to [22]. For further exploration of the graph-based approach, we generate a sequence of networks such that each network has 20 nodes with a pre-specified degree sequence (see below for details). A change in the sequence of networks happens at τ = 101. The first 50 networks in the sequence are treated as historical observations. Power is reported as the fraction of trials for which the change-point is detected within 50 observations after the change occurred. Among those trials where the change is detected within the first 50 observations, the expected detection delay and its standard deviation is also reported (EDD).
The networks are generated under two different settings: 1) A fixed degree change in the network: before the change, half of the nodes have out-degree and in-degree 1 and half the nodes have out-degree and in-degree 3. After the change, 5 nodes have out-degree and in-degree k 1 and 5 nodes have out-degree and in-degree k 2 . The remaining nodes remain unchanged. 2) A random degree change in the network: before the change, half of the nodes have out-degree and in-degree 1 and half the nodes have out-degree and in-degree 3. After the change, half of the nodes have out-degree and indegree randomly selected from k 3 to k 4 and the remaining half of the nodes have out-degree and in-degree 3.
Tables XIV and XV report the power and expected detection delay of the graph-based test statistics for similarity measure 1 and 2, respectively. For similarity measure 1, the performance of S and M dominate in almost all settings, with the exception of when k 1 = 2, k 2 = 2 ; in general S and M do well with respect to both power and expected detection delay. For similarity measure 2, the performance of W and Z improve substantially, while the performance of S and M remain stable.

V. A REAL DATA APPLICATION
We compare the new approaches to the method in [23] using the yellow taxi trip records data. The data set is publicly available on the NYC Taxi & Limousine Commission (TLC) website. It provides information on the taxi pickup and dropoff date/times, longitude and latitude coordinates of pickup and drop-off locations, trip distances, fares, rate types, payments types, and driver-reported passenger counts.
Based on this data set, a natural question to ask is: Can we detect a change in traffic patterns during peak travel seasons? Here, we focus on those trips that began at John F. Kennedy International Airport and we look at two different time periods: the months of June through August and November through December in 2015. The dataset has been completely collected at the time of analysis. However, we treat it as if the data were being observed in order to illustrate how the proposed method works. For simplicity, the boundary of JFK airport was set to be 40.63 to 40.66 latitude and −73.80 to −73.77 longitude. For those trips that began with a pickup at JFK, we extract information on their longitude and latitude drop-off coordinates. Using longitude/latitude coordinates, we create a 30 by 30 grid of New York City and count the number of taxi drop-offs that fall within each cell, where each cell represents a longitude, latitude coordinate range. Then for each day, we have a 30 by 30 matrix such that each element represents the number of taxi drop-offs in each location. We apply the new approaches, as well as Z, to detect changes in the months of June through August 2015. We use data from the month of May as historical data. Applying the offline change-point detection method in [32] and [25] on the observations in May, we find there is no change-point in the first 30 days, so we set L = 30, n 0 = 5, and n 1 = L − n 0 . We denote A i to be the 30 by 30 matrix on day i and v i to be the vector form of A i , which is now 900 by 1. The L 1 norm is used to construct the k-NN graph representing similarity between days. Here, the new test statistics (W , M , and S) all report a stopping time of July 3 and July 4 whereas Z is unable to detect any anomaly event (Table XVI). The change-point triggering these stopping times is estimated to be June 29. To perform a sanity check, we plot a heatmap of the L 1 distance matrix used to the construct the k-NN graph (see Figure 2, left panel). Based on the heatmap, we can see there is a clear signal happening around Day 30, which corresponds with the results from the new test statistics.
To detect changes in November and December 2015, we use data from the months of September and October 2015 as historical data. Applying the offline change-point detection method in [32] and [25] on the observations in September and October, we find there is no change-point in the first 50 days. Therefore, we treat the first 50 observations from Sept. 1 -Oct. 20 as historical observations and we begin the test at Oct. 21. We set L = 50, n 0 = 8, and n 1 = L − n 0 . The stopping times based on the new test statistics report back dates that seem to be quite reasonable (see Table XVII). We see that multiple stopping times are caused by the same anomaly event. When the signal is large enough, the new test statistics (W , M , and S) and Z perform similarly: all are able to detect a change in travel pattern close to Thanksgiving and the Christmas holidays. Again to check our results, we plot a heatmap of the L 1 distance matrix used to the construct the k-NN graph. We can see that there is a clear signal starting roughly around Day 30 and again around Day 60 and Day 67, which matches the results reported from the test statistics. In comparison with the heatmap from the months of June through August, the signal from the summer months is much weaker and in that case Z is unable to detect any anomaly event.

VI. CONCLUSION
We propose new graph-based test statistics under k-NN for detecting change-points sequentially as data are generated. We study the asymptotic properties of the stopping rules based on the new test statistics, and derived the analytic formulas to approximate the average run lengths of the new stopping rules. To accommodate finite samples, skewness corrected approximations were also derived for the weighted and maxtype edge-count statistic under k-NN. The skewness-corrected versions give much more accurate approximations to the average run lengths and can be used reliably in practice. The performance of the proposed test statistics are examined under various common scenarios. Simulation studies reveal that the new test statistics have shorter detection delays for a wider range of alternatives and exhibit power gains for scale change when compared to parametric tests and the test statistic proposed in [23]. Specifically, simulation results show that the weighted-edge count statistic (W ) is useful at quickly detecting mean changes. When a change in variance is also of interest, the generalized edge-count statistic (S) and maxtype edge-count statistic (M ) are more effective in detecting changes and obtain faster detection. Together with the fact that skewness corrected average run length approximations can be obtained for the max-type edge-count statistic, the stopping rule T M ξ is recommended for sequential detection of general changes.