Toeplitz Inverse Covariance-based Clustering of Multivariate Time Series Data

1 Introduction

Many applications, ranging from automobiles [Miyajima et al., 2007] to financial markets [Namaki et al., 2011] and wearable sensors [M¨orchen et al., 2005], generate large amounts of time series data. In most cases, this data is multivariate, where each timestamped observation consists of readings from multiple entities, or sensors. These long time series can often be broken down into a sequence of states, each defined by a simple "pattern", where the states can reoccur many times. For example, using automobile sensor data, a single driving session can be expressed as a sequential timeline of a few key states: turning, speeding up, slowing down, going straight, stopping at a red light, etc. This representation can be used to discover repeated patterns, understand trends, detect anomalies and more generally, better interpret large and high-dimensional datasets.

To achieve this representation, it is necessary to simultaneously segment and cluster the time series. This problem is more difficult than standard time series segmentation [Hallac et al., 2016; Himberg et al., 2001], since multiple segments can belong to the same cluster. However, it is also harder than subsequence clustering [Begum et al., 2015; Smyth, 1997] because each data point cannot be clustered independently (since neighboring points are encouraged to belong to the same cluster). Additionally, even if one is able to simultaneously segment and cluster the data, the question still arises as to how to interpret the different clusters. These clusters are rarely known a priori, and thus are best learned through data. However, without prior knowledge, it is difficult to understand what each of the clusters refers to. Traditional clustering methods are not particularly well-suited to discover interpretable structure in the data. This is because they typically rely on distance-based metrics, such as dynamic time warping [Berndt and Clifford, 1994]. These methods focus on matching the raw values, rather than looking for more nuanced structural similarities in the data.

In this paper, we propose a new method for multivariate time series clustering, which we call Toeplitz inverse covariance-based clustering (TICC). In our method, we define each cluster as a dependency network showing the relationships between the different sensors in a short (timeinvariant) subsequence (Figure 1). In these networks, known as Markov random fields (MRFs), an edge represents a partial correlation between two variables [Koller and Friedman, 2009; Rue and Held, 2005; Wytock and Kolter, 2013]. Partial correlations are used to control for the effect of other confounding variables, so the existence of an edge in an MRF implies that there is a direct dependency between two variables. Therefore, an MRF provides interpretable insights as to precisely what the key factors and relationships are that characterize each cluster. In TICC, we discover this structure by solving a constrained sparse inverse covariance estimation problem [Friedman et al., 2008; Yuan and Lin, 2006], which we call the Toeplitz graphical lasso, to learn the adjacency matrix of the MRF dependency network [Banerjee et al., 2008; Wainwright and Jordan, 2006].

To solve the TICC problem, we use an expectation maximization (EM)-like approach, based on alternating minimization, where we iteratively cluster the data and then update the cluster parameters. Even though TICC involves solving a highly non-convex maximum likelihood problem, our method is able to find a (locally) optimal solution very efficiently in practice. We then implement our TICC method and apply it to both real and synthetic datasets. We start by evaluating performance on several synthetic examples, where there are known ground truth clusters. We compare TICC with several state-of-the-art time series clustering methods, outperforming them all by at least 41% in terms of cluster assignment accuracy. We also quantify the amount of data needed for accurate cluster recovery for each method, and we see that TICC requires 3x fewer observations than the next best method to achieve similar performance. We then analyze an automobile sensor dataset to see an example of how TICC can be used to learn interpretable insights from real-world data. Applying our method, we discover that the automobile dataset has five true clusters, each corresponding to a "state" that cars are frequently in. We then validate our results by examining the latitude/longitude locations of the driving session, along with the resulting clustering assignments, to show how TICC can be a useful tool for unsupervised learning from multivariate time series.

2 Problem Setup

Consider a time series of T sequential observations, where xi ∈ Rn is the i-th multivariate observation. Our goal is to cluster these T observations into K clusters. However, instead of clustering each observation in isolation, we treat each point in the context of its predecessors in the time series. Thus, rather than just looking at xt, we instead cluster a short subsequence of size w ≪ T that ends at t. This consists of observations xt-w+1, . . . , xt, which we concatenate into an nw-dimensional vector Xt. Rather than clustering the observations directly, our approach consists of clustering these subsequences X1, . . . ,XT . We do so in such a way that encourages adjacent subsequences to belong to the same cluster, a goal called temporal consistency.

Toeplitz Inverse Covariance-Based Clustering. We define each cluster by a Gaussian inverse covariance Θi ∈ Rnw×nw . Recall that inverse covariances show the conditional independency structure between the variables [Koller and Friedman, 2009], so Θi defines a Markov random field encoding the structural representation of cluster i. In addition to providing interpretable results, sparse graphical representations are a useful way to prevent overfitting [Lauritzen, 1996]. Our objective is to solve for these K inverse covariances Θ = {Θ1, . . . ,ΘK}, one per cluster, and the resulting assignment sets P = {P1, . . . , PK}, where Pi ⊂ {1, 2, . . . , T}, and each point is assigned to exactly one cluster. Our overall optimization problem is

We call this the Toeplitz inverse covariance-based clustering (TICC) problem. Here, T is the set of symmetric block Toeplitz matrices and ‖λ ◦ Θi‖1 is an ℓ1-norm penalty of the Hadamard (element-wise) product to incentivize a sparse inverse covariance (where λ is a regularization parameter). Additionally, ℓℓ(Xt,Θi) is the log likelihood that Xt came from cluster i,

where µi is the empirical mean of cluster i. In Problem (1), β is a parameter that enforces temporal consistency, and 1{Xt-1 ̸∈ Pi} is an indicator function checking whether neighboring points are assigned to the same cluster. We constrain the Θi's to be block Toeplitz, to ensure time-invariance within each cluster.

3 TICCAlgorithm

Problem (1) is a mixed combinatorial and continuous optimization problem with two sets of variables, the cluster assignments P and inverse covariances Θ, coupled together to make the problem highly non-convex. As such, there is no tractable way to solve for the globally optimal solution. Instead, we use a variation of expectation maximization (EM) to alternate between assigning points to clusters and updating the cluster parameters.

Cluster Assignment. Given the model parameters (i.e., inverse covariances) for each of the K clusters, we assign each of the T subsequences, X1, . . . ,XT , to these K clusters in such a way that maximizes the likelihood of the data while also minimizing the number of times that the cluster assignment changes across the time series. This combinatorial optimization problem has KT possible assignments of points to clusters. However, we are able to solve for the globally optimal solution in only O(KT) operations. We do so through dynamic programming, since this is equivalent to finding the minimum cost Viterbi path [Viterbi, 1967] for the length-T sequence.

Solving the Toeplitz Graphical Lasso. Once we have the clustering assignments, we then update the inverse covariances, given the points assignments. Here, we solve for each Θi in parallel, as

where |Pi| is the number of points assigned to cluster i and Si is the empirical covariance of these points. We call this problem the Toeplitz graphical lasso, since it is a variation on the well-known graphical lasso problem [Friedman et al., 2008] where we add a block Toeplitz constraint on the inverse covariance. To solve it, we develop a scalable algorithm based on a distributed convex optimization approach know as the alternating direction method of multipliers (ADMM) [Boyd et al., 2011; Parikh and Boyd, 2014].

TICC. Our overall TICC algorithm iterates between assigning points to different clusters (E-step, dynamic programming) and updating the cluster parameters (M-step, ADMM), as outlined in Algorithm (1).

4 Experiments

We built a custom Python solver to run the TICC algorithm 1 . Our solver takes as inputs the original multivariate time series and the problem parameters. It then returns the clustering assignments of each point in the time series, along with the structural MRF representation of each cluster. We test TICC on several synthetic examples. We do so because there are known "ground truth" clusters to evaluate the accuracy of our method.

Generating the Datasets. We randomly generate synthetic multivariate data in R5 . The overall time series is generated by constructing a temporal sequence of cluster segments (for example, the sequence "1, 2, 1" with 200 samples in each of the three segments, coming from two inverse covariances Θ1 and Θ2). The data is then drawn one sample at a time, conditioned on the values of the previous w-1 samples. We run our experiments on four different temporal sequences. Each segment in each of the examples has 100K observations in R5 , where K is the number of clusters in that experiment. These examples were selected to convey various types of temporal sequences over various lengths of time.

1 Available at http://snap.stanford.edu/ticc/.

Clustering Accuracy. Since both TICC and the baseline approaches use very similar methods for selecting the appropriate number of clusters, we fix K to be the "true" number of clusters, for both TICC and for all the baselines. We measure the macro-F1 score for the four different temporal sequences in Table 1. Here, all eight methods are using the exact same synthetic data, to isolate each approach's effect on performance. As shown, TICC significantly outperforms the baselines. Our method achieves a F1 score between 0.90 and 0.98, averaging 0.95 across the four examples. This is 41% higher than the second best method (not counting TICC, β = 0).

Effect of the Total Number of Samples. We next focus on how many samples are required for each method to accurately cluster the time series. We take the "1,2,3,4,1,2,3,4" example from Table 1 and vary the number of samples. We plot the F1 score vs. number of samples per segment in Figure 2. As shown, when there are 100 samples, none of the methods are able to accurately cluster the data. However, as we observe more samples, both TICC and TICC, β = 0 improve rapidly. By the time there are 200 samples, TICC already has an F1 score above 0.9. Even when there is a limited amount of data, TICC is still able to accurately cluster the data. As the number of samples increases, TICC's F1 score goes to 1.0, but no other method (including TICC, β = 0) tops 0.9. We note that even with 500 samples, the distance-based methods struggle to identify the clusters.

Scalability of TICC. One iteration of the TICC algorithm consists of running the dynamic programming algorithm and then solving the Toeplitz graphical lasso problem for each cluster. These steps are repeated until convergence. The total number of iterations depends on the data, but typically is no more than a few tens of iterations. To evaluate the scalability of our algorithm, we vary the number of timestamps and compute the runtime of the algorithm over one iteration. We observe samples in R50 , estimate 5 clusters with a window size of 3, and vary T over several orders of magnitude. We plot the results in Figure 3. Note that our ADMM solver is independent of T , so this contributes to the constant offset in the plot. However, for large values of T , our algorithm scales linearly with the number of points. Our TICC solver can cluster 10 millions points, each in R50 , with a per-iteration runtime of approximately 25 minutes.

5 Case Study

Here, we apply our TICC method to a real-world example to demonstrate how this approach can be used to find meaningful insights from time series data in an unsupervised way. We analyze a dataset, provided by a large automobile company, containing sensor data from a real driving session. This session lasts for exactly 1 hour and occurs on real roads in the suburbs of a large European city. We observe 7 sensors every 0.1 seconds:

• Brake Pedal Position
• Forward (X-)Acceleration
• Lateral (Y-)Acceleration
• SteeringWheel Angle

We run TICC and segment the time series into 5 clusters. We then analyze the clusters to understand and interpret what "driving state" they each refer to. Each cluster has a multilayer MRF network defining its structure. To analyze the result, we use network analytics to determine the relative "importance" of each node in the cluster's network. We plot the betweenness centrality score [Brandes, 2001] of each node in Table 2. We see that each of the 5 clusters has a unique "signature", and that different sensors have different betweenness scores in each cluster. For example, the Y-Acceleration sensor has a non-zero score in only two of the five clusters: #2 and #5. Therefore, we expect these two clusters to refer to

• Vehicle Velocity
• Engine RPM
• Gas Pedal Position

states in which the car is turning, and the other three to refer to intervals where the car is going straight. As such, we can use these betweenness scores to interpret these clusters in a meaningful way. For example, from Table 2, a reasonable hypothesis might be that the clusters refer to 1) slowing down, 2) turning, 3) speeding up, 4) cruising straight, 5) driving on a curvy road segment.

Plotting the Resulting Clusters. To validate our hypotheses, we can plot the latitude/longitude locations of the drive, along with the resulting cluster assignments. Analyzing this data, we empirically discover that each of the five clusters has a clear real-world interpretation that aligns very closely with our estimates based on the betweenness scores in Table 2. Furthermore, we notice that many consistent and repeated patterns emerge in this one hour session. For example, whenever the driver is approaching a turn, he or she follows the same sequence of clusters: going straight, slowing down, turning, speeding up, then going straight again. We plot two typical turns in the dataset, coloring the timestamps according to their cluster assignments, in Figure 4.

6 Conclusion and FutureWork

In this paper, we have defined a method of clustering multivariate time series subsequences. Our method, Toeplitz Inverse Covariance-based Clustering (TICC), simultaneously segments and clusters the data, breaking down highdimensional time series into a clear sequential timeline. TICC's promising results on both synthetic and real-world data lead to many potential directions for future research. For example, our method could be extended to learn dependencies parameterized by any heterogeneous exponential family MRF. This would allow for a much broader class of datasets (such as boolean or categorical readings) to be incorporated into TICC.

Note

This is an abridged version of the full paper, originally published in ACM SIGKDD, 2017.