Clustering Large Datasets by Merging K-Means Solutions

Fast clustering of large datasets by hierarchically merging components of a K-means solution based on the pairwise overlap between the Gaussian mixture components implied by the K-means partition, as proposed by Melnykov and Michael (2020) . Implements the DEMP-K merging algorithm with single, Ward's, average, and complete linkages, the overlap map display for selecting the number of clusters, four K-means variants corresponding to Gaussian mixtures with spherical or elliptical, homoscedastic or heteroscedastic components, and a tool for selecting the number of K-means components.

MergeKmeans

Clustering large datasets by merging K-means solutions — an R implementation of the DEMP-K algorithm of Melnykov and Michael (2020), Journal of Classification, 37, 97–123 (doi:10.1007/s00357-019-09314-8).

A K-means solution with a deliberately large number of components K is interpreted as a Gaussian mixture fitted by classification EM; the pairwise overlap between all component pairs is computed in closed form; and components are merged by hierarchical clustering with 1 - overlap as the dissimilarity. The procedure keeps the speed of K-means (the post-processing operates on the K x K overlap matrix, not on the data) while recovering clusters of arbitrary shape, making it practical for millions of observations.

Installation

install.packages("MergeKmeans")

Quick start

library(MergeKmeans)

## two half-circular clusters: a classical K-means failure case
set.seed(7)
th <- runif(1000, 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(2000, sd = 0.1), 1000)

## fit 10 components, merge down to 2 clusters
fit <- MergeKmeans(x, K = 10, G = 2, linkage = "single", nstart = 25)
fit

plot(fit, what = "overlap")          # the overlap map (choose G visually)
plot(fit, what = "tree")             # dendrogram of the merging
plot(fit, what = "classification")   # final partition

recut(fit, G = 3)                    # re-cut the tree without refitting
predict(fit, rbind(c(0, 1)))         # classify new observations
chooseK(x, K = seq(2, 20, 2))        # SSB/SSTO elbow for selecting K

Features

• The four K-means variants of the paper: HoSC (classical K-means, default, recommended for large data), HoEC, HeSC, HeEC (Appendix A), with the variant-implied mixing proportions.
• Pairwise overlap computations for all four covariance structures, for the implied or for user-supplied mixing proportions; the general heteroscedastic elliptical case uses the Davies (1980) algorithm via the CompQuadForm package.
• Single, Ward's (ward.D), average, and complete linkage merging; tree cutting by number of clusters G or overlap threshold omega.star.
• The overlap map display of Section 3.3 for choosing the number of clusters.
• print, summary, plot, predict, as.hclust, as.dendrogram methods.

All overlap computations are validated in the test suite against MixSim::overlap() (Maitra and Melnykov, 2010); the heteroscedastic spherical overlap is additionally validated by Monte Carlo simulation in the companion paper's replication materials.

Guidance

• Pick K comfortably larger than the number of clusters you expect; the method is robust to this choice (see chooseK(), which reports a recommendation, or leave K = NULL for automatic selection).
• Use single linkage to find well-separated clusters of arbitrary shape, Ward's linkage for compact clusters.
• G = "auto" cuts at the largest gap, on the logarithmic overlap scale, in the single-linkage merge sequence of the overlap matrix; it targets well-separated clusters, is not available with Ward's linkage, and should be confirmed with the overlap map.

When to use MergeKmeans — and when not

Benchmarks against the main alternatives (half-moons, varying-density, and overlapping-Gaussian scenarios at n up to 10^6):

• Use MergeKmeans for arbitrary-shape, skewed, or varying-density clusters, especially at large n: it matches DBSCAN/HDBSCAN accuracy without a fragile eps radius or their memory ceilings (HDBSCAN exhausts 24 GB near n = 10^5; spectral clustering is impractical past n ~ 10^3), and it has a predict() for new data.
• Use mclust (or plain kmeans) when clusters genuinely are overlapping Gaussian ellipsoids: one component per cluster is then the right model, and over-segmentation followed by merging gives away accuracy at the cluster boundaries.

Massive datasets

The K-means step is the only material cost. Run it with any fast engine and hand the labels to MergeKmeans() via start=:

mb  <- ClusterR::MiniBatchKmeans(x, clusters = 30)
fit <- MergeKmeans(x, start = mb, G = "auto", linkage = "single")  # centroids suffice

License

GPL (>= 2)