K-means clustering using MUS and other pivotal algorithms

Leonardo Egidi


In this vignette we explore the K-means algorithm performed using the MUS algorithm and other pivotal methods through the function piv_KMeans of the pivmet package. First of all, we load the package:


Pivotal algorithms: how they works, and why

We present here a simulated case for applying our procedure. Given \(n\) units \(y_1,\ldots,y_n\):

We propose four alternative methods for achieving this task. Let \(j\) be the group containing units \(\mathcal J_j\), the user may choose \({i^*}\in\mathcal J_j\) that maximizes one of the quantities:

\[\begin{align*} & \sum_{p\in\mathcal J_j} c_{{i^*}p} \\ & \sum_{p\in\mathcal J_j} c_{{i^*}p} - \sum_{j\not\in\mathcal J_j} c_{{i^*}p}. \end{align*}\]

These methods give the unit that maximizes the global within similarity (maxsumint) and the unit that maximizes the difference between global within and between similarities (maxsumdiff), respectively. Alternatively, we may choose \(i^{*} \in\mathcal J_j\), which minimizes:

\[\sum_{p\not\in\mathcal J_j} c_{i^{*}p},\]

obtaining the most distant unit among the members that minimize the global dissimilarity between one group and all the others (minsumnoint).

MUS algorithm described in Egidi et al. (2018a) is a sequential procedure for extracting identity submatrices of small rank and pivotal units from large and sparse matrices. The procedure has already been satisfactorily applied for solving the label switching problem in Bayesian mixture models (Egidi et al. 2018b).

With the function MUS the user may detect pivotal units from a co-association matrix C, obtained through \(H\) different partitions, whose units may belong to \(k\) groups, expressed by the argument clusters. We remark here that MUS algorithm may be performed only when \(k <5\).

#generate some data

n  <- 620
centers  <- 3
n1 <- 20
n2 <- 100
n3 <- 500
x  <- matrix(NA, n,2)
truegroup <- c( rep(1,n1), rep(2, n2), rep(3, n3))

for (i in 1:n1){
 x[i,]=rmvnorm(1, c(1,5), sigma=diag(2))}
for (i in 1:n2){
 x[n1+i,]=rmvnorm(1, c(4,0), sigma=diag(2))}
for (i in 1:n3){
 x[n1+n2+i,]=rmvnorm(1, c(6,6), sigma=diag(2))}

H <- 1000
a <- matrix(NA, H, n)

  for (h in 1:H){
    a[h,] <- kmeans(x,centers)$cluster

#build the similarity matrix
sim_matr <- matrix(1, n,n)
 for (i in 1:(n-1)){
    for (j in (i+1):n){
      sim_matr[i,j] <- sum(a[,i]==a[,j])/H
      sim_matr[j,i] <- sim_matr[i,j]

cl <- KMeans(x, centers)$cluster
mus_alg <- MUS(C = sim_matr, clusters = cl, prec_par = 5)


Quite often, classical K-means fails in recognizing the true groups:

In such situations, we may need a more robust version of the classical K-means. The pivots may be used as initial seeds for a classical K-means algorithm. The function piv_KMeans works as the classical kmeans function, with some optional arguments

piv_res <- piv_KMeans(x, centers)

The function piv_KMeans has optional arguments:


Egidi, Leonardo, Roberta Pappadà, Francesco Pauli, and Nicola Torelli. 2018a. “Maxima Units Search (Mus) Algorithm: Methodology and Applications.” Studies in Theoretical and Applied Statistics: SIS 2016, Salerno, Italy, June 8-10. Perna, Cira and Pratesi, Monica and Ruiz-Gazen, Anne) 227. Springer: 71–81.

———. 2018b. “Relabelling in Bayesian Mixture Models by Pivotal Units.” Statistics and Computing 28 (4). Springer: 957–69.