Clustering piecewise stationary processes

TL;DR

This paper introduces a novel approach to clustering time-series data generated by piecewise stationary processes, relaxing the stationarity assumption and providing consistent, efficient algorithms for the first time.

Contribution

It formulates the clustering problem for piecewise stationary processes, introduces a notion of consistency, and proposes computationally efficient algorithms with proven consistency.

Findings

Algorithms are consistent without extra assumptions.

First application of clustering to piecewise stationary processes.

Provides a natural formulation for the problem.

Abstract

The problem of time-series clustering is considered in the case where each data-point is a sample generated by a piecewise stationary ergodic process. Stationary processes are perhaps the most general class of processes considered in non-parametric statistics and allow for arbitrary long-range dependence between variables. Piecewise stationary processes studied here for the first time in the context of clustering, relax the last remaining assumption in this model: stationarity. A natural formulation is proposed for this problem and a notion of consistency is introduced which requires the samples to be placed in the same cluster if and only if the piecewise stationary distributions that generate them have the same set of stationary distributions. Simple, computationally efficient algorithms are proposed and are shown to be consistent without any additional assumptions beyond piecewise stationarity.

Equations

$${\mu }_n(\mathbf{x},B)=\{\begin{array}{ll}\frac{1}{n-u+1}{\sum }_{i=1}^{n-u+1}\mathbb{I}\{X_{i..i+u}\in B\}& n\ge u\\ 0& \text{otherwise}\end{array}$$

$$\mu (X_{\mathrm{1..}u}\in B)=\mu (X_{1+j..u+j}\in B)$$

$$\underset{n\to \infty }{\lim }{\mu }_n(\mathbf{x},B)=\mu (B),\mu -\text{ a.s.}.$$

$${\mathcal{B}}^{\otimes m}:=\sigma (\{B_1\times \cdots \times B_m:B_i\in \mathcal{B},i\in \mathrm{1..}m\}).$$

$$\mathbf{X}:=\left(\begin{array}{cccc}{X}_{1,1}& X_{1,2}& X_{1,3}& \dots \\ ⋮& \dots & \ddots & \ddots \\ X_{m,1}& X_{m,2}& \dots & \dots \end{array}\right)$$

$$\mathcal{M}(\rho ):=\{{\mu }_i:i\in \mathrm{1..}m\}$$

$$d(\mu ,{\mu }^{\prime }):=\sum _{u,v\in \mathbb{N}}{w}_u{w}_v\sum _{B\in {\Delta }_{u,v}}|\mu (B)-{\mu }^{\prime }(B)|$$

$${\hat{d}}_n(\mathbf{x},{\mathbf{x}}^{\prime }):=\sum _{u,v\in \mathbb{N}}{w}_u{w}_v\sum _{B\in {\Delta }_{u,v}}|{\mu }_n(\mathbf{x},B)-{\mu }_n({\mathbf{x}}^{\prime },B)|,$$

$${\hat{d}}_n(\mathbf{x},\mu ):=\sum _{u,v\in \mathbb{N}}{w}_u{w}_v\sum _{B\in {\Delta }_{u,v}}|{\mu }_n(\mathbf{x},B)-\mu (B)|$$

$$\underset{n\to \infty }{\lim }{\hat{d}}_n({\mathbf{x}}_s,\mu )=d({\mu }_s,\mu )\rho -\text{ a.s.},$$

$$\underset{n\to \infty }{\lim }{\hat{d}}_n({\mathbf{x}}_s,{\mathbf{x}}_t)=d({\mu }_s,{\mu }_t)\rho -\text{ a.s.}$$

$$Y_1,\dots ,Y_{{\tau }_1},Y_{{\tau }_1+1},\dots ,Y_{{\tau }_2},\dots ,\dots ,Y_n,$$

$$c(j)\mapsto (t^{\ast }(j)+1,j-{\tau }_{t^{\ast }(j)})$$

$$\mathrm{Pr}(Y_{t_1}\in B_1,\dots ,Y_{t_{\ell }}\in B_{\ell })=\rho (X_{c(t_1)}\in B_1,\dots ,X_{c(t_{\ell })}\in B_{\ell }).$$

$$(\rho ,\mathbf{\tau })\sim ({\rho }^{\prime },{\mathbf{\tau }}^{\prime })\iff \mathcal{M}(\rho )=\mathcal{M}({\rho }^{\prime }).$$

$$[(\rho,\boldsymbol{\tau})]:=\big{\{}(\rho^{\prime},\boldsymbol{\tau}^{\prime})\in\mathcal{P}:(\rho^{\prime},\boldsymbol{\tau}^{\prime})\sim(\rho,\boldsymbol{\tau})\big{\}}$$

$$\mathcal{C}:=\{[(\rho ,\mathbf{\tau })]:(\rho ,\mathbf{\tau })\in \mathcal{P}\}$$

$$\mathcal{S}(n):=\{{\mathbf{y}}_1,\dots ,{\mathbf{y}}_N\}$$

$$\alpha :=\underset{i\in \mathrm{1..}N}{\min }\underset{k\in \mathrm{1..}{\kappa }_i+1}{\min }{\theta }_j^{(i)}-{\theta }_{j-1}^{(i)}.$$

$$\underset{n\to \infty }{\lim }f(\mathcal{S}(n),m)=\mathcal{G}.$$

$$\delta (C,C^{\prime })=\underset{\mu \in \mathcal{M}(\rho )}{\max }\underset{{\mu }^{\prime }\in \mathcal{M}({\rho }^{\prime })}{\min }d(\mu ,{\mu }^{\prime })+\underset{{\mu }^{\prime }\in \mathcal{M}({\rho }^{\prime })}{\max }\underset{\mu \in \mathcal{M}(\rho )}{\min }d({\mu }^{\prime },\mu )$$

$$\underset{n\to \infty }{\lim }\underset{k\in \mathrm{1..}\kappa }{\max }\frac{1}{n}|{\psi }_{i_k^{\ast }}(n)-{\tau }_k|=0.$$

$$\hat{\mathbf{\tau }}\leftarrow \mathcal{L}(\mathbf{y},\lambda )\text{and}$$

$$\mathcal{U}$$

$${\mathcal{U}}^{\prime }$$

$$n$$

$$\delta (\mathbf{y},{\mathbf{y}}^{\prime },\lambda )$$

$$\underset{n\to \infty }{\lim }{\hat{\delta }}_n(\mathbf{y},{\mathbf{y}}^{\prime })=\delta (C,C^{\prime }),P-\text{ a.s.}$$

$$\underset{n\to \infty }{\lim }f(\mathcal{S}(n),m)=\mathcal{G},P-\text{ a.s.}$$

$$i^{\ast }:=\min \underset{k\in 1..|\mathbf{\tau }|}{\mathrm{argmax}}\left|\{j\in u..v:(\underset{i\in 0..|\mathbf{\tau }|}{\max }{\tau }_i\le j)=k\}\right|$$

$$\underset{n\to \infty }{\lim }\underset{j=\mathrm{1..}\kappa }{\max }\frac{1}{n}|{\psi }_{[j]}(n)-{\tau }_j|=0.$$

$$\underset{n\to \infty }{\lim }\underset{\overline{\mathbf{y}}\in \mathcal{U}}{\max }{\hat{d}}_n(\overline{\mathbf{y}},{\mu }^{\ast }(\overline{\mathbf{y}}))=0.$$

Full text

1 Introduction

Clustering, in an informal sense, involves breaking a dataset into possibly disjoint subsets called clusters where the elements within the same cluster are somehow more similar to each other than to those in other clusters. This task, which is often a first step in data-analysis, is meant to help with the initial steps to making sense of the data that typically have complex structures and represent some unknown underlying phenomena to be inferred. Given the nature of the problem, it is desirable to make as little assumptions as possible about the underlying mechanisms generating the data. Moreover, the minimal assumptions made should ideally be qualitative and easily verifiable from an application’s perspective.

In this paper we consider a subclass of the clustering problem where each data-point is a time series. Indeed, such sequential data are ubiquitous in modern applications involving, for example, user-behaviour, social networks, as well as financial or biological data, where the observations are sequential by nature, and/or are collected over time. The common features in these real-world datasets are the absence of precise models as well as an abundance of data. From a mathematical perspective, a learning problem involving sequential data can be formulated as follows. Given sequences of the form $Y_{1},...,Y_{n}$ the aim is to make inference about the stochastic mechanisms that generate them. This task is typically done under the assumption that $Y_{i}$ are independently and identically distributed (i.i.d), or that their distribution belongs to a specific model class. However, such assumptions clearly undermine the possibly complex nature of data which may possess long-range dependencies.

To address this gap, one approach is to assume that the process distributions are stationary without requiring any conditions to hold on their memory. That is, $Y_{i}$ need not be independent nor do they need to possess any finite-memory or mixing properties. This allows for arbitrary long-range dependencies between the observations. Moreover, thanks to Birkhoff’s ergodic theorem, under this assumption alone, the frequency of each event almost surely converges to its underlying probability, even though there is no guarantee on the speed of this convergence. Intuitively, non-zero probability events that happen once, happen infinitely often, and their frequencies are meaningful. This characteristic is already enough to make some statistical inference. In particular, in [9] it was suggested to cluster time-series samples based on the distribution that generates them, putting together those and only those samples whose distribution is the same. Making use of the fact that in this setting the target clustering has the so-called strict separation property (e.g., [1]), it was shown that asymptotically consistent clustering is achievable under the assumption of stationarity alone; see also [8] and [11].

While a rather weak assumption, stationarity often breaks in applications: events can happen from which there is no coming back. Simple real-world examples concerning user-behaviour may include such events as changing a job, changing a mobile phone, or having a child. Events that were typical in the past stop happening and events of a new kind start occurring. In these cases, the time index itself bares information. It may still make sense to measure frequencies of events in-between the changes, even though we may not be able to tell from the data alone, when the change happened. Moreover, the changes themselves may be described as abrupt, or, more generally, such that the transition stages are negligible at least asymptotically. Thus, the model that we suggest for the data here is based on piecewise stationary processes. Each time series that is given can be broken into segments, such that within each segment the distribution of the data is stationary. The segments’ boundaries are arbitrary unknown, as are their corresponding distributions. Moreover, no assumptions beyond stationarity are made on the distributions of the stationary segments. In particular, the data within each segment are not assumed to be independent or possess any finite-memory or mixing properties.

We propose to identify a piecewise stationary distribution with a “bag of distributions” corresponding to the stationary segments. Thus, the locations where the change of distribution occurs, as well as the order of the stationary distributions, are disregarded. We consider a pair of piecewise stationary distributions equivalent if the set of stationary distributions of their segments coincide. Thus, the clustering objective is to put into the same cluster those and only those time-series samples whose distributions are equivalent in this sense. An algorithm to achieve this objective is said to be consistent. Our main result is an algorithm that, as we show, is asymptotically consistent under the only assumption that each time series is stationary ergodic. The algorithm needs to know the correct number of clusters, which, as is already the case for stationary time series, is a necessary assumption in this setting. In addition, the length of each segment is required to be linear in the sample size, and a lower-bound on the length of the segments is assumed to be provided. These conditions formalize the intuition that transition periods are negligible with respect to the segments, and each segment is sufficiently long to allow for making inference. The consistency result is obtained using a novel distance between equivalence classes of piecewise stationary distributions. This distance is based on minimax distances between the distributions that generate the segments. The main technical argument establishes that this distance can be estimated consistently based on samples. The proposed algorithm extends those on clustering stationary time series [9, 8]. It uses as sub-routines the algorithms for change-point analysis developed in [5, 6]. Note, however, that the definition of the clustering objective is more permissive here: the distributions of the time series within each cluster may be different, only the set of distributions of the segments has to be the same.

The results provided in this paper are theoretical, and their main appeal is in their generality. While the proposed clustering algorithm is shown to be computationally feasible, its detailed experimental investigation is left as future work. It should be noted that the literature on the related topic of change-point analysis is vast, but is mostly concerned with independent or mixing data, and also restricts the nature of the changes to single-dimensional marginals. The exceptions are the works cited above and references therein. To the best of our knowledge, there are no prior attempts to consider piecewise stationary distributions in the context of clustering. Outside of change-point analysis, related problems have been considered in the context of prediction, albeit with much more restrictive assumptions on the distributions between the changes. Thus, [13] considers the prediction problem with arbitrary change points i.i.d. segments. See also [3] and references therein.

2 Preliminaries

We sometimes use the abbreviation $u..v$ for $\{u,\dots,v\},~{}u\leq v\in\mathbb{N}$, and let $|\cdot|$ denote the cardinality of a set or the number of elements in a sequence; distinction will be apparent from the context. Let $(\mathcal{X},\mathcal{B}_{\mathcal{X}})$ be a measurable space; in this work we consider $\mathcal{X}=\mathbb{R}$, leaving extensions to more general spaces as future work. Denote by $\Delta_{u,v}:=\{[\frac{i_{1}}{2^{v}},\frac{i_{1}+1}{2^{v}})\times\dots\times[\frac{i_{u}}{2^{v}},\frac{i_{u}+1}{2^{v}}):i_{j}\in 0..2^{v}-1,~{}j\in 1..u\}$ the set of dyadic cubes in $\mathcal{X}^{u},~{}u\in\mathbb{N}$ of side-length $2^{-v}$, and let $\mathcal{B}^{(u)}:=\sigma(\{\Delta_{u,v},~{}v\in\mathbb{N}\})$ be the Borel subsets of $\mathcal{X}^{u},~{}u\in\mathbb{N}$. Let $\mathcal{X}^{\mathbb{N}}$ be the set of all $\mathcal{X}$-valued infinite sequences equipped with the Borel $\sigma$-algebra $\mathcal{B}:=\sigma(\{B\times\mathcal{X}^{\mathbb{N}}:B\in\Delta_{u,v},~{}u,v\in\mathbb{N}\})$. By a stochastic process we mean a probability measure on $(\mathcal{X}^{\mathbb{N}},\mathcal{B})$. Take a sequence of random variables $\mathbf{x}:=\langle X_{t}\rangle_{t\in\mathbb{N}}$ with joint distribution $\mu$ where for every $t\in\mathbb{N}$, $X_{t}:\mathcal{X}^{\mathbb{N}}\to\mathcal{X}$ is the coordinate projection of $\mathbf{a}:=\langle a_{t}\rangle_{t\in\mathbb{N}}\in\mathcal{X}^{\mathbb{N}}$ onto its $t^{\text{th}}$ element, i.e. $X_{t}(\mathbf{a})=a_{t}$. For each $n\in\mathbb{N}$ and $B\in\mathcal{B}^{(u)},~{}u\in\mathbb{N}$ we let $\mu_{n}(\mathbf{x},B):\mathcal{X}^{\mathbb{N}}\rightarrow[0,1],~{}n\in\mathbb{N}$ denote the empirical measure of $B$ where,

[TABLE]

and $\mathbb{I}$ is the indicator function. We often call $X_{1..n}$ a sample of length $n\in\mathbb{N}$ generated by a stochastic process $\mu$ with corresponding sequence of random variables $\langle X_{t}\rangle_{t\in\mathbb{N}}$.

Definition 1** (Stationary Ergodic Processes).**

A process $\mu$ is stationary if

[TABLE]

for every $B\in\mathcal{B}^{(u)},~{}u\in\mathbb{N}$ and every $j\in\mathbb{N}$. A stationary process $\mu$ with corresponding sequence of random variables $\mathbf{x}=\langle X_{t}\rangle_{t\in\mathbb{N}}$ is (stationary) ergodic if with for every $u\in\mathbb{N}$ and $B\in\mathcal{B}^{(u)}$ it holds that

[TABLE]

Remark 1**.**

The above definition can be shown to be equivalent to the standard definition involving triviality of shift-invariant measurable sets; see, e.g., [2].

Joint process distributions.

We may require to simultaneously consider multiple samples $X_{1..n},~{}n\in\mathbb{N}$ generated by different, possibly dependent stationary ergodic processes. To allow for this, we first define a distribution over a matrix of random variables, each row of which shall correspond to one of the samples. Next, we obtain each process as the marginal distribution of the corresponding row of the matrix. More specifically, we have the following formulation. For a fixed $m\in\mathbb{N}$, let $\rho$ be a measure on the space $(\mathcal{X}^{m\times\mathbb{N}},\mathcal{B}^{\otimes m})$ where,

[TABLE]

Define the matrix of $\mathcal{X}$-valued random variables

[TABLE]

where $X_{i,j}:\mathcal{X}^{\mathbb{N}\times\mathbb{N}}\rightarrow\mathcal{X},~{}i,j\in\mathbb{N}$ are jointly distributed according to ${\rho}$, so that for $B\in\mathcal{B}^{\otimes m}$ we have $\Pr({\bf X}\in B)=\rho(B)$. For each $i\in 1..m$, let $\mathbf{x}_{i}:=\langle X_{i,j}\rangle_{j\in\mathbb{N}}$ and define the projection map $\pi_{i}\mapsto\mathbf{x}_{i}$. The marginal distribution $\mu_{i}$ of $\mathbf{x}_{i}$ is then defined as the distribution induced by $\rho$ over the $i^{\text{th}}$ row, i.e. $\mu_{i}:=\rho\circ\pi_{i}^{-1}$. We denote by

[TABLE]

the set of marginal process distributions of $\rho$.

Definition 2** (Distributional Distance).**

A distributional distance between a pair of processes $\mu,\mu^{\prime}$ is defined as

[TABLE]

where $w_{j}=1/j(j+1),~{}j\in\mathbb{N}$ or any summable sequence of positive weights.

Definition 3** (Empirical Estimates of Distributional Distance).**

Consider a pair of process marginals $\mu,\mu^{\prime}\in\mathcal{M}(\rho)$ given by (4) with corresponding sequence of random variables $\mathbf{x}$ and $\mathbf{x}^{\prime}$ respectively, where $\mu:=\mu_{i},~{}\mu^{\prime}:=\mu_{j}$, and $\mathbf{x}:=\langle X_{i,t}\rangle_{t\in\mathbb{N}},~{}\mathbf{x}^{\prime}:=\langle X_{j,t}\rangle_{t\in\mathbb{N}}$ correspond to rows $i,j\in 1..m$ of $\bf{X}$ given by (3). Empirical estimates of $d(\mu,\mu^{\prime})$ are given by

[TABLE]

where $\mu_{n}(\cdot,\cdot)$ is given by (1) and $w_{j},~{}j\in\mathbb{N}$ is as in Definition 2.

Remark 2**.**

Note that (5) can be efficiently calculated with computational complexity $\mathcal{O}(n\operatorname{polylog}n)$ for $u_{n}:=\log n,~{}v_{n}:=-\log(s_{\min})$, where $s_{\min}$ is the minimal non-zero difference between the union of all the elements of the two sequences $\mathbf{x},\mathbf{x}^{\prime}$, see [8].

Proposition 1** ([8]).**

If the marginals in $\mathcal{M}(\rho)$ given by (4) are stationary ergodic, then for any $\mu\in\mathcal{M}(\rho)$ and any $s,~{}t\in 1..m$ it holds that,

[TABLE]

where $\mathbf{x}_{j}:=\langle X_{j,t}\rangle_{t\in\mathbb{N}}$ correspond to $j^{\text{th}}$ row of $\bf{X}$ given by (3) and $\mu_{j}\in\mathcal{M}(\rho)$, $j=s,~{}t$.

3 Problem formulation

Piecewise stationary ergodic processes.

In this paper, we shall be dealing with multiple samples of the form

[TABLE]

where the (stationary) segments $Y_{\tau_{i}},\dots,Y_{\tau_{i+1}}$ are generated by different, possibly dependent, stationary ergodic processes. To specify the distribution of the sample $\mathbf{y}$ we define a distribution over a matrix of random variables (3), each row of which shall correspond to a stationary segment of the sample. The stationary-segment distributions of (7) are then obtained by projecting the stationary segments onto the corresponding rows of the matrix. More formally, we specify a Piecewise Stationary Ergodic Process as follows. Consider a measure $\rho$ on $(\mathcal{X}^{\kappa\times\mathbb{N}},\mathcal{B}^{\otimes{\kappa+1}})$ for some fixed $\kappa\in\mathbb{N}$ with set of marginals $\mathcal{M}(\rho)=\{\mu_{i},~{}i\in 1..\kappa+1\}$ as specified by (4), where $\mu_{i}\neq\mu_{i+1},~{}i\in 1..\kappa$ are stationary ergodic. Fix some $n\in\mathbb{N}$ and a sequence $\boldsymbol{\tau}:=\langle\tau_{i}\rangle_{i\in 1..\kappa}$ with $\tau_{1}<\tau_{2}<\dots<\tau_{\kappa}\in 1..n$. Define the mapping $c:\mathbb{N}\rightarrow\mathbb{N}\times\mathbb{N}$ as

[TABLE]

where $t^{*}(j):=\max_{i\in 0..\kappa+1}\tau_{i}\leq j$ picks out the closest change point $\tau_{i}$ to $j\in\mathbb{N}$, with the convention that $\tau_{0}:=0$ and $\tau_{\kappa+1}:=n$. A Piecewise Stationary Ergodic Sample of the form (7) generated by $(\rho,\boldsymbol{\tau})$ can be specified as a sequence of coordinate projections $Y_{t}:\mathcal{X}^{n}\rightarrow\mathcal{X},~{}t\in 1..n$ such that for any $\ell\in 1..n,~{}t_{1},\dots,t_{\ell}\in 1..n$ and $B_{i}\in\mathcal{B}_{\mathcal{X}},~{}i\in 1..\ell$ it holds that

[TABLE]

Observe that by (9) the distribution of each segment $Y_{\tau_{i}+1..\tau_{i+1}}$ is given by a stationary ergodic process $\mu_{i},~{}i\in 1..\kappa+1$. Since it is assumed that $\mu_{i}\neq\mu_{i+1},~{}i\in 1..\kappa$, the indices $\tau_{i},~{}i\in 1..\kappa$ are referred to as change points. The pair $(\rho,\boldsymbol{\tau})$ composed of the measure $\rho$ and its corresponding sequence of change points $\boldsymbol{\tau}$ defines a piecewise stationary ergodic process.

Definition 4** (Equivalence of piecewise stationary ergodic processes).**

We consider a pair of piecewise stationary ergodic processes $(\rho,\boldsymbol{\tau})$ and $(\rho^{\prime},\boldsymbol{\tau}^{\prime})$ equivalent, if and only if they agree on their set of stationary-segment distributions, i.e.,

[TABLE]

Let $\mathcal{P}$ denote the set of all piecewise stationary ergodic processes. The equivalence relation defined above induces a partitioning of $\mathcal{P}$ into distinct classes $[(\rho,\boldsymbol{\tau})],~{}(\rho,\boldsymbol{\tau})\in\mathcal{P}$ where

[TABLE]

so that two piecewise stationary ergodic processes belong to the same class, if and only if they are equivalent in the sense of (10). We let

[TABLE]

denote the set of all such equivalence classes.

Clustering Problem.

Fix some $N\in\mathbb{N}$, which is the number of samples, and (unknown) sequence $\kappa_{i}\in\mathbb{N},~{}i\in 1..N$ corresponding to the number of change points in each sample. Moreover, define (unknown) sequences $\boldsymbol{\theta}_{i}:=\langle\theta^{(i)}_{j}\rangle_{j\in 1..\kappa_{i}+1},~{}i\in 1..N$ with $\theta^{(i)}_{j}\in(0,1)$. For any $n\in\mathbb{N}$, define the change points $\boldsymbol{\tau}_{i}(n):=\langle\tau^{i}_{j}(n)\rangle_{j\in 1..\kappa_{i}+1},~{}i\in\mathbb{N}$ where $\tau^{i}_{j}(n):=\lfloor n\theta_{j}^{(i)}\rfloor-\lfloor n\theta_{j-1}^{(i)}\rfloor$, with the convention that $\theta_{0}^{(i)}=0$ for all $i\in 1..N$. The problem is formulated as follows. For a fixed $n\in\mathbb{N}$, we are given a set

[TABLE]

of $N$ piecewise stationary ergodic samples of the form (7), each generated by an unknown piecewise stationary ergodic process $(\rho_{i},\boldsymbol{\tau}_{i}(n)),~{}i\in 1..N$. Thus, each sample is of length $n_{i}:=\lfloor n\theta_{\kappa_{i}+1}\rfloor$ and has $\kappa_{i}$ change points. It is assumed that each of $N$ piecewise stationary ergodic processes that generate the samples belongs to one of $m$ distinct classes $C_{1},\dots,C_{m}\in\mathcal{C}$, which are unknown. Define the normalized minimum separation between the change points as

[TABLE]

We assume that $\alpha>0$ so that each stationary segment in $\mathbf{y}_{i},~{}i\in 1..N$ is of length at least $n\alpha$.

Definition 5** (Ground Truth Clustering).**

Let $\mathcal{G}:=\{\mathcal{G}_{1},\dots,\mathcal{G}_{m}\}$ be a partitioning of $1..N$ where for any $i\in 1..N$, it holds that $i\in\mathcal{G}_{\ell}$ for some $\ell\in 1..m$ if and only if $(\rho_{i},\boldsymbol{\tau}_{i}(n))\in C_{\ell}$. We call $\mathcal{G}$ the ground-truth clustering.

Thus, in the ground-truth clustering two samples fall into the same cluster if and only if the piecewise-stationary distributions that generate them are equivalent, in the sense that they have the same set of stationary distributions of the segments.

A clustering function $f$ takes a set $\mathcal{S}$ of samples along with the number $m$ of target clusters to produce a partition $f(\mathcal{S},m)\mapsto\{J_{1},\dots,J_{m}\}$ of $1..N$, aiming to recover the ground-truth $\mathcal{G}$.

Definition 6** (Consistency).**

A clustering function $f$ is consistent for a set of samples $\mathcal{S}=\mathcal{S}(n),~{}n\in\mathbb{N}$ if $f(\mathcal{S},m)=\mathcal{G}$. Moreover, $f$ is called asymptotically consistent if with probability $1$ it holds that

[TABLE]

Joint distribution of piecewise stationary samples.

Observe that the problem requires us to simultaneously consider multiple samples, each generated by a piecewise stationary ergodic process. These samples can themselves be dependent. Formally, this is defined through the following construction. Consider the space $\mathcal{Y}:=\mathcal{X}^{\kappa_{1}\times\mathbb{N}}\times\dots\times\mathcal{X}^{\kappa_{N}\times\mathbb{N}}$. Denote by $\mathcal{F}:=\mathcal{B}_{1}\otimes\dots\otimes\mathcal{B}_{N}$ the product $\sigma$-algebra on $\mathcal{Y}$ where $\mathcal{B}_{i}:=\sigma(\{B_{1}\times\dots\times B_{\kappa_{i}}:B_{j}\in\mathcal{B}\}),~{}i\in 1..N$ is in turn the product $\sigma$-algebra on $\mathcal{X}^{\kappa_{i}\times\mathbb{N}}$. Let $P$ be a probability measure on $(\mathcal{Y},\mathcal{F})$. Consider a sequence ${\bf Z}:=\langle{\bf Y}_{i}\rangle_{i\in 1..N}$ of infinite matrices of $\mathcal{X}$-valued random variables $Y_{s,t}^{(i)}:\mathcal{X}^{\kappa_{i}\times\mathbb{N}}\rightarrow\mathcal{X},~{}s\in 1..\kappa_{i},~{}t\in\mathbb{N},i\in 1..N$, which can be easily shown to be $\mathcal{F}$-measurable. Suppose that $P$ is the distribution of ${\bf Z}$ so that $\Pr({\bf Z}\in F)=P(F)$ for all $F\in\mathcal{F}$. For each $i\in 1..N$, define the projection map $\widetilde{\pi}_{i}\mapsto\langle Y_{s,t}^{(i)}\rangle_{s\in 1..k_{i},~{}t\in\mathbb{N}}$. Then $\rho_{i}:=P\circ\widetilde{\pi}_{i}^{-1},~{}i\in 1..N$ is the measure of ${\bf Y}_{i}$. Our main probabilistic statements will be stated in terms of $P$.

The role of $n$.

The asymptotic results in this paper are all with respect to $n$ approaching infinity. While the formulation allows for all of the samples, and, of course, all of the stationary segments, to be of different lengths, $n$ parametrizes (via the unknown sequences ${\bf\theta}_{i}$) all these lengths. Thus, when $n$ goes to infinity, the lengths of all the samples and all the individual segments within the samples also approach infinity. Note, however, that the clustering protocol is not “online”, and thus the algorithm proposed below does not deal with sequences that grow over time, but only considers a fixed set of sequences of fixed lengths. The asymptotic results are thus to be interpreted as stating that if all the sequences are long enough then the algorithm is correct.

4 Main Results

In this section we outline the main results. We start by introducing a distance between equivalence classes of piecewise stationary distributions. Next, we show that this distance can be consistently estimated based on piecewise stationary samples. The distance estimates will then be used to construct an asymptotically consistent clustering algorithm. The proofs as well as auxiliary technical results are provided in Section 5.

Definition 7** (Distance between piecewise stationary ergodic classes).**

Let $C,C^{\prime}\in\mathcal{C}$ be two classes of piecewise stationary ergodic processes where $C=[(\rho,\boldsymbol{\tau})]$ and $C^{\prime}=[(\rho^{\prime},\boldsymbol{\tau}^{\prime})]$. We define a distance between $C,~{}C^{\prime}$ as follows.

[TABLE]

Proposition 2**.**

The distance $\delta$ induces a metric on the set of piecewise stationary ergodic equivalence classes.

We show that the distance $\delta$ can be estimated consistently. The proposed approach to estimate $\delta$ is outlined in Algorithm 1. This algorithm in turn relies on a procedure, introduced in [5] that, given a sample $\mathbf{y}$ of the form (7) with $\kappa$ change points generated by a piecewise stationary process outputs an exhaustive list of $K\geq\kappa$ candidate change-point estimates. While the algorithm from [5] does not attempt to estimate the number of change points $\kappa$, among the $K$ candidate change points that it outputs there are $\kappa$ consistent estimates of the unknown change points. It is worth noting that [5] establishes a much stronger property, namely, it sorts the list in such a way that its first $\kappa$ elements estimate the true change points of $\mathbf{y}$. However, we shall not require this feature here: it is enough to have a list of arbitrary order that includes a correct estimate for each change point; for simplicity, we assume that the $K$ candidate change-points are sorted left-to-right. More precisely we have the following definition.

Definition 8** (List-estimator [5]).**

A list-estimator is a function $\mathcal{L}$ that, given a sample of length $n\in\mathbb{N}$ with $\kappa\in\mathbb{N}$ change points generated by a piecewise stationary process $(\rho,\boldsymbol{\tau})$ produces a list $\psi_{1}(n)<\dots<\psi_{K}(n)\in\{1,\dots,n\}^{K}$ of some $K\geq\kappa$ candidate estimates. A list-estimator is said to be consistent if with probability $1$ there exists some $\mathcal{I}:=\{i^{*}_{1},\dots,i^{*}_{\kappa}\}\subseteq 1..K$ such that

[TABLE]

Proposition 3** ($\delta$ can be estimated consistently).**

Consider the samples $\mathbf{y},\mathbf{y}^{\prime}$ generated by a distribution $P$ with piecewise stationary ergodic marginals $(\rho,\boldsymbol{\tau})$ and $(\rho^{\prime},\boldsymbol{\tau^{\prime}})$; the lengths of the samples are parametrized by $n$. Let the estimate $\widehat{\delta}_{n}(\mathbf{y},\mathbf{y}^{\prime}):=\delta(\mathbf{y},\mathbf{y}^{\prime},\lambda)$ be obtained as the output of Algorithm 1 with $\mathbf{y},~{}\mathbf{y}^{\prime}$ and any $\lambda\in(0,\alpha]$ as input, where $\alpha$ is given by (13). Then

[TABLE]

where, $C:=[(\rho,\boldsymbol{\tau})]$ and $C^{\prime}:=[(\rho,\boldsymbol{\tau}^{\prime})]$ are the equivalence classes containing the piecewise stationary ergodic processes that generate $\mathbf{y}$ and $\mathbf{y}^{\prime}$ respectively.

Finally, the distance estimates given by Algorithm 1 can be used to construct a clustering algorithm. The clustering algorithm, i.e. Algorithm 2, starts by initializing the clusters using farthest-point initialisation, and then assigns the remaining samples to the nearest cluster. This is the same procedure as used in [9, 8] to cluster stationary samples, but using a different distance, $\widehat{\delta}_{n}$. The initialisation procedure is due to [4].

Theorem 1** (Algorithm 2 is asymptotically consistent).**

Let $f(\mathcal{S}(n),m)=\operatorname{Alg}\ref{alg:clust}(\mathcal{S}(n),m,\lambda)$ be the output of Algorithm 2 when provided with the set $\mathcal{S}(n)$ of piecewise stationary ergodic samples (12), along with the correct number $m$ of target clusters and some $\lambda\in(0,\alpha]$, where $\alpha$ is given by (13). It holds that,

[TABLE]

Moreover, the computational complexity of the algorithm is $\mathcal{O}(mN(n^{2}\operatorname{polylog}n+\lambda^{-2}n\operatorname{polylog}n))$.

5 Proofs

In this section we prove our main results. We start by providing a simple argument showing that $\delta$ is indeed a metric on $\mathcal{C}$. Next, we prove Proposition 3. To this end, we first introduce a notation via Definition 9 to identify the stationary-segment distribution which generates the largest connected portion of a sub-sequence of a piecewise stationary ergodic sample. The proof of the proposition relies on some technical results, namely, on [5, Theorem 1] and [6, Lemma 2] which are stated for the sake completeness, as well as on Lemma 1, stated and shown below.

Proof of Proposition 2.

It is easy to verify that, since $d(\cdot,\cdot)$ is a metric, $\delta$ is positive, symmetric and satisfies the triangle inequality. It remains to show that $\delta(C,C^{\prime})=0$ if and only if $C=C^{\prime}$ for all $C,~{}C^{\prime}\in\mathcal{C}$. Consider two piecewise stationary ergodic processes $(\rho,\boldsymbol{\tau})$ and $(\rho^{\prime},\boldsymbol{\tau}^{\prime})$ in $\mathcal{C}$. Let $\mathcal{M}(\rho):=\{\mu_{i}:i\in 1..k\}$ and $\mathcal{M}(\rho^{\prime}):=\{\mu^{\prime}_{i}:i\in 1..k^{\prime}\}$ denote the sets of distinct marginals corresponding to $\rho$ and $\rho^{\prime}$ respectively. If $\delta(\mathcal{C},\mathcal{C}^{\prime})=0$ it must hold that $\max_{i\in 1..k}\min_{j\in 1..k^{\prime}}d(\mu_{i},\mu^{\prime}_{j})=0$ and $\max_{i=1..k^{\prime}}\min_{j\in 1..k}d(\mu^{\prime}_{i},\mu_{j})=0$, leading to $\mathcal{M}(\rho)\subseteq\mathcal{M}(\rho^{\prime})$ and $\mathcal{M}(\rho^{\prime})\subseteq\mathcal{M}_{\rho}$ respectively, so that $\mathcal{M}(\rho)=\mathcal{M}(\rho^{\prime})$. On the other hand, suppose that $\mathcal{M}(\rho)=\mathcal{M}(\rho^{\prime})$. If $\mu\in\mathcal{M}(\rho)$ then $\mu\in\mathcal{M}(\rho^{\prime})$ and hence $\min_{\mu^{\prime}\in\mathcal{M}(\rho^{\prime})}d(\mu,\mu^{\prime})=0$, leading to $\max_{i\in 1..k}\min_{j\in 1..k^{\prime}}d(\mu_{i},\mu^{\prime}_{j})=0$. Similarly we obtain $\max_{i\in 1..k^{\prime}}\min_{j\in 1..k}d(\mu^{\prime}_{i},\mu_{j})=0$, so that $\delta(\rho,\rho^{\prime})=0$, and the result follows. ∎

Definition 9**.**

Consider a sample $\mathbf{y}:=\langle Y_{t}\rangle_{t\in 1..n},~{}n\in\mathbb{N}$ generated by a piecewise stationary ergodic process $(\rho,\boldsymbol{\tau})$. For any $u<v\in 1..n$ fix a sub-sequence $\overline{\mathbf{y}}:=\langle Y_{t}\rangle_{t\in u..v}$ of $\mathbf{y}$. We define $\mu^{*}(\overline{\mathbf{y}}):=\mu_{i^{*}}$ where

[TABLE]

where the $\min$-operator is used to consistently break ties.

** [5, Theorem 1].**

There exists an asymptotically consistent list-estimator $\mathcal{L}$ that, given a piecewise stationary ergodic sample of length $n\in\mathbb{N}$ with $\kappa\in\mathbb{N}$ change points generated by a piecewise stationary process $(\rho,\boldsymbol{\tau})$ produces a list $\psi_{1}(n),\dots,\psi_{K}(n)\in\{1,\dots,n\}^{K}$ of some $K\geq\kappa$ candidate estimates, that are at least $\lambda n$ apart. Let $\psi_{[1]}(n)\leq\psi_{[2]}(n)\leq\dots\leq\psi_{[\kappa]}(n)$ be the first $\kappa$ elements of $\psi_{j},~{}j\in 1..K$ sorted in increasing order of value. With probability $1$ it holds that

[TABLE]

[6, Lemma 2]****.

Let $\mathbf{y}$ be a piecewise stationary ergodic sample of length $n$ with $\kappa$ change points and at least $\alpha n$ apart for some $\alpha\in(0,1)$. Let $\mathcal{U}$ be the set of segments of the form (15). For all $\lambda\in(0,\alpha]$ with probability $1$ we have

[TABLE]

Lemma 1**.**

Consider the set $\mathcal{S}$ of piecewise stationary samples specified by (12). For any $j\in 1..|\mathcal{S}|$ suppose that $\mathbf{y}_{j}\in\mathcal{S}$ along with some $\lambda\in(0,\alpha]$ is provided as input to Algorithm 1, where $\alpha$ given by (13) denotes the minimum normalized separation between the change points. With probability $1$ it holds that, for every $\epsilon>0$ there exists some $N$ such that for all $n\geq N$ we have

[TABLE]

where $\mathcal{M}_{n}(\rho_{j}):=\{\mu^{*}(\overline{\mathbf{y}}):\overline{\mathbf{y}}\in\mathcal{U}\}$,  $\mathcal{U}$ is specified by (15) in Algorithm 1 and $\mu^{*}(\cdot)$ is given by Definition 9.

Proof.

For simplicity of notation let $\mathbf{y}:=\mathbf{y}_{j}$, $(\rho,\boldsymbol{\tau}(n))=(\rho_{j},\boldsymbol{\tau}_{j}(n)),~{}n\in\mathbb{N}$ and $\kappa:=\kappa_{j}$. First observe that by definition $\mathcal{M}_{n}(\rho)\subseteq\mathcal{M}(\rho)$ for all $n\in\mathbb{N}$. Therefore, to prove the statement, it suffices to show that for large enough $n$ we have $\mathcal{M}(\rho)\subseteq\mathcal{M}_{n}(\rho)$. Fix $\epsilon\in(0,\lambda/2)$. By [5, Theorem 1], there exists some $N$ such that for all $n\geq N$ the list of candidate change-points produced by $\mathcal{L}$ is at least $\kappa$ long and includes a subset of size $\kappa$ composed of consistent change point estimates. That is, there exists some $\mathcal{I}:=\{i^{*}_{1},\dots,i^{*}_{\kappa}\}\subseteq 1..|\boldsymbol{\widehat{\tau}}(n)|$ such that

[TABLE]

where $\psi_{0}:=0$ and $\psi_{|\boldsymbol{\widehat{\tau}}(n)|+1}:=n$. Note that (19) implies that the candidate estimates obtained by $\mathcal{L}$ are at least $n\lambda$ apart. We show that for every $\nu\in\mathcal{M}(\rho)$ there exists some $\overline{\mathbf{y}}\in\mathcal{U}$ such that $\nu=\mu^{*}(\overline{\mathbf{y}})$ so that $\nu\in\mathcal{M}_{n}(\rho)$ for all $n\geq N$. Take some $\nu\in\mathcal{M}(\rho)$. By construction, there exists some $k\in 0..\kappa$ such that $\mathbf{y}_{\tau_{k}+1..\tau_{k+1}}$ is generated by $\nu$. By (18) there exists some $i^{*}_{k}\in\mathcal{I}$ such that $\frac{1}{n}|\psi_{i^{*}_{k}}-\tau_{k}|\leq\epsilon$. We show that $[\psi_{i^{*}_{k}},\psi_{i^{*}_{k}+1}]\subseteq[\tau_{k}-\epsilon n,\tau_{k+1}+\epsilon n]$. We have two cases: either $i^{*}_{k}+1\in\mathcal{I}$, in which case by (18) it holds that $\frac{1}{n}|\psi_{i^{*}_{k}+1}-\tau_{k+1}|\leq\epsilon$ for all $n\geq N$, or $i^{*}_{k}+1\in\mathcal{I}^{c}:=\{1,\dots,|\boldsymbol{\widehat{\tau}}(n)|\}\setminus\mathcal{I}$. We argue that in the latter case $\psi_{i^{*}_{k}+1}<\tau_{k+1}$. To see this, assume by way of contradiction that $\psi_{i^{*}_{k}+1}>\tau_{k+1}$ where $\tau_{k+1}\neq n$; note that the statement trivially holds for $\tau_{k+1}=n$. By the consistency of $\mathcal{L}(\mathbf{y},\lambda)$ there exists some $j>i^{*}_{k}\in\mathcal{I}$ such that

[TABLE]

for all $n\geq N$. Moreover, by (18) and (19) for all $n\geq N$ the candidates indexed by the elements of $\mathcal{I}^{c}$ are linearly separated from the true change points, that is,

[TABLE]

Thus, from (18) and (20) we obtain that $\psi_{i^{*}_{k}}-\psi_{j}\geq\lambda-2\epsilon>0$. Since by construction $\psi_{i},~{}i\in 0..|\boldsymbol{\widehat{\tau}}(n)|+1$ are sorted in increasing order of value, this leads to a contradiction. Therefore, we have $[\psi_{i^{*}_{k}},\psi_{i^{*}_{k}+1}]\subseteq[\tau_{k}-\epsilon n,\tau_{k+1}+\epsilon n]$. It follows that $\mu^{*}(\overline{\mathbf{y}})=\nu$, where $\overline{\mathbf{y}}=\mathbf{y}_{\psi_{i^{*}_{k}}..\psi_{i^{*}_{k}+1}}$. This implies that $\nu\in\mathcal{M}_{n}(\rho)$, and the statement follows. ∎

Proof of Proposition 3.

Let

[TABLE]

denote the minimum non-zero distributional distance between the stationary-segment distributions in $\mathcal{M}(\rho)$ and $\mathcal{M}(\rho^{\prime})$. Fix some $\epsilon\in(0,\delta_{\min}/3)$ . Let $\mathcal{U}$ and $\mathcal{U}^{\prime}$, respectively specified by (15) and (16), be the stationary segments identified by Algorithm 1 in each sample. Note that by construction the maximum number of segments produced by $\mathcal{L}(\cdot,\lambda)$ is $1+\lambda^{-1}$ so that with $\lambda\leq\alpha$ we have

[TABLE]

where $\alpha$ given by (13) is the minimum normalized distance between the change points. Moreover, note that $n_{i}=n\theta^{(i)}_{\kappa_{i+1}}$ and $n_{j}=n\theta^{(j)}_{\kappa_{j+1}}$ where $\theta_{\kappa_{i+1}}^{(i)},~{}\theta^{(j)}_{\kappa_{j+1}}\in(0,1)$, we have $\max\{n_{i},n_{j}\}\leq n$ so that for any $\lambda\in(0,\alpha]$ it holds that

[TABLE]

where $\theta_{\min}:=\min\{\theta_{k}^{(i)},~{}\theta_{k^{\prime}}^{(j)},~{}k\in 1..\kappa_{i+1},~{}k^{\prime}\in 1..\kappa^{\prime}_{j+1}\}$ corresponds to the minimum length of the stationary segments in $\mathbf{y}$ and $\mathbf{y}^{\prime}$. By (23) the conditions for [6, Lemma 2] hold so that by this lemma there exists some $N_{0}\in\mathbb{N}$ such that for all $n\geq N_{0}$ we have

[TABLE]

Moreover, by Lemma 1 there exists some $N_{1}$ such that for all $n\geq N_{1}$ it holds that

[TABLE]

To prove the statement we proceed as follows. Take $n\geq\max\{N_{0},N_{1}\}$. We have,

[TABLE]

where (26) and (27) follow from the application of triangle inequality, (28) follows from (24), and (29) follows from the definition of $\mathcal{M}_{n}$ and (30) follows from(25). On the other hand,

[TABLE]

where (31) and (32) follow from the triangle inequality, (33) follows from (24), (34) follows from the definition of $\mathcal{M}_{n}$ and (35) follows from (25). By (29) and (35) be obtain

[TABLE]

Similarly, we have

[TABLE]

From (36) and (37) it follows that $|\delta(\mathbf{y},\mathbf{y}^{\prime},\lambda)-\delta(C,C^{\prime})|\leq 4\epsilon$. Since the choice of $\epsilon$ is arbitrary, the statement follows. ∎

Proof of Theorem 1.

The proof can be recovered from that of [8, Theorem 1] with the distributional distance $d$ replaced by $\delta$, since the algorithms only differ in this choice of distance. The consistency of the algorithm then follows from the consistency of $\widehat{\delta}_{n}$, namely, Proposition 3. For the sake of completeness, let us recall the main argument from the proof of [8, Theorem 1]. The ground-truth clustering $\mathcal{G}$ has the so-called strict separation property with respect to the distance $\delta$, meaning that any two distributions in the same cluster are closer to each other than to those in any other cluster: the $\delta$ distance is 0 if and only if the distributions are the same and thus in the same cluster. From the consistency of $\widehat{\delta}_{n}$ it follows that the same holds $P$- a.s.from some $n$ on for the $\widehat{\delta}_{n}$ distance. To calculate the computational complexity of the algorithm, first note that by [5, Remark 2] the computational complexity of the list-estimator is $\mathcal{O}(n^{2}\operatorname{polylog}n)$, and that of the distributional distance $\widehat{d}$ is $\mathcal{O}(n\operatorname{polylog}n)$. Moreover, to obtain an estimate of $\widehat{\delta}_{n}$, the algorithm needs two calls to the list estimator as well as an additional $1/\lambda^{2}$ (corresponding the maximal number of estimated segments in each sample) calculations of $\widehat{d}$. Given that the total number of calculations of $\widehat{\delta}_{n}$ of the algorithm is $mN$, this brings the overall computational complexity to $\mathcal{O}(mN(n^{2}\operatorname{polylog}n+\lambda^{-2}n\operatorname{polylog}n))$. ∎

6 Conclusion

In this paper, we have introduced a novel probabilistic framework for clustering time series, which is considerably more general than previously used models. At the same time, as shown here, it allows for provably consistent efficient algorithms. In this section we analyse the conditions of the main theorem, and briefly outline some possible extensions and generalizations of the results obtained.

Necessity of the conditions.

In Theorem 1 it is required that the correct number of clusters $m$, as well as a lower-bound on the minimal distance between change points, be provided. The former requirement is necessary. Indeed, as shown in [10] (see also [11]), for stationary ergodic distributions there is no asymptotically consistent algorithm that, given two samples, would answer whether they were generated by the same process or by different ones. As a corollary, without knowing the number of clusters, it is, in general, impossible to cluster even two samples generated by (single-piece and thus, of course, also piecewise) stationary distributions. As for the latter requirement, the lower-bound $\lambda$ on the minimal distance between change points, comes from the corresponding requirement of the change-point estimation algorithm of [6]. Currently, it is not known whether this requirement is necessary, although we would conjecture that it is. In contrast, if the number of change points is known, then these change points can be estimated consistently without $\lambda\in(0,\alpha]$ given, as is established in [7]. While of lesser practical interest, an analogous consistency result can be established for the clustering problem considered here. To do so, Algorithm 1 would obtain the list of change point estimates in each sample by applying the algorithm of [7], and then proceed without further modifications. Thus, the requirement of a known lower-bound on the minimum distance between change points can be traded for the requirement of a known number of change points in each sample.

Extensions and generalizations.

The time series studied in this paper are assumed to be generated by finitely many piecewise stationary distributions, with changes in distribution being abrupt. One can first ask whether stationary distributions can be generalized; a related question is whether the transitions can be made gradual. Both questions are answered by considering so-called asymptotic mean stationary distributions. Essentially, these are processes such that all frequencies converge almost surely; see [2] for a formal definition and results. Here, it is worth mentioning that, since all that we use is asymptotic convergence of frequencies, all the results can be directly extended to this case; this includes the corresponding results on change-point estimation, see [11]. What this means is that the changes between distributions do not have to be abrupt, as long as they happen over $o(n)$ time-steps. A different generalization could be achieved by considering the online setting of the problem. This means allowing the samples to grow with time, as well as adding new samples potentially at every time step. Such a setting allows one to accommodate a range of new applications that deal with growing bodies of data. The corresponding problem for stationary time series is addressed in [8]. Piecewise stationary samples present new challenges in this respect. Specifically, we would need to allow for an infinite number of change points. More importantly, it would become possible that two samples are generated by equivalent piecewise stationary distributions, but not all of the distributions of the samples have been revealed in each of the segments, and so the distributions look differently. A detailed investigation of what is possible to achieve in this setting is left for future work.

Finite-time guarantees.

In the framework of stationary (ergodic) time series, fundamental results establish the impossibility of obtaining any finite-time guarantees on the error of the resulting algorithms; indeed, already the speed of convergence of frequencies may be arbitrary slow [12]. Therefore, additional assumptions that go beyond stationarity and ergodicity are necessary if one wishes to obtain any finite-time guarantees on the performance. While this falls out of scope of the present paper, it would not be without interest to consider clustering piecewise stationary mixing or even i.i.d. time series: to the best of our knowledge, these questions would still be open. Making such assumptions would also lay grounds for the possibility of constructing consistent algorithms that do not require the knowledge of the correct number of clusters. The stationary mixing case of this problem is briefly considered in [8].