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   0 500
#>         2 500   0
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 = 8 (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] 1 6
#> 
#> $cluster
#> [1] 1 2

Cutting the tree

If neither G nor omega.star is supplied, the fit is returned uncut, together with the number of clusters suggested by the largest gap in the overlap structure. The tree can then be cut at a chosen number of clusters, at an overlap threshold, or at the automatic suggestion, without re-running K-means:

fit0 <- MergeKmeans(x, K = 10, linkage = "single", nstart = 25)
#> Merge tree built but not cut (no 'G' or 'omega.star'): the overlap structure suggests G = 2. Inspect plot(fit, what = "overlap") and call recut(fit, G = ...)
plot(fit0, what = "tree")

recut(fit0, G = "auto")$G
#> G = 2 selected by the largest gap in the overlap structure; confirm with plot(fit, what = "overlap")
#> [1] 2
recut(fit0, omega.star = 0.01)$G
#> [1] 2

The automatic suggestion (G = "auto") targets well-separated clusters, where the boundary between clusters shows up as a drop of orders of magnitude in the overlaps along the merge sequence; it is refused with Ward’s linkage, and it should always be confirmed against the overlap map.

Merging an external partition

The K-means stage is the only material cost of the procedure, and it can be delegated to any engine through the start argument, which accepts either labels (a kmeans object, a list with a cluster element, or a bare label vector) or centers (a K x p matrix, or a list with a centers/centroids element, as returned by ClusterR::MiniBatchKmeans()). Centers are converted to labels by nearest-center assignment:

km <- kmeans(x, 12, nstart = 10)
fit.lab <- MergeKmeans(x, start = km, G = 2, linkage = "single")
fit.ctr <- MergeKmeans(x, start = km$centers, G = 2, linkage = "single")
table(labels = fit.lab$cluster, centers = fit.ctr$cluster)
#>       centers
#> labels   1   2
#>      1 500   0
#>      2   0 500

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.