Llyod and Forgy algorithms

Several algorithms optimize the k-means criterion. For instance, the R function kmeans() provides four algorithms:

kmeans(x, centers, iter.max = 10, nstart = 1,
       algorithm = c("Hartigan-Wong", "Lloyd", "Forgy",
                     "MacQueen"), trace=FALSE)

In fact, Forgy and Lloyd algorithms are implemented the same way. We can see this in the source code of kmeans():

edit(kmeans)

It opens the source code in your favorite text editor. At lines 56 and 57, "Forgy" and "Lloyd" are assigned to the same number (2L) and are thus mapped to the same implementation:

nmeth <- switch(match.arg(algorithm), `Hartigan-Wong` = 1L,
    Lloyd = 2L, Forgy = 2L, MacQueen = 3L)

In the following, we will implement this algorithm. After the initilization, it iterates over two steps until convergence:

• an assignment step which assigns the points to the clusters;
• an update step which updates the centroids of the clusters.

Initialization

The initialization consists in selecting \(K\) points at random and consider them as the centroids of the clusters:

dataset = dataset %>% mutate(sample = row_number())
centroids = dataset %>% sample_n(2) %>% mutate(cluster = row_number()) %>% select(x, y, cluster)

Assignment step

The assignment step of k-means is equivalent to the E and C step of the CEM algorithm in the Gaussian mixture model. It assigns the points to the clusters according to the distances between the points and the centroids. Let's write \(z_k\) the set of points in the cluster \(k\):

We estimate \(z_k\) by:

A point \(x_i\) is set to be in the cluster \(k\) if the closest centroid (using the euclidean distance) is the centroid \(\mu_k\) of the cluster \(k\). This is done by the following R code:

assignmentStep = function(samplesDf, centroids) {
  d = samplesDf %>% select(x, y, sample)
  repCentroids = bind_rows(replicate(nrow(d), centroids, simplify = FALSE)) %>%
    transmute(xCentroid = x, yCentroid = y, cluster)
  d %>% slice(rep(1:n(), each=2)) %>%
    bind_cols(repCentroids) %>%
    mutate(s = (x-xCentroid)^2 + (y-yCentroid)^2) %>%
    group_by(sample) %>%
    top_n(1, -s) %>%
    select(cluster, x, y)
}

Update step

In the update step, the centroid of a cluster is computed by taking the mean of the points in the cluster, as defined in the previous step. It corresponds to the M step of the Gaussian mixture model and it is done in R with:

updateStep = function(samplesDf) {
  samplesDf %>% group_by(cluster) %>%
    summarise(x = mean(x), y = mean(y))
}

Iterations

Let's put together the steps defined above in a loop to complete the algorithm. We define a maximum number of iterations maxIter and iterate over the two steps until either convergence or maxIter is reached. It converges if the centroids are the same in two consecutive iterations:

maxIter = 10
d = data.frame(sample=c(), cluster=c(), x=c(), y=c(), step=c())
dCentroids = data.frame(cluster=c(), x=c(), y=c(), step=c())
for (i in 1:maxIter) {
  df = assignmentStep(dataset, centroids)
  updatedCentroids = updateStep(df)
  if (all(updatedCentroids == centroids )) {
    break
  }
  centroids = updatedCentroids
  d = bind_rows(d, df %>% mutate(step=i))
  dCentroids = bind_rows(dCentroids, centroids %>% mutate(step=i))
}

The above R code constructs two dataframes d and dCentroids which contain respectively the assignations of the points and the centroids. The column step indicates the iteration number and will be used to animate the plot.