Change Point Detection via Multivariate Singular Spectrum Analysis
Year: 2021

Abstract:
The objective of change point detection (CPD) is to detect significant and abrupt changes in the dynamics of the underlying system of interest through multivariate time series observations. In this work, we develop and analyze an algorithm for CPD that is inspired by a variant of the classical singular spectrum analysis (SSA) approach for time series by combining it with the classical cumulative sum (CUSUM) statistic from sequential hypothesis testing. In particular, we model the underlying dynamics of multivariate time series observations through the spatio-temporal model introduced recently in the multivariate SSA (mSSA) literature. The change point in such a setting corresponds to a change in the underlying spatio-temporal model. As the primary contributions of this work, we develop an algorithm based on CUSUM-statistic to detect such change points in an online fashion. We extend the analysis of CUSUM statistics, traditionally done for the setting of independent observations, to the dependent setting of (multivariate) time series under the spatio-temporal model. Specifically, for a given parameter h>0, our method achieves the following desirable trade-off: when a change happens, it detects it within O(h) time delay on average, while in the absence of change, it does not declare false detection for at least exp(⌦(h)) time length on average. We conduct empirical experiments using benchmark and synthetic datasets. We find that the proposed method performs competitively or outperforms the state-of-the-art change point detection methods across datasets.

heading: 1 Introduction
text: The task of change point detection (CPD) is concerned with detecting significant changes in the temporal evolution of a system from noisy observations. This task has attracted considerable attention in the statistics and machine learning communities as it has a broad range of applications including quality control [5], climate research [35], and speech recognition [36]. This interest resulted in many methods and algorithmswe refer to [3] for a thorough review.
Despite the relatively simple setup of the CPD problem, different variations of the problem have been studied and analyzed. We contrast two general settings characterized by the assumptions on the data generating process.The first is the classical view that aims to detect distribution changes from independent and identically distributed (i.i.d.) observations, cf. [34] and more recently [7, 23]. The second deals with the more intricate time series dynamics by assuming the data to be dependent and following a specific temporal structure. Therein, the aim is to detect changes in this temporal structure of the observations. Some of the time series models considered include auto-regressive models [4, 46] and state-space models [20].
An important model comes from the singular spectrum analysis (SSA) literature, where the data is assumed to follow an underlying function associated with a low-dimensional subspace. This time series model has inspired algorithms for different inference tasks in time series, including imputation, forecasting, and CPD [13]. The key idea behind using SSA for CPD is that the subspace spanned by the columns of the
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Hankel matrix constructed from an initial subsequence of the time series is approximately equivalent to the one spanned by the later subsequences of the time series. The theoretical analysis of this algorithm focused primarily on the probability of error in the asymptotic setting where the signal and noise are assumed to be separable. However, previous work on SSA for CPD does not provide finite-sample analysis on any metric on interest in CPD. Further, its extension to the case of multivariate time series (multiple SSA or mSSA) lacks theoretical understanding, despite being empirically effective [11]. In this work, we consider a variant of mSSA proposed and analyzed in [1] for the tasks of imputation and forecasting. We extend the use of this mSSA variant to devise an algorithm for CPD.
Another important distinction between CPD methods lies in whether the method is online versus offline. In the offline setting, the dataset is fixed and the detection is made retrospectively based on the knowledge of the entire dataset. In such settings, only the accuracy of the detection is of concern. In contrast, in the online setting, the data arrives sequentially and a decision must be made as quickly as possible. In addition to the detection accuracy, the delay in making the detection (i.e. the number of data points observed after the change until it is detected) is another important metric in the online setting.
A common overarching framework for the online CPD problem extends from sequential hypothesis testing, where the problem is formulated as a binary hypothesis test that is applied sequentially as new data comes in [40]. This framework provides a way to evaluate the performance of online CPD methods in terms of the trade-off between false alarm rate and detection delay. The classical setting of this framework uses parametric models and forms likelihood ratio-based test statistics. In more recent work, this framework has been utilized with different test statistics. For example, tests using statistics based on the kernel Fisher discriminant ratio [16], maximum mean discrepancy [25], and subspace distance [19] were proposed. However, rigorously characterizing the trade-off between false alarm rate and detection delay for time series is absent in the literature.
In this work, we are motivated to fill the gap in literature by providing an online method for CPD for the setting of multivariate time series data (not i.i.d.) to detect change points with provable performance guarantees and desirable empirical performance.
Summary of contributions. Below, we summarize the main contributions of this work. 1. We introduce an expressive spatio-temporal model as described in Section 2.1 to model the dynamics
and posit the question of detection as a change in these dynamics. 2. We develop a variant of the SSA algorithm for online change point detection in multivariate time series.
In particular, we propose an algorithm that utilizes the low-dimensional structure of the time series to construct a cumulative sum (CUSUM) statistic (a la [34]) based on subspace distance. We refer to Section 3 for a full description of the proposed algorithm. 3. We analyze the performance of this algorithm in terms of the average running length (ARL), a common metric used in CUSUM tests to measure both the delay in detection and the average running time until a false alarm. Our main result shows that in the case of no change, the expected running length is exponential in the selected threshold (Theorem 4.1). While in the presence of change, the expected detection delay is linear in the threshold (Theorem 4.2). The results are consistent with the behavior of the CUSUM tests in the i.i.d. case. 4. We conduct comparative empirical experiments against state-of-the-art change point detection methods using benchmark datasets. In Table 1, we provide a summary that suggests that our method is at least as good, if not better than, state-of-art methods (SOTA) across datasets.
Related work. Due to space constraints, we provide an overview of related topics in Appendix A.

heading: 2 Problem Setup
text: We consider the discrete time setting where for each time index t2 [T ] :={1,...T}, we observe a multivariate time series X(t) := [X1(t),...,XN(t)]2RN . For each n2 [N ], we denote the n-th latent time series by
fn :Z+!R, such that at each timestep t2 [T ], the observations take the formXn(t)=fn(t)+en(t),where en(t) is the per-step noise modeled as a mean-zero i.i.d. random variable. We define the latent multivariate time series as f :Z+!RN such that f(t):=[f1(t),...,fN(t)]. That is, the observation X(t) takes the form
X(t)=f(t)+e(t), (1)
where e(t):=[e1(t),...,en(t)]. We note that although the additive noise en(t) is i.i.d., the latent time series f(t) can have a complex dependence structure across both t and n, see Section 2.1.
Our objective is to detect changes in the latent time series f(t) as quickly as possible. Precisely, our goal is to detect a change point ⌧ such that
X(t)= ⇢ f0(t)+e(t) if t<⌧ f1(t)+e(t) if t ⌧.
where f0,f1 :Z+ !RN and f0(t) 6= f1(t). In this setting, we treat the change point ⌧ as an unknown deterministic quantity. We frame this CPD problem as a sequential hypothesis test where at each time t one of two hypotheses is accepted
H0 : IE[X(t)]=f0(t) (2) H1 : IE[X(t)]=f1(t).
Let H(t):Z+!{0,1} be the indicator random variable that indicates whether or not the null hypothesis is rejected at time t. We then estimate the change point ⌧ as
⌧̂=inf{t|H(t)=1}. (3)
Ideally, we wish to have ⌧̂ ⌧ , but also ⌧̂ ⌧ to be small.

heading: 2.1 Time Series Model
text: There are multiple settings for the CPD problem characterized by the probabilistic structure of the observations. In this work, we consider the spatio-temporal factor model for multivariate time series introduced in [1], which captures a wide variety of time series dynamics. Next, we describe the model in detail. The spatio-temporal factor model requires the underlying latent multivariate time series to satisfy two properties. The first property captures the “spatial” structure, i.e., the structure across the N time series (Property 2.1); and the second property captures the “temporal” structure (Property 2.2). Property 2.1 (Spatial structure). There exist R 2 N, where 1R⌧min(N,T), Wr : Z+ ! R, and ↵nr2R 8r2 [R], n2 [N ], such that for any n2 [N ], t2 [T ], fn(t)= PR r=1↵nr Wr(t), where |↵nr| ↵, |Wr(t)| W for constants ↵, W >0.
Property 2.1 implies that there are R “fundamental” time series, and each latent time series fn(·) 8n2 [N ] can be obtained through a weighted combination of these R times series. To capture the temporal structure, additional assumptions are imposed on these fundamental time series Wr(·) 8r2 [R]. To do so, we first introduce the Page matrix representation of a time series. Definition 2.1 (Page Matrix). Given a time series f :Z+ !R, its Page matrix representation over T observations with parameter 1LT is given by the matrix Z2RL⇥bT/Lc with Zij=f(i+(j 1)⇥L) for i2 [L], j2 [bT/Lc].
Property 2.2 (Temporal structure). For each r2 [R] and for any T >1, 1LT , let Z(r)W be RL⇥bT/Lc Page Matrix associated with time series Wr(t), t2 [T ]. Then, the rank of Z(r)W is at most G for any choice of 1LT .
Property 2.2 imposes a particular temporal structure in the time series Wr(·),8r2 [R]. While seemingly restrictive at first, it has been shown in the SSA literature that many standard functions that model time series dynamics satisfy this property. In particular, the Page matrix representation of any finite sum of products of harmonics, low-degree polynomials, and exponential functions is low-rank and satisfies Property 2.2 (refer to Proposition 2.1 in [1]). These functions can capture a rich class of time series dynamics which are typically modeled by three components: stationarity, periodicity, and trend. Periodicity is modeled as harmonics and trend as polynomials; the mixture of both are instances of this model. The inclusion of stationarity in this
model can be understood by considering the spectral representation of stationary processes (refer to Property 4.1 in [39]), which states that any sample-path of a stationary process can be decomposed into a sum of harmonics. Therefore, a finite sum of harmonics provides a good model representation for stationary processes.
For a given 1LT , let Z(n)f 2RL⇥bT/Lc denote the Page matrix induced by the n-th latent time series fn(t) 8n2 [N ]. Further, consider the stacked Page matrix Zf 2RL⇥(NbT/Lc), obtained by a column-wise concatenation of the Page matrices Z(1)f ,...,Z (N) f (See Figure 1). Specifically,
Zf = h Z(1)f Z (2) f ... Z (N) f i . (4)
Proposition 2.1 (Proposition 2.2 in [1]) establishes that under the spatio-temporal factor model satisfying Properties 2.1 and 2.2, the stacked Page matrix is low-rank. We refer to [1] for the proof of this property. Proposition 2.1 (Proposition 2.2 in [1]). Let Properties 2.1 and 2.2 hold. Then for any 1LT , the rank of Z(n)f for n2 [N ] is at most RxG. Moreover, the rank of the stacked Page matrix Zf is also at most RxG.

heading: 2.2 Changes in the Spatio-Temporal Factor Model
text: In this section, we give a characterization of the change points in view of the spatio-temporal model described in Section 2.1. Let us consider the following interpretation of Proposition 2.1: for a multivariate time series f(·) that satisfies Properties 2.1 and 2.2, and for a sufficiently largeL and number of observations T , the columns of its stacked Page matrix span the same linear subspace of RL independently of T . To formalize the notion of sufficiently large L and T , we introduce the definition of the order of the time series. Definition 2.2 (Order of the time series). Given a multivariate time series f :Z+ !RN , we define its order as the integer k>0 such that for all L k and T such that N⇥bT/Lc k, the stacked Page matrix representation of the time series with parameter L over T observations Zf 2RL⇥NbT/Lc has rank k.
Proposition 2.1 suggests that a time series that follows Properties 2.1 and 2.2 has a finite order kR⇥G. Further, any L, T that satisfy L k, N⇥bT/Lc k, are sufficient for this to hold. In lieu of the above, we can characterize changes in the latent time series using the distance from the subspace spanned by the columns of the stacked Page matrixZf . Before we do so, we introduce appropriate notations and definitions. Definition 2.3 (L-lagged vectors). Given a time series f : Z+ ! R and a lag parameter L > 1, the L-lagged vector at time t2{L,...,T} is given by f(t L+1:t):=[f(t L+1),...,f(t)]2RL.
For a multivariate time series f(t), we define the L-lagged matrix at time t2{L,...,T} as the column-wise concatenation of the L-lagged vectors fn(t L+1:t) for n2 [N ]. Specifically,
f(t L+1:t):=[f1(t L+1:t) ... fN(t L+1:t)]. (5) We analogously denote the L-lagged observation matrix at time t by X(t L+1:t). Let 1(f(t L+1:t)) be the largest singular values of the L-lagged matrix of f(·) at time t2{L,...,T}. We define the following quantities for the pre-change functions f0(·):
max,01 :=maxt<⌧ 21(f0(t L+1:t)),
min,0 1 :=mint<⌧ 21(f0(t L+1:t), (6)
with max,11 and min,1 1 analogously defined for f1(·) over the time range t ⌧ + L 1. Now let L0 := span(Zf0)⇢RL, where Zf0 is the stacked Page matrix induced by f0(·) with parameter L over T0 observations (assuming no change occurs at any tT0). If the latent time series continues to follow the same function for t>T0, then the columns of its L-lagged matrix belong to the subspace L01. The
1Since Proposition 2.1 is stated for the Page matrix representation, this is true for L-lagged vectors that are non-overlapping and spaced by L observations. However, an analogous result to Proposition 2.1 can be proven for the Hankel matrix representation making this valid for any segment of L observations. For simplicity, we omit discussion of the Hankel matrix but we refer to [1] for such results.
occurrence of the change point at t= ⌧ forces some or all of the columns of the the L-lagged matrix to leave the pre-change subspace L0. Eventually, i.e., for t ⌧ +L 1, the columns of the L-lagged matrix will be in the post change subspace L1 which corresponds to the function f1(·). Following this characterization of change points, we revisit our hypotheses in (2) and reformulate them as
H0 : span(IE[X(t L+1:t)])⇢L0 (7) H1 : span(IE[X(t L+1:t)]) 6⇢L0.
The occurrence of a change point can be captured by a test statistic that monitors the distance between the lagged vectors and the pre-change subspace. This motivates the CPD algorithm we propose next.

heading: 3 Algorithm
text: The core steps in the proposed CPD algorithm are to first estimate the pre-change subspace from the stacked Page matrix of noisy observations X(t), then repeatedly estimate the distance between this subspace and future lagged vectors. If our estimation of the subspace is accurate, the distance will remain small as long as the observations continue to follow the same latent time series. Once a change occurs, the lagged vectors are forced to leave the subspace causing the distance to become larger. The algorithm uses the subspace distance as a detection score to construct a CUSUM statistic to perform the hypothesis test in (7). In what follows, we review the CUSUM procedure for sequential hypothesis testing then we introduce the full CPD algorithm.

heading: 3.1 CUSUM Procedure
text: Definition 3.1 (CUSUM Statistic). Given an observation X(t) we assign a detection score D(t), and define the CUSUM statistic of the observations {X(1),...,X(t)} as y(t):=max1it Pt j=iD(j)
It is easy to verify that the CUSUM statistic in Definition 3.1 can be, more conveniently, computed recursively (also known as Lindley’s recursion in queueing theory [26]) as
y(t)=max{y(t 1)+D(t),0}, y(0)=0. (8) Using the CUSUM statistic, and for a threshold h>0, we define the decision rule at each time t as
H(t)=H1y(t) h,
and the estimated change point in (3) can be written as ⌧̂=inf{t|y(t) h}. The CUSUM procedure for change point detection, as introduced by Page [34] requires the assigned detection score D(t) to satisfy the following Property 3.1. Property 3.1. A valid CUSUM detection score D(t) has a negative expectation when no change is present, and a positive expectation when a change occurs, that is E[D(t) |H0]<0 and E[D(t) |H1]>0.

heading: 3.2 Algorithm Description
text: Our mSSA algorithm for CPD has five parameters: (1) the number of initial observations T0 used to estimate the subspace L0, (2) the lag parameter 1<L T0, (3) the estimated order of the time series k̂ >0, (4) the CUSUM test threshold h>0, and (5) the shift-downwards constant c 0. For simplicity and without loss of generality2, let T0 be an integer multiple of L. We introduce two additional integer quantities: M :=T0/L and M̄ :=N⇥M . In what follows, we detail the four main steps of the algorithm. Step I: Base Matrix Construction. Transform the observationsXn(1:T0):=[Xn(1),...,Xn(T0)], n2 [N ] into a Page matrix Z(n)X 2RL⇥M as per Definition 2.1. Then, form the stacked Page matrix ZX2RL⇥M̄ as in (4). We refer to ZX as the base matrix.
Step II: Subspace Estimation. Let ui for i2 [L] be the left singular vectors of the base matrix ZX sorted as per associated singular values in the decreasing order. We group the singular vectors into two matrices Û0=[u1,...,uk̂] and Û?=[uk̂+1,...,uL]. Then, we estimate the pre-change subspace L0 as:
L̂0=span(Û0). (9)
Let L̂?=span(Û?) be the orthogonal complement of the subspace L̂0 in RL. 2If T0 is not integer multiple of L we can introduce a new T 00=L⇥bT0/Lc.
Step III: Detection Scores. For t>T0, construct the L-lagged matrix X(t L+1 : t) analogous to (5). We define the subspace detection score, which measure the Euclidean distance between the columns of the L-lagged matrix and the estimated subspace L̂0, as
D(t)=kÛT?X(t L+1:t)k2F c. (10)
Step IV: CUSUM Test. For t>T0, compute the CUSUM statistic y(t) using (8) with y(T0)=0. Steps III and IV are repeated until a change point is detected when y(t)>h. That is,
⌧̂= inf t>T0
{t |y(t) h} (11)
To adapt this algorithm for multiple change point detection, when a change point is detected at ⌧̂ , the algorithm restarts with a new base matrix being constructed using the segment X(⌧̂ : ⌧̂+T0 1). Refer to Figure 2 for a visual depiction of these steps, and Algorithm 1 in Appendix C for a summary of the algorithm.

heading: 3.3 Choice of Parameters
text: Below, we provide general guidelines for selecting the parameters T0 and L, k̂. Further, we detail a data-driven procedure for choosing the CUSUM parameters c and h in Appendix D.1. A general recommendation is to choose T0 to be as large as possible and set L= b p min(N,T0)⇥T0c to achieve a good estimation of the pre-change subspace (see Proposition 4.1). However, if T0 is too large we might potentially miss or smooth out a change point. In practice, one heuristic used to choose T0 is as follows: given a starting T0, construct the stacked Page matrix and compute its “effective rank” (e.g. top k singular values whose sum of squares is at least 90% of the squared Frobenius norm). Then keep increasing T0 till the “effective rank” of ZX stops changing. Once it does, this indicates that we have captured the principal subspace L0. Additionally, the “effective rank” of matrix ZX can be used as k̂, the estimated order of time series. Indeed, it is well understood that the largest singular value of a random A⇥B matrix with independent zero mean sub-Gaussian entries scales as O( p A+ p B). That is, the L⇥M̄ matrix ZX, which can be viewed as a summation of the ‘signal’ matrix (corresponding to the stacked time series f0(·)) and ‘noise’ matrix (corresponding to the stacked observation noise e(·)), is such that the top k singular values scale as ⇥( p L⇥M̄) while other singular values scale as O( p L+ p M̄). Therefore, with large enough L and M̄ , a clear separation is observed and is the reason for this heuristic to work. See [8] and [10] for use of such reasoning and heuristic in the context of matrix estimation literature.

heading: 4 Theoretical Analysis
text: In this section, we establish theoretical performance guarantees of the described algorithm. To that end, we formally introduce the performance objective. Recall that if a change point is present at time t=⌧ , then
f(t)= ⇢ f0(t) if t<⌧ f1(t) if t ⌧,
where f0(t) and f1(t) each satisfying Properties 2.1 and 2.2 with parameters Ri,Gi, ↵i, and Wi for i2{0,1}. The algorithm declares ⌧̂ as the estimate of ⌧ . Therefore, ideally if ⌧<1 then we wish to have ⌧̂ ⌧ , and ⌧̂ ⌧ as small as possible. Further, if ⌧=1, i.e. no change point is present, then we want ⌧̂ to be as large as possible. This is captured through average running length of type 0 and type 1 defined next.
Average running length. This metric is typically considered in the setting of sequential hypothesis testing [34]. Let T1=T0+1, i.e., one time-step after the time horizon used for the base matrix. Then, define
ARL0=E[⌧̂ T0|⌧=1]⌘E1[⌧̂ ], ARL1=E[⌧̂ T0|⌧=T1]⌘ET1[⌧̂ ]. The desired behavior is to have a large ARL0 and small ARL1. In what follows, we provide bounds on these two quantities for the specific choice of the detection score D(t) as a function of the threshold h.
Subspaces distance and affinity. We introduce the definitions of two quantities that will be used in the main results. The first is ✏ which describes how well is our estimation of the pre-change subspace. The other is , which describes the similarity between the pre and post change subspaces. Recall the subsapces L0, L1, L̂0, and L̂?, and their corresponding basis matrices. Then define
✏⌘dist(L0,L̂0):=kÛT?U0kop, ⌘affinity(L0,L1):=kUT0U1kop (12)
Operating assumptions. We make certain operating assumptions to establish our results stated next.
Assumption 4.1 (Parameter Selection). We assume that the parameters of the algorithm are chosen such that T0⌧ , L=b p min(N,T0)⇥T0c, LM̄ , and L k.
Assumption 4.2 (Gaussian Noise). We shall assume i.i.d. Gaussian noise en(t)⇠N(0, 2).3
Assumption 4.3 (Balanced spectra). Let Zf be the stacked Page matrix associated with the time series f(·) that satisfies properties 2.1 and 2.2 and has order k, where under the setup of Proposition 2.1, kRxG. Then for L=b p min(N,T0)⇥T0c k, the stacked page matrix is such that k(Zf) K p NT0/ p k for some absolute constant K>0, where k(·) is the k-th largest singular value of its argument.
Based on the these assumptions, we can quantify the bound on the estimation error ✏, as follows. Proposition 4.1. Let assumptions 4.1-4.3 hold. Let ✏ be the estimation error defined in (12). Let
E :={✏q}, with q := 4 p R0⇥G0
K(min(T0,N)⇥T0)1/4 . (13)
Then P(E) 1 2exp( 2 p NT0). Note that the choice of L=b p min(N,T0)⇥T0c plays a key role in minimizing the subspace estimation error. Precisely, the bound on the estimation error (q) scales inversely with the square root of the minimum dimension of the Page matrix. Thus, to minimize the error, we construct a nearly square base matrix.
Main results. Now we state the main result providing a lower bound on ARL0 and an upper bound on ARL1 for a given parameter h>0. In particular, it shows thatARL0 scales exponentially in hwhileARL1 scales linearly in h. Proofs of Proposition 4.1, Theorems 4.1, 4.2, and 4.3 can be found in Appendix B. Theorem 4.1 (Bound on ARL0). Let Assumptions 4.1-4.3 hold. Let ⌧̂ be the estimate of ⌧ as per (11). Given h>0, k̂=k, and a choice of c such that (recall notation in (6))
c>N(L k) 2+qR0( max,01 )2. (14) Then, if ⌧=1,
ARL0=E1[⌧̂ ] C0 exp(C1h), where C0,C1 are positive constants that depend on the parameters ,R0,G0,N,L and T0. Theorem 4.2 (Bound on ARL1). Let Assumptions 4.1-4.3 hold. Let ⌧̂ be the estimate of ⌧ as per (11). Conditioned on E and given h>0, k̂=k, and a choice of c such that (recall notation in (6))
c<N(L k) 2+(1 q)( min,11 )2, (15) 3The analysis and guarantees obtained will extend easily for generic sub-Gaussian distribution.
if ⌧=T0+1, then
ARL1C̃0+C̃1h+C̃2 exp ⇣ C̃3h ⌘ , (16)
where C̃0, C̃1, C̃2, and C̃3 are positive constants that depend on the parameters ,c,N,L,R1, ↵1, W 1.
Next, we show that there exists a feasible choice of c for both theorems to hold, provided that the pre and post change subspaces are sufficiently “different". Theorem 4.3. Let be the affinity between the pre and post change subspaces defined as per (12). Then there exists a feasible choice of c that satisfies both equations (14) and (15) if is such that
<1 q
0
@1+R0 max,01 min,11
!21
A. (17)
Interpretation. Recall that is a measure of "similarity" between subspaces which takes a value in [0,1]: it is 1 if the two subspaces intersect, otherwise, it is <1. Theorem 4.3 implies that as q gets smaller, which may be achieved as more data points (NT0) become available, we can tolerate larger without compromising the feasibility of c. Under perfect estimation (e.g. no noise in data or infinite data availability), i.e. q = 0, any <1 will guarantee the feasibility of c.

heading: 5 Evaluation on Benchmark Datasets
text: We present a comparative evaluation of two variants of the proposed mSSA algorithm and four other standard methods for CPD: BinSeg [21], KL-CPD [7], Microsoft SSA [29], and BOCPDMS [23].
Datasets. The evaluation is conducted on four benchmark datasets (see Table 2). Details of these datasets, their sources, and physical interpretations of change points are given in Appendix D.2.
Evaluation metric. For evaluation we use the F1-score metric, which is defined as: F1= TPTP+12 (FP+FN) , where TP, FP, and FN denote the number of true positives, false positives, and false negatives, respectively. Following standard practice cf. [42, 27, 46], we use the following soft true positive (TP) rule: for a constant ⌘>0, if the algorithm detects a change point ⌧̂ such that |⌧̂ ⌧ |⌘ where ⌧ is a labeled change point, a TP is recorded, otherwise a FP is recorded. Here we report the results for ⌘=10. To avoid double-counting, we ensure that each labeled and detected change point is associated with at most one TP. Further, we use the convention of counting the first index of the time series (t = 1) as a TP.
Note that we choose the F1-score as our metric over the average running length (i.e. ARL0 and ARL1) for two practical reasons: (i) it is a more commonly used metric in the change point detection literature; (ii) it can be used for both online and offline change point detection methods, both are used in our experiments.
Practical implementation (mSSA-MW). The Algorithm described in Section 3 is meant for a single change point detection. In practice, there can be multiple change points in a time series and hence we use a natural Moving Window variant. This is described in details in Algorithm 2 in Appendix C. In this section, we run the evaluations for both the original variant (mSSA) and the Moving Window variant (mSSA-MW).
Parameter Selection. When applying a CPD algorithm as a statistical tool for the automatic detection of change points, a common challenge that arises in practice is tuning the algorithm’s parameters. This
is particularly challenging in the task of CPD as it is often applied in an unsupervised manner due to the scarcity of labels (change points are anomalous events). To overcome this, we follow the evaluation setup proposed in [42], where for each dataset we report (1) the F1-score using default parameters, and (2) the best F1-score over a grid search of parameter configurations. We provide details about the default parameters and the grid search for all algorithms in Appendix D.1.
Results. We compare with retrospective change point detection methods (BinSeg and KL-CPD) and real-time methods (BOCPDMS and Microsoft SSA). Two of these algorithms (BinSeg and Microsoft SSA) only support univariate time series. For datasets of higher dimension, we report the results of these algorithms on a univariate time series constructed from the L2 norm of the observation vector (in Appendix D.5, we explore other approaches of applying univariate algorithms to multivariate data).
In Table 3, the best and default performance of each algorithm are shown as the average and standard deviation of the F1-scores across time series in each dataset4. mSSA-MW algorithm shows consistently superior performance over other methods, including the retrospective ones, across all datasets in the best parameters setting, and all but one dataset using the default parameters.
For all datasets, we notice a significant improvement of mSSA-MW over the mSSA variant. We describe one possible explanation for this behavior. If the order of the time series is large compared to the dimensions of the base matrix, then the estimation of the subspace is inaccurate, leading to a large detection score even when no change occurs. While this is also true for the moving window approach, therein the base matrix estimates a new subspace in every time step, which captures the local structure of the time series at that time. While it doesn’t necessarily capture the full signal, its estimated subspace has more proximity to the test lagged vectors. Another plausible explanation is the existence of gradual and small changes in the data that aren’t labeled as change points. In such a case, constructing a fixed base matrix makes the algorithm more sensitive to such changes, while having a moving window smooths out such changes.

heading: 6 Evaluation on Synthetic Data
text: To further explore its performance, we evaluate mSSA on synthetic datasets, using the same setup of Section 5. We follow the practice of [7] and [27] and create synthetic datasets using our proposed time series model, each with a different representative mode of change. Specifically, we generate four sets
4For the occupancy, we only have one time series, and hence the standard deviation is not reported.
of synthetic time series, each is a mixture of harmonics and a polynomial trend with the following modes of change: (i) jumping mean (Mean), a univariate time series where at each change point, a constant offset is added; (ii) scaling signal energy (Energy), a univariate time series where at each change point, the amplitude of the harmonics is changed; (iii) changing frequency (Frequency), a univariate time series where at each change point, the frequencies of the harmonics are changed; (iv) mixed change (Mixed), a multivariate time series where at each change point, a mixture of the above changes is applied. Refer to Table 2 for more details and Appendix D.3 for details about the data generation process.
Results. In Table 3, the best and default performance of each algorithm are shown as the average and standard deviation of the F1-scores across time series in each synthetic dataset. Both mSSA and mSSA-MW implementations show consistently superior performance over other methods across all datasets in the best parameters setting, and they compete with Microsoft SSA in the default parameters setting. Notice that the results on the simulated data show a similar performance of the two variants of the algorithm (mSSA and mSSA-MW). Here the data is constructed to follow the proposed spatio-temporal model with a low-rank Page matrix, and enough time steps are available to construct a sufficiently large base matrix. As a result, the subspaces of the fixed base matrix (in mSSA) and the moving window base matrices (in mSSA-MW) become almost identical with small variations due to the noise.
For the univariate datasets, another algorithm that performs well is the Microsoft SSA, which is proposed for a similar time series model. However, we notice a significant degradation in the performance in the case of the multivariate time series. This highlights the importance of extending SSA to accommodate multivariate time series. The BinSeg algorithm performs relatively well on the Energy and Mean datasets. However, the algorithm is not capable of detecting the Frequency mode of change, nor the Mixed mode. We note that BinSeg has designated methods for detecting mean and variance changes, which might be the reason for its good performance in the Energy and Mean datasets. KL-CPD and BPCPDMS algorithms are built with the assumption of i.i.d. samples, so the highly dependant structure and the time-varying mean of our data might explain the bad performance of these algorithms on the simulated datasets.
Multivariate data helps. Here, we explore the benefits of utilizing information across time series in mSSA. Specifically, we study the trade-off between ARL0 and ARL1 as N (number of time series) increases. We start by generating 25 synthetic time series as described in details in Appendix D.4. Then for each value ofN2 [5,10,15,20,25] we perform the following experiment. First, we assume that there is no change point (⌧=1) and run the algorithm for 10 trials. In each trial, we vary the threshold h and record ⌧̂ which represents the time of false alarm. Next, we introduce a change point at T0+1, and run the algorithm for another 10 trials. Again we vary the threshold h and record ⌧̂ which represents the detection delay.
Figure 3 shows the average results of the 10 trials for each pair (N,h). Each line of connected dots corresponds to the same value of N with different values of h. As shown in the figure, as we fix the value of ARL0 at a certain acceptable rate, we can achieve better ARL1 (i.e. detection delay) using a time series of a larger dimension. This shows that mSSA is effectively utilizing information across multiple time series to improve the detection performance over the univariate SSA.

heading: 7 Limitations
text: Our method achieves encouraging theoretical and empirical results for change point detection. However, like any such research program, there are various limitations suggesting exciting directions for future research. A prominent one being extending the theoretical analysis to the setting of multiple change points. In a sense, this would be important step forward for the field of change point or anomaly detection.
Acknowledgements and Funding Disclosure
This work was supported in parts by the MIT-IBM project on time series anomaly detection, NSF TRIPODS Phase II grant towards Foundations of Data Science Institute, KACST project on Towards Foundations of Reinforcement Learning, and scholarship from KACST (for Arwa Alanqary and Abdullah Alomar).

Review summaries:
The article presents a change point detection (CPD) method for multivariate data. The proposed methodology operates by monitoring the space generated by segments of the data and then assesses how this space varies. A sound theoretical treatment and convincing experiments are reported as well.  

This paper presents a multivariate singular spectrum analysis to perform change point detection. A spatio-temporal model is introduced to model the dynamics and CUSUM statistic was used to perform online change point detection. The experiment results showed the effectiveness of the proposed method.

This paper studies the online change-point detection for multivariate time series and proposes a new cusum-like detection statistic based on singular spectrum analysis. Theoretical approximations for the ARL and EDD are provided. This paper also validates the performance of the proposed method using a variety of synthetic and real data sets.

The paper studies change-point detection in time-series using a combination multi-variate singular spectrum analysis (SSA) and CUSUM.   While multi-variate SSA and CUSUM are established techniques, the contribution of the paper is to (a) analyze the online CPD setting with a dependent (non-iid) timeseries. (b) Conduct rigorous analysis of the tradeoff of detection delay vs. false-alarm within the MSSA model.  (c) The exact version of MSSA is a bit different from the classical one.