The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Clustering by Merging K-Means Solutions with MergeKmeans

The K-means algorithm is fast and intuitive, but in its traditional form it can only succeed when clusters are spherical with approximately equal volumes. MergeKmeans implements the DEMP-K procedure of Melnykov and Michael (2020), which removes this restriction while preserving the speed of K-means:

  1. run K-means with a number of components \(K\) deliberately larger than the number of clusters sought, so that small spherical components fill the clusters whatever their shapes;
  2. interpret the solution as a Gaussian mixture fitted by classification EM and compute the pairwise overlap \(\omega_{kk'}\) (the sum of the two misclassification probabilities) between all pairs of components in closed form — for classical K-means simply \(\omega_{kk'} = 2\Phi\!\left(-\|\mu_k - \mu_{k'}\| / (2\sigma)\right)\);
  3. merge components by hierarchical clustering with \(1-\omega_{kk'}\) as the dissimilarity, and cut the tree at the requested number of clusters \(G\) or at an overlap threshold \(\omega^*\).

Because steps 2 and 3 operate on the \(K \times K\) overlap matrix rather than on the \(n\) observations, the whole post-processing cost is negligible and the procedure scales to millions of observations.

Two half-circular clusters

A classical failure case for K-means, taken from Section 4.1 of the paper: two well-separated half-circles.

library(MergeKmeans)

n <- 1000
th <- runif(n, 0, pi)
x <- rbind(
  cbind(cos(th[1:500]), sin(th[1:500])),
  cbind(1 - cos(th[501:1000]), 0.45 - sin(th[501:1000])))
x <- x + matrix(rnorm(2 * n, sd = 0.1), n)
truth <- rep(1:2, each = 500)
plot(x, col = truth, pch = 20, asp = 1, xlab = "x1", ylab = "x2")

Plain K-means with two centers cuts both half-circles in half. Instead we fit ten components and merge. Since the goal is to detect well-separated clusters of arbitrary shape, single linkage is the appropriate choice:

fit <- MergeKmeans(x, K = 10, G = 2, linkage = "single", nstart = 25)
fit
#> Merging K-means solutions (DEMP-K)
#>   variant: HoSC, K = 10 components, linkage: single
#>   n = 1000 observations, p = 2 variables
#>   G = 2 clusters of sizes 500, 500
table(estimated = fit$cluster, truth = truth)
#>          truth
#> estimated   1   2
#>         1 499   1
#>         2   1 499
plot(fit, what = "classification", asp = 1)

The overlap map

When the number of clusters is unknown, the overlap map displays all pairwise overlaps with components ordered so that strongly overlapping ones are adjacent. Groups of dark cells correspond to components modeling a common cluster; pale cells in the bottom summary row mark the gaps between well-separated clusters — here, one pale cell splits the ten components into two groups of five.

plot(fit, what = "overlap")

The merge tree offers a complementary view, and recut() lets you cut it at a different number of clusters (or at an overlap threshold \(\omega^*\)) without refitting:

plot(fit, what = "tree")

recut(fit, omega.star = 0.01)
#> Merging K-means solutions (DEMP-K)
#>   variant: HoSC, K = 10 components, linkage: single
#>   n = 1000 observations, p = 2 variables
#>   G = 2 clusters of sizes 500, 500

An inscribed cluster

The second illustration of the paper: a circular cluster inscribed in a ring-shaped one. No rigid-motion of two spherical prototypes can separate these, yet merging a 12-component solution recovers the structure:

m <- 500
th2 <- runif(m, 0, 2 * pi)
y <- rbind(cbind(rnorm(m, 3, 0.2) * cos(th2), rnorm(m, 3, 0.2) * sin(th2)),
           matrix(rnorm(2 * m, sd = 0.4), m))
truth2 <- rep(1:2, each = m)
fit2 <- MergeKmeans(y, K = 12, G = 2, linkage = "single", nstart = 25)
table(estimated = fit2$cluster, truth = truth2)
#>          truth
#> estimated   1   2
#>         1 500   0
#>         2   0 500
plot(fit2, what = "classification", asp = 1)

Choosing K and predicting new data

The merging procedure is robust to the choice of \(K\) as long as \(K\) is large enough for components to fill the clusters. A practical guide is the elbow of the proportion of variability explained by the between-component variability:

chooseK(x, K = seq(2, 20, by = 2), nstart = 10)

#> Recommended K = 10 (one step beyond the elbow; K must exceed the number of clusters sought)

New observations are assigned to the nearest component and then mapped to its merged cluster:

predict(fit, rbind(c(0, 1), c(1, -0.6)))
#> $component
#> [1] 2 3
#> 
#> $cluster
#> [1] 1 2

Compact overlapping clusters

For compact, possibly overlapping clusters, Ward’s linkage is the recommended choice (it is also the package default). Here three Gaussian blobs with unequal spreads are recovered from a 15-component solution:

z <- rbind(matrix(rnorm(600), ncol = 2),
           matrix(rnorm(600, mean = 4, sd = 1.4), ncol = 2),
           cbind(rnorm(300, 8), rnorm(300)))
fit3 <- MergeKmeans(z, K = 15, G = 3, linkage = "ward.D", nstart = 25)
plot(fit3, what = "classification", asp = 1)

Beyond classical K-means

MergeKmeans() also provides the three generalized K-means variants of the paper’s Appendix A, each corresponding to a Gaussian mixture with a different covariance restriction: "HoEC" (common elliptical covariance), "HeSC" (component-specific spherical variances), and "HeEC" (unrestricted covariances; the overlap computation then uses the Davies algorithm from the CompQuadForm package). The spherical "HoSC" default remains the recommendation for large datasets.

Reference

Melnykov, V. and Michael, S. (2020). Clustering large datasets by merging K-means solutions. Journal of Classification, 37, 97–123. doi:10.1007/s00357-019-09314-8

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.