k-medoids聚类算法实现
Yanjun发表于 2015-12-11 13:31:09
k-medoids聚类算法，即k-中心聚类算法，它是基于k-means聚类算法的改进。我们知道，k-means算法执行过程，首先需要随机选择初始质心，只有第一次随机选择的初始质心才是实际待聚类点集中的点，而后续将非质心点指派到对应的质心点后，重新计算得到的质心并非是待聚类点集中的点，而且如果某些非质心点是离群点的话，导致重新计算得到的质心可能偏离整个簇，为了解决这个问题，提出了改进的k-medoids聚类算法。
k-medoids聚类算法也是通过划分的方式来计算得到聚类结果，它使用绝对差值和（Sum of Absolute Differences，SAD）来衡量聚类结果的优劣，在n维欧几里德空间中，计算SAD的公式如下所示：
围绕中心点划分（Partitioning Around Medoids，PAM）的方法是比较常用的，使用PAM方法进行处理，可以指定一个最大迭代次数的参数，在迭代过程中基于贪心策略来选择使得聚类的质量最高的划分。PAM的方法，每次交换一个中心点和非质心点，进行指派以后，计算得到的SAD值越小，则聚类质量越好，如此不断地迭代，直到找到一个最好的划分。
在维基百科上，它给出的基于PAM方法计算聚类的过程，描述如下：
- 从待聚类的数据点集中随机选择k个点，作为初始中心点；
- 将待聚类的数据点集中的点，指派到最近的中心点；
- 进入迭代，直到聚类的质量满足指定的阈值（可以通过计算SAD），使总代价减少：
- 对每一个中心点o，对每一个非中心点p，执行如下计算步骤：
- 交换点o和p，重新计算交换后的该划分所生成的代价值；
- 如果本次交换造成代价增加，则取消交换。
上面算法描述，应该是按照某种顺序的取中心点集合中的点，也从非中心点集合中取点，分别计算生成新的划分的代价，我们可以考虑，每次取点的时候，采用随机取点的策略，随机性越强越好，只要满足最终迭代终止的条件即可。通常，如果能够迭代所有情况，那么最终得到的划分一定是最优的划分，即聚类结果最好，这通常适用于比较小的点的集合，但是如果待聚类的点比较大，则需要通过限制迭代次数来终止迭代计算，从而得到一个能够满足实际需要的聚类结果。
我们在下面的实现中，选择随机取中心点和非中心点进行交换，通过设置允许最大迭代次数这个参数值，来使聚类迭代计算停止。
聚类算法实现
首先，便于理解后面的实现过程，我们描述一下代码实现聚类过程的基本逻辑，如下所示：
- 输入待聚类点集，以及参数k、maxIterations、parallism；
- 同k-means算法一样，随机选择初始中心点集合；
- 启动parallism个线程，用来将非中心点指派给最近的中心点；
- 开始执行迭代，使得聚类结果的SAD值最小：
- 将非中心点，基于Round-Robin策略，分配给多个线程，并行指派：将非中心点指派给距离最近的中心点；
- 将多个线程指派的局部结果进行合并，得到一个全局的指派结果；
- 根据指派结果计算SAD值：如果是第一次进行指派，直接计算其SAD值，保存在previousSAD变量中，该变量保存的是最小的SAD值，第一次初始化第一次指派结果计算得到的SAD值；如果不是第一次进行指派，也计算SAD值，将SAD值保存在变量currentSAD中，继续执行步骤8；
- 随机选择一个非中心点；
- 创建一个ClusterHolder对象，该对象保存了该轮迭代指派结果，根据随机选择的非中心点修改ClusterHolder对象中的结果，将随机选择非中心点和对应的中心点进行交换，为下一轮指派过程准备数据；
- 最后，判断是否达到指定的最大迭代次数，如果达到则终止计算，处理最终聚类结果，否则执行下一轮迭代计算。
我们实现的k-medoids聚类算法，需要指定2个聚类相关参数，另外一个参数是程序计算并行度，可以通过构造方法看到，代码如下所示：
public KMedoidsClustering(int k, int maxIterations, int parallism) {
super(k, maxIterations, parallism);
distanceCache = new DistanceCache(Integer.MAX_VALUE);
executorService = Executors.newCachedThreadPool(new NamedThreadFactory("SEEKER"));
latch = new CountDownLatch(parallism);
}
上面代码中的参数含义如下：
- k：聚类最终想要得到的簇的个数
- maxIterations：因为k-medoids聚类算法的最终目标是最小化SAD的值，所以聚类算法执行迭代的次数越大，最终的结果可能越接近最优，如果是对一个不大的点集进行聚类，可以设置该参数的值大一些
- parallism：每一次迭代过程中，我们都需要将非中心点（Non-medoid Point）指派到最近的中心点，所以将原待聚类点集划分成多组，有多个处理线程并行处理可能速度会更快，该参数就是并行度
聚类实现的核心代码如下所示：
@Override
public void clustering() {
// parse sample files
FileUtils.read2DPointsFromFiles(allPoints, "[\t,;\\s]+", inputFiles); // 从文件读取点数据，加入到集合allPoints中
LOG.info("Total points: count=" + allPoints.size());
ClusterHolder currentHolder = new ClusterHolder(); // 每一次迭代过程中的需要的数据结构都封装到ClusterHolde对象中
ClusterHolder previousHolder = null;
currentHolder.medoids = initialCentroidsSelectionPolicy.select(k, allPoints); // 随机选择初始中心
LOG.info("Initial selected medoids: " + currentHolder.medoids);
// start seeker threads
for (int i = 0; i < parallism; i++) { // 启动parallism个线程，执行非中心点到中心点的指派
final NearestMedoidSeeker seeker = new NearestMedoidSeeker(seekerQueueSize);
executorService.execute(seeker);
seekers.add(seeker);
}
// /////////////////
// make iterations
// /////////////////
boolean firstTimeToAssign = true;
int numIterations = 0;
double previousSAD = 0.0;
double currentSAD = 0.0;
try {
while(!finallyCompleted) {
try {
LOG.debug("Current medoid set: " + currentHolder.medoids);
if(firstTimeToAssign) {
// 第一次处理时，只是根据随机选择的初始中心集合，和全部点的集合，指派给多个线程处理
assignNearestMedoids(currentHolder, true);
firstTimeToAssign = false;
} else {
// 非第一次处理时，每次迭代得到的聚类结果，都是基于中心点进行分组的，处理逻辑稍微不同
assignNearestMedoids(currentHolder, false);
}
// merge result
mergeMedoidAssignedResult(currentHolder); // 每个线程处理一部分，最后要合并多个线程分别处理的结果
LOG.debug("Merged result: " + currentHolder.medoidWithNearestPointSet);
// compare cost for 2 iterations, we use SAD (sum of absolute differences)
if(previousSAD == 0.0) {
// first time compute SAD
previousSAD = currentSAD;
currentSAD = computeSAD(currentHolder); // 第一次计算SAD
} else {
RandomPoint randomPoint = selectNonCenterPointRandomly(currentHolder); // 随机选择一个非中心点
LOG.debug("Randomly selected: " + randomPoint);
// compute current cost when using random point to substitute for the medoid
currentSAD = computeSAD(currentHolder); // // 计算用随机选择非中心点替换一个中心点得到的SAD值
// compare SADs
if(currentSAD - previousSAD < 0.0) { // 如果此次迭代得到的SAD值，比上次迭代计算得到SAD小，替换previousHolder和previousSAD，以保证最终算法终止后，该最小值对应的划分能够保留下来
previousHolder = currentHolder;
previousSAD = currentSAD;
}
// construct new cluster holder
currentHolder = constructNewHolder(currentHolder, randomPoint); // 根据随机选择的中心点，创建一个新的 ClusterHolde对象，用于下次迭代
}
LOG.info("Iteration #" + (++numIterations) + ": previousSAD=" + previousSAD + ", currentSAD=" + currentSAD);
if(numIterations > maxIterations) { // 如果达到指定的最大迭代次数，则终止
finallyCompleted = true;
}
} catch(Exception e) {
Throwables.propagate(e);
} finally {
try {
if(!finallyCompleted) {
latch = new CountDownLatch(parallism);
completeToAssignTask = false;
}
Thread.sleep(10);
synchronized(signalLock) {
signalLock.notifyAll();
}
} catch (InterruptedException e) {}
}
}
} finally {
LOG.info("Shutdown executor service: " + executorService);
executorService.shutdown();
}
// finally result
centerPointSet.addAll(previousHolder.medoids); // 处理最终的聚类结果
Iterator<Entry<CenterPoint, List<Point2D>>> iter = previousHolder.medoidWithNearestPointSet.entrySet().iterator();
while(iter.hasNext()) {
Entry<CenterPoint, List<Point2D>> entry = iter.next();
int clusterId = entry.getKey().getId();
Set<ClusterPoint<Point2D>> set = Sets.newHashSet();
for(Point2D p : entry.getValue()) {
set.add(new ClusterPoint2D(p, clusterId));
}
clusteredPoints.put(clusterId, set);
}
}
通过上面代码及其注释，我们可以了解到基本的处理流程，首先看一下几个工具类ClusterHolder和RandomPoint：
private class ClusterHolder {
/** snapshot of clustering result: medoids of clustering result, as well as non-medoid points */
private TreeMap<CenterPoint, List<Point2D>> medoidWithNearestPointSet;
/** center point set represented by Point2D */
private Set<Point2D> centerPoints;
/** center point set represented by CenterPoint */
private TreeSet<CenterPoint> medoids;
public ClusterHolder() {
super();
}
}
private class RandomPoint {
/** medoid which the random point belongs to */
private final CenterPoint medoid; // 随机选择的中心点
/** a non-medoid point selected randomly */
private final Point2D point; // 随机选择的非中心点，该点被指派给上面的中心点medoid
public RandomPoint(CenterPoint medoid, Point2D point) {
super();
this.medoid = medoid;
this.point = point;
}
@Override
public String toString() {
return "RandomPoint[medoid=" + medoid + ", point=" + point + "]";
}
}
下面我们看一下上面代码调用的比较重要的方法的实现逻辑。
- 并行将非中心点指派到最近的中心点
将非中心点指派到最近的中心点的计算，是调用assignNearestMedoids方法，该方法的代码实现，如下所示：
private void assignNearestMedoids(final ClusterHolder holder, boolean firstTimeToAssign) {
LOG.debug("firstTimeToAssign=" + firstTimeToAssign);
try {
// assign tasks to seeker threads
if(firstTimeToAssign) { // 第一次进行指派，因为还没有进行指派过，所以只有随机选择的一组中心点，和全部待聚类的点的集合
holder.centerPoints = Sets.newHashSet();
for(CenterPoint medoid : holder.medoids) {
holder.centerPoints.add(medoid.toPoint()); // 构造ClusterHolder对象，将中心点加入到集合中
}
LOG.debug("holder.centerPoints: " + holder.centerPoints);
for(Point2D p : allPoints) { // 对全部待聚类的点作为任务，加入到每个线程的队列中，但是要排除已经被选择为中心点的点
LOG.debug("Assign point: " + p);
if(!holder.centerPoints.contains(p)) {
selectSeeker().q.put(new Task(holder.medoids, p));
}
}
} else {
for(List<Point2D> points : holder.medoidWithNearestPointSet.values()) { // 如果笔试第一次进行指派，已经在构造ClusterHolder对象的时候，将随机选择的中心点和非中心点进行了交换，这里直接进行指派即可
for(Point2D p : points) {
selectSeeker().q.put(new Task(holder.medoids, p));
}
}
}
} catch(Exception e) {
Throwables.propagate(e);
} finally {
try {
completeToAssignTask = true;
latch.await();
} catch (InterruptedException e) { }
}
}
上面代码调用selectSeeker()方法，获取到一个NearestMedoidSeeker线程，将待指派的点加入到其队列中，selectSeeker()方法实现代码如下所示：
private NearestMedoidSeeker selectSeeker() {
int index = taskIndex++ % parallism;
return seekers.get(index);
}
下面，我们看一下NearestMedoidSeeker线程的实现，它也比较简单，实现了将队列q中的点取出，计算到该点最近的中心，然后指派给该中心，线程实现代码如下所示：
private class NearestMedoidSeeker implements Runnable {
private final Log LOG = LogFactory.getLog(NearestMedoidSeeker.class);
private final BlockingQueue<Task> q;
private Map<CenterPoint, List<Point2D>> clusteringNearestPoints = Maps.newHashMap();
private int processedTasks = 0;
public NearestMedoidSeeker(int qsize) {
q = new LinkedBlockingQueue<Task>(qsize);
}
@Override
public void run() {
while(!finallyCompleted) { // 每一轮迭代，调用一次assign方法
try {
assign();
Thread.sleep(200);
} catch (Exception e) {
e.printStackTrace();
}
}
}
private void assign() throws InterruptedException {
try {
LOG.debug("Q size: " + q.size());
while(!(q.isEmpty() && completeToAssignTask)) {
processedTasks++;
final Task task = q.poll();
if(task != null) {
final Point2D p1 = task.point;
double minDistance = Double.MAX_VALUE;
CenterPoint nearestMedoid = null;
for(CenterPoint medoid : task.medoids) {
final Point2D p2 = medoid.toPoint();
Double distance = distanceCache.computeDistance(p1, p2); // 计算非中心点p1到中心点p2的欧几里德距离
if(distance < minDistance) {
minDistance = distance;
nearestMedoid = medoid;
}
}
LOG.debug("Nearest medoid seeked: point=" + p1 + ", medoid=" + nearestMedoid);
List<Point2D> points = clusteringNearestPoints.get(nearestMedoid);
if(points == null) {
points = Lists.newArrayList();
clusteringNearestPoints.put(nearestMedoid, points);
}
points.add(p1); // 将非中心点p1，指派给到中心点的欧几里德距离最近的点
} else {
Thread.sleep(150);
}
}
} catch (Exception e) {
e.printStackTrace();
} finally {
latch.countDown();
LOG.debug("Point processed: processedTasks=" + processedTasks);
synchronized(signalLock) {
signalLock.wait();
}
clusteringNearestPoints = Maps.newHashMap();
processedTasks = 0;
}
}
}
每一轮指派，多个线程都计算得到一个非中心点指派的子集，最后还要将这些子集合并为一个全局的指派结果，合并的实现在mergeMedoidAssignedResult()方法中，代码如下所示：
private void mergeMedoidAssignedResult(ClusterHolder currentHolder) {
currentHolder.medoidWithNearestPointSet = Maps.newTreeMap();
for(NearestMedoidSeeker seeker : seekers) {
LOG.debug("seeker.clusteringNearestPoints: " + seeker.clusteringNearestPoints);
Iterator<Entry<CenterPoint, List<Point2D>>> iter = seeker.clusteringNearestPoints.entrySet().iterator();
while(iter.hasNext()) {
Entry<CenterPoint, List<Point2D>> entry = iter.next();
List<Point2D> set = currentHolder.medoidWithNearestPointSet.get(entry.getKey());
if(set == null) {
set = Lists.newArrayList();
currentHolder.medoidWithNearestPointSet.put(entry.getKey(), set);
}
set.addAll(entry.getValue());
}
}
}
合并后的指派结果，都存放在ClusterHolder对象中，为下一轮迭代准备了数据。
- 随机选择中心点和非中心点
随机选择一个中心点和非中心点，并交换，实现代码如下所示：
private RandomPoint selectNonCenterPointRandomly(ClusterHolder holder) {
List<CenterPoint> medoids = new ArrayList<CenterPoint>(holder.medoidWithNearestPointSet.keySet());
CenterPoint selectedMedoid = medoids.get(random.nextInt(medoids.size())); // 随机选择一个中心点
List<Point2D> belongingPoints = holder.medoidWithNearestPointSet.get(selectedMedoid);
Point2D point = belongingPoints.get(random.nextInt(belongingPoints.size())); // 随机选择一个非中心点
return new RandomPoint(selectedMedoid, point); // 返回这2个点
}
因为每一次迭代，我们都得到一个非中心点指派到最近的中心点的聚类结果集合，所以在设计随机选择中心点和非中心点进行交换时，我们首先从中心点集合中选择一个中心点，然后再从该中心点对应的非中心点的簇的集合中选择一个非中心点，当然也可以考虑用其他的方法，比如，待交换的中心点和非中心点不在同一个簇中。
- 创建ClusterHolder对象，交换非中心点和中心点
我们处理的策略是，事后处理，也就是每次先实现非中心点和中心点的交换，再进行指派，计算SAD值，如果此轮迭代得到的SAD值比上一轮的大，则直接丢弃结果，将上一轮的指派结果作为最终候选结果，直到最后，保留着具有最小SAD值的指派结果。创建ClusterHolder对象，交换非中心点和中心点的实现逻辑，如下所示：
private ClusterHolder constructNewHolder(final ClusterHolder holder, RandomPoint randomPoint) {
ClusterHolder newHolder = new ClusterHolder();
// collect center points with type Point2D for a holder object
// from previous result of clustering procedure
newHolder.centerPoints = Sets.newHashSet();
for(CenterPoint c : holder.medoidWithNearestPointSet.keySet()) {
newHolder.centerPoints.add(c.toPoint());
}
Point2D newPoint = randomPoint.point;
CenterPoint oldMedoid = randomPoint.medoid;
// create a new center point with type CenterPoint based on the randomly selected non-medoid point
// and it's id is equal to the old medoid's
CenterPoint newMedoid = new CenterPoint(oldMedoid.getId(), newPoint);
// use new medoid above to substitute the old medoid
newHolder.centerPoints.remove(oldMedoid.toPoint());
newHolder.centerPoints.add(newPoint);
newHolder.medoids = Sets.newTreeSet();
newHolder.medoids.addAll(holder.medoidWithNearestPointSet.keySet());
newHolder.medoids.remove(oldMedoid); // remove old medoid from center point set of new holder object
newHolder.medoids.add(newMedoid);
// copy the holder's medoidWithNearestPointSet, and modify it
newHolder.medoidWithNearestPointSet = Maps.newTreeMap();
newHolder.medoidWithNearestPointSet.putAll(holder.medoidWithNearestPointSet);
List<Point2D> oldPoints = newHolder.medoidWithNearestPointSet.get(oldMedoid);
oldPoints.remove(newPoint); // remove new randomly selected non-medoid point from previous result set of clustering
oldPoints.add(oldMedoid.toPoint()); // add old medoid point to the non-medoid set
newHolder.medoidWithNearestPointSet.put(newMedoid, oldPoints);
return newHolder;
}
为了保留上一次迭代指派的结果，这里不要修改holder对应的结果的集合，而是拷贝出一份，在拷贝的结果上交换中心点和非中心点。
聚类效果对比
我们将使用k-medoids算法与k-means算法聚类，分别对结果进行比较。其中k-means算法由于随机选择初始质心，每次执行聚类结果不同，而k-medoids算法聚类依赖于与迭代次数，所以我们选择maxIterations分别为300、1000、3000时，对比效果，如下图所示：
通过上图可以看出，k-medoids聚类算法，当maxIterations越大的时候，可能更加靠近最优解，聚类结果的质量越高，也就是SAD值越小。
参考链接