# --- Simulate signal // groups A & B ---
n <- 600
signal <- rep(c("A", "B"), each = n / 2)
Z <- matrix(0, nrow = n, ncol = 2)
Z[signal == "A", ] <- cbind(rnorm(n / 2, -0.8, 0.8), rnorm(n / 2, 0.0, 0.5))
Z[signal == "B", ] <- cbind(rnorm(n / 2, 0.8, 0.8), rnorm(n / 2, 0.0, 0.5))
# --- Assign batches independently ---
batch <- sample(c("Batch 1", "Batch 2"), n, replace = TRUE)
# --- Add strong batch effect ---
shift <- c(0., 2.)
Z[batch == "Batch 2", ] <- sweep(Z[batch == "Batch 2", , drop = FALSE], 2, shift, "+")Maximum Diversity Clustering
1 Introduction
The mdclust package provides an independent implementation of the maximum diversity clustering algorithm introduced by Korsunsky et al. (2019) as the clustering component of the Harmony framework for single-cell data integration. Harmony popularized maximum diversity clustering as an effective way of constructing clusters that are robust to batch effects. However, the clustering algorithm is tightly integrated into Harmony’s embedding correction workflow and is not readily reusable in other contexts. mdclust provides an independent implementation of the clustering algorithm, exposing all optimization parameters and intermediate quantities for standalone use and integration into custom workflows. mdclust supports user-defined centroid initialization, making it suitable for iterative optimization procedures requiring repeated diversity clustering.
This vignette introduces the basic workflow for using mdclust as both a standalone clustering method and as a reusable component in custom analysis pipelines.
2 Background
The objective optimized by mdclust consists of three components:
- a spherical clustering objective that encourages observations to be close to their assigned centroids,
- an entropy term that allows observations to belong probabilistically to multiple clusters, and
- a diversity penalty that promotes balanced batch composition within each cluster.
Setting the diversity penalty (theta = 0) recovers ordinary soft spherical k-means, while increasing theta progressively favours clusters with improved batch mixing.
Throughout this vignette we use a small synthetic data set to illustrate the behaviour of the algorithm. The same workflow can be applied directly to embeddings obtained from common single-cell analysis pipelines, including principal component analysis, latent semantic indexing, variational autoencoders, or other dimensionality reduction methods.
3 A complete example
We simulate two signal groups A and B separated but overlapping Gaussians in 2D, and randomly assign them to two batches Batch 1 and Batch 2. A strong batch effect is added by offsetting one of the batches along the y-axis.
We can plot the simulated data to confirm that, indeed, we have a strong batch separation of the biological signal. Ordinary clustering (e.g. k-means) tends to struggle with such data and overwhelmingly aligns with undesired batch variation.
Code
library(ggplot2)
library(patchwork)
df <- data.frame(
Dim1 = Z[, 1],
Dim2 = Z[, 2],
signal = signal,
batch = batch
)
p_truth <- ggplot(df, aes(Dim1, Dim2)) +
geom_point(
aes(colour = signal),
size = 2,
alpha = 0.8
) +
coord_equal() +
labs(
title = "Ground truth",
colour = "Signal",
x = "Dimension 1",
y = "Dimension 2"
) +
scale_colour_brewer(palette = "Dark2") +
theme_classic(base_size = 12)
p_batch <- ggplot(df, aes(Dim1, Dim2)) +
geom_point(
aes(colour = batch),
size = 2,
alpha = 0.8
) +
coord_equal() +
labs(
title = "Batch effect",
colour = "Batch",
x = "Dimension 1",
y = "Dimension 2"
) +
scale_colour_brewer(palette = "Set2") +
theme_classic(base_size = 12)
p_truth + p_batch +
plot_layout(widths = c(1, 1))
A minimal call to mdclust on this data looks as follows:
library(mdclust)
set.seed(111)
mdc <- mdclust(t(Z), batch, K = 2)
mdc
<mdclust>
Elements : 600
Clusters : 2
Embedding dims : 2
Converged : yes
Iterations : 23
Final objective : 308.935
Components:
R soft assignments
cluster hard assignments
Y centroids
objective optimization trace
O, E, logOE batch diagnosticsPrinting the fitted object confirms that the optimization converged. The resulting clustering is shown below.
df$cluster <- factor(mdc$cluster)
p_mdclust <- ggplot(df, aes(Dim1, Dim2)) +
geom_point(
aes(colour = cluster),
size = 2,
alpha = 0.8
) +
coord_equal() +
labs(
title = "Maximum diversity clustering",
colour = "Cluster",
x = "Dimension 1",
y = "Dimension 2"
) +
scale_colour_brewer(palette = "Set1") +
theme_classic(base_size = 12)
p_truth + p_mdclust +
plot_layout(widths = c(1, 1))
The recovered clusters closely follow the underlying biological groups despite the pronounced batch effect, illustrating the purpose of the diversity penalty.
Any run of mdclust requires at least these three arguments. Here, t(Z) is the low-dimensional embedding used for clustering. Observations are in columns, so be mindful of orientation of your embedding. t(Z) should be of shape d × N, and an error will be raised if there is a mismatch. Columns are internally L2-normalized per column for spherical k-means. batch encodes the label that is penalized during the optimization. It can be either factor, character or integer and only one label per observation is permitted. K is simply the number of clusters and completely analogous to ordinary k-means.
However, mdclust exposes most internals to the user and the algorithm is fully configurable. The following sections describe the available initialization methods and optimization parameters.
3.1 Initialization
Maximum diversity clustering requires an initial matrix of centroid positions Y_init, and allows for three different initializations.
mdclust(..., Y_init = "kmeans++")
mdclust(..., Y_init = "kmeans")
mdclust(..., Y_init = Y)The default is "kmeans++", which is an ordinary k-means optimization with an enhanced starting point selection. "kmeans" performs a single k-means optimisation from a random starting point. In addition, Y_init can be assigned any valid centroid matrix Y of shape d × K. This allows for custom initialisations, and especially warm starts or use of mdclust as a diversity clustering engine in iterative procedures.
3.2 Clustering parameters
The behaviour of maximum diversity clustering is primarily controlled by four parameters. These are in principle analogous to their Harmony counterparts, hence more information on them is available in the Harmony documentation and the respective publication.
thetacontrols the strength of the diversity penalty. Settingtheta = 0disables the diversity objective entirely, reducing the algorithm to ordinary soft spherical k-means. Larger values increasingly favour clusters with balanced batch composition. The default (theta = 2) matches the Harmony implementation.sigmacontrols the softness of cluster assignments. Smaller values produce assignments closer to hard clustering, whereas larger values distribute observations across multiple clusters. The default issigma = 0.1.tauprotects very small batches from being over-penalized by reducing the effective diversity penalty for batches containing fewer than approximatelyK × tauobservations. The default (tau = 0) disables this behaviour.alphais a pseudocount used when computing observed-to-expected enrichment ratios for the diversity penalty. It improves numerical stability and rarely needs adjustment.
3.3 Optimization parameters
The remaining parameters control the optimization procedure.
block_sizespecifies the fraction of observations updated simultaneously. Smaller values provide more accurate online updates at the expense of runtime.max_itersets the maximum number of optimization iterations.epsilonspecifies the convergence threshold on the objective.
The k-means implementations have their own convergence controls. Convergence is declared once subsequent runs fall under the threshold abs(new - old)/abs(old) < epsilon_kmeans. If this condition is not met, k-means terminates after max_iter_kmeans. Both parameters are accessible with defaults mdclust(..., epsilon_kmeans = 1e-3, max_iter_kmeans = 10L).
3.4 Output
An additional concise overview to printing of the finished clustering can be accessed via
summary(mdc)
<mdclust summary>
------------------
Elements : 600
Clusters : 2
Embedding dims : 2
Converged : yes
Iterations : 23
Final objective
Total : 308.935
Distance : 311.769
Entropy : -2.061
Diversity : -0.774
Cluster sizes
1 2
270 330 The complete object contains almost all internal parameters of interest.
str(mdc, max.level = 1)List of 9
$ R : num [1:2, 1:600] 8.43e-15 1.00 5.73e-13 1.00 1.98e-04 ...
..- attr(*, "dimnames")=List of 2
$ Y : num [1:2, 1:2] 0.882 0.471 -0.707 0.707
..- attr(*, "dimnames")=List of 2
$ O : num [1:2, 1:2] 148 138 122 192
..- attr(*, "dimnames")=List of 2
$ E : num [1:2, 1:2] 129 157 141 173
..- attr(*, "dimnames")=List of 2
$ logOE : num [1:2, 1:2] -0.146 0.128 0.143 -0.11
..- attr(*, "dimnames")=List of 2
$ objective :List of 4
$ converged : logi TRUE
$ cluster : int [1:600] 2 2 2 2 2 1 1 2 2 2 ...
$ batch_levels: chr [1:2] "Batch 1" "Batch 2"
- attr(*, "class")= chr "mdclust"Here
R: Soft-assignment matrix of shape K × N. Columns sum to one.Y: Final centroids.cluster: Hard-assignment viamax.col(t(mdc$R)).objective: List containing the partial and full objective traces of the run.O,E,logOE: Expected and observed batch compositions. A useful diagnostic of the batch diversity objective.converged: Convergence flag.batch_levels: Batch levels used internally during optimization.
mdclust also supports rudimentary plots of the objective trace for quick visual diagnostics.
plot(mdc)
4 Relationship to Harmony
This package is an independently developed, standalone implementation of Algorithm 2 from Korsunsky et al. (2019). Compared with the current Harmony implementation (2.0.x), mdclust differs in a few key aspects.
Most notably, mdclust exposes all parameters of the underlying algorithm, making maximum diversity clustering available as a standalone method and facilitating its use in custom workflows.
The initialization procedure has also been extended. Whereas Harmony initializes diversity clustering using a k-means solution, mdclust additionally supports random initialization, k-means++, and user-supplied centroid matrices. Custom centroid initialization enables warm starts and simplifies the use of mdclust as a diversity clustering engine within iterative algorithms.
Finally, mdclust provides independent convergence criteria and optimization controls for both the initialization procedure and the diversity clustering optimization.
5 Session information
sessionInfo()R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 26.04 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: Etc/UTC
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] mdclust_0.99.2 patchwork_1.3.2 ggplot2_4.0.3
loaded via a namespace (and not attached):
[1] vctrs_0.7.3 cli_3.6.6 knitr_1.51 rlang_1.3.0
[5] xfun_0.60 otel_0.2.0 generics_0.1.4 S7_0.2.2
[9] jsonlite_2.0.0 labeling_0.4.3 glue_1.8.1 buildtools_1.0.0
[13] htmltools_0.5.9 maketools_1.3.2 sys_3.4.3 scales_1.4.0
[17] rmarkdown_2.31 grid_4.6.1 tibble_3.3.1 evaluate_1.0.5
[21] fastmap_1.2.0 yaml_2.3.12 lifecycle_1.0.5 compiler_4.6.1
[25] dplyr_1.2.1 RColorBrewer_1.1-3 Rcpp_1.1.2 pkgconfig_2.0.3
[29] farver_2.1.2 digest_0.6.39 R6_2.6.1 tidyselect_1.2.1
[33] pillar_1.11.1 magrittr_2.0.5 withr_3.0.3 tools_4.6.1
[37] gtable_0.3.6