k-MLE, k-Bregman, k-VARs:
Theory, Convergence, Computation

Given independently distributed (i.d.) data vectors (of dimension $d$ ) $x_{n},n=1,\dots,N$, each belonging to one of $K$ clusters, $C_{1},\dots,C_{K}$, we introduce the label matrix or binary cluster membership array $\tau=[\tau_{n,k}]_{N\times K}$

The clustering problem is to estimate $\tau_{N\times K}$ from the $d\times N$ data matrix $X_{1}^{N}=[x_{1},\cdots,x_{N}]$. Note that there are only finitely many $\tau$’s. We denote the $k^{th}$ column of $\tau$ by $\tau_{k}$.

Clustering algorithms are often categorised in various ways: partition methods versus hierarchical methods; hard clustering, where each data vector is assigned to only one cluster, versus soft (a.k.a. fuzzy) clustering, where each data vector may belong to more than one cluster; similarity based versus model based. Here we introduce a new division: between deterministic label clustering (DL-C) and stochastic label clustering (SL-C). DL-C treats the membership variables as deterministic, whereas SL-C treats them as random. We note that SL-C always produces soft clustering whereas DL-C can produce either hard clustering or soft clustering.

The classic DL-C method is k-means[1], where clustering is achieved by minimizing the criterion

w.r.t. $\tau,\mu$. Here $\mu_{k}$ is a cluster centre and $d(x,\mu)$ is a distance or similarity measure, i.e. it has two properties: (i) positivity $d(x,\mu)\geq 0$; (ii) uniqueness $d(x,\mu)=0$ iff $x=\mu$. Sometimes symmetry is useful: (iii) $d(x,\mu)=d(\mu,x)$. k-means uses Euclidean distance but other distance measures are also used, e.g. k-medoids [1].

The classic SL-C method is based on a mixture model [1] in which the $\tau_{nk}$ are independently distributed with $P(\tau_{nk}=1)=\pi_{k}$.

Historically almost all hard clustering methods have been based on similarity measures whereas almost all soft clustering methods are based on models. This seems to have fostered the impression that hard clustering methods could not be model based. But we will see that k-MLE is a DL-C method but is model based. This suggests that a deeper division is between DL and SL rather than between similarity based and model based methods. This paper concentrates on DL methods.

Suppose that $x_{n}$ are i.d. with those in cluster $k$ having density $p(x_{n};\theta_{k})$ where $\theta_{k}\in\Omega$ is $q$-dimensional. Then the joint density or likelihood of the data in cluster $k$ is $\prod_{n=1}^{N}\left[p(x_{n};\theta_{k})\right]^{\tau_{nk}}$ with corresponding log-likelihood

where $\ell(x_{n},\theta_{k})=\mathrm{ln}\,p(x_{n};\theta_{k})$. Denoting $\Theta=[\theta_{1},\dots,\theta_{K}]$, the joint log-likelihood of all the data is then

We call this the deterministic label clustering likelihood (DL-CL). The corresponding clustering likelihood for a SL approach would be called a stochastic label clustering likelihood (SL-CL).

We now note two important properties of the DL-CL which now follow immediately.

• •

  Property P1: $\mathcal{L}(\tau,\Theta)$ is linear in $\tau$.

• •

  Property P2: $\mathcal{L}(\tau,\Theta)$ is separable in $(\tau,\Theta)$.

Separability means that e.g. maximizing over $\Theta$ can be done by separately maximizing each cluster-specific likelihood.

The DL-CL depends on binary and analog (i.e. continuous valued) parameters. This leads to a hybrid maximum likelihood estimation (MLE) problem as follows.

For this optimization to be well defined we need to introduce some assumptions.

Note that (a) and (b) ensure that $\ell(x_{n},\theta)$ is bounded. Since $\mathcal{B}$ is a discrete, bounded set, then under Assumption 1, the k-MLE problem has at least one solution.

The k-MLE algorithm solves the hybrid k-MLE problem by cyclic ascent¹¹1Cyclic optimization is such a simple idea that it is repeatedly reinvented and renamed as if it was new. (a.k.a. coordinate ascent[27]).

Notice there is no requirement that the log-density is a similarity measure. So this provides a considerable generalization of DL algorithms such as k-means.

After the second author had developed this method we discovered it is not new. In [14] the DL-CL likelihood appears under the name ‘classification likelihood’. But it is only presented in the context of a multivariate Gaussian density. When the densities are multivariate spherical Gaussian then k-MLE reduces to k-means. No algorithm details are given in [14].

The DL-CL has also been developed (in a different direction) by [28] under the name ‘model-based k-means’. It is derived heuristically from a SL-CL and an algorithm is given without recognizing that it is cyclic ascent. There is no convergence analysis and no tuning parameter selection. [28] attributes the method to [29],[30]. The first reference only deals with two clusters and has no convergence analysis or tuning parameter selection. The second reference fits a mixture model with an EM algorithm and so is not a DL-C method.

In a seminal piece of work [11] showed ‘there exists a unique Bregman divergence corresponding to every regular exponential family’. Then, using a notion of Bregman information, they developed a hard clustering algorithm based on Bregman divergence. It is not noted in [11] but their algorithm is a coordinate descent algorithm [27]. Further [11] provided a convergence result, which we will see below leaves major questions unanswered.

We now show that the Bregman hard clustering algorithm for exponential families is a special case of k-MLE; we henceforth refer to it as k-Bregman. It is shown in [11] that if $p(x;\theta)$ is a regular exponential family²²2For expontial families in general, the same result is remarked in [31]. it has a unique decomposition

where $d_{\phi}(x;\mu(\theta))\geq 0$ is the Bregman divergence, $\phi$ is the conjugate function of $\psi$, $\mu(\theta)=\nabla\psi(\theta)$ is the expectation parameter corresponding to $\theta$, and the reminder term $b_{\phi}(x)$ does not depend on the parameter $\theta$ and is a normalization factor for the probability density. Thus the log-likelihood for the whole data can be written as

However, by the property of $\tau_{n,k}$, the remainder term $\zeta(X_{1}^{n})$ reduces to $\sum_{n=1}^{N}\ln(b_{\phi}(x_{n}))$ and can be dropped. Therefore for an exponential family k-MLE reduces to the k-Bregman similarity based method.

The possibilities for k-MLE are extensive. Below we develop a new algorithm for clustering vector time series which is not similarity based; we also prove its convergence. We previously developed a scalar version [32] by taking a limit in an AR mixture model. That limiting argument can be extended to the vector case but yields a different algorithm which will be discussed elsewhere.

More generally, k-MLE can handle any data type e.g. multi-modal data, hybrid data with both continuous valued and discrete valued measurements; spatio-temporal data; temporal data sampled on mixed time-scales; multi-subject data and so on.

It turns out to be convenient to consider a ‘relaxed’ or purely analog version of the k-MLE problem in which the membership variables are analog lying in $[0,1]$. This is formalised by introducing an analog constraint set.

Note that $\mathcal{A}$ is bounded and $\mathcal{L}(\tau,\Theta)$ is well defined and bounded for $\tau\in\mathcal{A},\Theta\in\Omega^{\times}$. We can thus introduce the:

We can now introduce the concentrated k-MLE problem:

While Theorem 4 is encouraging, the problem is that, a partial maximum may not be a local maximum; [15] gives an example. We discuss this further in detail in the next section where we develop conditions for local maxima.

The papers [11] and [28] claim to show convergence to a local optimum, but what they show is that the algorithm stops in a finite number of steps. They do not obtain Theorem 4 because they do not have the notion of partial maximum which is a prerequisite for a local maximum.

We introduce the

Since there are finitely many $\tau$’s, there are finitely many $M(\tau)$’s. Next, recalling Property P2, for fixed $\tau$, $\mathcal{L}(\tau,\Theta)$ is separable. Thus

where the operator “$\times$” is the Cartesian product. We say that $M(\tau)$ is a singleton set if it consists of one point. This singleton property turns out to be crucial for convergence.

We need to determine the local maxima of $F(\tau)$ and so need to characterize its gradient. The one-sided directional derivative of $F(\tau)$ in direction $d$ (so that $d$ is a unit vector) is defined by

Then we have the following characterization and a well-known optimality condition.

Now we are ready to present our main result on the local maximum of k-MLE problem.

To establish the singleton property of $M(\tau^{*})$, we need conditions to ensure the uniqueness of MLE. We first introduce the following classic result.

Remark. We state the result for the log-likelihood whereas [34] state it for the likelihood. But the statements are equivalent since the Hessians at a stationary point are equivalent.

Remark. In Proposition 10 rather than applying the conditions of Theorem 9 to establish MLE uniqueness, it is sometimes possible to establish MLE uniqueness directly or by other means. We will do this for k-VARs below.

We can now state the main result.

Remark. As indicated in the previous remark we can alternatively directly show MLE uniqueness when possible.

We could now apply Proposition 10, Theorem 11 to k-Bregman. However it is simpler instead to replace Theorem 11 with yet another result from [34].

So we have to check the Theorem 12 conditions. In view of (7) we need only check Theorem 12 for exponential family densities. The density of an exponential family $\ln p(x;\theta)$ in (7) has the equivalent form [11, Section 4]

where $\psi(\theta)$ is the log partition function

and is convex [11], and $p_{0}(\cdot)$ is an arbitrary function endowing a measure (see [11, Section 4.1]). This means that the log-density is concave and so the Hessian is negative-semi definite at each point of $\Omega$. But we need negative definiteness and to see what that entails we need to explicitly calculate the Hessian.

We find the score function $\frac{\partial\ln p(x;\theta)}{\partial\theta}=x-\frac{\partial\psi}{\partial\theta}$ and

hence $\displaystyle\frac{\partial\ln p}{\partial\theta}=x-\mu(\theta).$ Then we have

in which $\mathbb{E}(X)$ and $\mathrm{Var}(X)$ denote the mean and variance of random variable $X$. Now we can state the Bregman result.

We model a time series in the $k^{th}$ cluster by a Gaussian vector autoregression (VAR), of block order $p_{k}$. The conditional likelihood function of the $k^{th}$ VAR model is

where for each $k$

Also $A_{k0}$ is a vector, while for $r\geq 1$, $A_{k,r}$ are the VAR coefficient matrices and $\Sigma_{k}$ the driving noise variance matrix. The set-up now fits in the k-MLE framework.

The classification log-likelihood is given by

where $\theta_{k}$ denotes the collection of parameters $A_{k}$ and $\Sigma_{k}$. Applying the coordinate ascent solver for k-MLE is straightforward, since we have a classic vector or multivariate regression which leads to the following algorithm:

• •

  label update:

  	$\displaystyle\tau_{n,k}$	$\displaystyle=\begin{cases}1&\text{if }k=\operatorname{arg\,min}_{k^{\prime}}D_{n,k^{\prime}},\\ 0&\text{otherwise},\end{cases}$		(19)
  	$\displaystyle D_{n,k^{\prime}}$	$\displaystyle=(T-p)\log|\Sigma_{k^{\prime}}|+\textstyle\sum_{t=p+1}^{T}\mathbf{e}_{nk^{\prime}t}^{T}\Sigma_{k^{\prime}}^{-1}\mathbf{e}_{nk^{\prime}t},$		(20)

• •

  parameter update:

  	$\displaystyle I_{k}$	$\displaystyle=\big{\{}n:\tau_{n,k}=1;\ n=1,\dots,N\big{\}},$		(21)
  	$\displaystyle A_{k}^{T}$	$\displaystyle=\textstyle\begin{array}[t]{@{}l@{}l}&\Big{(}\sum_{n\in I_{k}}\big{(}\sum_{t=p+1}^{T}X_{kt}^{(n)}X_{kt}^{(n)T}\big{)}\Big{)}^{-1}\\ &\Big{(}\sum_{n\in I_{k}}\big{(}\sum_{t=p+1}^{T}X_{kt}^{(n)}Y_{t}^{(n)T}\big{)}\Big{)},\end{array}$		(24)
  	$\displaystyle\Sigma_{k}$	$\displaystyle=\Big{(}\textstyle\sum_{n\in I_{k}}\big{(}\sum_{t=p+1}^{T}\mathbf{e}_{nkt}\mathbf{e}_{nkt}^{T}\big{)}\Big{)}/\Big{(}(T-p)|I_{k}|\Big{)},$		(25)

  where $|I_{k}|$ denotes the number of elements of $I_{k}$ (i.e. the cardinality of set $I_{k}$), and $p\triangleq\max\{p_{1},\dots,p_{K}\}$.

We now consider: initialisation, stopping conditions, and fast matrix computation.

Convergence analysis of k-VARs algorithm is provided by applying the general k-MLE results. As indicated in the remarks following Proposition 9 and Theorem 10 we will directly show the solution to the MLE equations is unique.

Note that the model is not from a (vector) exponential family. This is because the exponent can only be represented as an inner product of parameters with data if $\Sigma_{k}$ is known. So we cannot apply (a vector version of) k-Bregman.

Differentiating a cluster-specific likelihood gives the least squares equations shown in (21),(22). From these it is implicit that uniqueness requires that, for each $k$: (i) the inverse of the ${\bf X}_{nk}^{T}{\bf X}_{nk}$ matrix in (21) exists; (ii) $\Sigma_{k}$ in (21) has full rank. This will follow if (i),(ii) hold for every cluster of size 1. But we are now in a classic situation of requiring uniqueness of the least squares estimate of a VAR and of its noise variance from a single time series record. For cluster size 1, the ${\bf X}_{nk}^{T}{\bf X}_{nk}$ matrix in (22) is very nearly a block Toeplitz matrix. So long as each time series is generated by (i) a stable VAR with (ii) full rank noise covariance matrix and $T>p+1$, then the ${\bf X}_{nk}^{T}{\bf X}_{nk}$ matrix will have full rank with probability 1 [35]. Further the estimated noise variance matrix will have full rank with probability 1 [35]. So uniqueness of the MLE is established and convergence of the k-VARs algorithm now follows.

The clustering algorithms, presented in section IV-A and IV-B, require the number of clusters $K$ and cluster-specific model orders $(p_{1},\dots,p_{K})$ to be specified. This section develops a Bayesian information criterion (BIC) to resolve the model selection problem. We only discuss the special case with equal $p_{k}$’s. The general case is computationally challenging and may not offer big improvements on large data sets.

For the case of equal $p_{k}$’s, we choose $K$ and $p$ as the joint minimizer of BIC on a $(K,p)$ grid as follows.

Since the $N$ time series are assumed to be sampled independently and the $K$ models are independent then the likelihood is given by (18) with all $p_{k}$ in $D_{n,k}(\cdot,\cdot)$ being set to $p$.

A faster way to find the minimum of BIC is to use cyclic descent. Here we alternate between the following two steps until a convergence criterion is met.
(i) Given $p$, minimize over $K$;
(ii) Given $K$, minimize over $p$.
Of course this may get stuck in a local minimum. But that could be managed by trying a number of random starts.

k-VARs is implemented in MATLAB. All simulations are executed on 2.4 GHz Intel Core i9 machines with 32GB RAM and macOS Catalina. The ‘state of the art’ comparator methods are as follows.

• •

  k-DBA: k-Means with DTW barycenter averaging (DBA) [19], the implementation TimeSeriesKMeans with option metric="dtw" from an official Python package [36], tslearn, version 0.5.0.5.

• •

  k-Shape: the well-known algorithm proposed in [22, 23], the official implementation KShape from the same Python package tslearn.

• •

  k-GAK: Kernel k-Means [37] with Global Alignment Kernel [38] for time series clustering, the implementation KernelKMeans with option kernel="gak" from the Python package tslearn.

• •

  k-SC: k-Spectral Centroid proposed in [21], the MATLAB implementation given on the companion website of SNAP software [39]. (Note: this version does not include the Haar-Wavelet based incremental implementation.)

To avoid the drawbacks of Rand Index (RI) and the squared version of Normalised Mutual Information ($\text{NMI}_{\text{sqrt}}$) discussed in [40], we used the improved measures, adjusted Rand Index (ARI) and Normalized Information Distance (NID) to evaluate clustering performance

where $N_{11}$ is the number of pairs that are in the same cluster in both $U$ and $V$, $N_{00}$ is the number of pairs that are in different clusters in both $U$ and $V$, $N_{01}$ is the number of pairs that are in the same cluster in $U$ but in different clusters in $V$, $N_{10}$ is the number of pairs that are in different clusters in $U$ but in the same cluster in $V$; and $H(U)$ is the average amount of information in $U$, $I(U,V)$ is the mutual information, given as

in which $n_{ij}$ denotes the number of objects that are common to clusters $U_{i}$ and $V_{j}$, $a_{i}$ is the number of objects in cluster $U_{i}$, and $b_{j}$ is that in cluster $V_{j}$. For both ARI and $\mathrm{NMI}_{\text{max}}$, values close to 1 indicate high clustering performance.

To evaluate performance in a statistical manner, we randomly generate $40$ independent datasets. Due to scalability limits on some of the competing methods, the problem size had to be restricted. Each time-series dataset is generated by simulating a class of VAR models with the following specifications:

• •

  time series dimensions: $m=2,4,8$;

• •

  model order: $p=5$;

• •

  time series length: $T=80$;

• •

  number of groups: $K=8$;

• •

  number of time series per cluster: $N_{c}=30$.

The parameters of each VAR model are generated randomly and roots scaled to ensure stability. The noise covariance matrix is computed from randomly generated Cholesky factors.

The results in Figure 1, show two versions of k-VARs (defined in the caption). k-VARs (oracle), is superior to the state-of-the-art methods. And k-VARs (rnd) is competitive with them. Note also that most of the other methods require a number of tuning parameters to be chosen. We used default values suggested in the various implementations. k-VARs has no free tuning parameters and is competitive because the other methods ignore the autocorrelation feature.

Next we illustrate BIC. The following set-up applied.

• •

  cluster-specific VAR models: $m=4$, $p=5$;

• •

  time series: $T=200$, $K=10$, $N_{c}=20$.

The grid for $K,p$ is $2\!:\!2\!:\!20\times 2\!:\!1\!:\!8$. log-BIC is plotted as a heatmap in Figure 2.

We see that BIC attains its minimum at the ground truth $K=10$. However there is a relatively flat region in the vicinity of the minimum. We found in the simulations that good clustering performance was tolerant to some flatness in BIC near its minimum.

Here we study two things. The ability of k-VARS to deal with non-Gaussian driving noise. And the effect of signal to noise ratio (SNR).

We use t-distributed driving noise with degrees of freedom (dof): $[5,10,50,100,10000]$. Larger dof is closer to Gaussianity. The setup is otherwise the same as for Figure 1, but with $m=3$. The results in Figure 3 illustrate the robustness of k-VARs with little loss of performance for low dof.

The next simulation studies the effect of SNR for Gaussian driving noise. We modify a definition of SNR due to [41]. Denote $\Pi$ as the zero-lag auto-covariance of a stationary VAR(p) time series $Y_{t}$ whose white noise variance is $\Sigma$. Then [41] uses SNR$=10\log\mbox{VNSR},\mbox{VSNR}=\frac{\lambda_{\mathrm{max}}(\Pi)}{\lambda_{\mathrm{max}}(\Sigma)}$. This has a scaling problem. If we multiply $Y_{t}$ by a diagonal scaling matrix $L\neq cI$ (where $c$ is a constant) then this SNR changes. So instead we use VSNR$=\frac{1}{2}\lambda_{max}(\Sigma^{-1/2}\Pi\Sigma^{-1/2})$ which is unaffected by matrix scaling. The divisor of $2$ is due to a second adjustment as follows. In steady state we can write $Y_{t}=\hat{Y}_{t}+\epsilon_{t}$ where $\hat{Y}_{t}$ is the one step ahead predictor and $\epsilon_{t}$ is the one step ahead prediction error and so has variance $\Sigma$. Also it is uncorrelated with $\hat{Y}_{t}$ which thus has variance $M=\Pi-\Sigma$. In the scalar case the when $M=\Sigma$ i.e. ‘signal’ and noise variance are equal this would be taken to yield a SNR of 0 dB. So to reproduce that in the vector case we divide by $2$. The model setup is the same as that for Figure 3, except we fixed $m=3$. We have to modify $A_{k}$ to achieve a given SNR. The trick is to scale the roots of the VAR polynomial $A_{k_{1}},\cdots,A_{k_{p}}$. The SNR values are [0,5,10,15,20] (dB), The results in Figure 4, show k-VARS to be reliable.