Akaike Information criterion for k-means

I am trying use the AIC & BIC for selecting the number of clusters for k-means. I found a post on Stackoverflow where some R code is posted, but I could not find any valid explanation for its correctness.

The code:

kmeansAIC = function(fit){
m = ncol(fit$centers) n = length(fit$cluster)
k = nrow(fit$centers) D = fit$tot.withinss
return(D + 2*m*k)
}


basically states

$$AIC = -2 \mathcal{L}(\hat{\theta}) + 2p = (\sum_{l=1}^k\sum_{i=1}^{n_l}(\vec{x_{i,l}}-\vec{\mu_l})^2)+2 d k$$ where $k$ is the number of clusters, $n_l$ is the number of point in cluster $l$, $\mu_l$ is the cluster center of cluster $l$ and $d$ is the dimensionality of each point.

Now my problem is with the derivation of the log-likelihood function. If I consider k-means as a form of an EM algorithm with spherical clusters I get a gaussian cluster $cl_i$ with density $f_{cl_i}$ like

$$f_{cl_i}(\vec{x}, \vec{\mu_i}, \mathbf{\Sigma}) = \frac{1}{\sqrt{(2 \pi)^d det(\mathbf{\Sigma}})} \exp(-\frac{1}{2} (\vec{x}-\vec{\mu})^T \mathbf{\Sigma}^{-1} (\vec{x}-\vec{\mu}))$$ with $\bf{\Sigma}$ as covariance of the cluster $cl_i$.

As the cluster is spherical, $\bf{\Sigma}$ can be written as $\bf{\Sigma}=\sigma^2 \bf{I}$ where $\bf{I}$ is the identiy matrix. This gives $f_{cl_i}$ as $$f_{cl_i}(\vec{x}, \vec{\mu_i}, \Sigma) = \frac{1}{\sqrt{(2 \pi \sigma^2)^d}} \exp(-\frac{(\vec{x}-\vec{\mu})^2}{2 \sigma^2})$$ and therefore a log-likelihood function $\mathcal{L}$ for a cluster \begin{aligned} \mathcal{L}(\Theta;\vec{x}) = \sum_{i=1}^{n} log( \frac{1}{\sqrt{(2 \pi \sigma^2)^d}} \exp(-\frac{(\vec{x_i}-\vec{\mu})^2}{2 \sigma^2}) ) = \\ - n \log( \sqrt{(2 \pi \sigma^2)^d}) - \frac{1}{2 \sigma^2} \sum_{i=1}^{n} (\vec{x_i}-\vec{\mu})^2 \end{aligned}

I understand, that the first term can be omitted, as I am only interested in differences of the AIC for different $k$ and that term is constant over $k$.

This leads me to the following formula for the AIC: $$AIC = \frac{1}{\sigma^2} (\sum_{l=1}^k\sum_{i=1}^{n_l}(\vec{x_{i,l}}-\vec{\mu_l})^2) + 2 k d$$

which differs in the factor $\frac{1}{\sigma^2}$ compared to the formula from Stackoverflow.

I know, that $\sigma^2$ does not matter for k-means itself but it is obvious, that a greater $\sigma^2$ penalizes the model complexity in the AIC. Thus it seems strange to me, to just assume $\sigma^2 = 1$.

Is there an error in my line of argumentation or am I missing something?

I also used the same concept for calculation of the Bayesian Information Criterion and observed, that I end up with a suspiciously large number of clusters if I select the minimum of the BIC. The BIC drops very fast in the beginning and then decreases very slowly.

Is it common to select an AIC/BIC value where the scoring function doesn't show big decreases anymore (similar to the elbow technique) or should it always be the minimum of the function?

• I'm confused. Isn't the $\sigma$ dependent on the cluster? In other words, should they be $\sigma_i$? Aug 1, 2019 at 15:32

I think that the answer you refer to in the link in Stackoverflow, is in turn based on another link: http://sherrytowers.com/2013/10/24/k-means-clustering/.

In that post the data used is standardized the variance is equal to 1 and the formula makes sense. In any case I think your line of reasoning is absolutely correct, and one should refrain to use that expression if the underlying assumption is not correct.

I know this is an old (and answered) question, but I wanted to jump in because I ran into this same question, and was surprised by how little information there was available.

Before diving in to my answer, I wanted to echo what was said by FinanceGuyThatCantCode: you have a small error in your proposal. The $$\sigma$$ in your question really should be $$\sigma_i$$, since it depends on each cluster. So what you're asking is not why total variance is assumed equal to 1, but rather why each cluster variance is assumed equal to 1.

The answer is that the $$\sigma_i=1$$ assumption is actually an assumption of the k-means algorithm. Lets extend your analogy that k-means is really a latent Gaussian problem, now look at the k-means loss function: $$\min_{\mu_i} \sum_i \sum_{x_i \in S_i} |x_i - \mu_i|^2$$
Where $$S_i$$ is the set of $$x_i$$ assigned to center $$\mu_i$$. If we imagine that this loss function is the negative log-likelihood, and that what we're really doing is a likelihood maximization of the latent Gaussian means, we'd get the following equivalent likelihood maximization: $$\max_{\mu_i} \prod_i \prod_{x_i \in S_i} \exp \left[ (x_i - \mu_i)^2 \right]$$
Playing a little fast and loose with scaling constants, this is the same as: $$\max_{\mu_i} \prod_i \prod_{x_i \in S_i} N(x_i;\mu_i,\sigma_i^2 = 1)$$
Where $$N(x_i;\mu_i,\sigma_i^2 = 1)$$ is the normal PDf.

So we see that k-means implicitly assumes that the latent Gaussians have unit variance. This is because k-means isn't really meant to be a probabilistic model of the data; it's not a "generative" model which is intended to predict future datapoints. It's just an algorithm which uses a reasonable heuristic to produce label assignments.

Speculation incoming: I believe that if you wanted to relax the assumption that $$\sigma_i=1$$, then your best estimate would probably be something like $$\frac{1}{|S_i|-1} \sum_{x_i \in S_i} |x_i-\mu_i|^2$$, and so the AIC you propose (with the correction of $$\sigma$$ to $$\sigma_i$$) would collapse into something like: $$\sum_i |S_i| + 2 k d = n + 2 k d$$
Which is just a direct penalty on $$k$$.

• Is there a minus sign inside exp[(xi-ui)^2] for probability<1? Sep 21, 2021 at 15:13

The model for which k-means is maximum likelihood is not a mixture model (there are no mixture proportions, and the classification of observations to clusters is crisp whereas in mixture models there are posterior probabilities for this), but a model for which every observation has its own parameter indicating to which cluster it belongs. This is nonstandard as the number of parameters converges to infinity, and standard maximum likelihood theory doesn't hold (it is well known that the estimated k-means are not consistent for the means of the corresponding model, which is contrary to standard ML theory). This means that also any theory behind AIC and BIC doesn't hold, and these methods are theoretically invalid. (Obviously it may be that they do a good job in some situations anyway; note also that it doesn't matter whether these strange discrete observation parameters are counted into the general number of parameters or not, because for a given dataset their number is constant. One would need to deal with them though in any attempt to derive asymptotic theory.)

The model does not assume $$\sigma^2=1$$, and it is invalid as well to assume $$\sigma^2$$ to be dependent of the data. In fact it makes some sense to have a penalty that depends on $$\sigma^2$$, because if the data have a certain variation, more clusters with a low variance are needed to cover them in a Voronoi tesselation fashion (which is what k-means does) than if the clusters have a larger variance. In fact in this way using the AIC you could also choose a variance estimator, which is connected to the number of clusters in that fashion - smaller variance, more clusters needed. This looks surprising because the variance in k-means is usually ignored because people are only interested in the clustering for given k.

Despite this, AIC and BIC are invalid anyway, for the reason stated above, or rather, nobody has yet demonstrated how they may be valid, and the standard arguments don't apply.

An alternative to assume $$\sigma^2=1$$, is to assume \begin{aligned} \sigma^2 = \text{average residual sum of square} = \\ \frac{RSS}{\text{number of cell}} = \\ \frac{\sum_{l=1}^k\sum_{i=1}^{n_l}(\vec{x_{i,l}}-\vec{\mu_l})^2}{nd} \end{aligned}

The resulting formula seems to agree with the BIC in the Gaussian special case, where

The modified code:

kmeansAIC = function(fit){
d = ncol(fit$$centers) n = length(fit$$cluster)
k = nrow(fit$$centers) D = fit$$tot.withinss

# return(D + 2*d*k) # original
# return(n*d*log(D/n/d) + 2*d*k) # for individual cell
return(n*log(D/n) + 2*k) # for individual row
}


Following your log-likelihood function $$\mathcal{L}$$ for a cluster \begin{aligned} \mathcal{L}(\Theta;\vec{x}) = \sum_{i=1}^{n} log( \frac{1}{\sqrt{(2 \pi \sigma^2)^d}} \exp(-\frac{(\vec{x_i}-\vec{\mu})^2}{2 \sigma^2}) ) = \\ - n \log( \sqrt{(2 \pi \sigma^2)^d}) - \frac{1}{2 \sigma^2} \sum_{i=1}^{n} (\vec{x_i}-\vec{\mu})^2 \end{aligned} Instead of assuming $$\sigma^2=1$$, we get \begin{aligned} =-\frac{nd}{2}(log(2 \pi) + log( \sigma^2)) - \frac{1}{2 \sigma^2} \sum_{i=1}^{n} (\vec{x_i}-\vec{\mu})^2 = \\ -\frac{nd}{2}log(2 \pi) -\frac{nd}{2} log( \frac{RSS}{nd}) - \frac{1}{2 \frac{RSS}{nd}}\text{RSS} = \\ -\frac{nd}{2} log( \frac{RSS}{nd}) + Constant(nd) \end{aligned} which brings us back to the special case BIC above except that nd refers to number of cells and n to number of rows.

It doesn't matter anyway because you typically compare which AIC is better, in which case the constant doesn't matter (if $$A > B$$, $$cA > cB$$ for positive $$c$$).

• The problem with this is that one often is interested in how different the two values might be. It helps to use a universal unit of measure for AIC.
– whuber
Sep 8 at 15:49