A Nested Two-Stage Clustering Method for Structured Temporal Sequence Data

Abstract

Mining patterns of temporal sequence data is an important problem across many disciplines. Under appropriate preprocessing procedures, a structured temporal sequence can be organized into a probability measure or a time series representation, which grants a potential to reveal distinctive temporal pattern characteristics. In this paper, we propose a nested two-stage clustering method that integrates optimal transport and the dynamic time warping distances to learn the distributional and dynamic shape-based dissimilarity at the respective stage. The proposed clustering algorithm preserves both the distribution and shape patterns present in the data, which are critical for the datasets composed of structured temporal sequences. The effectiveness of the method is tested against existing agglomerative and K-shape-based clustering algorithms on Monte Carlo simulated synthetic datasets, and the performance is compared through various cluster validation metrics. Furthermore, we apply the developed method to real-world datasets from three domains: temporal dietary records, online retail sales, and smart meter energy profiles. The expressiveness of the cluster and subcluster centroid patterns shows significant promise of our method for structured temporal sequence data mining.

Appendices

Appendix

Remark 1

The Wasserstein barycenter ${\overline{\varPhi }}_k$ of $n_k$ continuous distributions $\{\varPhi _1,\ldots ,\varPhi _{n_k}\}$ of cluster k under the objective of Definition (13) satisfies

$$\begin{aligned} {\overline{\varPhi }}_k^{-1}(w) = \frac{1}{n_k} \sum _{i:g(i)=k} \varPhi _i^{-1}(w),\quad \forall ~ w \in [0,1]. \end{aligned}$$

(12)

Remark 2

(Wasserstein Barycenter, [2]) A Wasserstein barycenter of N measures $\{\nu _i: i=1,\ldots ,N\}$ in ${\mathbb {P}} \subset P(\varOmega )$ is a minimizer of f over ${\mathbb {P}}$, where

$$\begin{aligned} \mu ^* :=\mathrm{arg\,min}_{\mu } f(\mu ) = \mathrm{arg\,min}_{\mu } \sum _{i=1}^N \lambda _i W_2^2(\mu ,\nu _i). \end{aligned}$$

(13)

Fig. 18

Cluster validity index DB and CH for experiments with K ranging from 2 to 10, to determine the optimal choice of K for the first-stage OT-means on temporal dietary dataset

Fig. 19

Cluster validity index DB and CH for experiments with K ranging from 2 to 10, to determine the optimal choice of K for the first-stage OT-means on Online Retail Dataset

Fig. 20

Mean energy ratio weighted average time for each cluster of Online Retail Dataset. Distinct first-stage clusters are colored in distinct colors

Fig. 21

Two-stage OT–DTW cluster output chart for Online Retail Dataset. The piechart in the middle shows the sample size of the six cluster outputs of the first stage. The six bar charts show the subcluster output sample size distribution of the second stage

Fig. 22

Mean energy ratio weighted average time for each cluster of Smart Meter Energy Consumption Dataset. Distinct first-stage clusters are colored in distinct colors

Fig. 23

Two-stage OT–DTW cluster output chart for Smart Meter Energy Consumption Dataset. The piechart in the middle shows the sample size of the five cluster outputs of the first stage. The five bar charts show the subcluster output sample size distribution of the second stage

Fig. 24

The left is a 15 oscillator network with shown connection topology, and phase difference-based synchronization clusters are colored in different colors. The right is an illustration of the cluster result using out OT–DTW method

Remark 3

(DTW Barycenter) A DTW barycenter of N time series $P=\{{\mathbf {p}}_1, \ldots ,{\mathbf {p}}_N\}$ in a space ${\mathbb {E}}$ induced by DTW metric is a minimizer of the sum of squared distance to the set P, where

$$\begin{aligned} \eta ^* :=\mathrm{arg\,min}_{\eta } \frac{1}{N} \sum _{i=1}^N E^2(\eta ,{\mathbf {p}}_i). \end{aligned}$$

(14)

A. Results

Based on the definition of DB and CH indices, we seek to find the local minimum of DB index and the local maximum of CH index. From Fig. 18b, the CH index strictly decreases with increasing K and there is no clear kink point toward plateau, which provides little information for the optimal choice of K. From Fig. 18a, due to the relative smaller DB index and clearer separation of cluster centroids, we set $K=4$ in the current experiment (Example of temporal dietary dataset). From Fig. 19b, the CH index also strictly decreases with increasing K and provides little information for the optimal choice of K. But from Fig. 19a, $K=6$ becomes a good candidate for the number of clusters since the DB index achieves local minimum then (Figs. 20, 21, 22, 23).

B. Applications

Apart from the temporal pattern discovery in the applications discussed in Sect. 6, the proposed clustering algorithm appears to posses some desirable properties which would extend its use in synchronization detection application in an oscillator network [23]. The synchronization detection problem is defined as follows: In an oscillator network, each oscillator can be treated as a node in the network, and the coupling between oscillators is the edges. Each oscillator’s dynamics consists of two parts—its own intrinsic dynamics and the coupling functions from other oscillators. The network starts from an arbitrary initial condition and evolves over time (according to the oscillator dynamical equations). Given the time series measurement corresponding to the output of each oscillator, we aim to determine which of the oscillators (nodes) are phase synchronized. Traditionally, this problem requires preprocessing of the data by peak-finding or Hilbert transform (to extract phase information from the measured data) and further clustering according to the oscillator phase model [23]. Our method saves the expensive phase processing step, and can directly work with the recordings. For example, Fig. 24 shows an illustration of a synthetic oscillator network with 15 oscillators and cluster results from our OT–DTW method. The colored nodes in the left network plot provide the synchronization clusters based on phase difference calculation. On the right is our two-stage cluster outputs, and except oscillator 14, our cluster results match very well with the phase-based synchronization clusters (our results also separate oscillator 7, 12, and 13 into a separate cluster from oscillator 2, 3, 6, and 8). This leads to our conjecture that the distributional difference and the dynamic shape difference in the time domain have some intrinsic correlation with the phase synchronization and we plan to pursue this direction in a future study.