EconBase
← Back to paper

Coresets for Time Series Clustering

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

91,526 characters · 17 sections · 90 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Coresets for Time Series Clustering

abstractWe 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 $\mathbb{R}^d$. 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/\varepsilon$, where $\varepsilon$ is the error parameter. We empirically assess the performance of our coresets with synthetic data.

Introduction

A multivariate time series dataset, represented as ${\mathcal{X}} = \left\{X_i=\left(x_{i1}, \ldots, x_{i,T_i}\right)\subset \mathbb{R}^{T_i\times d}\mid i\in [N]\right\}$, where $N$ is the number of entities, $T_i$ 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 baltagi2008econometric 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 aghabozorgi2015time.

\sloppy Computational problems of interest with time series data include pattern discovery marlin2012unsupervised, regression/prediction hung2015identifying, forecasting bandara2020forecasting, and clustering liao2005clustering, fu2011review,aghabozorgi2015time which arises in applications such as anomaly detection li2021clustering. 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 hennig2015handbook,blum2020foundations. Time series clustering introduces a number of additional challenges and is an active area of research, with many types of clustering methods proposed caiado2015time; 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 liao2005clustering,aghabozorgi2015time.

We focus on model-based time series clustering using a likelihood framework that naturally extends clustering with static data blum2020foundations to time series data liao2005clustering, where each multivariate time series is generated by one of $k$ different parametrically specified models. \textcolor{black}{For a survey of the importance of this sub-literature and its various real-world applications, see fruhwirth2011panel.} A prevalent approach is to assume a finite mixture of data generating models, with Gaussian mixture models being a common specification ouyang2004gaussian,reddy2017using. Roughly, the problem is to partition ${\mathcal{X}}$ probabilistically into $k$ clusters, grouping those series generated by the same time series model into one cluster. Generically, the model-based $k$-clustering problem can be formulated as: $\arg\max_{\alpha,\theta^{(1)},\ldots,\theta^{(k)}} \sum_{i\in [N]} \ln \sum_{l\in [k]} \alpha_l \cdot p_{i}({\cal{X}}|\theta^{(l)})$, where $\alpha_l$ is the probability of a time series belonging to cluster $l$ and $p_i({\cal{X}}|\theta^{(l)})$ is the likelihood of the data of entity $i$, given that the data is generated from model $l$ (with parameters $\theta^{(l)}$). Temporal relationships in time series are commonly modeled as autocorrelations (AR)/moving averages (MA) d2009autocorrelation, xiong2004time that account for correlations across observations, and hidden Markov models (HMM) where the underlying data generating process is allowed to switch periodically ernst2005clustering, yang2014hmm. This paper focuses on the case when the data generation model for each segment $l \in [k]$ is a multivariate Gaussian with autocorrelations. More formally, the generative model for each cluster is from a Gaussian mixture: $x_{it} := \mu^{(l)} + e_{it}$, where $e_{it} = N(0,\Sigma^{(l)})+ \Lambda^{(l)} e_{i, t-1}$, where $N(0,\Sigma^{(l)})$ captures the mixture of Gaussian distributions from which entity level observations are drawn, and $\Lambda$ captures the correlation between two successive observations through an AR(1) process haddad1998simple,lesage1999theory. Overall, this model can be represented with cluster level model parameters $\theta^{(l)} = \{\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)}\}$ and cluster probability $\alpha_l$.

Time series datasets are also much larger than static datasets liao2005clustering. For instance, as noted in keogh2006decade, 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 feldman2009private, ke2020privacy. 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 harpeled2004on. 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 bachem2017practical, feldman2020core. 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 huang2020coresets. However, as we explain in Section (ref), there are challenges in extending coresets for static data clustering and time series (panel) data regressions for time series clustering.

\paragraph{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 (ref)). 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 (ref)). To estimate this objective using a coreset, we introduce an analogous clustering objective on a subset of data (Definition (ref)). Our main result is an algorithm to construct coresets for this objective for the aforementioned mixture model (Theorem (ref)). Our algorithm is efficient (linear on both $N$ and $\sum_{i\in [N]}T_i$) and, assuming that the condition number of the covariance matrices of the Gaussians ($\Sigma^{(l)}$) are bounded and the eigenvalues of the autocorrelation matrices ($\Lambda^{(l)}$) are bounded away from $1$ (Assumption (ref)), the size of our coreset is independent of the number of entities $N$ and the number of observations $T_i$s; it only polynomially depends on $k$, $d$ and $1/\varepsilon$, where $\varepsilon$ is the error parameter (Theorem (ref)).

Our coreset construction (Algorithm (ref)) leverages the Feldman-Langberg (FL) framework feldman2011unified. 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 (ref)) 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 $\Sigma$ and $\Lambda$ 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 (ref)), and compare the performance with two benchmarks: uniform sampling and a static coreset benchmark (LFKF lucic2017training) 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 ($T_i \gg N$).

Related work

Time series or panel data analysis are important in many fields-- biological cohen2014analyzing, engineering akaike2012practice, economics pesaran2015time and the social sciences wei2006time. 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 bock2007clustering. While most standard treatments of clustering in machine learning tend to be on static data hennig2015handbook,blum2020foundations, there is by now a significant and rapidly growing literature on time series clustering; recent surveys include liao2005clustering,aghabozorgi2015time. 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).

\paragraph{Coresets for clustering.} Coresets for various clustering objectives on static data have been well studied, including $k$-median harpeled2004on,chen2009on,sohler2018strong,huang2020coresets_a,Cohen-Addad2021new, $k$-means harpeled2004on,chen2009on,feldman2013turning,braverman2016new,becchetti2019oblivious,sohler2018strong,huang2020coresets_a,Cohen-Addad2021new, and $k$-center agarwal2002exact,har2004clustering. For surveys of static coresets, see phillips2016coresets,feldman2020introduction. lucic2017training,feldman2019coresets present coreset constructions for clustering problems on static data generated using Gaussian mixture models. lucic2017training constructs a coreset of size $\ensuremath{\mathsf{poly}}(k,d, 1/\varepsilon)$ under the same boundedness assumption as ours. feldman2019coresets remove this assumption and construct a coreset whose size additionally depends on $\ensuremath{\mathsf{poly}}\log N$, but require another assumption that all $x_{i}$s 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 (ref)) to a $k$-means clustering problem whose points are linear combinations of those in $X_i$ (Lemma (ref)), and show that the effect of autocorrelation parameters $\Lambda$ on the entity objectives can be upper bounded (Lemma (ref)).

\paragraph{Regression with time series data.} huang2020coresets construct coresets for a regression problem, called GLSE$_k$, for time-series data including autocorrelations and a $k$-partition of entities. This is a supervised learning task that includes an additional label $y_{it}$ for each entity-time pair. It aims to find a $k$-partition of entities and additional regression vectors $\beta^{(l)}\in \mathbb{R}^d$ for each partition $l\in [k]$ such that a certain linear combination of terms $\|y_{it}-(\beta^{(l)})^\top x_{it}\|_2^2$ is minimized. In contrast, our problem is an unsupervised learning task, estimates points $x_{it}$ by Gaussian means $\mu$ instead of $\beta$, and includes covariances $\Sigma$ between different features. The definitions of coresets in huang2020coresets 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 (ref)) 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. huang2020coresets'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 $\mu^{(l)}$ are far away from observations of this entity. To bypass this, we reduce our problem to a $k$-means clustering problem (Definition (ref)), whose total sensitivity is provably $O(k)$ (Lemma (ref)), by upper bounding the affect of covariance and autocorrelation parameters on the clustering objectives of entities (Lemma (ref)) under the assumption that the condition number of covariance matrix is bounded (Assumption (ref)), which provides an upper bound for the total sensitivity of entities on our problem (Lemma (ref)).

Another related direction is to consider coreset construction for $k$-segmentation with time series data rosman2014coresets,feldman2015idiary, which aims to estimate the trace of an entity by a $k$-piecewise linear function ($k$-segment). Note that the case of $k=1$ is equivalent to the linear regression problem. rosman2014coresets,feldman2015idiary 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.

Clustering model and coreset definition

Given ${\mathcal{X}} = \left\{X_i=\left(x_{i,1}, \ldots, x_{i,T_i}\right)\in \mathbb{R}^{d\times T_i}\mid i\in [N]\right\}$, we first model a general class of time series clustering problems. Then we specify our setting with Gaussian mixture time series data (Problem (ref)), and define coresets accordingly (Definition (ref)).

\paragraph{Clustering with time series data.} Given an integer $k\geq 1$, let $\Delta_k\subset \mathbb{R}^k$ denote a probability simplex satisfying that for any $\alpha\in \Delta_k$, we have $\alpha_i\in [0,1]$ for $i\in [k]$ and $\sum_{l\in [k]} \alpha_i = 1$. Let ${\mathcal{P}}$ denote a parameter space where each $\theta\in {\mathcal{P}}$ represents a specific generative model of time series data. For each entity $i\in [N]$ and model $\theta\in {\mathcal{P}}$, define $p_i({\mathcal{X}}\mid \theta) := \Pr\left[X_i\mid \theta\right]^{1/T_i}$ to be the average likelihood/realized probability of $X_i$ from model $\theta$. The ratio $1/T_i$ is used to normalize the contribution to the objective of entity $i$ due to different lengths $T_i$. The time series clustering problem is to compute $\alpha\in \Delta_k$ and $\theta = (\theta^{(1)}, \ldots, \theta^{(k)})\in {\mathcal{P}}^k$ that minimizes the negative log-likelihood, i.e., $ -\sum_{i\in [N]} \ln \sum_{l\in [k]} \alpha_l \cdot p_{i}({\mathcal{X}}\mid \theta^{(l)}). $ Here, for each $l\in [k]$, $\alpha_l$ represents the probability that each time series is generated from model $\theta^{(l)}$. Note that the clustering objective depends on the choice of model family ${\mathcal{P}}$. In this paper, we consider the following specific model family ${\mathcal{P}}$ that each time series is generated from Gaussian mixtures.

\sloppy \paragraph{GMM clustering with time series data.} Let $\alpha\in \Delta_k$ be a given probability vector. Let ${\mathcal{S}}^d$ denote the collection of all symmetric positive definite matrices in $\mathbb{R}^{d\times d}$. For each $i\in [N]$, with probability $\alpha_l$ ($l\in [k]$), let $x_{it} := \mu^{(l)} + e_{it}$ for each $t\in [T_i]$ where $\mu^{(l)}\in \mathbb{R}^d$ represents a Gaussian mean and $e_{it}\in \mathbb{R}^d$ represents the error vector drawn from the following normal distribution: $e_{it} := \Lambda^{(l)} e_{i, t-1} + N(0,\Sigma^{(l)})$, where $\Lambda^{(l)}\in {\mathcal{S}}^d$ is an AR(1) autocorrelation matrix and $\Sigma^{(l)}\in {\mathcal{S}}^d$ is the covariance matrix of a multivariate Gaussian distribution. Now we let ${\mathcal{P}} := \mathbb{R}^d\times {\mathcal{S}}^d\times {\mathcal{S}}^{d}$ and note that each Gaussian generative model can be represented by a tuple $\theta = (\mu,\Sigma, \Lambda)\in {\mathcal{P}}$. Moreover, we have that the realized probability of each entity $i\in [N]$ is

eqnarray[eqnarray omitted — 209 chars of source]

where $|\Sigma|$ is the determinant of $\Sigma$ and $\psi_i(\mu, \Sigma, \Lambda) := \sum_{t\in [T_i]} \psi_{it}(\Sigma, \Lambda)$ with \[ \textstyle \psi_{i1}(\Sigma, \Lambda) := (x_{i,1}-\mu)^\top \Sigma^{-1} (x_{i,1}-\mu)- \left(\Lambda(x_{i,1}-\mu)\right)^\top \Sigma^{-1} \left(\Lambda(x_{i,1}-\mu)\right), \text{ and } \] \[ \psi_{it}(\Sigma, \Lambda) := \left((x_{it}-\mu)- \Lambda(x_{i,t-1}-\mu)\right)^\top \Sigma^{-1} \left((x_{it}-\mu)-\Lambda(x_{i,t-1}-\mu)\right), \forall \; \; 2\leq t\leq T_i. \] We note that each sub-function $\psi_{it}$ contains at most two entity-time pairs: $x_{it}$ and $x_{i,t-1}$. Our Gaussian mixture time series model gives rise in the following clustering problem.

problem[\bf{Clustering with GMM time series data}] Given a time series dataset ${\mathcal{X}} = \left\{X_i=\left(x_{i,1}, \ldots, x_{i,T_i}\right)\in \mathbb{R}^{d\times T_i}\mid i\in [N]\right\}$ and integer $k\geq 1$, the GMM time series clustering problem is to compute $\alpha\in \Delta_k$ and $\theta=(\mu^{(l)},\Sigma^{(l)}, \Lambda^{(l)})_{l\in [k]}\in {\mathcal{P}}^k$ that minimize \[ \textstyle f(\alpha,\theta):= \sum_{i\in [N]} f_i(\alpha,\theta) = -\sum_{i\in [N]} \ln \sum_{l\in [k]} \alpha_l \cdot p_i({\mathcal{X}}\mid\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)}). \]

From Equation (ref) it follows that the coefficient of each Gaussian before the $\exp$ operation is $\frac{\alpha_l }{(2\pi)^{d/2} |\Sigma^{(l)}|^{1/2}}$, whose summation $Z(\alpha, \theta) := \sum_{l\in [k]} \frac{\alpha_l }{(2\pi)^{d/2} |\Sigma^{(l)}|^{1/2}}$ is not a fixed value and depends on $\Sigma^{(l)}$s. To remove this dependence on the coefficients, we define $\alpha'_l(\theta) := \frac{\alpha_l }{(2\pi)^{d/2} |\Sigma^{(l)}|^{1/2} Z(\alpha, \theta)}$ to be the normalized coefficient for $l\in [k]$, define offset function $\phi: \Delta_k\times {\mathcal{P}}^k\rightarrow \mathbb{R}$ to be $\phi(\alpha, \theta) := -N \ln Z(\alpha, \theta)$, and define $ f'(\alpha, \theta):= -\sum_{i\in [N]} \ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i}\psi_i(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)}))$ whose summation of coefficients before the $\exp$ operations is 1. This idea of introducing $f'$ also appears in feldman2011scalable,feldman2019coresets, which is useful for coreset construction. This leads to the following observation.

observationFor any $\alpha\in \Delta_k$ and $\theta\in {\mathcal{P}}^k$, $ f(\alpha, \theta) = f'(\alpha'(\theta), \theta) + \phi(\alpha,\theta). $

\paragraph{Our coreset definition.} The clustering objective $f$ can only be decomposed into the summation of $f_i$s 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 $T_i$s and fail to be independent of both $N$ and $T_i$s. To address this problem, we define our coreset as a collection of weighted entity-time pairs, which is similar to huang2020coresets. Let $P_{\mathcal{X}}:=\left\{(i,t)\mid i\in [N], t\in [T_i]\right\}$ denote the collection of indices of all $x_{it}$. Given $S\subseteq P_{\mathcal{X}}$, we let $I_S := \left\{i\in [N]\mid \exists t\in [T_i], s.t., (i,t)\in S \right\}$ denote the collection of entities that appear in $S$. Moreover, for each $i\in I_S$, we let $J_{S,i}:=\left\{t\in [T_i]: (i,t)\in S\right\}$ denote the collection of observations for entity $i$ in $S$.

definition[\bf{Coresets for GMM time series clustering}] Given a time series dataset ${\mathcal{X}} = \left\{X_i=\left(x_{i,1}, \ldots, x_{i,T_i}\right)\in \mathbb{R}^{d\times T_i}\mid i\in [N]\right\}$, constant $\varepsilon \in (0,1)$, integer $k\geq 1$, and parameter space $\Delta_k\times {\mathcal{P}}^k$, an $\varepsilon$-coreset for GMM time series clustering is a weighted set $S\subseteq P_{\mathcal{X}}$ together with weight functions $w: I_S\rightarrow \mathbb{R}_{\geq 0}$ and $w^{(i)}: J_{S,i}\rightarrow \mathbb{R}_{\geq 0}$ for $i\in I_S$ such that for any $\alpha\in \Delta_k$ and $\theta=(\mu^{(l)},\Sigma^{(l)}, \Lambda^{(l)})_{l\in [k]}\in {\mathcal{P}}^k$, $ f'_S(\alpha,\theta) := -\sum_{i\in I_S} w(i) \cdot \ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i} \sum_{t\in J_{S,i}} w^{(i)}(t) \cdot \psi_{it}(\mu^{(l)},\Sigma^{(l)}), \Lambda^{(l)})) \in (1\pm \varepsilon)\cdot f'(\alpha,\theta). $

Combining the above definition with Observation (ref), we note that for any $\alpha\in \Delta_k$ and $\theta\in {\mathcal{P}}^k$,

align[align omitted — 246 chars of source]

As $\varepsilon$ tends to 0, $f'_S(\alpha'(\theta), \theta) + \phi(\alpha, \theta)$ converges to $f(\alpha, \theta)$. Moreover, given such a coreset $S$, if we additionally have that $\phi(\alpha, \theta)\lesssim f(\alpha, \theta)$, we can minimize $f'_S(\alpha'(\theta),\theta) + \phi(\alpha,\theta)$ to solve Problem (ref). Thus, if $\phi(\alpha, \theta)\lesssim f(\alpha, \theta)$, we conclude that $f'_S(\alpha'(\theta), \theta)+ \phi(\alpha,\theta) \in (1\pm 3\varepsilon)\cdot f(\alpha, \theta)$. Note that we introduce multiple weight functions ($w$ for entities and $w^{(i)}$ for time periods of selected entity $i$) in Definition (ref) unlike the coreset for time series regression in huang2020coresets which uses only one weight function. This is because $f'_i$ contains $\ln$ and $\exp$ operators instead of linear combinations of $\psi_{it}$; it is unclear how to use a single weight function to capture both the entity and time levels.

Theoretical results

We present our coreset algorithm (Algorithm (ref)) and the main result on its performance (Theorem (ref)). Theorem (ref) needs the following assumptions on covariance and autocorrelation parameters.

assumptionAssume (1) that there exists constant $D\geq 1$ such that $\frac{\max_{l\in [k]} \lambda_{\max}(\Sigma^{(l)})}{\min_{l\in [k]} \lambda_{\min}(\Sigma^{(l)})}\leq {D}$ where $\lambda_{\max}(\cdot)$ represents the largest eigenvalue and $\lambda_{\min}(\cdot)$ represents the smallest eigenvalue, and (2) for $l\in [k]$, $\Lambda^{(l)}\in {\mathcal{D}}^d_{\lambda}$ for some constant $\lambda\in (0,1)$. Here ${\mathcal{D}}^d_{\lambda}$ is a collection of all diagonal matrix in $\mathbb{R}^{d\times d}$ whose diagonal elements are at most $1-\sqrt{\lambda}$.

The first assumption requires that the condition number of each covariance matrix is upper-bounded, which also appears in won2013condition,li2020efficient,lucic2017training 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 huang2020coresets. Note that both the eigenvalues of $\Sigma^{(i)}$ and the positions of means $\mu^{(l)}$s affect cost function $\psi_i$s. For instance, consider the case that all eigenvalues are 1 and all autocorrelations are 0, i.e., for all $l\in [k]$, $\Sigma^{(l)}= I_d$ and $\Lambda^{(l)}=0_d$. In this case, it is easy to see that $\frac{\max_{\mu\in \mathcal{R}^d}\psi_i(\mu, I_d, 0_d)}{\min_{\mu\in \mathcal{R}^d}\psi_i(\mu, I_d, 0_d)} = \infty$, i.e., the value of $\psi_i(\mu, I_d, 0_d)$ is unbounded as $\mu$ changes. Differently, the component's cost functions are bounded in huang2020coresets since $\mu^{(l)}$s do not appear in huang2020coresets that consider regression problems.

Let ${\mathcal{P}}_{\lambda}^k := (\mathbb{R}^d\times {\mathcal{S}}^d\times {\mathcal{D}}^d_{\lambda})^k$ denote the parameter space. For preparation, we propose the following $k$-means clustering problem.

definition[\bf{$k$-means clustering of entities}] Given an integer $k \geq 1$, the goal of the $k$-means clustering problem of entities of $X$ is to find a set $C^\star = \left\{c^\star_1,\ldots, c^\star_k \right\}\subset \mathbb{R}^d$ of $k$ centers that minimizes $ \sum_{i\in [N]} \min_{l\in [k]} \|\frac{\sum_{t\in [T_i]}x_{it}}{T_i} - c_l\|_2^2 $ over all $k$ center sets $C=\left\{c_1, \ldots, c_k\right\}\subset \mathbb{R}^d$. Let $\ensuremath{\mathsf{OPT}}^{(O)}$ denote the optimal $k$-means value of $C^\star$.

Note that there exists an $O(2^{\ensuremath{\mathsf{poly}}(k)}Nd)$ time algorithm to compute a nearly optimal solution for both $\ensuremath{\mathsf{OPT}}^{(O)}$ and $C^\star$ kumarSS2004simple. Another widely used algorithm is called $k$-means++ arthur2007kmeans, which provides an $O(\ln k)$-approximation in $O(Ndk \ln N \ln k)$ time but performs well in practice. We also observe that $\psi_i(\mu, I_d, 0_d) = \|\frac{\sum_{t\in [T_i]}x_{it}}{T_i} - \mu\|_2^2 + \frac{1}{T_i}\sum_{t\in [T_i]} \|x_{it}\|_2^2 - \frac{\|\sum_{t\in [T_i]} x_{it}\|_2^2}{T_i^2}$ for any $\mu\in \mathbb{R}^d$. This observation motivates us to consider a reduction from Problem (ref) to Definition (ref), which is useful for constructing $I_S$.

Our coreset algorithm

We first give a summary of our algorithm (Algorithm (ref)). The input to Algorithm (ref) is a time series dataset ${\mathcal{X}}$, an error parameter of coreset $\varepsilon\in (0,1)$, an integer $k\geq 1$ of clusters, $\lambda\in (0,1)$ (vector norm bound), and a covariance eigenvalue gap $D\geq 1$. We develop a two-staged importance sampling framework which first samples a subset $I_S$ of entities (Lines 1-8) and then samples a subset of time periods $J_{S,i}$ for each selected entity $i\in I_S$ (Lines 9-14).

In the first stage, we first set $M$ as the number of selected entities $|I_S|$ (Line 1). Then we solve the $k$-means clustering problem over entity means $b_i$ (Definition (ref)), e.g., by $k$-means++ arthur2007kmeans (Lines 2-3), and obtain an (approximate) optimal clustering value $\ensuremath{\mathsf{OPT}}^{(O)}$, a $k$-center set $C^\star:=\left\{c^\star_1, \ldots, c^\star_k\right\}\subset \mathbb{R}^d$, and an offset value $a_i$ for each $i\in [N]$. Next, based on the distances of points $b_i$ to $C^\star$, we partition $[N]$ into $k$ clusters where entity $i$ belongs to cluster $c^\star_{p(i)}$ (Line 4). Based on $b_i$, $c^\star_{p(i)}$ and $\ensuremath{\mathsf{OPT}}^{(O)}$, we compute an upper bound $s(i)$ for the sensitivity w.r.t. entity $i\in [N]$ (Lines 5-6). Finally, we sample a weighted subset $I_S$ 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 (ref)).

In the second stage, we first set $L$ as the number of time periods $|J_{S,i}|$ for each selected entity $i\in I_S$ (Line 9). Then for each selected entity $i\in I_S$, we compute $\ensuremath{\mathsf{OPT}}_i^{(O)}$ as the optimal $1$-means clustering objective of $x_{it}$s (Line 10). Next, we compute an upper bound $s_i(t)$ for the sensitivity w.r.t. time period $t\in [T_i]$ (Lines 11-12) based on $b_i$ and $\ensuremath{\mathsf{OPT}}^{(O)}_i$. Finally, we sample a weighted subset $J_{S,i}$ of time periods by importance sampling based on $s_i(t)$ (Lines 13-14).

algorithm[algorithm omitted — 3,325 chars of source]

\footnotetext{Here, we slightly abuse the notation by using $\ensuremath{\mathsf{OPT}}^{(O)}$ and $C^\star$ to represent the obtained approximate solution.}

Our main theorem

Our main theorem is as follows, which indicates that Algorithm (ref) provides an efficient coreset construction algorithm for GMM time series clustering.

theorem[\bf{Main result}] Under Assumption (ref), with probability at least $0.9$, Algorithm (ref) outputs an $\varepsilon$-coreset of size $O\left(\frac{D^2 k^5d^{16}\ln^2\frac{k}{\lambda}}{\lambda^2 \varepsilon^4}\right)$ for GMM time-series clustering in $O(d\sum_{i\in [N]}T_i + Ndk\ln N \ln k)$ time.

Note that the success probability in the theorem can be made $1-\delta$ for any $\delta \in (0,1)$ at the expense of an additional factor $\log^2 1/\delta$ in the coreset size. The coreset size guaranteed by Theorem (ref) has a polynomial dependence on factors $k,d,1/\varepsilon,D$, and $1/\lambda$ and, importantly, does not depend on $N$ or $T_i$s. Compared to the clustering problem with GMM static data lucic2017training, the coreset has an additional dependence on $1/\lambda$ due to autocorrelation parameters. Compared to the regression problem (GLSE$_k$) with time series data huang2020coresets, 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 $\sum_{i\in [N]}T_i$, which is efficient. The proof of Theorem (ref) can be found in Section (ref).

\paragraph{Proof overview of Theorem (ref).} 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.

{\em Parts 1 & 2: Coreset size and correctness guarantee of Theorem (ref).} We first note that the coreset size in Theorem (ref) is $|S| = ML$. Here, $M$ is the number of sampled entities $|I_S|$ that guarantee that $I_S$ is a coreset for $f'$ at the entity level (Lemma (ref)). $L$ is the number of sampled observations $|J_{S,i}|$ for each $i\in I_S$ that guarantees that $J_{S,i}$ is a coreset for $\psi_i$ at the time level (Lemma (ref)).

lemma[\bf{The 1st stage outputs an entity-level coreset}] With probability at least $0.95$, the output $I_S$ of Algorithm (ref) with $M:=|I_S|=O\left(\frac{kD}{\lambda \varepsilon^2}(k^4d^4+k^3d^8)\ln\frac{k}{\lambda}\right)$ satisfies that for any $\alpha\in \Delta_k$ and $\theta\in {\mathcal{P}}_\lambda^k$, $ -\sum_{i\in I_S} w(i)\cdot \ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i}\psi_i(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)})) \in (1\pm \varepsilon)\cdot f'(\alpha, \theta). $
lemma[\bf{The 2nd stage outputs time level coresets for each $i\in I_S$}] With probability at least $0.95$, for all $i\in I_S$, the output $J_{S,i}$ of Algorithm (ref) with $L:=|J_{S,i}|=O\left(\frac{D d^8 \ln\frac{1}{\lambda}}{\lambda\varepsilon^2}\right)$ satisfies that for any $(\mu, \Sigma, \Lambda)\in {\mathcal{P}}_{\lambda}$, $ \sum_{t\in J_{S,i}} w^{(i)}(t) \cdot \psi_{it}(\mu, \Sigma,\Lambda) \in (1\pm \varepsilon)\cdot \psi_i(\mu, \Sigma,\Lambda). $

From Lemmas (ref) and (ref), we can conclude that $S$ is an $O(\varepsilon)$-coreset for GMM time-series clustering, since the errors $\varepsilon$ for the entity level and the time level are additively accumulated, and that the size of the coreset $ML$ as claimed; see Section (ref) for a complete proof.

Both lemmas rely on the Feldman-Langberg\ framework feldman2011unified,braverman2016new (Theorem (ref)). The key is to upper bound both the “pseudo-dimension” $\dim$ (Definition (ref)) that measures the complexity of the parameter space, and the total sensitivity ${\mathcal{G}}$ (Definition (ref)) that measures the sum of maximum influences of all sub-functions ($f'_i$ for Lemma (ref) and $\psi_{it}$ for Lemma (ref)). The coreset size guaranteed by the Feldman-Langberg\ framework then is $O(\varepsilon^{-2} {\mathcal{G}} \dim \ln {\mathcal{G}})$.

For Lemma (ref), we can upper bound the pseudo-dimension by $\dim = O(d^8)$ (Lemma (ref)). For the total sensitivity, a key step is to verify that $s_i$ (Line 12) is a sensitivity function of $\psi_i$. This step uses a similar idea as in huang2020coresets and leads to showing a ${\mathcal{G}} = O(D/\lambda)$ bound for the total sensitivity (Lemma (ref)). Then we verify that $L=O(\varepsilon^{-2} {\mathcal{G}} \dim \ln {\mathcal{G}})=O\left(\frac{D d^8 \ln\frac{1}{\lambda}}{\lambda\varepsilon^2}\right)$ is enough for Lemma (ref) by the Feldman-Langberg\ framework (Theorem (ref)). The proof can be found in Section (ref).

Lemma (ref) is the most technical. The proof can be found in Section (ref) and we provide a proof sketch.

proofsketchBy anthony2009neural,vidyasagar2002theory, the pseudo-dimension is determined by the number of parameters that is upper bounded by $O(kd^2)$ and the number of operations on parameters that is upper bounded by $O(kd^6)$ including $O(k)$ exponential functions. Then we can upper bound the pseudo-dimension by $\dim = O(k^4 d^4 + k^3 d^8)$ (Lemma (ref) in Section (ref)). 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 (ref) below. Then the total sensitivity is at most ${\mathcal{G}}=(16D+12Dk)/\lambda$. This completes the proof of Lemma (ref) since $M=O(\varepsilon^{-2} {\mathcal{G}} \dim \ln {\mathcal{G}})=O\left(\frac{kD}{\lambda \varepsilon^2}(k^4d^4+k^3d^8)\ln\frac{k}{\lambda}\right)$ is enough by the Feldman-Langberg\ framework (Theorem (ref)). \begin{lemma}[\bf{$s$ is a sensitivity function w.r.t. $f'$}] For each $i\in [N]$, we have $ s(i) \geq \max_{\alpha\in \Delta_k, \theta \in {\mathcal{P}}_\lambda^k} \frac{f'_i(\alpha, \theta)}{f'(\alpha,\theta)}. $ Moreover, $\sum_{i\in [N]} s(i) \leq (16D+12Dk)/\lambda$. \end{lemma} {\em Proof sketch of Lemma (ref).} We first define $\psi^{(O)}_i(\mu) := \sum_{t\in [T_i]} \|x_{it}-\mu\|_2^2$ for any $\mu\in \mathbb{R}^d$, and define $f_i^{(O)}(\alpha, \theta^{(O)})= - \ln \sum_{l\in [k]} \alpha_l \cdot \exp\left(-\frac{1}{2T_i\cdot \min_{l'\in [k]}\lambda_{\min}(\Sigma^{(l')})}\psi^{(O)}_i(\mu^{(l)})\right)$. Based on $f_i^{(O)}$, we define another sensitivity function $ s^{(O)}(i) := \max_{\alpha\in \Delta_k, \theta^{(O)}\in \mathbb{R}^{d\times k}} \frac{f^{(O)}_i(\alpha, \theta^{(O)})}{f^{(O)}(\alpha, \theta^{(O)})}$ and want to reduce $s$ to $s^{(O)}_i$. Note that $\psi^{(O)}_i(\mu) = \psi_i(\mu, I_d, 0_d)$, i.e., $\psi^{(O)}_i$ removes the covariance and autocorrelation parameters from $\psi_i$. Due to this observation, we have that $s^{(O)}(i) \leq \max_{\alpha\in \Delta_k, \theta \in {\mathcal{P}}_\lambda^k} \frac{f'_i(\alpha, \theta)}{f'(\alpha,\theta)} \leq 4D\cdot s^{(O)}(i)/\lambda$ by upper bounding the affect of covariance parameters (controlled by factor $D$) and autocorrelation parameters (controlled by factor $1/\lambda$), summarized as Lemma (ref). Then to prove Lemma (ref), it suffices to prove that $s^{(O)}(i) \leq \frac{4 \|b_i - c^\star_{p(i)}\|_2^2}{\ensuremath{\mathsf{OPT}}^{(O)}+A} + 3 s^c(i)$ which implies that the total sensitivity of $s^{(O)}$ is upper bounded by $4+ 3k$ (Lemma (ref)), due to fact that $\sum_{i\in [N]} s^c(i) \leq k$ (Lemma (ref)). Finally, the proof of Lemma (ref) is based on a reduction from $\psi^{(O)}$ to the $k$-means clustering problem of entities (Definition (ref)) by rewriting $\frac{1}{T_i}\psi^{(O)}_i(\mu) = \|b_i - \mu\|_2^2 + a_i$ and then projecting each $b_i$ to its closest center in $C^\star$. The full proof of this lemma can be found in Section (ref). {\em Part 3: Running time in Theorem (ref).} The first stage costs $O(d\sum_{i\in [N]}T_i + Ndk\ln N \ln k)$ time. The dominating steps are 1) to compute $b_i$ in Line 2 which costs $O(d\sum_{i\in [N]}T_i)$ time; 2) to solve the $k$-means clustering problem of entities (Line 3) which costs $Ndk\ln N \ln k)$ time. The second stage costs at most $O(d\sum_{i\in [N]}T_i)$ time. The dominating step is to compute $\ensuremath{\mathsf{OPT}}_i^{(O)}$ for all $i\in I_S$ in Line 10, where it costs $O(d T_i)$ time to compute each $\ensuremath{\mathsf{OPT}}_i^{(O)}$.
remark[\bf{Technical comparison with prior works}] Note that our entity coreset construction uses some ideas from feldman2011scalable,lucic2017training and also develops novel technical ideas. A generalization of feldman2011scalable,lucic2017training to time series data is to treat all observations $x_{it}$ independent and compute a sensitivity for each $x_{it}$ directly for importance sampling. However, this idea cannot capture the property that multiple observations $x_{i,1},\ldots,x_{i,T_i}$ are drawn from the same Gaussian model (certain $l\in [k]$). To handle multiple observations, we show that although each $\psi_i$ consists of $T_i$ sub-functions $\psi_{it}$, it can be approximated by a single function on the average observation $b_i = \frac{\sum_{t\in [T_i]}x_{it}}{T_i}$, specifically, we have $\psi_i(\mu, I_d, 0_d) = d(b_i, \mu)^2 + O(1)$. This property enables us to "reduce" the representative complexity of $\psi_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 $b_i = \frac{\sum_{t\in [T_i]}x_{it}}{T_i}$ of entities (Definition 4.2) by upper bounding the maximum effect of covariances and autocorrelations and applying the fact that $\psi_i(\mu, I_d, 0_d) = d(b_i, \mu)^2 + O(1)$. 2) For the pseudo-dimension, we prove that there are only $\mathrm{poly}(k,d)$ intristic operators in $\psi_i$ between parameters $\theta$ and observations $x_{it}$, based on the reduction of the representative complexity of $\psi_i$. Generalizing our results in the same manner as feldman2019coresets does over feldman2011scalable,lucic2017training would be an interesting future direction since we may get rid of Assumption (ref). Currently, it is unclear how to generalize the approach of feldman2019coresets to time series data since feldman2019coresets assumes that each point is an integral point within a bounded box, while we consider a continuous GMM generative model (ref), and hence, each coordinate of an arbitrary observation $x_{itr}$ is drawn from a certain continuous GMM distribution which is not integral with probability $\approx 1$ and can be unbounded.
remark[Discussion on lower bounds] There is no provable lower bound result for our GMM coreset with time series data. We conjecture that without the first condition in Assumption (ref), the coreset size should depend exponentially in $k$ and logarithmic in $n$. The motivation is from a simple setting that all $T_i=1$ (GMM with static data), in which feldman2019coresets reduces the problem to projective clustering whose coreset size depends exponentially in $k$ and logarithmic in $n$. Moreover, feldman2019coresets believes that these dependencies are unavoidable for GMM coreset.

Empirical results

We compare the performance of our coreset algorithm $\ensuremath{\mathsf{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.

\paragraph{Datasets.} We generate two sets of synthetic data all with 250K observations with different number of entities $N$ and observations per individual $T_i$: (i) $N=T_i = 500$ for all $i\in [N]$ and (ii) $N=200$, $T_i = 1250$ for all $i\in [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 $T_i$).

We fix $d=2$, $k=3$, $\lambda = 0.01$ and generate multiple datasets with model parameters randomly varied as follows: (i) Draw $\alpha$ from a uniform distribution over $\Delta_k$. (ii) For each $l\in [k]$, draw $\mu^{(l)}\in \mathbb{R}^d$ from $N(0, I_d)$. (iii) Draw $\Sigma^{(l)}$ by first generating a random matrix $A\in [0,1]^{d\times d}$ and then let $\Sigma^{(l)} = (AA^\top)^{-1}$, and (iv) draw all diagonal elements of $\Lambda^{(l)}\in {\mathcal{D}}^d$ from a uniform distribution over [0, $1-\sqrt{\lambda}$]. Given these draws of parameters, we generate a GMM time-series dataset as follows: For each $i\in [N]$, draw $l\in [k]$ given $\alpha$. Then, for all $t\in [T_i]$ draw $e_{it}\in \mathbb{R}^d$ with covariance matrix $\Sigma^{(l)}$ and autocorrelation matrix $\Lambda^{(l)}$ and compute $x_{it}=\mu^{(l)} + e_{it}\in \mathbb{R}^{d}$.

\paragraph{Baseline and metrics.} We use coresets based on uniform sampling (Uni) and based on lucic2017training (LFKF) that constructs a coreset for Gaussian mixture model with static data as the baseline. Given an integer $\Gamma$, uniformly sample a collection $S$ of $\Gamma$ entity-time pairs $(i, t) \in P_{\mathcal{X}}$, let $w(i) = \frac{N}{|I_S|}$ for each $i\in I_S$, and let $w^{(i)}(t) = \frac{T_i}{|J_{S,i}|}$ for each $t\in J_{S,i}$; and LFKF regard all entity-time pairs as independent points and sample a collection $S$ of $\Gamma$ entity-time pairs $(i, t) \in P_{\mathcal{X}}$ via importance sampling. For comparability, we set $\Gamma$ 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^\star$ (negative log-likelihood) denote the objective of Problem (ref) with full data. Given a weighted subset $S\subseteq [N]\times [T]$, we first compute $(\alpha^\star,\theta^\star):=\arg\min_{\alpha\in \Delta_k, \theta\in {\mathcal{P}}_\lambda^k} f'_S(\alpha'(\theta), \theta) + \phi(\alpha, \theta)$ as the maximum likelihood solution over $S$.\footnote{We solve this optimization problem using an EM algorithm similar to arcidiacono2003finite. The M step in each iteration is based on IRLS jorgensen2006iteratively and the E-step involves an individual level Bayesian update for $\alpha$.} Then $V^\star_S = f(\alpha^\star, \theta^\star)$ 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., $\gamma_S := 2(V^\star_S - V^\star)$ as a measure of the quality of $S$ buse1982likelihood. The running time for GMM clustering with full data ($V^\star$) and coreset ($V^\star_S$) are $T_{\mathcal{X}}$ and $T_S$ respectively. $T_C$ is the coreset $S$ construction time.

\paragraph{Empirical setup.} We vary $\varepsilon = 0.1, 0.2, 0.3, 0.4, 0.5$. For each $\varepsilon$, we run $\ensuremath{\mathsf{CRGMM}}$ and Uni to generate 5 coresets each. We estimate the model with the full dataset ${\mathcal{X}}$ and the coresets $S$ and record $V^\star$, $V^\star_S$ and the run times.

\FloatBarrier

table[table omitted — 2,622 chars of source]

\paragraph{Results.} Table (ref) summarizes the quality-size-time trade-offs of our coresets for different error guarantees $\varepsilon$. The $V^\star_S$ achieved by our coreset is always close to the optimal value $V^\star$ and is always smaller than that of Uni and LFKF. We next assess coreset size. For the $N=T_i=500$ case, the coreset achieves a likelihood ratio $\gamma_S=684$ with only 72 entity-time pairs (0.03%), but Uni/LFKF achieves the closest (but worse) likelihood ratio of $\gamma_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 (ref), we see that not only is the quality of our coreset ($\gamma_S$) superior, it has lower standard deviations suggesting lower variance in performance. In terms of total computation time relative to the full data ($\frac{T_{\mathcal{X}}}{T_C+T_S}$), our coresets speed up by 14x-171x. Also, the computation time ($T_C+T_S$) of Uni and \textbf{LFKF} is always larger than that of our coreset (3x-22x).\footnote{A possible explanation is that our coreset selects $I_S$ of entities that are more representative that \textbf{Uni}, and hence, achieves a better convergence speed.} Finally, $V^\star_S$ is lower for Synthetic 2 given that there are fewer entities and cross-entity variance is greater than within entity variance. From Figure (ref), 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 ($T_i \gg N$)---a feature of sensor based measurement that motivates this paper.

figure[figure omitted — 586 chars of source]

Proof of Lemma (ref): The first stage of Algorithm (ref) outputs an entity-level coreset

For preparation, we first introduce an importance sampling framework for coreset construction, called the Feldman-Langberg\ framework feldman2011unified,braverman2016new.

The Feldman-Langberg\ framework

We first give the definition of query space and the corresponding coresets.

definition[\bf{Query space feldman2011unified,braverman2016new}] Let ${\mathcal{X}}$ be a finite set together with a weight function $u: {\mathcal{X}}\rightarrow \mathbb{R}_{\geq 0}$. Let ${\mathcal{P}}$ be a set called queries, and $f_x:{\mathcal{P}}\rightarrow \mathbb{R}_{\geq 0}$ be a given loss function w.r.t. $x\in {\mathcal{X}}$. The total cost of ${\mathcal{X}}$ with respect to a query $\theta\in {\mathcal{P}}$ is $ f(\theta) := \sum_{x\in {\mathcal{X}}} u(x)\cdot f_x(\theta). $ The tuple $({\mathcal{X}},u,{\mathcal{P}},f)$ is called a query space. Specifically, if $u(x)=1$ for all $x\in {\mathcal{X}}$, we use $({\mathcal{X}},{\mathcal{P}},f)$ for simplicity.

Intuitively, $f$ represents a linear combination of weighted functions indexed by ${\mathcal{X}}$, and ${\mathcal{P}}$ represents the ground set of $f$. Due to the separability of $f$, we have the following coreset definition.

definition[\bf{Coresets of a query space feldman2011unified,braverman2016new}] Let $({\mathcal{X}},u,{\mathcal{P}},f)$ be a query space and $\varepsilon \in (0,1)$ be an error parameter. An $\varepsilon$-coreset of $({\mathcal{X}},u,{\mathcal{P}},f)$ is a weighted set $S\subseteq {\mathcal{X}}$ together with a weight function $w: S\rightarrow \mathbb{R}_{\geq 0}$ such that for any $\theta\in {\mathcal{P}}$, $ \sum_{x\in S} w(x)\cdot f_x(\theta) \in (1\pm \varepsilon)\cdot f(\theta). $
remarkFor instance, we set ${\mathcal{X}} = [N]$, $u=1$, ${\mathcal{P}} = \Delta_k\times {\mathcal{P}}_\lambda^k$ and $f=f'$ in the GMM time-series clustering problem. Then by Definition (ref), Lemma (ref) represents that $(I_S, w)$ is an $\varepsilon$-coreset for the query space $([N], 1, \Delta_k\times {\mathcal{P}}_\lambda^k, f')$. Another example is to set ${\mathcal{X}} = [T_i]$, $u = 1$, ${\mathcal{P}} = {\mathcal{P}}_\lambda$ and $f = \psi_i$. Then Lemma (ref) represents that $(J_{S,i}, w^{(i)})$ is an $\varepsilon$-coreset for the query space $([T_i], 1, {\mathcal{P}}_\lambda, \psi_i)$.

Now we are ready to give the Feldman-Langberg\ framework.

\paragraph{The Feldman-Langberg framework.} Feldman and Langberg feldman2011unified show how to construct coresets by importance sampling and the coreset size has been improved by braverman2016new. For preparation, we first give the notion of sensitivity which measures the maximum influence for each point $x\in {\mathcal{X}}$.

definition[\bf{Sensitivity feldman2011unified,braverman2016new}] Given a query space $({\mathcal{X}},u,{\mathcal{P}},f)$, the sensitivity of a point $x\in {\mathcal{X}}$ is $ s(x):= \sup_{\theta\in {\mathcal{P}}} \frac{u(x)\cdot f_x(\theta)}{f({\mathcal{X}},u,\theta)}. $ The total sensitivity of the query space is $\sum_{x\in {\mathcal{X}}} s(x)$.

We also introduce a notion which measures the combinatorial complexity of a query space.

definition[\bf{Pseudo-dimension feldman2011unified,braverman2016new}] For a query space $({\mathcal{X}},u,{\mathcal{P}},f)$, we define $ \ensuremath{\mathsf{range}}(\theta, r) = \left\{x\in {\mathcal{X}}: u(x)\cdot f_x(\theta)\leq r\right\} $ for every $\theta\in {\mathcal{P}}$ and $r\geq 0$. The (pseudo-)dimension of $({\mathcal{X}},u,{\mathcal{P}},f)$ is the largest integer $t$ such that there exists a subset $A\subseteq {\mathcal{X}}$ of size $t$ satisfying that $ |\left\{A\cap \ensuremath{\mathsf{range}}(\theta,r): \theta\in {\mathcal{P}}, r\geq 0\right\}| = 2^{|A|}. $

Pseudo-dimension plays the same role as VC-dimension vapnik1971uniform. Specifically, if the range of $f$ is $\left\{0,1\right\}$ and $u=1$, pseudo-dimension can be regarded as a generalization of VC-dimension to function spaces. Now we are ready to describe the Feldman-Langberg framework.

theorem[\bf{Feldman-Langberg framework feldman2011unified,braverman2016new}] Let $({\mathcal{X}},u,{\mathcal{P}},f)$ be a given query space and $\varepsilon,\delta\in (0,1)$. Let $\dim$ be an upper bound of the pseudo-dimension of every query space $({\mathcal{X}},u',{\mathcal{P}},f)$ over $u'$. Suppose $s:{\mathcal{X}}\rightarrow \mathbb{R}_{\geq 0}$ is a function satisfying that for any $x\in {\mathcal{X}}$, $ s(x) \geq \sup_{\theta\in {\mathcal{P}}} \frac{u(x)\cdot f_x(\theta)}{f({\mathcal{X}},u,\theta)}, $ and define ${\mathcal{G}} := \sum_{x\in {\mathcal{X}}} s(x)$ to be the total sensitivity. Let $S\subseteq {\mathcal{X}}$ be constructed by taking $ O\left(\varepsilon^{-2} {\mathcal{G}} (\dim \cdot \ln {\mathcal{G}} +\ln(1/\delta))\right) $ samples, where each sample $x\in {\mathcal{X}}$ is selected with probability $\frac{s(x)}{{\mathcal{G}}}$ and has weight $w(x):= \frac{{\mathcal{G}}}{|S|\cdot s(x)}$. Then, with probability at least $1-\delta$, $S$ is an $\varepsilon$-coreset of $({\mathcal{X}},u,{\mathcal{P}},f)$.

Bounding the pseudo-dimension of $f'$

Our proof idea is similar to that in lucic2017training. For preparation, we need the following lemma which is proposed to bound the pseudo-dimension of feed-forward neural networks.

lemma[\bf{Restatement of Theorem 8.14 of anthony2009neural}] Let $({\mathcal{X}},u,{\mathcal{P}},f)$ be a given query space where $f_x(\theta)\in \left\{0,1\right\}$ for any $x\in {\mathcal{X}}$ and $\theta\in {\mathcal{P}}$, and ${\mathcal{P}}\subseteq \mathbb{R}^{m}$. Suppose that $f$ can be computed by an algorithm that takes as input the pair $(x,\theta)\in {\mathcal{X}}\times {\mathcal{P}}$ and returns $f_x(\theta)$ after no more than $l$ of the following operations: \begin{itemize} • the exponent function $a \rightarrow e^a$ on real numbers. • the arithmetic operations $+,-,\times$, and $/$ on real numbers. • jumps conditioned on $>,\geq,<,\leq,=$, and $\neq$ comparisons of real numbers, and • output 0,1. \end{itemize} If the $l$ operations include no more than $q$ in which the exponential function is evaluated, then the pseudo-dimension of $({\mathcal{X}},u,{\mathcal{P}},f)$ is at most $O(m^2q^2 + mq(l+\ln mq))$.

Note that the above lemma requires that the range of functions $f_x$ is $[0,1]$. We have the following lemma which can help extend this range to $\mathbb{R}$.

lemma[\bf{Restatement of Lemma 4.1 of vidyasagar2002theory}] Let $({\mathcal{X}},u,{\mathcal{P}},f)$ be a given query space. Let $g_x:{\mathcal{P}}\times \mathbb{R}\rightarrow \left\{0,1\right\}$ be the indicator function satisfying that for any $x\in {\mathcal{X}}$, $\theta\in {\mathcal{P}}$ and $r\in \mathbb{R}$, \[ g_x(\theta,r) = I\left[u(x)\cdot f(x,\theta)\geq r\right]. \] Then the pseudo-dimension of $({\mathcal{X}},u,{\mathcal{P}},f)$ is precisely the pseudo-dimension of the query space $({\mathcal{X}},u,{\mathcal{P}}\times \mathbb{R},g_f)$.

Now we are ready to prove bound the pseudo-dimension of $f'$ by the following lemma.

lemma[\bf{Pseudo-dimension of $f'$}] The pseudo-dimension of $([N], u, \Delta_k\times {\mathcal{P}}_\lambda^k, f')$ over weight functions $u: [N]\rightarrow \mathbb{R}_{\geq 0}$ is at most $O(k^4 d^4 + k^3 d^8)$.
proofOur argument is similar to that in huang2020coresets. Fix a weight function $u: [N]\rightarrow \mathbb{R}_{\geq 0}$. We only need to consider the following indicator function $g_i: \Delta_k\times{\mathcal{P}}_\lambda^k\times \mathbb{R}_{\geq 0}\rightarrow \left\{0,1\right\}$ where for any $\alpha\in \Delta_k$, $\theta\in {\mathcal{P}}_\lambda^k$ and $r\in \mathbb{R}_{\geq 0}$, \[ g_i(\alpha,\theta,r):= I\left[\sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i}\psi_i(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)})) \geq r\right]. \] Note that the parameter space is $\Delta_k\times {\mathcal{P}}_\lambda^k$ which consists of at most $m=O(kd^2)$ parameters. For any $(\mu,\Sigma,\Lambda)\in {\mathcal{P}}_\lambda$, function $\psi_i(\mu, \Sigma, \Lambda)$ can be represented as a multivariate polynomial that consists of $O(d^6)$ terms $\mu_{c_1}^{b_1} \mu_{c_2}^{b_2} \Lambda_{c_3,c_3}^{b_3}\Lambda_{c_4,c_4}^{b_4}\left(\Sigma^{-1}\right)_{c_5, c_6}^{b_5}$ where $c_1,c_2,c_3,c_4,c_5,c_6\in [d]$, and $b_1,b_2,b_3,b_4,b_5\in \left\{0,1\right\}$. Thus, $g_i$ consists of $l=O(kd^6)$ arithmetic operations, $q=k$ exponential functions, and $k$ jumps. By Lemmas (ref) and (ref), we complete the proof.

Bounding the total sensitivity of $f'$

Next, we prove that function $s$ (Line 6 of Algorithm (ref)) is a sensitivity function w.r.t. $f'$; summarized as follows.

lemma[\bf{$s$ is a sensitivity function w.r.t. $f'$}] For each $i\in [N]$, we have \[ s(i) \geq \max_{\alpha\in \Delta_k, \theta \in {\mathcal{P}}_\lambda^k} \frac{f'_i(\alpha, \theta)}{f'(\alpha,\theta)}. \] Moreover, $\sum_{i\in [N]} s(i) \leq (16D+12Dk)/\lambda$.

To prove the lemma, we will use a reduction from general $\Sigma$ to $I_d$ and from $\Lambda$ to $0_d$ (without both covariances and autocorrelations), which upper bounds the affect of the covariance matrix and autocorrelation matrices. We define $\psi^{(O)}_i: \mathbb{R}^d \rightarrow \mathbb{R}_{\geq 0}$ to be

eqnarray[eqnarray omitted — 89 chars of source]

for any $\mu\in \mathbb{R}^d$, and define $f^{(O)}: \Delta_k\times \mathbb{R}^{d\times k} \rightarrow \mathbb{R}_{\geq 0}$ to be

eqnarray*[eqnarray* omitted — 291 chars of source]

for any $\alpha\in \Delta_k$ and $\theta^{(O)} = (\mu^{(l)})_{l\in [k]}\in \mathbb{R}^{d\times k}$. Compared to $f'$, we note that $f^{(O)}$ does not contain covariance and autocorrelation matrices.

\paragraph{Clustering cost of entities.} By the definition of $\psi_i^{(O)}$, we have that for any $\mu\in \mathbb{R}^d$,

eqnarray[eqnarray omitted — 223 chars of source]

Next, we introduce another function $s^{(O)}: [N]\rightarrow \mathbb{R}_{\geq 0}$ as a sensitivity function w.r.t. $f^{(O)}$, i.e., for any $i\in [N]$, \[ s^{(O)}(i) := \max_{\alpha\in \Delta_k, \theta^{(O)}\in \mathbb{R}^{d\times k}} \frac{f^{(O)}_i(\alpha, \theta^{(O)})}{f^{(O)}(\alpha, \theta^{(O)})}. \] Define ${\mathcal{G}}^{(O)} := \sum_{i\in [N]} s^{(O)}(i)$ to be the total sensitivity w.r.t. $f^{(O)}$. We first have the following lemma.

lemma[\bf{Relation between sensitivities w.r.t. $f_i^{(O)}$ and $f'_i$}] For each $i\in [N]$, we have \[ s^{(O)}(i) \leq \max_{\alpha\in \Delta_k, \theta \in {\mathcal{P}}_\lambda^k} \frac{f'_i(\alpha, \theta)}{f'(\alpha,\theta)} \leq 4D\cdot s^{(O)}(i)/\lambda. \]
proofIt is easy to verify $s^{(O)}(i) \leq \max_{\alpha\in \Delta_k, \theta \in {\mathcal{P}}_\lambda^k} \frac{f'_i(\alpha, \theta)}{f'(\alpha,\theta)}$ since \begin{align*} & && \max_{\alpha\in \Delta_k, \theta \in {\mathcal{P}}_\lambda^k} \frac{f'_i(\alpha, \theta)}{f'(\alpha,\theta)} & \\ & \geq && \max_{\alpha\in \Delta_k, \theta \in (\mathbb{R}^d \times I_d\times 0_d)^k} \frac{f'_i(\alpha, \theta)}{f'(\alpha,\theta)} & (0_d \in {\mathcal{D}}^d_{\lambda}, I_d\in {\mathcal{S}}^d) \\ & = && \max_{\alpha\in \Delta_k, \theta^{(O)}\in \mathbb{R}^{d\times k}} \frac{f^{(O)}_i(\alpha, \theta^{(O)})}{f^{(O)}(\alpha, \theta^{(O)})} & (Defn. of $f_i^{(O)}$)\\ & = && s^{(O)}(i). \end{align*} For the other side, we have the following claim that for any $i\in [N]$ and $\theta = (\mu,\Sigma, \Lambda)\in {\mathcal{P}}_\lambda$, \begin{eqnarray} \frac{\lambda}{\lambda_{\max}(\Sigma)}\cdot \psi^{(O)}_i(\mu)\leq \psi_i(\mu, \Sigma,\Lambda) \leq \frac{4}{\lambda_{\min}(\Sigma)}\cdot \psi^{(O)}_i(\mu). \end{eqnarray} Then for any $\alpha\in \Delta_k$ and $\theta=(\mu^{(l)},\Sigma^{(l)},\Lambda^{(l)})_{l\in [k]}\in {\mathcal{P}}_\lambda^k$, letting $\theta^{(O)} = (\mu^{(l)})_{l\in [k]}\in \mathbb{R}^{d\times k}$ and $\beta = \min_{l'\in [k]}\lambda_{\min}(\Sigma^{(l')})^2$, we have \begin{align*} \frac{f'_i(\alpha, \theta)}{f'(\alpha,\theta)} & = && \frac{-\ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i}\psi_i(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)}))}{-\sum_{j\in [N]} \ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i}\psi_j(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)}))} & (by definition)\\ & \leq && \frac{-\ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{2}{T_i \lambda_{\min}(\Sigma^{(l)})}\psi^{(O)}_i(\mu^{(l)}))}{-\sum_{j\in [N]} \ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{\lambda}{2T_i \lambda_{\max}(\Sigma^{(l)})}\cdot \psi^{(O)}_j(\mu^{(l)}))} & (Ineq. (ref)) \\ & \leq && \frac{4f^{(O)}_i(\alpha, \theta^{(O)})}{-\sum_{j\in [N]} \ln \left(\sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i\cdot \beta}\cdot \psi^{(O)}_j(\mu^{(l)}))\right)^{\frac{\lambda\beta}{ \lambda_{\max}(\Sigma^{(l)})}}} & (Defn. of $f_i^{(O)}$) \\ & \leq && \frac{4f^{(O)}_i(\alpha, \theta^{(O)})}{-\sum_{j\in [N]} \ln \left(\sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i\cdot \beta}\cdot \psi^{(O)}_j(\mu^{(l)}))\right)^{\lambda/D}} & (Assumption (ref)) \\ & = && \frac{4D\cdot f^{(O)}_i(\alpha, \theta^{(O)})}{\lambda\cdot f^{(O)}(\alpha, \theta^{(O)})}. & (by definition) \end{align*} Consequently, we have $\max_{\alpha\in \Delta_k, \theta \in {\mathcal{P}}_\lambda^k} \frac{f'_i(\alpha, \theta)}{f'(\alpha,\theta)}\leq 4D\cdot s^{(O)}(i)/\lambda$, which completes the proof. \paragraph{Proof of Claim (ref).} It remains to prove Claim (ref). \sloppy For any $i\in [N]$ and $\theta = (\mu,\Sigma, \Lambda)\in {\mathcal{P}}_\lambda$, we have \begin{eqnarray*} \begin{split} & && \psi_i(\mu, \Sigma, \Lambda) \\ & = && (x_{i,1}-\mu)^\top \Sigma^{-1} (x_{i,1}-\mu)- \left(\Lambda(x_{i,1}-\mu)\right)^\top \Sigma^{-1} \left(\Lambda(x_{i,1}-\mu)\right) \\ & && + \sum_{t=2}^{T_i} \left((x_{it}-\mu)- \Lambda(x_{i,t-1}-\mu)\right)^\top \Sigma^{-1} \left((x_{it}-\mu)-\Lambda(x_{i,t-1}-\mu)\right) \\ & \in && [\frac{1}{\lambda_{\max}(\Sigma)},\frac{1}{\lambda_{\min}(\Sigma)}]\cdot \left((x_{i,1}-\mu)^\top (x_{i,1}-\mu)- \left(\Lambda(x_{i,1}-\mu)\right)^\top \left(\Lambda(x_{i,1}-\mu)\right)\right) \\ & && + [\frac{1}{\lambda_{\max}(\Sigma)},\frac{1}{\lambda_{\min}(\Sigma)}]\cdot \sum_{t=2}^{T_i} \left((x_{it}-\mu)- \Lambda(x_{i,t-1}-\mu)\right)^\top \left((x_{it}-\mu)-\Lambda(x_{i,t-1}-\mu)\right). \end{split} \end{eqnarray*} Hence, it suffices to prove that \begin{eqnarray*} \begin{split} & && (x_{i,1}-\mu)^\top \Sigma^{-1} (x_{i,1}-\mu)- \left(\Lambda(x_{i,1}-\mu)\right)^\top \Sigma^{-1} \left(\Lambda(x_{i,1}-\mu)\right) \\ & && + \sum_{t=2}^{T_i} \left((x_{it}-\mu)- \Lambda(x_{i,t-1}-\mu)\right)^\top \left((x_{it}-\mu)-\Lambda(x_{i,t-1}-\mu)\right) \\ &\in && [\lambda, 4]\cdot \psi_i^{(O)}(\mu) \end{split} \end{eqnarray*} Since $\Lambda\in {\mathcal{D}}^d_\lambda$, we suppose $\Lambda = \left(\Delta_1, \ldots, \Delta_d\right)$. Then we have \begin{eqnarray*} \begin{split} & && (x_{i,1}-\mu)^\top (x_{i,1}-\mu)- \left(\Lambda(x_{i,1}-\mu)\right)^\top \left(\Lambda(x_{i,1}-\mu)\right) \\ & && + \sum_{t=2}^{T_i} \left((x_{it}-\mu)- \Lambda(x_{i,t-1}-\mu)\right)^\top \left((x_{it}-\mu)-\Lambda(x_{i,t-1}-\mu)\right) \\ & = && \sum_{r\in [d]} (1-\Delta_r^2)(x_{i1r}-\mu_r)^2 \\ & && + \sum_{t=2}^{T_i} \left((x_{itr}-\mu_r)- \Delta_r(x_{i,t-1,r}-\mu_r)\right)^\top \left((x_{itr}-\mu_r)-\Delta_r(x_{i,t-1,r}-\mu_r)\right). \end{split} \end{eqnarray*} On one hand, we have \begin{eqnarray*} \begin{split} & && (1-\Delta_r^2)(x_{i1r}-\mu_r)^2 + \sum_{t=2}^{T_i} \left((x_{itr}-\mu_r)- \Delta_r(x_{i,t-1,r}-\mu_r)\right)^2 & \\ & \leq && (1-\Delta_r^2)(x_{i1r}-\mu_r)^2 + \sum_{t=2}^{T_i} 2\left((x_{itr}-\mu_r)^2 + \Delta_r^2(x_{i,t-1,r}-\mu_r)^2\right) & (\text{Arithmetic Ineq.}) \\ & \leq && 4 \sum_{t\in [T_i]} (x_{itr}-\mu_r)^2. & (\Delta_r\leq 1) \end{split} \end{eqnarray*} On the other hand, we have \begin{eqnarray*} \begin{split} & && (1-\Delta_r^2)(x_{i1r}-\mu_r)^2 + \sum_{t=2}^{T_i} \left((x_{itr}-\mu_r)- \Delta_r(x_{i,t-1,r}-\mu_r)\right)^2 &\\ & = && (x_{i1r}-\mu_r)^2 + \sum_{t=2}^{T_i} (1+\Delta_r^2)(x_{itr}-\mu_r)^2 - 2\Delta_r(x_{itr}-\mu_r)(x_{i,t-1,r}-\mu_r) &\\ & \geq && (1-\Delta_r)^2 \sum_{t\in [T_i]} (x_{itr}-\mu_r)^2 & (\text{Arithmetic Ineq.}) \\ & \geq && \lambda\sum_{t\in [T_i]} (x_{itr}-\mu_r)^2. & (\text{Assumption (ref)}) \end{split} \end{eqnarray*} This completes the proof.

By the definition of $s$ and the above lemma, it suffices to prove the following lemma that provides an upper bound for $s^{(O)}$.

lemma[\bf{Sensitivities w.r.t. $f^{(O)}$}] The following holds: \begin{enumerate} • For each $i\in [N]$, we have \[ s^{(O)}(i) \leq \frac{4 \|b_i - c^\star_{p(i)}\|_2^2}{\ensuremath{\mathsf{OPT}}^{(O)}+A} + 3 s^c(i) \] • ${\mathcal{G}}^{(O)} \leq 4+3k$. \end{enumerate}

For preparation, we introduce some notations related to the clustering problem (Definition (ref)). For any $\mu\in \mathbb{R}^d$, we define \[ h^c_i(\mu) := \|c^\star_{p(i)}-\mu\|_2^2, \] and for any $\alpha\in \Delta_k, \theta^{(O)}\in \mathbb{R}^{d\times k}$, \[ f^c_i(\alpha,\theta^{(O)}) := -\ln \sum_{l\in [k]} \alpha_k \cdot \exp(-\frac{1}{2T_i\cdot \beta}h^c_i(\mu^{(l)})), \] where $\beta = \min_{l'\in [k]}\lambda_{\min}(\Sigma^{(l')})^2$. Let $f^c := \sum_{i\in [N]} f^c_i$. Then similarly, we can prove that $s^c$ (Line 5 of Algorithm (ref)) is a sensitivity function w.r.t. $f^c$; summarized as follows.

lemma[\bf{$s^c$ is a sensitivity function w.r.t. $f^c$}] For each $i\in [N]$, \[ s^c(i) \geq \max_{\alpha\in \Delta_k, \theta^{(O)}} \frac{f^c_i(\alpha, \theta^{(O)})}{f^c(\alpha, \theta^{(O)})}. \] Moreover, $\sum_{i\in [N]} s^c(i) \leq k$.
proofThis lemma is a direct corollary by the fact that there are only $k$ different centers $C^\star_i$, which implies that there are at most $k$ different functions $f^c_i$ accordingly. We partition $[N]$ into at most $k$ groups $A_l$ where each element $i\in A_l$ satisfies that $c^\star_{p(i)}= l$. Then we observe that $f^c_i = f^c_j$ if $i,j\in A_l$. Then for any $\alpha, \theta^{(O)}\in \Delta_k\times \mathbb{R}^{d\times k}$, \[ \frac{f^c_i(\alpha, \theta^{(O)})}{f^c(\alpha, \theta^{(O)})} \leq \frac{f^c_i(\alpha, \theta^{(O)})}{\sum_{j\in A_l} f^c_j(\alpha, \theta^{(O)})} = \frac{1}{|A_l|} \leq s^c(i), \] which implies the lemma.

To prove Lemma (ref), the main idea is to relate $s^{(O)}(i)$ to $s^c(i)$. The idea is similar to lucic2017training. For preparation, we also need the following key observation.

lemma[\bf{Upper bounding the projection cost}] For a fixed number $L>0$ and a fixed $\theta = (\alpha, \theta^{(O)})\in \Delta_k\times \mathbb{R}^{d\times k}$ and a fixed value $a\geq 0$, define $\pi_{a,\theta}: \mathbb{R}^d \rightarrow \mathbb{R}_{\geq 0}$ as \[ \pi_{a,\theta, L}(y) = -\ln \sum_{l\in [k]} \alpha_l \cdot \exp\left(-\frac{1}{2L}\left(\|y - \mu^{(l)}\|_2^2 + a\right)\right). \] Then, for every $y, y'\in \mathbb{R}^d$ it holds that \[ \pi_{a,\theta}(y) \leq \frac{2}{L}\|y-y'\|_2^2 + \pi_{a,\theta, L}(y'). \]
proofUse the relaxed triangle inequality for $l_2^2$-norm, we have \begin{align*} \pi_{a,\theta, L}(y) & = && -\ln \sum_{l\in [k]} \alpha_l \cdot \exp\left(-\frac{1}{2L}\left(\|y-\mu^{(l)}\|_2^2 + a\right)\right) \\ & \leq && -\ln \sum_{l\in [k]} \alpha_l \cdot \exp\left(-\frac{2}{L}\|y'-\mu^{(l)}\|_2^2+ \frac{1}{2L}\left(\|y-y'\|_2^2 + a\right) \right) \\ & && (relaxed triangle ineq.) \\ & \leq && -\ln \left(\exp\left(-\frac{2}{L}\|y-y'\|_2^2\right)\cdot \sum_{l\in [k]} \alpha_l \cdot \exp\left(-\frac{1}{2L}\left(\|y'-\mu^{(l)}\|_2^2 +a \right) \right) \right) \\ & \leq && \frac{2}{L} \|y-y'\|_2^2 -\ln \sum_{l\in [k]} \alpha_l \cdot \exp\left(-\frac{1}{2L}\left(\|y'-\mu^{(l)}\|_2^2 +a \right) \right) \\ & \leq && \frac{2}{L}\|y-y'\|_2^2 -\ln \left(\sum_{l\in [k]} \alpha_l \cdot \exp\left(-\frac{1}{2L}\left(\|y'-\mu^{(l)}\|_2^2 + a\right)\right) \right)^2 \\ & = && \frac{2}{L}\|y-y'\|_2^2 + \pi_{a,\theta, L}(y'). \end{align*}

Recall that $b_i\leftarrow\frac{\sum_{t\in [T_i]}x_{it}}{T_i}$. Now we are ready to prove Lemma (ref).

proof[Proof of Lemma (ref)] For each $i\in [N]$ and $\theta = (\alpha, \theta^{(O)})\in \Delta_k\times \mathbb{R}^{d\times k}$, letting $\beta = \min_{l'\in [k]}\lambda_{\min}(\Sigma^{(l')})^2$, we have \begin{eqnarray} \begin{split} & && f^{(O)}_i(\theta) &\\ & = && \pi_{a_i,\theta, \beta}(b_i) &\\ & \leq && \frac{2}{\beta} \|b_i-c^\star_{p(i)}\|_2^2 + \pi_{a_i,\theta, \beta}(c^\star_{p(i)}) & (Lemma (ref))\\ & \leq && \frac{2}{\beta} \|b_i-c^\star_{p(i)}\|_2^2 + s^c(i) \cdot \sum_{j\in [N]} \pi_{a_i,\theta, \beta}(c^\star_{p(j)}) & (Defns. of $s^c$) \\ & \leq && \frac{2}{\beta} \|b_i-c^\star_{p(i)}\|_2^2 + s^c(i) \cdot \sum_{j\in [N]} \left(\frac{2}{\beta}\cdot \|b_j-c^\star_{p(j)}\|_2^2 + \pi_{a_i,\theta, \beta}(b_j)\right) & (Lemma (ref)) &\\ & \leq && \frac{2}{\beta} \|b_i-c^\star_{p(i)}\|_2^2 + s^c(i) \cdot ( \frac{2\cdot \ensuremath{\mathsf{OPT}}^{(O)}}{\beta} + f^{(O)}(\theta)).& (Defns. of $\ensuremath{\mathsf{OPT}}^{(O)}$) \end{split} \end{eqnarray} Let $\ensuremath{\mathsf{OPT}}:=\min_{\theta} \sum_{i\in [N]} f_i^{(O)}(\theta)$. We have that \begin{eqnarray} \begin{split} & && \ensuremath{\mathsf{OPT}} \\ & = && \min_{\theta} \sum_{i\in [N]} f_i^{(O)}(\theta) &\\ & = && \min_{\theta} \sum_{i\in [N]} -\ln \sum_{l\in [k]} \alpha_l \cdot \exp\left(-\frac{1}{2 \beta}\left(\|b_i - \mu^{(l)}\|_2^2 +a_i\right)\right) &\\ & \geq && \min_{\theta} \sum_{i\in [N]} -\ln \sum_{l\in [k]} \alpha_l \cdot \exp\left(-\frac{1}{2 \beta} \left(\min_{l'\in [k]}\|b_i - \mu^{(l')}\|_2^2 +a_i\right)\right) &\\ & \geq && \min_{\theta} \sum_{i\in [N]} -\ln \sum_{l\in [k]} \alpha_l \cdot \exp\left(-\frac{1}{2 \beta} \left(\min_{l'\in [k]}\|b_i - \mu^{(l')}\|_2^2 + a_i\right)\right) &\\ & = && \frac{1}{2\beta} \left(\min_{\theta} \sum_{i\in [N]} \min_{l'\in [k]}\|b_i - \mu^{(l')}\|_2^2 +a_i\right) &\\ & \geq && \frac{1}{2\beta}\cdot \left(\ensuremath{\mathsf{OPT}}^{(O)} + A\right). & (Defns. of $\ensuremath{\mathsf{OPT}}^{(O)}$ and $A$) \end{split} \end{eqnarray} Hence, we have \begin{align*} s^{(O)}(i) & = && \max_{\theta \in \Delta_k\times \mathbb{R}^{d\times k}} \frac{f^{(O)}_i(\theta)}{f^{(O)}(\theta)} & \\ & \leq && \frac{ 2\cdot \|b_i-c^\star_{p(i)}\|_2^2}{\beta \cdot \ensuremath{\mathsf{OPT}}} + \frac{ s^c(i)\cdot \ensuremath{\mathsf{OPT}}^{(O)}}{\beta \cdot \ensuremath{\mathsf{OPT}}} + s^c(i) & (Ineq. (ref))\\ & \leq && \frac{4 \cdot \|b_i-c^\star_{p(i)}\|_2^2}{\ensuremath{\mathsf{OPT}}^{(O)}+A} + 3 s^c(i). & (\text{Ineq. (ref)}) \end{align*} The second property is a direct conclusion.

Now we are ready to prove Lemma (ref).

proof[Proof of Lemma (ref)] Lemma (ref) is a direct corollary of Lemmas (ref) and (ref) since for each $i\in [N]$, \begin{eqnarray*} \begin{split} \max_{\alpha\in \Delta_k, \theta \in {\mathcal{P}}_\lambda^k} \frac{f'_i(\alpha, \theta)}{f'(\alpha,\theta)} & \leq && s^{(O)}(i)/\lambda & (Lemma (ref)) \\ & \leq && 4D\left(\frac{4 \cdot \|b_i-c^\star_{p(i)}\|_2^2}{\ensuremath{\mathsf{OPT}}^{(O)} + A} + 3 s^c(i)\right)/\lambda & (Lemma (ref)) \\ & = && s(i). & (Line 6 of Algorithm (ref)) \end{split} \end{eqnarray*}

By the Feldman-Langberg\ framework (Theorem (ref)), we note that Lemma (ref) is a direct corollary of Lemmas (ref) and (ref).

Proof of Lemma (ref): The second stage of Algorithm (ref) outputs a time-level coreset

The proof idea is similar to that in Lemma (ref), i.e., to bound the pseudo-dimension and the total sensitivity for the query space $([T_i], 1, {\mathcal{P}}_\lambda, \psi_i)$.

Bounding the pseudo-dimension of $\psi_i$

We have the following lemma.

lemma[\bf{Pseudo-dimension of $\psi_i$}] The pseudo-dimension of $([T_i], 1, {\mathcal{P}}_\lambda, \psi_i)$ over weight functions $u: [N]\rightarrow \mathbb{R}_{\geq 0}$ is at most $O(d^8)$.
proofThe argument is almost the same as in Lemma (ref). The parameter space of $\psi_i$ consists of at most $m=O(d^2)$ parameters and $\psi_i$ can be represented by at most $l=O(d^6)$ arithmetic operations. By Lemmas (ref) and (ref), it completes the proof.

Bounding the total sensitivity of $\psi_i$

Next, we again focus on proving $s_i$ (Line 12 of Algorithm (ref)) is a sensitivity function w.r.t. $\psi_i$; summarized as follows.

lemma[\bf{$s_i$ is a sensitivity function for $\psi_i$}] For each $i\in [N]$, we have that for each $t\in [T_i]$ \[ s_i(t) \geq \max_{\mu\in \mathbb{R}^d, \Lambda \in {\mathcal{P}}_{\tau,\lambda}} \frac{\psi_{it}(\mu, \Sigma,\Lambda)}{\psi_i(\mu,\Sigma,\Lambda)}. \] Moreover, $\sum_{t\in [T_i]} s_i(t) = O(D/\lambda)$.

Similar to Section (ref), we introduce another function $s^{(O)}_i: [T_i]\rightarrow \mathbb{R}_{\geq 0}$ as a sensitivity function w.r.t. $\psi^{(O)}_i$, i.e., for any $t\in [T_i]$, \[ s^{(O)}_i(t) := \max_{\mu\in \mathbb{R}^d} \frac{\psi^{(O)}_{it}(\mu)}{\psi^{(O)}_i(\mu)}. \] Define ${\mathcal{G}}^{(O)}_i := \sum_{t\in [T_i]} s^{(O)}_i(t)$ to be the total sensitivity w.r.t. $\psi^{(O)}_i$. We first have the following lemma, whose proof idea is simply from Lemma (ref) and huang2020coresets.

lemma[\bf{Relation between sensitivities w.r.t. $\psi_{it}^{(O)}$ and $\psi_{it}$}] For each $t\in [T_i]$, we have \[ s_i(t) \leq 4D\lambda^{-1}\cdot \left(s^{(O)}_i(t)+\sum_{j=1}^{\min\left\{t-1,1\right\}}s^{(O)}_i(t-j)\right). \]
proofBy the same argument as in Lemma (ref), we have that \[ s_i(t) \leq D\lambda^{-1}\cdot \max_{\mu\in \mathbb{R}^d, \Lambda\in {\mathcal{D}}^d_\lambda} \frac{\psi_{it}(\mu, I_d, \Lambda)}{\psi^{(O)}_i(\mu)}. \] By a similar argument as in huang2020coresets, we have that \[ \psi_{it}(\mu, I_d, \Lambda) \leq \psi^{(O)}_{it}(\mu) +\sum_{j=1}^{\min\left\{t-1,1\right\}}\psi^{(O)}_{i,t-q}(\mu). \] Combining the above two inequalities, we complete the proof.

Then we have the following lemma that relates $s^{(O)}_i$ and $s^c_i$ (Line 11 of Algorithm (ref)), whose proof follows from varadarajan2012sensitivity for the case that $k=1$.

lemma[\bf{Sensitivities w.r.t. $\psi^{(O)}_i$}] For each $i\in [N]$, the following holds: \begin{enumerate} • For each $t\in [T_i]$, we have \[ s^{(O)}_i(t) \leq s^c_i(t). \] • ${\mathcal{G}}^{(O)}_i \leq 8$. \end{enumerate}

Note that Lemma (ref) is a direct corollary of Lemmas (ref) and (ref).

Proof of Theorem (ref)

proofNote that the coreset size $|S|= ML$ matches the bound in Theorem (ref). We first prove the correctness. For any $i\in I_S$, we have \begin{align*} & && -\ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i}\sum_{t\in J_{S,i}} w^{(i)}(t)\cdot \psi_{it}(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)})) & \\ & \leq && -\ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i} (1 + \varepsilon)\cdot \psi_{i}(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)})) & (Lemma (ref))\\ & \leq && -\ln \sum_{l\in [k]} \left(\alpha_l \cdot \exp(-\frac{1}{2T_i} \psi_{i}(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)}))\right)^{1+\varepsilon} & ($x^{1+\varepsilon}$ is convex when $x\in [0,1]$)\\ & \leq && (1+ \varepsilon)\cdot \left(-\ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i}\psi_i(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)}))\right). & \end{align*} Symmetrically, we can also verify that \begin{align*} & && -\ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i}\sum_{t\in J_{S,i}} w^{(i)}(t)\cdot \psi_{it}(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)})) \\ &\geq && (1- \varepsilon)\cdot \left(-\ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i}\psi_i(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)}))\right). \end{align*} Consequently, we have \begin{eqnarray} \begin{split} & && -\ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i}\sum_{t\in J_{S,i}} w^{(i)}(t)\cdot \psi_{it}(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)})) \\ &\in && (1 \pm \varepsilon)\cdot \left(-\ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i}\psi_i(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)}))\right). \end{split} \end{eqnarray} Combining with Lemma (ref), we have that \begin{align*} & && f'_S(\alpha,\theta) & \\ & = && -\sum_{i\in I_S} w(i)\cdot \ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i} \sum_{t\in J_{S,i}} w^{(i)}(t) \cdot \psi_{it}(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)})) & \\ & \in && (1\pm \varepsilon)\cdot \left(-\sum_{i\in I_S} w(i)\cdot \ln \sum_{l\in [k]} \alpha_l \cdot \exp(-\frac{1}{2T_i} \psi_{i}(\mu^{(l)}, \Sigma^{(l)}, \Lambda^{(l)})) \right) & (Ineq. (ref)) \\ & \in && (1\pm \varepsilon)^2 \cdot f'(\alpha,\theta). & (Lemma (ref)) \end{align*} By replacing $\varepsilon$ with $O(\varepsilon)$, we prove the correctness. For the computation time, the computation in Line 2 costs $O(d\sum_{i\in [N]}T_i)$ time since each $b_i$ and $a_i$ can be computed in $O(d T_i)$ time, and $A$ can be computed in $O(N)$ time. In Line 3, it costs $O(Ndk\ln N \ln k)$ time to solve the $k$-means clustering problem by $k$-means++. Line 4 costs $O(Ndk)$ time since each $p(i)$ costs $O(dk)$ time to compute. Lines 5-6 cost $O(Nd)$ time for computing sensitivity function $s$. Lines 7-8 cost $O(N)$ time for constructing $I_S$. Overall, it costs $O(d\sum_{i\in [N]}T_i + Ndk\ln N \ln k)$ at the first stage. Line 10 costs at most $O(d T_i)$ time to compute $\ensuremath{\mathsf{OPT}}_i^{(O)}$. Lines 11-12 cost $O(T_i)$ time to compute $s_i$. Lines 13-14 cost $O(T_i)$ time to construct $J_{S,i}$. Since $|I_S|\leq N$, we have that it costs at most $O(d\sum_{i\in [N]}T_i)$ time at the second stage. We complete the proof.

Limitations, conclusion, and future work

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 celis2018fair,huang2019coresets; future research needs to extend such techniques for time series coresets.

Acknowledgments

This research was supported in part by an NSF CCF-1908347 grant.