Coresets for Time Series Clustering

Abstract:
We study the problem of constructing coresets for clustering problems with time series data. This problem has gained importance across many fields including biology, medicine, and economics due to the proliferation of sensors for real-time measurement and rapid drop in storage costs. In particular, we consider the setting where the time series data on N entities is generated from a Gaussian mixture model with autocorrelations over k clusters in R. Our main contribution is an algorithm to construct coresets for the maximum likelihood objective for this mixture model. Our algorithm is efficient, and, under a mild assumption on the covariance matrices of the Gaussians, the size of the coreset is independent of the number of entities N and the number of observations for each entity, and depends only polynomially on k, d and 1/ε, where ε is the error parameter. We empirically assess the performance of our coresets with synthetic data.

heading: 1 Introduction
text: A multivariate time series dataset, represented as X = { Xi = (xi1, . . . , xi,Ti) ⊂ RTi×d | i ∈ [N ] } , where N is the number of entities, Ti is the number of time periods corresponding to entity i and d is the number of features, tracks features of a cross-section of entities longitudinally over time. Such data is also referred to as panel data [8] and has seen rapid growth due to proliferation of sensors, IOT and wearables that facilitate real time measurement of various features associated with entities and the rapid drop in storage costs. Specific examples of time series data include biomedical measurements (e.g., blood pressure and electrocardiogram), epidemiology and diffusion through social networks, weather, consumer search, content browsing and purchase behaviors, mobility through cell phone and GPS locations, stock prices and exchange rates in finance [2].
Computational problems of interest with time series data include pattern discovery [49], regression/prediction [39], forecasting [9], and clustering [47, 31, 2] which arises in applications such as anomaly detection [46]. Though clustering is a central and well-studied problem in unsupervised learning, most standard treatments of clustering tend to be on static data with one observation per entity [35, 11]. Time series clustering introduces a number of additional challenges and is an active area of research, with many types of clustering methods proposed [15]; direct methods on raw data, indirect methods based on features generated from the raw data, and model based clustering where the data are assumed to be generated from a model. For surveys of time series clustering, see [47, 2].
We focus on model-based time series clustering using a likelihood framework that naturally extends clustering with static data [11] to time series data [47], where each multivariate time series is generated by one of k different parametrically specified models. For a survey of the importance of this subliterature and its various real-world applications, see [30]. A prevalent approach is to assume a finite mixture of data generating models, with Gaussian mixture models being a common specification [50, 53]. Roughly, the problem is to partition X probabilistically into k clusters, grouping those series
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
generated by the same time series model into one cluster. Generically, the model-based k-clustering problem can be formulated as: arg maxα,θ(1),...,θ(k) ∑ i∈[N ] ln ∑ l∈[k] αl ·pi(X|θ(l)), where αl is the probability of a time series belonging to cluster l and pi(X|θ(l)) is the likelihood of the data of entity i, given that the data is generated from model l (with parameters θ(l)). Temporal relationships in time series are commonly modeled as autocorrelations (AR)/moving averages (MA) [20, 61] that account for correlations across observations, and hidden Markov models (HMM) where the underlying data generating process is allowed to switch periodically [21, 62]. This paper focuses on the case when the data generation model for each segment l ∈ [k] is a multivariate Gaussian with autocorrelations. More formally, the generative model for each cluster is from a Gaussian mixture: xit := µ(l) + eit, where eit = N(0,Σ
(l)) + Λ(l)ei,t−1, where N(0,Σ(l)) captures the mixture of Gaussian distributions from which entity level observations are drawn, and Λ captures the correlation between two successive observations through an AR(1) process [32, 44]. Overall, this model can be represented with cluster level model parameters θ(l) = {µ(l),Σ(l),Λ(l)} and cluster probability αl. Time series datasets are also much larger than static datasets [47]. For instance, as noted in [42], ECG (electrocardiogram) data requires 1 GB/hour, a typical weblog requires 5 GB/week. Such large sizes lead to high storage costs and also make it necessary to work with a subset of the data to conduct the clustering analysis to be computationally practical. Further, long time series data on entities also entails significant data security and privacy risks around the entity as it requires longer histories of entity behaviors to be stored and maintained [25, 41]. Thus, there is significant value for techniques that can produce approximately similar clustering results with smaller samples as with complete data.
Coresets have emerged as an effective tool to both speed up and reduce data storage by taking into account the objective and carefully sampling from the full dataset in a way that any algorithm run on the sampled set returns an approximate solution for the original objective with guaranteed bounds on approximation error [34]. Coresets have been developed for both unsupervised clustering (e.g., k-means, Gaussian mixture models) and supervised machine learning (e.g., regression) methods for static data; for surveys on coresets for static data, see [7, 22]. A natural question is whether coresets can be also useful for addressing the time series clustering problem. Recently, there is evidence that small coresets can be efficiently constructed for time series data, albeit for regressions with time series (panel) data [37]. However, as we explain in Section 2, there are challenges in extending coresets for static data clustering and time series (panel) data regressions for time series clustering.
Our contributions. We study coresets for a general class of time series clustering in which all entities are drawn from Gaussian mixtures with autocorrelations (as described above; see Problem 3.1). We first present a definition of coresets for the log-likelihood k-clustering objective for this mixture model. One issue is that this objective cannot be decomposed into a summation of entity-time objectives due to the interaction of the ln term from log-likelihood and the exp term from Gaussians (Problem 3.1). To estimate this objective using a coreset, we introduce an analogous clustering objective on a subset of data (Definition 3.3). Our main result is an algorithm to construct coresets for this objective for the aforementioned mixture model (Theorem 4.3). Our algorithm is efficient (linear on both N and ∑ i∈[N ] Ti) and, assuming that the condition number of the covariance matrices of the Gaussians (Σ(l)) are bounded and the eigenvalues of the autocorrelation matrices (Λ(l)) are bounded away from 1 (Assumption 4.1), the size of our coreset is independent of the number of entities N and the number of observations Tis; it only polynomially depends on k, d and 1/ε, where ε is the error parameter (Theorem 4.3).
Our coreset construction (Algorithm 1) leverages the Feldman-Langberg (FL) framework [27]. Due to multiple observations for each entity and associated autocorrelations, the objective is non-convex and quite complicated in contrast to the static setting where the objective only contains one observation drawn from Gaussian mixtures for each entity. Thus, bounding the “sensitivity” (Lemma 4.6) of each entity, which is a key step in coreset construction using the FL-framework, becomes challenging. We handle this technical difficulty by 1) upper-bounding the maximum effect of covariances and autocorrelations by using the observation that the gap between the clustering objectives with and without Σ and Λ is always constant, and 2) reducing the Gaussian mixture time series clustering problem to a certain k-means clustering problem whose total sensitivity is upper bounded.
Empirically, we assess the performance of our coreset on synthetic data for a three cluster problem (Section 5), and compare the performance with two benchmarks: uniform sampling and a static coreset benchmark (LFKF [48]) in which we do not consider autocorrelations and regard time series
data as static. We find that our coreset performs better relative to uniform sampling and LFKF on both data storage and computation speed for a range of accuracy guarantees: To achieve a similar fit with the full data, our coreset needs fewer entity-time observations (<40%); the computation time for a given level of accuracy reduces by a 3x-22x factor when compared to uniform sampling and LFKF. Moreover, our coreset speeds up the computation time relative to the full data by 14x-171x. We note that the performance advantage of our coreset is greater when there is more time series data relative to entities (Ti N ).

heading: 2 Related work
text: Time series or panel data analysis are important in many fields– biological [18], engineering [3], economics [51] and the social sciences [59]. Clustering has been a central problem in unsupervised learning and also has a long history of use across many fields; for a historical overview of applications across fields, see [12]. While most standard treatments of clustering in machine learning tend to be on static data [35, 11], there is by now a significant and rapidly growing literature on time series clustering; recent surveys include [47, 2]. While there are many coresets for algorithms on static data (including for clustering), there has been little work on coresets on time series data (see below).
Coresets for clustering. Coresets for various clustering objectives on static data have been well studied, including k-median [34, 17, 55, 38, 19], k-means [34, 17, 28, 13, 10, 55, 38, 19], and k-center [1, 33]. For surveys of static coresets, see [52, 23]. [48, 26] present coreset constructions for clustering problems on static data generated using Gaussian mixture models. [48] constructs a coreset of size poly(k, d, 1/ε) under the same boundedness assumption as ours. [26] remove this assumption and construct a coreset whose size additionally depends on poly logN , but require another assumption that all xis are integral and in a bounded range. In contrast, we consider a generalized problem where each entity has multiple observations over time and accounts for autocorrelations across successive observations. The new idea is to reduce the GMM time series clustering problem (Problem 3.1) to a k-means clustering problem whose points are linear combinations of those in Xi (Lemma A.12), and show that the effect of autocorrelation parameters Λ on the entity objectives can be upper bounded (Lemma A.11).
Regression with time series data. [37] construct coresets for a regression problem, called GLSEk, for time-series data including autocorrelations and a k-partition of entities. This is a supervised learning task that includes an additional label yit for each entity-time pair. It aims to find a kpartition of entities and additional regression vectors β(l) ∈ Rd for each partition l ∈ [k] such that a certain linear combination of terms ‖yit − (β(l))>xit‖22 is minimized. In contrast, our problem is an unsupervised learning task, estimates points xit by Gaussian means µ instead of β, and includes covariances Σ between different features. The definitions of coresets in [37] and ours consist of entity-time pairs but our coresets need to include multiple weight functions (w for entities and w(i) for time periods of selected entity i; see Definition 3.3) instead of one weight function as in the regression problem. This is to deal with the interaction of the ln term from the log-likelihood and the exp term from the Gaussians. [37]’s result relies on an assumption that the regression objective of each entity is in a bounded range over the parameter space. This may not be satisfied in our setting since the clustering objective of an entity may be unbounded when all Gaussian means µ(l) are far away from observations of this entity. To bypass this, we reduce our problem to a k-means clustering problem (Definition 4.2), whose total sensitivity is provably O(k) (Lemma A.12), by upper bounding the affect of covariance and autocorrelation parameters on the clustering objectives of entities (Lemma A.11) under the assumption that the condition number of covariance matrix is bounded (Assumption 4.1), which provides an upper bound for the total sensitivity of entities on our problem (Lemma 4.6).
Another related direction is to consider coreset construction for k-segmentation with time series data [54, 29], which aims to estimate the trace of an entity by a k-piecewise linear function (ksegment). Note that the case of k = 1 is equivalent to the linear regression problem. [54, 29] proposed efficient coreset construction algorithms which accelerate the computation time of the k-segmentation problem. The main difference from our setting is that they consider a single entity observed at multiple time periods, and the objective is an additive function on all time periods. This enables them to relate the k-segmentation problem to the static setting. It is interesting to investigate coreset construction for more optimization problems with time series data.

heading: 3 Clustering model and coreset definition
text: Given X = { Xi = (xi,1, . . . , xi,Ti) ∈ Rd×Ti | i ∈ [N ] } , we first model a general class of time series clustering problems. Then we specify our setting with Gaussian mixture time series data (Problem 3.1), and define coresets accordingly (Definition 3.3).
Clustering with time series data. Given an integer k ≥ 1, let ∆k ⊂ Rk denote a probability simplex satisfying that for any α ∈ ∆k, we have αi ∈ [0, 1] for i ∈ [k] and ∑ l∈[k] αi = 1. Let P denote a parameter space where each θ ∈ P represents a specific generative model of time series data. For each entity i ∈ [N ] and model θ ∈ P , define pi(X | θ) := Pr [Xi | θ]1/Ti to be the average likelihood/realized probability of Xi from model θ. The ratio 1/Ti is used to normalize the contribution to the objective of entity i due to different lengths Ti. The time series clustering problem is to compute α ∈ ∆k and θ = (θ(1), . . . , θ(k)) ∈ Pk that minimizes the negative log-likelihood, i.e., − ∑ i∈[N ] ln ∑ l∈[k] αl · pi(X | θ(l)). Here, for each l ∈ [k], αl represents the probability that each time series is generated from model θ(l). Note that the clustering objective depends on the choice of model family P . In this paper, we consider the following specific model family P that each time series is generated from Gaussian mixtures.
GMM clustering with time series data. Let α ∈ ∆k be a given probability vector. Let Sd denote the collection of all symmetric positive definite matrices in Rd×d. For each i ∈ [N ], with probability αl (l ∈ [k]), let xit := µ(l) + eit for each t ∈ [Ti] where µ(l) ∈ Rd represents a Gaussian mean and eit ∈ Rd represents the error vector drawn from the following normal distribution: eit := Λ (l)ei,t−1 +N(0,Σ (l)), where Λ(l) ∈ Sd is an AR(1) autocorrelation matrix and Σ(l) ∈ Sd is the covariance matrix of a multivariate Gaussian distribution. Now we let P := Rd × Sd × Sd and note that each Gaussian generative model can be represented by a tuple θ = (µ,Σ,Λ) ∈ P . Moreover, we have that the realized probability of each entity i ∈ [N ] is
pi(X | θ) = pi(X | µ,Σ,Λ) := exp(− 12Ti ψi(µ,Σ,Λ))
(2π)d/2|Σ|1/2 , (1) where |Σ| is the determinant of Σ and ψi(µ,Σ,Λ) := ∑ t∈[Ti] ψit(Σ,Λ) with
ψi1(Σ,Λ) := (xi,1 − µ)>Σ−1(xi,1 − µ)− (Λ(xi,1 − µ))>Σ−1 (Λ(xi,1 − µ)) , and
ψit(Σ,Λ) := ((xit − µ)− Λ(xi,t−1 − µ))>Σ−1 ((xit − µ)− Λ(xi,t−1 − µ)) ,∀ 2 ≤ t ≤ Ti. We note that each sub-function ψit contains at most two entity-time pairs: xit and xi,t−1. Our Gaussian mixture time series model gives rise in the following clustering problem. Problem 3.1 (Clustering with GMM time series data). Given a time series dataset X ={ Xi = (xi,1, . . . , xi,Ti) ∈ Rd×Ti | i ∈ [N ] } and integer k ≥ 1, the GMM time series clustering problem is to compute α ∈ ∆k and θ = (µ(l),Σ(l),Λ(l))l∈[k] ∈ Pk that minimize f(α, θ) := ∑ i∈[N ] fi(α, θ) = − ∑ i∈[N ] ln ∑ l∈[k] αl · pi(X | µ(l),Σ(l),Λ(l)).
From Equation (1) it follows that the coefficient of each Gaussian before the exp operation is αl (2π)d/2|Σ(l)|1/2 , whose summation Z(α, θ) := ∑ l∈[k] αl (2π)d/2|Σ(l)|1/2 is not a fixed value and depends on Σ(l)s. To remove this dependence on the coefficients, we define α′l(θ) := αl (2π)d/2|Σ(l)|1/2Z(α,θ) to be the normalized coefficient for l ∈ [k], define offset function φ : ∆k ×Pk → R to be φ(α, θ) := −N lnZ(α, θ), and define f ′(α, θ) := − ∑ i∈[N ] ln ∑ l∈[k] αl ·exp(− 1 2Ti ψi(µ (l),Σ(l),Λ(l))) whose summation of coefficients before the exp operations is 1. This idea of introducing f ′ also appears in [24, 26], which is useful for coreset construction. This leads to the following observation. Observation 3.2. For any α ∈ ∆k and θ ∈ Pk, f(α, θ) = f ′(α′(θ), θ) + φ(α, θ).
Our coreset definition. The clustering objective f can only be decomposed into the summation of fis instead of sub-functions w.r.t. entity-time pairs. A simple idea is to select a collection of weighted entities as a coreset. However, in doing so, the coreset size will depend on Tis and fail to be independent of both N and Tis. To address this problem, we define our coreset as a collection of weighted entity-time pairs, which is similar to [37]. Let PX := {(i, t) | i ∈ [N ], t ∈ [Ti]} denote the collection of indices of all xit. Given S ⊆ PX , we let IS := {i ∈ [N ] | ∃t ∈ [Ti], s.t., (i, t) ∈ S} denote the collection of entities that appear in S. Moreover, for each i ∈ IS , we let JS,i := {t ∈ [Ti] : (i, t) ∈ S} denote the collection of observations for entity i in S.
Definition 3.3 (Coresets for GMM time series clustering). Given a time series dataset X ={ Xi = (xi,1, . . . , xi,Ti) ∈ Rd×Ti | i ∈ [N ] } , constant ε ∈ (0, 1), integer k ≥ 1, and parameter space ∆k ×Pk, an ε-coreset for GMM time series clustering is a weighted set S ⊆ PX together with weight functions w : IS → R≥0 and w(i) : JS,i → R≥0 for i ∈ IS such that for any α ∈ ∆k and θ = (µ(l),Σ(l),Λ(l))l∈[k] ∈ Pk, f ′S(α, θ) := − ∑ i∈IS w(i) · ln ∑ l∈[k] αl · exp(− 1 2Ti ∑ t∈JS,i w
(i)(t) · ψit(µ (l),Σ(l)),Λ(l))) ∈ (1± ε) · f ′(α, θ).
Combining the above definition with Observation 3.2, we note that for any α ∈ ∆k and θ ∈ Pk,
f ′S(α ′(θ), θ) + φ(α, θ) ∈ (1± ε) · f ′(α, θ) + φ(α, θ) ∈ (1± ε) · f(α, θ)± 2εφ(α, θ). (2)
As ε tends to 0, f ′S(α ′(θ), θ) + φ(α, θ) converges to f(α, θ). Moreover, given such a coreset S, if we additionally have that φ(α, θ) . f(α, θ), we can minimize f ′S(α ′(θ), θ) + φ(α, θ) to solve Problem 3.1. Thus, if φ(α, θ) . f(α, θ), we conclude that f ′S(α ′(θ), θ)+φ(α, θ) ∈ (1±3ε)·f(α, θ). Note that we introduce multiple weight functions (w for entities and w(i) for time periods of selected entity i) in Definition 3.3 unlike the coreset for time series regression in [37] which uses only one weight function. This is because f ′i contains ln and exp operators instead of linear combinations of ψit; it is unclear how to use a single weight function to capture both the entity and time levels.

heading: 4 Theoretical results
text: We present our coreset algorithm (Algorithm 1) and the main result on its performance (Theorem 4.3). Theorem 4.3 needs the following assumptions on covariance and autocorrelation parameters.
Assumption 4.1. Assume (1) that there exists constant D ≥ 1 such that maxl∈[k] λmax(Σ (l))
minl∈[k] λmin(Σ(l)) ≤ D
where λmax(·) represents the largest eigenvalue and λmin(·) represents the smallest eigenvalue, and (2) for l ∈ [k], Λ(l) ∈ Ddλ for some constant λ ∈ (0, 1). Here Ddλ is a collection of all diagonal matrix in Rd×d whose diagonal elements are at most 1− √ λ.
The first assumption requires that the condition number of each covariance matrix is upper-bounded, which also appears in [60, 45, 48] that consider Gaussian mixture models with static data. The second assumption, roughly, requires that there exist autocorrelations only between the same features. The upper bound for diagonal elements ensures that the autocorrelation degrades as time period increases, which is also assumed in [37]. Note that both the eigenvalues of Σ(i) and the positions of means µ(l)s affect cost function ψis. For instance, consider the case that all eigenvalues are 1 and all autocorrelations are 0, i.e., for all l ∈ [k], Σ(l) = Id and Λ(l) = 0d. In this case, it is easy to see that max µ∈Rd ψi(µ,Id,0d)
min µ∈Rd ψi(µ,Id,0d) =∞, i.e., the value of ψi(µ, Id, 0d) is unbounded as µ changes. Differently, the component’s cost functions are bounded in [37] since µ(l)s do not appear in [37] that consider regression problems.
Let Pkλ := (Rd × Sd ×Ddλ)k denote the parameter space. For preparation, we propose the following k-means clustering problem. Definition 4.2 (k-means clustering of entities). Given an integer k ≥ 1, the goal of the k-means clustering problem of entities of X is to find a set C? = {c?1, . . . , c?k} ⊂ Rd of k centers that minimizes ∑ i∈[N ] minl∈[k] ‖ ∑ t∈[Ti] xit Ti − cl‖22 over all k center sets C = {c1, . . . , ck} ⊂ Rd. Let OPT(O) denote the optimal k-means value of C?.
Note that there exists an O(2poly(k)Nd) time algorithm to compute a nearly optimal solution for both OPT(O) and C? [43]. Another widely used algorithm is called k-means++ [6], which provides an O(ln k)-approximation in O(Ndk lnN ln k) time but performs well in practice. We also observe
that ψi(µ, Id, 0d) = ‖ ∑ t∈[Ti] xit
Ti − µ‖22 + 1Ti ∑ t∈[Ti] ‖xit‖ 2 2 − ‖ ∑ t∈[Ti] xit‖22 T 2i
for any µ ∈ Rd. This observation motivates us to consider a reduction from Problem 3.1 to Definition 4.2, which is useful for constructing IS .
Our coreset algorithm. We first give a summary of our algorithm (Algorithm 1). The input to Algorithm 1 is a time series dataset X , an error parameter of coreset ε ∈ (0, 1), an integer k ≥ 1 of
clusters, λ ∈ (0, 1) (vector norm bound), and a covariance eigenvalue gap D ≥ 1. We develop a two-staged importance sampling framework which first samples a subset IS of entities (Lines 1-8) and then samples a subset of time periods JS,i for each selected entity i ∈ IS (Lines 9-14). In the first stage, we first set M as the number of selected entities |IS | (Line 1). Then we solve the k-means clustering problem over entity means bi (Definition 4.2), e.g., by k-means++ [6] (Lines 2-3), and obtain an (approximate) optimal clustering value OPT(O), a k-center set C? := {c?1, . . . , c?k} ⊂ Rd, and an offset value ai for each i ∈ [N ]. Next, based on the distances of points bi to C?, we partition [N ] into k clusters where entity i belongs to cluster c?p(i) (Line 4). Based on bi, c ? p(i) and OPT(O), we compute an upper bound s(i) for the sensitivity w.r.t. entity i ∈ [N ] (Lines 5-6). Finally, we sample a weighted subset IS of entities as our coreset for the entity level by importance sampling based on s(i) (Lines 7-8), which follows from the Feldman-Langberg framework (Theorem A.6).
In the second stage, we first set L as the number of time periods |JS,i| for each selected entity i ∈ IS (Line 9). Then for each selected entity i ∈ IS , we compute OPT(O)i as the optimal 1-means clustering objective of xits (Line 10). Next, we compute an upper bound si(t) for the sensitivity w.r.t. time period t ∈ [Ti] (Lines 11-12) based on bi and OPT(O)i . Finally, we sample a weighted subset JS,i of time periods by importance sampling based on si(t) (Lines 13-14).
Algorithm 1: CRGMM: Coreset construction for GMM time series clustering Input: X = { Xi = (xi,1, . . . , xi,Ti) ∈ Rd×Ti | i ∈ [N ] } , constant ε, λ ∈ (0, 1), number of clusters k ≥ 1,
dimension of autocorrelation vectors q ≥ 0, constant of the variance gap D ≥ 1, and parameter space Pkλ = (Rd × Sd ×Ddλ)k.
Output: a subset S ⊆ P together with weight functions w : IS → R≥0 and w(i) : JS,i → R≥0 for i ∈ IS . % Constructing a subset of entities
1: M ← O ( kD λε2 (k4d4 + k3d8) ln k λ ) .
2: For each i ∈ [N ], compute point bi ∈ Rd by bi ← ∑ t∈[Ti] xit
Ti and value ai ≥ 0 by ai ← 1Ti ∑ t∈[Ti] ‖xit‖ 2 2 − ‖ ∑ t∈[Ti] xit‖22 T2i . Let A← ∑ i∈[N ] ai.
3: Compute (an O(1)-approximate for) OPT(O) and C? = {c?1, . . . , c?k} ⊂ Rd for the k-means clustering problem over points bi (Definition 4.2),1by [43, 6]. 4: For each i ∈ [N ], compute index p(i) ∈ [k] by p(i)← arg minl∈[k] ‖bi − c?l ‖22. 5: For each i ∈ [N ], sc(i)← 1
| { i′∈[N ]:p(i′)=c?
p(i)
} | .
6: For each i ∈ [N ], s(i)← min { 1, 4D ( 4‖bi−c?p(i)‖ 2 2
OPT(O)+A + 3sc(i)
) /λ } .
7: Pick a random sample IS ⊆ [N ] of size M , where each i ∈ IS is selected w.p. s(i)∑ i′∈[N] s(i ′)
8: For each i ∈ IS , w(i)← ∑ i′∈[N] s(i ′)
M·s(i) . % Constructing a subset of time periods for each selected entity
9: L← O ( Dd8 ln 1
λ λε2
) .
10: For each i ∈ IS , compute OPT(O)i ← ∑ t∈[Ti] ‖xit − bi‖ 2 2. 11: For each i ∈ IS and t ∈ [Ti], sci (t)← 2‖xit−bi‖22
OPT (O) i
+ 6 Ti . 12: For each i ∈ IS and t ∈ [Ti], si(t)← min { 1, 4Dλ−1 ( sci (t) + ∑min{t−1,1} j=1 s c i (t− j) )} . 13: For each i ∈ IS , pick a random sample JS,i of L points, where each t ∈ JS,i is selected w.p. si(t)∑
t′∈[Ti] si(t′)
.
14: For each i ∈ IS and t ∈ JS,i, w(i)(t)← ∑ t′∈[Ti] si(t ′)
L·si(t) .
% Output coreset S of entity-time pairs 15: Let S ← {(i, t) ∈ PX : i ∈ IS , t ∈ JS,i}. 16: Output ( S,w, { w(i) } i∈IS ) .
1Here, we slightly abuse the notation by using OPT(O) and C? to represent the obtained approximate solution.
Our main theorem. Our main theorem is as follows, which indicates that Algorithm 1 provides an efficient coreset construction algorithm for GMM time series clustering. Theorem 4.3 (Main result). Under Assumption 4.1, with probability at least 0.9, Algorithm 1 outputs an ε-coreset of size O ( D2k5d16 ln2 kλ
λ2ε4
) for GMM time-series clustering in O(d ∑ i∈[N ] Ti +
Ndk lnN ln k) time.
Note that the success probability in the theorem can be made 1− δ for any δ ∈ (0, 1) at the expense of an additional factor log2 1/δ in the coreset size. The coreset size guaranteed by Theorem 4.3 has a polynomial dependence on factors k, d, 1/ε,D, and 1/λ and, importantly, does not depend on N or Tis. Compared to the clustering problem with GMM static data [48], the coreset has an additional dependence on 1/λ due to autocorrelation parameters. Compared to the regression problem (GLSEk) with time series data [37, Theorem 5.2], the coreset does not contain the factor “M” that upper bounds the gap between the maximum and the minimum entity objective. The construction time linearly depends on the total number of observations ∑ i∈[N ] Ti, which is efficient. The proof of Theorem 4.3 can be found in Section C.
Proof overview of Theorem 4.3. The proof consists of three parts: 1) Bounding the size of the coreset, 2) proving the approximation guarantee of the coreset, and 3) bounding the running time.
Parts 1 & 2: Coreset size and correctness guarantee of Theorem 4.3. We first note that the coreset size in Theorem 4.3 is |S| = ML. Here, M is the number of sampled entities |IS | that guarantee that IS is a coreset for f ′ at the entity level (Lemma 4.4). L is the number of sampled observations |JS,i| for each i ∈ IS that guarantees that JS,i is a coreset for ψi at the time level (Lemma 4.5). Lemma 4.4 (The 1st stage outputs an entity-level coreset). With probability at least 0.95, the output IS of Algorithm 1 with M := |IS | = O ( kD λε2 (k 4d4 + k3d8) ln kλ )
satisfies that for any α ∈ ∆k and θ ∈ Pkλ ,− ∑ i∈IS w(i)·ln ∑ l∈[k] αl ·exp(− 1 2Ti ψi(µ (l),Σ(l),Λ(l))) ∈ (1±ε)·f ′(α, θ). Lemma 4.5 (The 2nd stage outputs time level coresets for each i ∈ IS). With probability at least 0.95, for all i ∈ IS , the output JS,i of Algorithm 1 with L := |JS,i| = O ( Dd8 ln 1λ λε2 ) satisfies that for
any (µ,Σ,Λ) ∈ Pλ, ∑ t∈JS,i w (i)(t) · ψit(µ,Σ,Λ) ∈ (1± ε) · ψi(µ,Σ,Λ).
From Lemmas 4.4 and 4.5, we can conclude that S is anO(ε)-coreset for GMM time-series clustering, since the errors ε for the entity level and the time level are additively accumulated, and that the size of the coreset ML as claimed; see Section C for a complete proof.
Both lemmas rely on the Feldman-Langberg framework [27, 13] (Theorem A.6). The key is to upper bound both the “pseudo-dimension” dim (Definition A.5) that measures the complexity of the parameter space, and the total sensitivity G (Definition A.4) that measures the sum of maximum influences of all sub-functions (f ′i for Lemma 4.4 and ψit for Lemma 4.5). The coreset size guaranteed by the Feldman-Langberg framework then is O(ε−2G dim lnG). For Lemma 4.5, we can upper bound the pseudo-dimension by dim = O(d8) (Lemma B.1). For the total sensitivity, a key step is to verify that si (Line 12) is a sensitivity function of ψi. This step uses a similar idea as in [37] and leads to showing a G = O(D/λ) bound for the total sensitivity (Lemma B.2). Then we verify that L = O(ε−2G dim lnG) = O ( Dd8 ln 1λ λε2 ) is enough for Lemma 4.5 by the Feldman-Langberg framework (Theorem A.6). The proof can be found in Section B.
Lemma 4.4 is the most technical. The proof can be found in Section A and we provide a proof sketch.
Proof sketch of Lemma 4.4. By [4, 58], the pseudo-dimension is determined by the number of parameters that is upper bounded by O(kd2) and the number of operations on parameters that is upper bounded by O(kd6) including O(k) exponential functions. Then we can upper bound the pseudo-dimension by dim = O(k4d4 + k3d8) (Lemma A.9 in Section A.2).
For the total sensitivity, the key is to verify that s (Line 6) is a sensitivity function w.r.t. f ′; this is established by Lemma 4.6 below. Then the total sensitivity is at most G = (16D + 12Dk)/λ. This completes the proof of Lemma 4.4 since M = O(ε−2G dim lnG) = O ( kD λε2 (k 4d4 + k3d8) ln kλ )
is enough by the Feldman-Langberg framework (Theorem A.6). Lemma 4.6 (s is a sensitivity function w.r.t. f ′). For each i ∈ [N ], we have s(i) ≥ maxα∈∆k,θ∈Pkλ f ′i(α,θ) f ′(α,θ) . Moreover, ∑ i∈[N ] s(i) ≤ (16D + 12Dk)/λ.
Proof sketch of Lemma 4.6. We first define ψ(O)i (µ) := ∑ t∈[Ti] ‖xit−µ‖ 2 2 for any µ ∈ Rd, and define f (O) i (α, θ (O)) = − ln ∑ l∈[k] αl · exp ( − 1
2Ti·minl′∈[k] λmin(Σ(l ′)) ψ
(O) i (µ
(l)) )
. Based on f (O)i , we
define another sensitivity function s(O)(i) := maxα∈∆k,θ(O)∈Rd×k f
(O) i (α,θ (O))
f(O)(α,θ(O)) and want to reduce s
to s(O)i . Note that ψ (O) i (µ) = ψi(µ, Id, 0d), i.e., ψ (O) i removes the covariance and autocorrelation parameters from ψi. Due to this observation, we have that s(O)(i) ≤ maxα∈∆k,θ∈Pkλ f ′i(α,θ) f ′(α,θ) ≤ 4D · s(O)(i)/λ by upper bounding the affect of covariance parameters (controlled by factor D) and autocorrelation parameters (controlled by factor 1/λ), summarized as Lemma A.11. Then to prove Lemma 4.6, it suffices to prove that s(O)(i) ≤ 4‖bi−c ? p(i)‖ 2 2
OPT(O)+A + 3sc(i) which implies that the total sensitivity of s(O) is upper bounded by 4 + 3k (Lemma A.12), due to fact that ∑ i∈[N ] s
c(i) ≤ k (Lemma A.13). Finally, the proof of Lemma A.12 is based on a reduction from ψ(O) to the k-means clustering problem of entities (Definition 4.2) by rewriting 1Tiψ (O) i (µ) = ‖bi − µ‖22 + ai and then projecting each bi to its closest center in C?. The full proof of this lemma can be found in Section A.3.
Part 3: Running time in Theorem 4.3. The first stage costs O(d ∑ i∈[N ] Ti + Ndk lnN ln k) time.
The dominating steps are 1) to compute bi in Line 2 which costs O(d ∑ i∈[N ] Ti) time; 2) to solve the k-means clustering problem of entities (Line 3) which costs Ndk lnN ln k) time. The second stage costs at most O(d ∑ i∈[N ] Ti) time. The dominating step is to compute OPT (O) i for all i ∈ IS in Line 10, where it costs O(dTi) time to compute each OPT (O) i . Remark 4.7 (Technical comparison with prior works). Note that our entity coreset construction uses some ideas from [24, 48] and also develops novel technical ideas. A generalization of [24, 48] to time series data is to treat all observations xit independent and compute a sensitivity for each xit directly for importance sampling. However, this idea cannot capture the property that multiple observations xi,1, . . . , xi,Ti are drawn from the same Gaussian model (certain l ∈ [k]). To handle multiple observations, we show that although each ψi consists of Ti sub-functions ψit, it can be
approximated by a single function on the average observation bi = ∑ t∈[Ti] xit
Ti , specifically, we
have ψi(µ, Id, 0d) = d(bi, µ)2 + O(1). This property enables us to "reduce" the representative complexity of ψi, and motivate two key steps of our construction: 1) For the sensitivity function, we
give a reduction to a certain k-means clustering problem on average observations bi = ∑ t∈[Ti] xit
Ti of
entities (Definition 4.2) by upper bounding the maximum effect of covariances and autocorrelations and applying the fact that ψi(µ, Id, 0d) = d(bi, µ)2 +O(1). 2) For the pseudo-dimension, we prove that there are only poly(k, d) intristic operators in ψi between parameters θ and observations xit, based on the reduction of the representative complexity of ψi.
Generalizing our results in the same manner as [26] does over [24, 48] would be an interesting future direction since we may get rid of Assumption 4.1. Currently, it is unclear how to generalize the approach of [26] to time series data since [26] assumes that each point is an integral point within a bounded box, while we consider a continuous GMM generative model (1), and hence, each coordinate of an arbitrary observation xitr is drawn from a certain continuous GMM distribution which is not integral with probability ≈ 1 and can be unbounded.

heading: 5 Empirical results
text: We compare the performance of our coreset algorithm CRGMM, with uniform sampling as a benchmark using synthetic data. The experiments are conducted with PyCharm IDE on a computer with 8-core CPU and 32 GB RAM.
Datasets. We generate two sets of synthetic data all with 250K observations with different number of entities N and observations per individual Ti: (i) N = Ti = 500 for all i ∈ [N ] and (ii) N = 200, Ti = 1250 for all i ∈ [N ]. As cross-entity variance is greater than within entity variance, this helps to assess coreset performance sensitivity when there are more entities (high N ) vs. more frequent measurement/ longer measurement period (high Ti).
We fix d = 2, k = 3, λ = 0.01 and generate multiple datasets with model parameters randomly varied as follows: (i) Draw α from a uniform distribution over ∆k. (ii) For each l ∈ [k], draw
µ(l) ∈ Rd from N(0, Id). (iii) Draw Σ(l) by first generating a random matrix A ∈ [0, 1]d×d and then let Σ(l) = (AA>)−1, and (iv) draw all diagonal elements of Λ(l) ∈ Dd from a uniform distribution over [0, 1 − √ λ]. Given these draws of parameters, we generate a GMM time-series dataset as follows: For each i ∈ [N ], draw l ∈ [k] given α. Then, for all t ∈ [Ti] draw eit ∈ Rd with covariance matrix Σ(l) and autocorrelation matrix Λ(l) and compute xit = µ(l) + eit ∈ Rd. Baseline and metrics. We use coresets based on uniform sampling (Uni) and based on [48] (LFKF) that constructs a coreset for Gaussian mixture model with static data as the baseline. Given an integer Γ, uniformly sample a collection S of Γ entity-time pairs (i, t) ∈ PX , let w(i) = N|IS | for each i ∈ IS , and let w(i)(t) = Ti|JS,i| for each t ∈ JS,i; and LFKF regard all entity-time pairs as independent points and sample a collection S of Γ entity-time pairs (i, t) ∈ PX via importance sampling. For comparability, we set Γ to be the same as our coreset size. The goal of selecting LFKF as a baseline is to see the effect of autocorrelations to the objective and show the difference between static clustering and time series clustering.
Let V ? (negative log-likelihood) denote the objective of Problem 3.1 with full data. Given a weighted subset S ⊆ [N ] × [T ], we first compute (α?, θ?) := arg minα∈∆k,θ∈Pkλ f ′ S(α
′(θ), θ) + φ(α, θ) as the maximum likelihood solution over S.2 Then V ?S = f(α
?, θ?) over the full dataset serves as a metric of model fit given model estimates from a weighted subset S. We use the likelihood ratio, i.e., γS := 2(V ? S −V ?) as a measure of the quality of S [14]. The running time for GMM clustering with full data (V ?) and coreset (V ?S ) are TX and TS respectively. TC is the coreset S construction time.
Empirical setup. We vary ε = 0.1, 0.2, 0.3, 0.4, 0.5. For each ε, we run CRGMM and Uni to generate 5 coresets each. We estimate the model with the full dataset X and the coresets S and record V ?, V ?S and the run times.

heading: 1 0.1 2050 2058 2069 2041 22 34 56 1514 109 2416 1355 3436

heading: 2 0.1 811 825 841 812 2 26 58 1718 694 1859 1687 9787
text: Results. Table 1 summarizes the quality-size-time trade-offs of our coresets for different error guarantees ε. The V ?S achieved by our coreset is always close to the optimal value V
? and is always smaller than that of Uni and LFKF. We next assess coreset size. For the N = Ti = 500 case, the coreset achieves a likelihood ratio γS = 684 with only 72 entity-time pairs (0.03%), but Uni/LFKF achieves the closest (but worse) likelihood ratio of γS = 714/686, with 191 entity-time pairs (0.08%). Thus, our coreset achieves a better likelihood than uniform and LFKF with less than 40% of observations. Further, from Figure 1, we see that not only is the quality of our coreset (γS) superior, it has lower standard deviations suggesting lower variance in performance. In terms of total computation time relative to the full data ( TXTC+TS ), our coresets speed up by 14x-171x. Also, the
2We solve this optimization problem using an EM algorithm similar to [5]. The M step in each iteration is based on IRLS [40] and the E-step involves an individual level Bayesian update for α.
computation time (TC + TS) of Uni and LFKF is always larger than that of our coreset (3x-22x).3 Finally, V ?S is lower for Synthetic 2 given that there are fewer entities and cross-entity variance is greater than within entity variance. From Figure 1, we also see that both the fit and std in performance for coresets is much lower for Synthetic 2. Thus, overall our coreset performs better when there is more time series data relative to entities (Ti N )—a feature of sensor based measurement that motivates this paper.

heading: 6 Limitations, conclusion, and future work
text: In this paper, we study the problem of constructing coresets for clustering problems with time series data; in particular, we address the problem of constructing coresets for time series data generated from Gaussian mixture models with auto-correlations across time. Our coreset construction algorithm is efficient under a mild boundedness assumption on the covariance matrices of the underlying Gaussians, and the size of the coreset is independent of the number of entities and the number of observations and depends only polynomially on the number of clusters, the number of variables and an error parameter. Through empirical analysis on synthetic data, we demonstrate that the coreset sampling is superior to uniform sampling in computation time and accuracy.
Our work leaves several interesting directions for future work on time series clustering. While our current specification with autocorrelations assumes a stable time series pattern over time, future work should extend it to a hidden Markov process, where the time series process may switch over time. Further, while our focus here is on model-based clustering, it would be useful to consider how coresets should be constructed for other clustering methods such as direct clustering of raw data and indirect clustering of features.
Overall, we hope the paper stimulates more research on coreset construction for time series data on a variety of unsupervised and supervised machine learning algorithms. The savings in storage and computational cost without sacrificing accuracy is not only financially valuable, but also can have sustainability benefits through reduced energy consumption. Finally, recent research has shown that summaries (such as coresets) for static data need to include fairness constraints to avoid biased outcomes for under-represented groups based on gender and race when using the summary [16, 36]; future research needs to extend such techniques for time series coresets.

heading: Acknowledgments
text: This research was supported in part by an NSF CCF-1908347 grant.
3A possible explanation is that our coreset selects IS of entities that are more representative that Uni, and hence, achieves a better convergence speed.

Review summaries:
This paper uses the Feldman Langberg sensitivity framework to compute coresets for time series clustering under a Gaussian Mixture Model.  This imposes additional challenges since there is an instance domain as well as a time domain to approximate. The mixture model introduces covariances and to reflect the time dependence there is an additional autocorrelation introduced on the time domain.  Those covariance and correlation structures are being bounded in order to bound their impact on the sensitivities and after getting rid of all those structures the problem is reduced to a simpler clustering model for which sensitivity bounds are obtained more easily.   The resulting coresets are assessed empirically on synthetic data.  The writing is very good and the main body provides all information to understand even on a technical level, although the details are in the appendix.  In principle I would tend to accept this paper but knowing the previous literature on GMM as well as a recent coreset work on "panel data" my main impression is that the submission is just a dejavu of those two works.

In this paper the authors provide the definition of a coreset for time series clustering for entities generated from a gaussian mixture model and then build a coreset for the same problem.  The coreset is independent of both the number of entities N and the the number of Time stamps T_i for each entity i. They demonstarte the accuracy and efficiency of their coreset by superior empirical perfromance over synthetic data as compared to uniform sampling.

The authors propose a novel coreset construction framework for the problem of clustering multivariate time series. To my knowledge, this is the first construction of coresets in this setting. The framework, as well as the proof scheme, is close to the work in [45], in which Lucic et al. studied the problem of coreset constructions for the MLE task under the GMM model.   Here, the authors add to the GMM model a temporality and thus an additional autocorrelation in the underlying model. The proof's scheme is similar to that of [45] in the sense that they also reduce the GMM time series clustering problem to a $k$-means clustering problem (and then use more classical results for $k$-means from the state-of-the-art).

The paper considers the problem of clustering time series data. The authors provide a coreset for this problem, under some assumptions. The coreset size depends polynomially on the number of clusters k, the dimension d, and 1/\varepsilon and the parameters related to the assumption on the data, and requires near-linear time to compute.  The authors validate their coreset empirically on synthetic data, as compared to a uniform random sample. They show that the coreset: (i) achieves better results can the random sample, and (i) can help accelerate the computational time, relative to the full data, by 1-2 orders of magnitude, with only a small compromise in the accuracy.