---
title: "Clustering by Merging K-Means Solutions with MergeKmeans"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Clustering by Merging K-Means Solutions with MergeKmeans}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 6, fig.height = 5, fig.align = "center")
set.seed(2026)
```

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.

```{r half-circles-data}
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:

```{r half-circles-fit}
fit <- MergeKmeans(x, K = 10, G = 2, linkage = "single", nstart = 25)
fit
table(estimated = fit$cluster, truth = truth)
```

```{r half-circles-plot}
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.

```{r overlap-map}
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:

```{r tree}
plot(fit, what = "tree")
recut(fit, omega.star = 0.01)
```

## 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:

```{r inscribed}
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)
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:

```{r chooseK}
chooseK(x, K = seq(2, 20, by = 2), nstart = 10)
```

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

```{r predict}
predict(fit, rbind(c(0, 1), c(1, -0.6)))
```

## 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:

```{r ward}
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](https://doi.org/10.1007/s00357-019-09314-8)