Dynamic Interpretable Change Point Detection

A CP can be caused by a change in the correlation between features. This results in a change in the covariance of the joint distribution that can be identified as a state change in the feature network. The evolving dynamics of features can be modeled using graphical models, i.e. at every time point, the interactions of features can be modeled as a graph network, with nodes and links corresponding to each feature and correlation between sets of variables, respectively (Figure 2(a)). However, in a time-series setting, estimating the graphical network at every time step is computationally challenging. To estimate these networks to detect CPs, we use Time-Varying Graphical Lasso (TVGL) (Hallac et al., 2017), an efficient algorithm for estimating the inverse covariance matrix (precision matrix) of multivariate time-series with time-varying structures. TVGL infers the structure of graph networks by estimating a sparse time-varying inverse covariance matrix, $P_{t}=\Sigma^{-1}_{t}$ ($\Sigma_{t}$ indicating the covariance matrix at time $t$), of the variables for all $t\in[T]$. TVGL extends the graphical lasso problem to dynamic networks by allowing covariance $\Sigma$ to vary over time, taking into account how the relationships between signals evolve. TVGL method enforces sparsity in the covariance matrix for less computational cost and easier interpretability of the graphical network structure. A scalable message-passing algorithm called Alternating Direction Method of Multipliers (Banerjee et al., 2008) is employed to estimate the sparse inverse covariance. Additional details on the TVGL method can be found in the Supplementary Material.

Learning the inverse covariance estimates from windows of data reveals the underlying evolutionary patterns present in the time-series. We employ TVGL technique to identify points of change in feature interactions between adjacent local windows of time. To promote the identification of shifts in the covariance pattern of features, we integrated an L2-norm penalty function into the estimation of the matrix inverse. Moreover, to ensure the invertibility of $\Sigma_{t}$, we applied feature standardization and removed highly correlated features before feeding them to TVGL. The partial correlation of two features can be estimated from the joint precision matrix entries as $-\frac{P(i,j)}{\sqrt{P(i,i)P(j,j)}}$, which means contrasting consecutive precision matrix entries over time quantifies the change in the features’ correlation. By taking the difference between the absolute values of adjacent precision matrices, we quantify these changes in the distribution. A value close to 0 in this matrix indicates nearly identical estimations of the network and therefore no CP in that feature interaction (see Algorithm 1). A negative value indicates an increase in absolute correlation and a positive value indicates a decrease.

Assuming in a time-series sample $X$, each $X_{t}$ is independently generated from a joint probability distribution $p_{t}(\cdot)$, a CP occurs at time $t^{*}$ if observations after $t^{*}$ are generated from a different distribution. To compare the probability distributions of adjacent windows, we employ a non-parametric two-sample testing procedure called MMD Aggregate (MMDAgg), introduced in Schrab et al. (2021). Let $\Delta^{-}$ represent the initial size of the window of $X$ before a query point $t$, and let this prior window be denoted by $X_{\Delta^{-}}^{t}=\{X_{t-\Delta^{-}},...,X_{t-1},X_{t}\}$. Similarly, the window of future observations can be denoted by $X_{\Delta^{+}}^{t}=$ $\{X_{t+1},...,X_{t+\Delta^{+}-1},X_{t+\Delta^{+}}\}$, where ${\Delta^{+}}$ represents the length of a future window. Kernel-based MMD tests serve as a measure between two probability distributions. With the statistical test threshold $\alpha$, if the null hypothesis $H_{0}$, is rejected, the time-series may be partitioned by a CP at $t^{*}$, signifying that measurements in the $X_{{t^{*}}-\Delta^{-}:t^{*}}$ windows come from a different distribution than measurements in $X_{t^{*}:t^{*}+\Delta^{+}}$. The performance of a single kernel-based MMD test typically depends on the choice of kernel and kernel bandwidths. Since we compare adjacent windows with restricted number of samples, any loss of data to kernel bandwidth selection can be detrimental to our method’s performance. To overcome this problem, MMDAgg aggregates multiple MMD tests using different kernel bandwidths, ensuring maximized test power over the collection of kernels used and eliminating the need for data splitting or arbitrary kernel selection. The method considers a finite collection of bandwidths where the aggregated test is defined as a test that rejects $H_{0}$ if one of the tests over a given bandwidth rejects $H_{0}$.

We propose to dynamically establish the window size based on the presence of CPs (Algorithm 2). Let $\Delta^{-}$ represent the size of the dynamic window of data points from the last estimated CP ($\widetilde{t}$), up until the current time point ($t$). Starting with a constant $\Delta^{+}$ and a small $\Delta^{-}$ window, the length of the running window $\Delta$ increases with each new observation until a new CP occurs, according to the MMD test. If a significant change in distribution isn’t detected by the MMD test, i.e. the MMD score is smaller than a pre-defined threshold $\epsilon$, the two sub-sequences are combined and compared against the next sub-sequence in the series. This process is also explained in Figure 2(b). Our dynamic windowing method eliminates the need for repetitive fixed-window comparisons and utilizes a growing sample set for the MMD test.

For determining the final CP score, we need to meaningfully ensemble the DistScore with the CovScore, which is challenging because the covariance score is bounded while MMD is a positive unbounded score. Hence, TiVaCPD incorporates kernel normalization in the MMDAgg algorithm. we use a generalization of Cosine normalization (Ah-Pine, 2010) to normalize our kernels so as to have a similarity index. For a given kernel function, $K^{z=1}(x,y)$ represents the normalized kernel of order $z$. We use the generalized mean with exponent $z=1$ (arithmetic mean), which means $K^{z=1}(x,y)=\frac{K(x,y)}{M^{z=1}(K(x,x),K(y,y))}$, where $M^{z=1}(a_{i=1}^{p})=\frac{1}{p}\sum^{p}_{i=1}(a^{z}_{i})$. This normalization technique projects the objects from the feature space to a unit hypersphere and guarantees $|K^{z=1}(x,y)|\leq 1$.

An effective method for combining unsupervised CPD scores is poorly studied. Here, we use an ensemble method that utilizes the score differences to highlight the scores that contribute the most to representing the CPs. Algorithm 3 indicates our ensemble technique for generating a unified score $S$ by combing CovScore and DistScore. We use four score variants derived from CovScore and DistScore to generate dynamic dissimilarity weight $W$ to effectively aggregate the scores into a unified one (As shown in Figure 1). The four scores used are as follows: a) b) Normalize(CovScore) and Normalize(DistScore): standardized CP scores using z-score normalization to bring them into the same scale, c) SG(CovScore): Since CovScore is sensitive to small distributional changes that can lead to false change point detection, we mitigate the risk of detecting spurious CPs by applying Savitzky-Golay (SG) smoothing filter. SG filtering technique is a widely used method for smoothing patterns and reducing noise in time series data (Press & Teukolsky, 1990). d) SG_combined: We also apply the filter to the sum of filtered CovScore and DistScore, which effectively reduces noise and improves the performance of CPD. The ensemble approach utilizes a weighted average of aforementioned four scores. The importance weight $W$ is determined by computing the mean absolute difference between scores. By assigning greater weight to scores that are more dissimilar, we can effectively identify CPs that may have been missed by other scoring methods. The weights are calculated for each 21 time-points window, as the distribution of scores’ importance may vary over time. To locate the exact time of the CPs, we look for peaks in the ensemble scores by searching for local maxima with a threshold that is used to remove number of false positive CPs created by noise. Using a threshold for peak detection is a common practice that requires careful tuning based on circumstances.

TiVaCPD offers valuable insights into the underlying nature of the observed change points. Detecting both changes in correlation and data distribution, such as an inverse correlation between heart rate variability and heart rate, is crucial in clinical practice, as it can indicate an adverse health event. The DistScore detects changes in the underlying distribution of the time-series, while the CovScore detects changes in the correlation structure of feature pairs, providing a detailed analysis of the feature dynamics at each time step. This information allows for the categorization of CPs, thereby enhancing interpretability, as shown in Figure 3. This figure presents a multivariate time series sample showcasing TiVaCPD results, featuring ConvScore, DistScore, and a weighted ensemble score. In addition, CovScore’s heatmap illustrates the feature pairs that caused a CP, and TiVaCPD identifies the direction of correlation change at each CP. It also includes a comparative analysis with other CPD methods, such as KLCPD, Roerich, GraphTime, and TIRE, providing a comprehensive evaluation of the TiVaCPD results against alternative CPD methods. The first two CPs ($CP_{1}$ and $CP_{2}$) are caused by changes in the correlation between features 0 and 1, where $CP_{1}$ corresponds to a negative change in correlation and $CP_{2}$ corresponds to a positive change in correlation. CP 3 is caused by changes in the mean between features 0 and 1 as the DistScore is high. CP 4 is caused by changes in a combination of variance, correlation, and mean, which are simulated to resemble real-world data.

We demonstrate the performance of our method compared to multiple baselines on four simulated and two real-life multivariate datasets commonly used in the CPD literature. Different simulations test the functionality of CPD methods on a variety of potential CP scenarios. All hyper-parameters are determined based on random search over 10% of the datasets (more details on best parameters and sensitivity to hyperparameter change are provided in the supplementary material).

We compare the performance of TiVaCPD with SOTA CPD methods on various types of CPs. ¹¹1The implementation of the baselines was done using publicly released source code by the authors. The selected SOTA approaches include those that measure a change in distribution and those that focus on the graphical structure of features over time. In addition, we conducted an ablation study on various parts of the TiVaCPD score to analyse each component impact.

Kernel Change Point Detection (KLCPD) (Chang et al., 2019) is a kernel learning framework for CPD that uses a two-sample test and optimizes a lower bound of test power via an auxiliary generative model. For this method, we used window sizes $w\in[10,25]$ for all experiments and trained the model for $25$ epochs, unless more training led to improved results. Consistent with our own post-processing steps, we performed peak detection to detect the exact time of change.

Roerich (Hushchyn & Ustyuzhanin, 2021) is a CPD method based on direct density ratio estimation to detect the change in distribution. We set all parameters to default and use window sizes $w\in[10,25]$ for all experiments.

Group Fused Graph Lasso (GraphTime) (Gibberd & Nelson, 2015) is a time-varying graphical model based on the group fused-lasso. Similar to TiVaCPD  it uses the graphical model to model the dependencies of variables in time series. GraphTime models the temporal dependencies between variables while TiVaCPD models the pairwise dependencies to identify CPs.

Time-Invariant Representation (TIRE) De Ryck et al. (2021) is an autoencoder-based CPD method, we used window sizes $w$ between ${10,25}$ for all experiments. For the $domain$ parameter, we chose $both$ to include both time, and frequency domains. We used $200$ epochs to train the model.

Tables 1-3 show the performance of TiVaCPD and all baselines on 4 simulated and 2 real-world data sets. Results are reported as $F_{1}$ scores, as well as precision and recall, measuring how well it detects CP locations in time. To estimate performance metrics, we define a margin of error $M$ for the exact location of CP, which is common practice in the CPD literature (van den Burg & Williams, 2020; Deldari et al., 2021b). Given a user-defined margin of error, $M>0$, an estimated CP is a True Positive (TP) if the distance between the ground truth ($t^{*}$) and the estimated CP ($\widetilde{t}$) is smaller than the margin, i.e. $|t^{*}-\widetilde{t}|\leq M$. As explained in Figure 4, if an estimated CP falls outside the margin, then it is considered False Positive (FP), i.e. $\widetilde{t}\notin[t^{*}-M/2,t^{*}+M/2]$. We show the impact of margin values $M=\{5,10\}$ on the performance of all baselines in the supplementary material.

Our method, TiVaCPD, effectively detects various types of CPs in simulated datasets by using different components of the score to capture specific CP’s types. The CovScore component excels at detecting CPs caused by changes in correlation, the DistScore component detects CPs caused by changes in distribution, and the ensemble score combines the strengths of both components for improved overall performance. Additionally, each component of the score individually outperforms baseline methods, indicating that the performance boost is not just the result of ensembling the scores, but each component of the score itself does a better job at detecting the relevant CPs compared to baseline methods. Although Roerich performs well in detecting simple CPs; however, it fails to address more complex CPs in both synthetic and real-world data. Our method is the only one that performs well in detecting complex CPs. The performance results on the real-life datasets are reported in Table 3, for all baselines with a margin value of 5. We show that our method outperforms all baselines in detecting the exact time of CP in real-world and simulation datasets. Figure 3 shows a graphical representation and comparison of the different CPD methods for a time-series sample (row 1), and demonstrates how different methods generate scores for CPs. The different components of TiVaCPD are shown in rows 2 and 3, and show what category of CPs they identify. The heatmaps of CovScore can be used to identify in which pair of features the change in correlation has occurred, and also in which direction the change was.