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:

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.