Kmeans聚类算法

Kmeans是最流行的，以及最简单的用于挖掘数据潜在结构的机器学习算法之一。Kmeans的目标很简单：根据数据的均值，将数据划分为若干个簇。假定每个簇的均值可以很好地代表簇内的每一个观察值。

Kmeans算法

假定我们想把数据划分为k个簇，那么我们就需要找出这k个簇的k个中心。该如何定义，以及寻找这些中心呢？

我们只需要求解方程：$min\sum_i^{N}\sum_j^KO_i^j||x_i-u_j||^2$，其中当观察点i为簇j的中心时$O_i^j=1$，否则为0.

我们正在寻找k个中心，使得各个簇内的点到簇中心的距离最小。这是一个最优化问题，但是上面的目标函数是非凸的，而且有一些变量是二元的，无法用传统的梯度下降法求解。

解决这个问题的方法如下：

1. 随机初始化k个中心。
2. 更新每个中心。新的中心为相应簇中的所有观察值的平均值。
3. 更新收敛准则。

用R语言实现K-means算法

既然，我们已经有了伪代码算法，接下来我们可以用R语言实现Kmeans算法。首先，我们创建5类数据，它们都服从2维高斯分布。

require(MASS)
require(ggplot2)
set.seed(1234)
set1 <- mvrnorm(n = 300, mu = c(-4, 10), Sigma = matrix(c(1.5, 1, 1, 1.5), 2))
set2 <- mvrnorm(n = 300, mu = c(5, 7), Sigma = matrix(c(1, 2, 2, 6), 2))
set3 <- mvrnorm(n = 300, mu = c(-1, 1), Sigma = matrix(c(4, 0, 0, 4), 2))
set4 <- mvrnorm(n = 300, mu = c(10, -10), Sigma = matrix(c(4, 0, 0, 4), 2))
set5 <- mvrnorm(n = 300, mu = c(3, -3), Sigma = matrix(c(4, 0, 0, 4), 2))
DF <- data.frame(rbind(set1, set2, set3, set4, set5), cluster=as.factor(rep(1:5, each = 300)))

ggplot(DF, aes(x=X1, y=X2, color=cluster)) + geom_point()

在这个数据集中，Kmeans算法将会得到一个很好的结果，因为每个分布都呈现圆形，如上图所示。

初始化簇中心

簇中心的初始化非常重要，它能够影响算法。因此，我们从数据集中随机选取K的点。

# 中心初始化
centroids <- data[sample.int(nrow(data), K), ]
# 停止准则初始化
current_stop_crit <- 10e10
# 初始化每个点的指定中心
cluster <- rep(0, nrow(data))
# 算法是否收敛
converged <- FALSE
iter <- 1

将每个点归为指定的簇

在每一次迭代中，每个点将会归为离它最近的簇。我们可以使用欧几里得距离来计算每个点到每个中心的距离，并保存各个点到各个中心的最近距离以及相应的簇中心。

# 在观察值上进行迭代
for (i in 1:nrow(data)){
  # 设置最小距离
  min_dist <- 10e10
  # 在中心上进行迭代
  for (centroid in 1:nrow(centroids)){
    # 计算欧式距离
    distance_to_centroid <- sum((centroids[centroid, ] - data[i, ]) ^ 2)
    # 距离点最近的中心
    if (distance_to_centroid <= min_dist){
      # 这个点将会被归为这个簇
      cluster[i] <- centroid
      min_dist <- distance_to_centroid
    }
  }
}

簇中心更新

一旦每个观察值都被归为距离最近的簇，每个簇的中心坐标就会被更新。簇中心的新坐标为相应簇中所有点的平均值。收敛准则为：当各个簇的中心停止变动时，算法就应该终止。

while (current_stop_crit >= stop_crit & converged == FALSE){
  iter <- iter + 1
  if (current_stop_crit <= stop_crit){
    converged <- TRUE
  }
  old_centroids <- centroids
  # 更新停止准则
  current_stop_crit <- mean((old_centroids - centroids) ^ 2))
}

完整Kmeans函数

kmeans <- function(data, K=4, stop_crit=10e-5){
  # 初始化簇中心
  centroids <- data[sample.int(nrow(data), K), ]
  current_stop_crit <- 1000
  cluster <- rep(0, nrow(data))
  converged <- FALSE
  iter <- 1
  while (current_stop_crit >= stop_crit & converged == FALSE){
    iter <- iter + 1
    if (current_stop_crit <= stop_crit){
      converged <- TRUE
    }
    old_centroids <- centroids
    # 把每个点归入相应的簇
    for (i in 1:nrow(data)){
      min_dist <- 10e10
      for (centroid in 1:nrow(centroids)){
        distance_to_centroid <- sum((centroids[centroid, ] - data[i, ]) ^ 2)
        if (distance_to_centroid <= min_dist){
          cluster[i] <- centroid
          min_dist <- distance_to_centroid
        }
      }
    }
    # 更新簇中心
    for (i in 1:nrow(centroids)){
      centroids[i, ] <- apply(data[cluster == i, ], 2, mean)
    }
    current_stop_crit <- mean((old_centroids - centroids) ^ 2)
  }
  return(list(data=data.frame(data, cluster), centroids=centroids))
}

我们可以使用上述代码来观察我们的数据：

res <- kmeans(DF[1:2], K=5)
res$centroids$cluster <- 1:5
res$data$isCentroid <- FALSE
res$centroids$isCentroid <- TRUE
data_plot <- rbind(res$centroids, res$data)
ggplot(data_plot,aes(x=X1,y=X2,color=as.factor(cluster),size=isCentroid,alpha=isCentroid)) +
  geom_point()