Determining the number of clusters in a data set

X-means clustering

In statistics and data mining, X-means clustering is a variation of k-means clustering that refines cluster assignments by repeatedly attempting subdivision, and keeping the best resulting splits, until a criterion such as the Akaike information criterion (AIC) or Bayesian information criterion (BIC) is reached.

Information criterion approach

Another set of methods for determining the number of clusters are information criteria, such as the Akaike information criterion (AIC), Bayesian information criterion (BIC), or the deviance information criterion (DIC) — if it is possible to make a likelihood function for the clustering model. For example: The k-means model is "almost" a Gaussian mixture model and one can construct a likelihood for the Gaussian mixture model and thus also determine information criterion values.

Information–theoretic approach

Rate distortion theory has been applied to choosing k called the "jump" method, which determines the number of clusters that maximizes efficiency while minimizing error by information-theoretic standards. The strategy of the algorithm is to generate a distortion curve for the input data by running a standard clustering algorithm such as k-means for all values of k between 1 and n, and computing the distortion (described below) of the resulting clustering. The distortion curve is then transformed by a negative power chosen based on the dimensionality of the data. Jumps in the resulting values then signify reasonable choices for k, with the largest jump representing the best choice.

The distortion of a clustering of some input data is formally defined as follows: Let the data set be modeled as a p-dimensional random variable, X, consisting of a mixture distribution of G components with common covariance, Γ. If we let ${\displaystyle c_1\dots c_K}$ be a set of K cluster centers, with ${\displaystyle c_X}$ the closest center to a given sample of X, then the minimum average distortion per dimension when fitting the K centers to the data is:

${\displaystyle d_K=\frac{1}{p}\underset{c_1\dots c_K}{min}E[(X-c_X{)}^T{\mathrm{\Gamma }}^{-1}(X-c_X)]}$

This is also the average Mahalanobis distance per dimension between X and the closest cluster center ${\displaystyle c_X}$. Because the minimization over all possible sets of cluster centers is prohibitively complex, the distortion is computed in practice by generating a set of cluster centers using a standard clustering algorithm and computing the distortion using the result. The pseudo-code for the jump method with an input set of p-dimensional data points X is:

JumpMethod(X): <b id="mwdg">Let</b> Y = (p/2) <b id="mwdw">Init</b> a list D, of size n+1 <b id="mweA">Let</b> D[0] = 0 <b id="mweQ">For</b> k = 1 ... n: Cluster X with k clusters (e.g., with k-means) <b id="mweg">Let</b> d = Distortion of the resulting clustering D[k] = d^(-Y) <b id="mwew">Define</b> J(i) = D[i] - D[i-1] <b id="mwfA">Return</b> the k between 1 and n that maximizes J(k)

The choice of the transform power ${\displaystyle Y=(p/2)}$ is motivated by asymptotic reasoning using results from rate distortion theory. Let the data X have a single, arbitrarily p-dimensional Gaussian distribution, and let fixed ${\displaystyle K=\lfloor {\alpha }^p\rfloor }$, for some α greater than zero. Then the distortion of a clustering of K clusters in the limit as p goes to infinity is ${\displaystyle {\alpha }^{-2}}$. It can be seen that asymptotically, the distortion of a clustering to the power ${\displaystyle (-p/2)}$ is proportional to ${\displaystyle {\alpha }^p}$, which by definition is approximately the number of clusters K. In other words, for a single Gaussian distribution, increasing K beyond the true number of clusters, which should be one, causes a linear growth in distortion. This behavior is important in the general case of a mixture of multiple distribution components.

Let X be a mixture of G p-dimensional Gaussian distributions with common covariance. Then for any fixed K less than G, the distortion of a clustering as p goes to infinity is infinite. Intuitively, this means that a clustering of less than the correct number of clusters is unable to describe asymptotically high-dimensional data, causing the distortion to increase without limit. If, as described above, K is made an increasing function of p, namely, ${\displaystyle K=\lfloor {\alpha }^p\rfloor }$, the same result as above is achieved, with the value of the distortion in the limit as p goes to infinity being equal to ${\displaystyle {\alpha }^{-2}}$. Correspondingly, there is the same proportional relationship between the transformed distortion and the number of clusters, K.

Putting the results above together, it can be seen that for sufficiently high values of p, the transformed distortion ${\displaystyle d_K^{-p/2}}$ is approximately zero for K < G, then jumps suddenly and begins increasing linearly for K ≥ G. The jump algorithm for choosing K makes use of these behaviors to identify the most likely value for the true number of clusters.

Although the mathematical support for the method is given in terms of asymptotic results, the algorithm has been empirically verified to work well in a variety of data sets with reasonable dimensionality. In addition to the localized jump method described above, there exists a second algorithm for choosing K using the same transformed distortion values known as the broken line method. The broken line method identifies the jump point in the graph of the transformed distortion by doing a simple least squares error line fit of two line segments, which in theory will fall along the x-axis for K < G, and along the linearly increasing phase of the transformed distortion plot for K ≥ G. The broken line method is more robust than the jump method in that its decision is global rather than local, but it also relies on the assumption of Gaussian mixture components, whereas the jump method is fully non-parametric and has been shown to be viable for general mixture distributions.

Silhouette method

The average silhouette of the data is another useful criterion for assessing the natural number of clusters. The silhouette of a data instance is a measure of how closely it is matched to data within its cluster and how loosely it is matched to data of the neighboring cluster, i.e., the cluster whose average distance from the datum is lowest. A silhouette close to 1 implies the datum is in an appropriate cluster, while a silhouette close to −1 implies the datum is in the wrong cluster. Optimization techniques such as genetic algorithms are useful in determining the number of clusters that gives rise to the largest silhouette. It is also possible to re-scale the data in such a way that the silhouette is more likely to be maximized at the correct number of clusters.

Cross-validation

One can also use the process of cross-validation to analyze the number of clusters. In this process, the data is partitioned into v parts. Each of the parts is then set aside at turn as a test set, a clustering model computed on the other v − 1 training sets, and the value of the objective function (for example, the sum of the squared distances to the centroids for k-means) calculated for the test set. These v values are calculated and averaged for each alternative number of clusters, and the cluster number selected such that further increase in number of clusters leads to only a small reduction in the objective function.

Finding number of clusters in text databases

When clustering text databases with the cover coefficient on a document collection defined by a document by term D matrix (of size m×n, where m is the number of documents and n is the number of terms), the number of clusters can roughly be estimated by the formula ${\displaystyle \frac{mn}{t}}$ where t is the number of non-zero entries in D. Note that in D each row and each column must contain at least one non-zero element.

Analyzing the kernel matrix

Kernel matrix defines the proximity of the input information. For example, in Gaussian radial basis function, it determines the dot product of the inputs in a higher-dimensional space, called feature space. It is believed that the data become more linearly separable in the feature space, and hence, linear algorithms can be applied on the data with a higher success.

The kernel matrix can thus be analyzed in order to find the optimal number of clusters. The method proceeds by the eigenvalue decomposition of the kernel matrix. It will then analyze the eigenvalues and eigenvectors to obtain a measure of the compactness of the input distribution. Finally, a plot will be drawn, where the elbow of that plot indicates the optimal number of clusters in the data set. Unlike previous methods, this technique does not need to perform any clustering a-priori. It directly finds the number of clusters from the data.

The gap statistics

Robert Tibshirani, Guenther Walther, and Trevor Hastie proposed estimating the number of clusters in a data set via the gap statistic. The gap statistics, based on theoretical grounds, measures how far is the pooled within-cluster sum of squares around the cluster centers from the sum of squares expected under the null reference distribution of data. The expected value is estimated by simulating null reference data of characteristics of the original data, but lacking any clusters in it. The optimal number of clusters is then estimated as the value of k for which the observed sum of squares falls farthest below the null reference.

Unlike many previous methods, the gap statistics can tell us that there is no value of k for which there is a good clustering, but the reliability depends on how plausible the assumed null distribution (e.g., a uniform distribution) is on the given data. This tends to work well in synthetic settings, but cannot handle difficult data sets with, e.g., uninformative attributes well because it assumes all attributes to be equally important.

The gap statistics is implemented as the clusGap function in the cluster package in R.