k-means（k-均值）算法是一种基于距离的聚类算法，它用质心（Centroid）到属于该质心的点距离这个度量来实现聚类，通常可以用于N维空间中对象。下面，我们以二维空间为例，概要地总结一下k-means聚类算法的一些要点：
- 除了随机选择的初始质心，后续迭代质心是根据给定的待聚类的集合S中点计算均值得到的，所以质心一般不是S中的点，但是标识的是一簇点的中心。
- 基本k-means算法，开始需要随机选择指定的k个质心，因为初始k个质心是随机选择的，所以每次执行k-means聚类的结果可能都不相同。如果初始随机选择的质心位置不好，可能造成k-means聚类的结果非常不理想。
- 计算质心：假设k-means聚类过程中，得到某一个簇的集合Ci={p(x1,y1), p(x2,y2), …,p(xn,yn)}，则簇Ci的质心，质心x坐标为(x1+x2+ …+xn)/n，质心y坐标为(y1+y2+ …+yn)/n。
- k-means算法的终止条件：质心在每一轮迭代中会发生变化，然后需要重新将非质心点指派给最近的质心而形成新的簇，如果只有很少的一部分点在迭代过程中，还在改变簇（如，更新一次质心，有些点从一个簇移动到另一个簇），那么满足这样一个收敛条件，可以提前结束迭代过程。
- k-means算法的框架是：首先随机选择k个初始质心点，然后执行聚类处理迭代，不断更新质心，直到满足算法收敛条件。由于该算法收敛于局部最优，所以多次执行聚类算法，通过比较，选择聚类效果最好的结果作为最终的结果。
- k-means算法聚类完成后，没有离群点，所有的点都会被指派到对应的簇中。
由于k-means算法比较简单，对于算法的实现过程，我们概要地描述如下：
- 随机选择k个初始质心；
- 如果没有满足聚类算法终止条件，则继续执行步骤3，否则转步骤5；
- 计算每个非质心点p到k个质心的欧几里德距离，将p指派给距离最近的质心；
- 根据上一步的k个质心及其对应的非质心点集，重新计算新的质心点，然后转步骤2；
- 输出聚类结果，算法可以执行多次，使用散点图比较不同的聚类结果。
下面，我们详细说明上述步骤：
随机选择初始质心
由于随机选择初始质心，每次执行聚类选择的初始质心都不相同，这也导致k-means算法聚类后，没有确定的结果，或者说，可能两次聚类的结果完全不同。该过程的实现，比较简单，只要随机选择给定待聚类点集合中的点即可，初始质心是实际存在的点，代码如下所示：
@Override
public TreeSet<Centroid> select(int k, List<Point2D> points) {
TreeSet<Centroid> centroids = Sets.newTreeSet();
Set<Point2D> selectedPoints = Sets.newHashSet();
while(selectedPoints.size() < k) { // 先随机选择k个点
int index = random.nextInt(points.size());
Point2D p = points.get(index);
selectedPoints.add(p);
}
Iterator<Point2D> iter = selectedPoints.iterator();
int id = 0;
while(iter.hasNext()) { // 构造Centroid质心对象，分配一个id作为簇的唯一标识
centroids.add(new Centroid(id++, iter.next()));
}
return centroids;
}
有一些方法，可以在这一步中，解决初始质心选择的随机性，可以将选择初始质心作为选择策略的设计，根据需要选择不同的策略，比如，可以这样设计策略接口：
public interface SelectInitialCentroidsPolicy {
TreeSet<Centroid> select(int k, List<Point2D> points);
}
我们这里只给了简单地随机选择策略，也是基本k-means算法最基础的策略。其他方法，可以查阅相关资料。
计算欧几里德距离，指派点到质心所在簇
计算每个非质心点到全部k个质心点的距离，将该非质心点指派给距离最小的质心点所在的簇。如果输出的数据量比较大，可以将数据集合进行分割，基于多线程去并行处理，最后再合并结果。我们的实现思路是：每个线程都共享k个质心的集合，然后将非质心点均匀分发到多个线程的队列中，然后每个线程从队列取出非质心点，计算非质心点到k个质心的距离，并计算出距离最短的质心，将该非质心点指派给该质心所在的簇。实现代码如下所示：
Task task = q.poll();
Point2D p1 = task.point;
// assign points to a nearest centroid
Distance minDistance = null;
for(Centroid centroid : task.centroids) { // 计算一个非质心点，到k个质心的距离，并计算距离最短的
double distance = MetricUtils.euclideanDistance(p1, centroid);
if(minDistance != null) {
if(distance < minDistance.distance) {
minDistance = new Distance(p1, centroid, distance);
}
} else {
minDistance = new Distance(p1, centroid, distance);
}
}
LOG.debug("Assign Point2D[" + p1 + "] to Centroid[" + minDistance.centroid + "]");
Multiset<Point2D> pointsBelongingToCentroid = localClusteredPoints.get(minDistance.centroid);
if(pointsBelongingToCentroid == null) {
pointsBelongingToCentroid = HashMultiset.create();
localClusteredPoints.put(minDistance.centroid, pointsBelongingToCentroid); // localClusteredPoints是局部的，key为质心，value为属于该质心的非质心点的集合
}
pointsBelongingToCentroid.add(p1);
这样，经过一轮的迭代计算，每个线程都处理完，得到一个局部的指派的簇的集合，然后对每个局部集合进行合并，得到一个全局的、质心到属于该质心的点的簇的集合，作为下一次迭代的输入，也比较容易处理。
迭代终止条件计算
这一步应该算是k-means算法聚类过程中比较核心的步骤。我们考虑了如下3个终止条件：
- 比较相邻的2轮迭代结果，在2轮过程中移动的非质心点的个数，设置移动非质心点占比全部点数的最小比例值，如果达到则算法终止
- 为了防止k-means聚类过程长时间不收敛，设置最大迭代次数，如果达到最大迭代次数还没有达到上述条件，则也终止计算
- 如果相邻2次迭代过程，质心没有发生变化，则算法终止，这是最强的终止约束条件。能够满足这种条件，几乎是不可能的，除非两次迭代过程中没有非质心点重新指派给到另一个不同的质心。
我们计算k-means聚类的核心代码框架，如下所示：
@Override
public void clustering() {
// start centroid calculators
for (int i = 0; i < parallism; i++) { // 启动parallism个线程计算距离并指派簇
CentroidCalculator calculator = new CentroidCalculator(calculatorQueueSize);
calculators.add(calculator);
executorService.execute(calculator);
LOG.info("Centroid calculator started: " + calculator);
}
// sort by centroid id ASC
TreeSet<Centroid> centroids = selectInitialCentroidsPolicy.select(k, allPoints); // 随机选择初始质心
LOG.info("Initial selected centroids: " + centroids);
// 下面进入迭代过程
int iterations = 0;
boolean stopped = false;
CentroidSetWithClusteringPoints lastClusteringResult = null; // 上一轮聚类结果
CentroidSetWithClusteringPoints currentClusteringResult = null; // 当前轮聚类结果
int totalPointCount = allPoints.size();
float currentClusterMovingPointRate = 1.0f;
try {
// enter clustering iteration procedure
while(currentClusterMovingPointRate > maxMovingPointRate
&& !stopped
&& iterations < maxIterations) { // 3个终止条件约束
LOG.info("Start iterate: #" + (++iterations));
currentClusteringResult = computeCentroids(centroids); // 每一轮重新计算质心点
LOG.info("Re-computed centroids: " + centroids);
// compute centroid convergence status
int numMovingPoints = 0;
if(lastClusteringResult == null) {
numMovingPoints = totalPointCount;
} else {
// compare 2 iterations' result for centroid computation
numMovingPoints = analyzeMovingPoints(lastClusteringResult.clusteringPoints, currentClusteringResult.clusteringPoints); // 分析两轮聚类结果：在簇之间移动的非质心点的集合
// check iteration stop condition
boolean isIdentical = (currentClusteringResult.centroids.size() ==
Multisets.intersection(HashMultiset.create(lastClusteringResult.centroids), HashMultiset.create(currentClusteringResult.centroids)).size()); // 检测终止最强约束条件：两轮迭代是否没有非质心点发生重新指派，即质心完全没变
if(iterations > 1 && isIdentical) {
stopped = true;
}
}
lastClusteringResult = currentClusteringResult;
centroids = currentClusteringResult.centroids;
currentClusterMovingPointRate = (float) numMovingPoints / totalPointCount; // 计算非质心点移动比例
LOG.info("Clustering meta: k=" + k +
", numMovingPoints=" + numMovingPoints +
", totalPointCount=" + totalPointCount +
", stopped=" + stopped +
", currentClusterMovingPointRate=" + currentClusterMovingPointRate );
// reset some structures
reset();
for(CentroidCalculator calculator : calculators) {
calculator.reset();
}
LOG.info("Finish iterate: #" + iterations);
}
} finally {
// notify all calculators to exit normally
clusteringCompletedFinally = true;
LOG.info("Shutdown executor service: " + executorService);
executorService.shutdown();
// process final clustering result
LOG.info("Final clustering result: ");
Iterator<Entry<Centroid, Multiset<Point2D>>> iter = currentClusteringResult.clusteringPoints.entrySet().iterator();
while(iter.hasNext()) { // 达到终止条件后，处理最终的结果
Entry<Centroid, Multiset<Point2D>> entry = iter.next();
int id = entry.getKey().getId();
Set<ClusterPoint<Point2D>> set = Sets.newHashSet();
for(Point2D p : entry.getValue()) {
set.add(new ClusterPoint2D(p, id));
}
clusteredPoints.put(id, set);
id++;
}
centroidSet = currentClusteringResult.clusteringPoints.keySet();
}
}
下面，我们讨论一下，如何根据两次聚类迭代结果，计算在簇之间移动的点的个数。如果把两轮聚类迭代结果中的k个簇分别从整体上来比较，得出在前后两轮迭代结果中在簇之间移动的非质心点的个数，可能比较麻烦，也容易陷入混乱的计算逻辑中。
我们可以这么思考：假设a、b两轮迭代结束，a轮中生成k个簇的集合Ca={C(a1),C(a2), …,C(ak)}，b轮中生成k个簇的集合Cb={C(b1),C(b2), …,C(bk)}，我们假设生成的簇是有编号的，而且，a轮生成的簇C(ai)，在b轮重新计算质心后生成的新簇为C(bi)，这样一一对应起来，分别计算在簇C(ai)与簇C(bi)之间移动的点的个数，首先计算簇C(ai)与簇C(bi)的交集S：
S = C(ai) ∩ C(bi)
然后，分别计算簇C(ai)、簇C(bi)与S的差集Dai、Dbi：
Dai = Ca - S = Ca - (C(ai) ∩ C(bi) ) Dbi = Cb - S = Cb - (C(ai) ∩ C(bi) )
这样，差集Dai和Dbi中的点都是在两轮聚类中移动的非质心点，由于一个簇中的点可能移动到另一个簇中，如某非质心点p，从C(ai)移动到C(bj)，其中i不等于j，那么在计算差集Dai与Dbi时，发现C(ai)中少了点p，点p被放入差集Dai；在计算簇C(aj)与簇C(bj)时，发现C(bj)中多了一个点p，则点p又被放入差集Dbj。可见，点p被放入到两个差集Dai和Dbj中，所以我们需要对最终得到的k个差集先做并计算：
D = Σ(Dai ∪ dbi), i=1,2, ...k
然后再对集合D做一个去重操作，得到的点的集合就是两轮迭代过程中，在簇之间移动的点的集合。
我们基于上述计算思路实现的代码，对应上面代码中的analyzeMovingPoints方法，代码实现如下所示：
private int analyzeMovingPoints(TreeMap<Centroid, Multiset<Point2D>> lastClusteringPoints,
TreeMap<Centroid, Multiset<Point2D>> currentClusteringPoints) {
// Map<current, Map<last, intersected point count>>
Set<Point2D> movingPoints = Sets.newHashSet(); // 用来收集移动的点，使用Set集合类去重
Iterator<Entry<Centroid, Multiset<Point2D>>> lastIter = lastClusteringPoints.entrySet().iterator();
Iterator<Entry<Centroid, Multiset<Point2D>>> currentIter = currentClusteringPoints.entrySet().iterator();
while(lastIter.hasNext() && currentIter.hasNext()) {
Entry<Centroid, Multiset<Point2D>> last = lastIter.next();
Entry<Centroid, Multiset<Point2D>> current = currentIter.next();
Multiset<Point2D> intersection = Multisets.intersection(last.getValue(), current.getValue()); // 计算交集S = C(ai) ∩ C(bi)
movingPoints.addAll(Multisets.difference(last.getValue(), intersection)); // 计算差集Dai = Ca - S = Ca - (C(ai) ∩ C(bi) )
movingPoints.addAll(Multisets.difference(current.getValue(), intersection)); // 计算差集Dbi = Cb - S = Cb - (C(ai) ∩ C(bi) )
}
return movingPoints.size();
}
通过上面的计算逻辑，就能够计算出两轮聚类过程中，在簇之间移动的点的集合和个数。
聚类效果
每次执行k-means聚类，得到的结果都不相同，我们可以执行两次，取k=10，看一下聚类结果的散点图，如下图所示：
图中，标号为9999的点为质心点，上面两图对比可以看出，聚类结果中簇的形状是不同的，其中红色值满足迭代停止条件的质心的坐标位置。
下面，我们选择不同的k值：5、10、20、50，分别执行k-means聚类，然后对比聚类结果，如下图所示：
总结
通过上面的实现，我们知道基本k-means聚类算法的实现过程比较简单，很容易实现。另外，该聚类算法适用于处理具有中心的球形簇，而且运行相当有效。但是，该聚类算法的结果受随机选择的质心的影响，每次计算都得到不同的结果，而且当待聚数据的具有不同的尺寸，或者密度非常不均匀，聚类结果非常的差。为了解决k-means聚类随机算法选择初始质心的问题，会有很多处理方法，可以查阅相关资料，其中bisecting k-means算法（二分k-均值）就是基于基本k-means得到的一种变体，能够比较好地处理，不受随机选择初始质心的影响，后续我们会实现并详细讨论。
本文基于署名-非商业性使用-相同方式共享 4.0许可协议发布，欢迎转载、使用、重新发布，但务必保留文章署名时延军（包含链接：http://shiyanjun.cn），不得用于商业目的，基于本文修改后的作品务必以相同的许可发布。如有任何疑问，请与我联系。