谱聚类、Chameleon聚类、PCCA、SOM、Affinity Propagation

层次聚类和K-means聚类，可以说是聚类算法里面最基本的两种方法（wiki的cluster analysis页面都把它们排前两位）。这次要探讨的，则是两个相对“高级”一点的方法：谱聚类和chameleon聚类。

1、谱聚类

一般说到谱聚类，都是从降维(Dimensionality Reduction)或者是图分割(Graph Cut)的角度来理解。但是实际上，从物理学的简正模式的角度，可以更为直观地理解这个算法的本质。

这里先把基本的算法步骤写出来，然后再讨论算法的原理。

谱聚类基本步骤：

1、给出N个数据点两两之间的相似性。也就是一个N*N的相似性矩阵A，A(i,j)代表i和j两个数据点的相似度，数值越大则表示越相似。注意A(i,j)=A(j,i)，A(i,i)=0。

2、计算矩阵D，使它的对角元是A矩阵的对应的那一列（或行）的值之和，其余地方为0。也就是使得

3、令B=D-A

4、求B矩阵的前k个本征值和本征矢，将数据点投影到一个k维空间。第i本征矢的第j个值，就表示第j个数据点在k维空间中第i维的投影。就是说如果把k个特征矢量并成一个N*k的矩阵，则每一行代表一个数据点在k维空间的坐标。

5、根据每个数据点的k维空间坐标，使用K-means或者其它聚类算法在k维空间对数据进行聚类。

从算法的第4、5步就可以看出，谱聚类的本质实际上就类似于PCA，先将数据点投影到一个更能反映数据特征的空间，然后再用其它办法进行聚类。这也就是一种降维的思想（实际上也可能是升维）。那么问题的关键就在于，它把数据点投影到什么空间去了？为什么这个空间更能反映数据特征？这个问题可以从图分割的角度来理解，不过我这里要从简谐振动的角度来讨论这个问题，这也是一个更为直观的理解。

Python的几行代码：

[python] view plaincopy

1. #获取聚类中心
2. def spectralProject(M):
3. #计算矩阵D，使它的对角元是A矩阵的对应的那一列（或行）的值之和，其余地方为0。
4. #也就是使得 D[i,i] = A[i,j]
5. D =np.zeros((Len(M),Len(M) ),dtype=float32 );
6. for i in range(Len(M) ):
7. for j in range(Len(M) ):
8. D[i,i] += M[i,j];
9. #令B =D-A
10. B = D - M;#M[1,2]
11. #求B矩阵的前k个本征值和本征矢，将数据点投影到一个k维空间。
12. #第i本征矢的第j个值，就表示第j个数据点在k维空间中第i维的投影。
13. eVals, eVecs = np.linalg.eig(A);
14. #根据每个数据点的k维空间坐标，使用K-means或者其它聚类算法在k维空间对数据进行聚类。
15. ClusterCenters= eVecs;
16. return ClusterCenters;
18. #谱聚类的本质实际上就类似于PCA，先将数据点投影到一个更能反映数据特征的空间，
19. #然后再用其它办法进行聚类。这也就是一种降维的思想（实际上也可能是升维）。
20. def spectralClustering(M ,dataSet, k):
21. ClusterCenters= spectralProject(M);
22. #根据每个数据点的k维空间坐标，使用K-means或者其它聚类算法在k维空间对数据进行聚类。
23. centroids, clusterAssment = km.kMeans(dataSet, k, distMeas=distEclud, createCent=ClusterCenters);
25. return centroids, clusterAssment;

简正模式

说起简谐振动，学过高中物理的童鞋都不会陌生：两个小球连上一根弹簧，就是最简单的简谐振动模型。为了简单起见，写成一维的形式，而且弹簧的平衡距离设为0，于是，当小球的坐标给定时，弹性势能就是

我们把上面那个算法套用在这个例子上试试，两个小球的“相似度”就看成是它们之间弹簧的弹性系数k，k越大，小球之间的关系自然就越紧密了。这样上面要求的矩阵就是

给出整个系统的坐标矢量

，容易证明

B只有两个本征矢量，分别是

这两个本征向量就代表了体系的两个简正运动模式，向量中的值表示对应的小球在这个运动模式中的运动方向。比如p1之中两个小球往同一个方向运动，这其实是系统的整体平动，对应的本征值为0；p2则表示两个小球往相反方向运动，这就符合我们想象中这两个小球的简谐振动了。

究竟什么是简正运动模式？为什么用上面的方法就能得到系统的简正运动模式呢？其实所谓的简正模式，就类似于傅立叶分析里面，用来将原函数展开的那组相互正交的基函数组。这里所使用的，就是简谐振动这样的一种基组，将整个系统的复杂运动表示为这些简正模式的叠加。

无论我们有多少个小球，只要小球之间是以弹簧相连的，那么根据它们之间的连接方式，总是可以将系统的势能表示为

但是，我们希望的是将运动方式去耦合，写成多个简正模式之和，也就是

所以需要对原来的B矩阵对角化，而对角化过程中得到的本征矢量和本征值，也就是所要求的简正模式以及它们的频率的平方值。

上面说了那么多关于简正模的东西，可是到底为什么要求简正模呢？这是因为谱聚类的目的是要找到一个能很好地反映数据点特征的空间，然后在新空间中进行聚类。试着想象一下，如果两个数据之间相似性很大，那么也就是说它们之间的“弹簧”弹性系数很大，就跟用一个棒子连起来一样，那么自然在运动的时候，它们就会倾向于往同一个方向运动。类似地，如果一堆数据点之间很相似，那么它们就会形成一个rigid的整体，就像一个刚体一样，刚体内部的小球喜欢一起动。而两个刚体之间，则会产生简谐运动，倾向于往不同的方向运动。

用一个简单的例子来说明这个现象，我们可以想象6个小球，分成对称的两组123和456，组内小球两两之间连在一起，两组之间则在3和4间有一根弹簧相连。这样一个结构，很明显应该是分为123和456两组。

本征值最小的前4个本征矢量：第一个一样是整体平动，没有什么意义；第二个表示两组小球之间的相对运动，两组小球往不同方向运动，这是我们想要得到的结果；第三、四个表示每组小球内部的相对运动。

在这个结构里，组内部的相对运动相比起组间运动是很弱的，这可以从它们对应的本征值看出来。根据能均分定理，能量应该在每个简正模式之间均分，所以模式的振幅反比于它们的频率，也就是跟本征值的开方的倒数成正比，这里A2:A3:A4~3:1:1。这其实也就告诉了我们，对传统的谱聚类算法可以根据它的物理意义进行改进，根据本征值对本征矢量进行加权，而不是同一对待。这样，不重要的模式即便被考虑进来，因为振幅很小，所以也不会对结果产生什么影响。这样，我们在算法的第4步，考虑k个本征矢量来进行投影时，就不用担心会多取了多余的本征矢了，而且也可以根据本征值谱的变化来判断k的合理取值，就像在层次聚类中那样。

对PCA比较熟悉的童鞋，会发现这个方法在形式上跟主元分析有类似的地方。其实简正模分析和PCA是两个相反的思路，简正模是根据系统的性质来推断系统的特征运动模式，而PCA则是根据系统的运动结果来得到特征方向。一个是从原理来进行推断，一个则是从结果来进行反推。

聚类结果

使用谱聚类对样品1进行聚类，可以得到下图。两个结果分别对应聚类数目k取为3和8的情况，可以看到并不会像K-means那样把大的cluster分离，只会把少量异常点分离出去，总体的聚类结果十分稳定。因为算法最后还是使用了K-means进行聚类，所以我们可以想象谱聚类在投影到新空间的时候，应该是很好地把不同的cluster远远地分离了开来。

可惜，谱聚类对特殊形状的cluster的聚类效果依然不尽如人意。不过相比起K-means这样的算法，谱聚类已经辨认出一些形状信息了（有成环状的cluster，而不是都是球型的）。

谱聚类聚类的效果比较好，性能也比较稳定。算法需要的输入只是相似矩阵，不需要数据点的坐标矩阵，适用性也较广。一个潜在的问题是，如果数据量很大的话，对大矩阵的对角化可能会导致算法效率低下。但是如果是稀疏矩阵的情况，只计算前k个本征矢量和本征值的效率还是很高的。所以谱聚类算法总体而言是一个不错的选择。

2、Chameleon聚类

之前介绍的三个算法都没办法分辨出样品2的甜甜圈，而这次介绍的chameleon算法则可以说是专门用来干这活的。下面是算法的提出者在他们的文献中给出的一些测试组，可以看到这个算法就是对各类奇葩形状都应对自如……

我们不妨来想当然地考虑一下，怎样才能识别出甜甜圈结构的cluster。简单起见，考虑最极端的情况，假设数据的噪声不是像样品2那么大，而是分界很明显的三个环带。想象我们把一只甲虫放在了其中一个环带上，甲虫的视野很小，而且它会随机地走到它能看到的数据点上。如果环带之间的间距足够大，那么甲虫就不会走到其它环带上。最后，甲虫能走到的区域就是一个环形的cluster。

上面这个问题的关键就在于，要主动地把甲虫的视野变小，也就是根据近邻数据来进行聚类，然后不断延伸。这其实也就类似于层次聚类中的single-linkage，实际上single-linkage也确实可以识别出样品2。

但是这样会带来新的问题，比如对于下面的情况。在fig 4中，如果只看最近邻的连接，算法会倾向于合并c、d而不是a、b，又或者说，如果甲虫的视野足够大到会合并a和b，那么c和d也就一定会被看作一个cluster。但事实上，a、b的邻接区域较大，距离也不远（相对于a、b内部），所以是应该被认为是属于同一个cluster，而c、d则显然不应该被看作一类。fig 5则表示另外一种情况，就是过分强调邻接区域的大小，这样就会倾向于把a与c进行合并，而不是与b合并。

Chameleon算法就是努力在这两种情况之间保持平衡，既考虑Closeness，即近邻节点的靠近程度，也考虑Inter-Connectivity，即邻接区域的大小。

Chameleon本质上也是一个从下而上的层次聚类算法，不过它只考虑每个节点邻近的K个节点（K由用户给定），也就是说，只有最接近节点的其它K个节点会被认为与节点存在连接。两个节点越“接近”，连接强度越强。两个cluster的“距离”由两个参量决定：relative inter-connectivity和relative closeness。

Relative Inter-Connectivity

cluster i和j的relative inter-connectivity表征了他们之间邻接区域的大小，其中EC{Ci,Cj}表示跨越i和j的所有连接的强度之和，EC{Ci}表示将cluster i分割为大小接近的两部分所需要切断的连接的强度值和。

Relative Closeness

cluster i和j的relative closeness表征了他们之间近邻连接的强度，其中SEC{Ci,Cj}表示跨越i和j的所有连接的强度平均值，SEC{Ci}表示将cluster i分割为大小接近的两部分所需要切断的连接的强度平均值。

最后，两个cluster之间的距离为

a是用来调节两个参量的比重的参数。

聚类结果

因为原算法需要使用一些额外的算法来进行图分割，我这里只是使用了一个简化版本的Chameleon算法，使用了绝对的inter-connectivity和closeness，没有对cluster的大小进行normalization（也就是没有考虑上面两条式子中的分母）。

对样品1进行聚类，分别取聚类数目k=5和8。类似于谱聚类，Chameleon算法也可以稳定地对数据进行聚类，不会对k的选择过分敏感。

经过多次的调整参数，我也终于把样品2的cluster给识别了出来……

从算法的角度来说，Chameleon可以用于识别形状特别的cluster，但是实际上调整参数殊为不易（当然这也跟我使用简化版算法有关），而且关键是层次聚类的效率终归不高。所以Chameleon可以在一些特殊的场合使用，个人认为不是一个十分通用的算法。

3、PCCA

PCCA算法的全称是Perron Cluster Cluster Analysis，名称里有两个cluster是因为这样简写就可以和PCA区分开来（无语……）。PCCA并不是设计来处理传统的聚类问题的，而是专门用于得到马尔科夫链中的cluster。当然，对于一般的聚类问题，只要根据系统特点构造出一个概率转移矩阵，也可以使用PCCA算法。

要解释马尔科夫模型中的cluster，让我们想象有一只跳蚤在了数据点间跳跃转移。它下一个时刻跳到特定数据点上的概率，只跟它当前落在哪个数据点上有关，这显然是一个经典的马尔可夫过程。再让我们假定，跳蚤在点与点之间的跃迁概率跟数据点的“距离”成反比。如果数据点可以分成几个分界明显的cluster，跳蚤大多数时间就只会在某个cluster内部的数据点间转移，在cluster之间的跳跃则相对罕见。

先解释一下马尔科夫模型的一些性质。马尔科夫模型需要的是一个转移矩阵A，元素A(i,j)表示系统从状态i转移到状态j的概率。矩阵的每一列元素之和必须为1，这是因为转移概率总和必须为1。转移矩阵有一个本征值为1的本征矢量，对应着系统的稳态，亦即系统到达这个状态后，它在各个状态的概率分布就不会再发生变化。

为了说明PCCA的原理，我们直接来考虑最为极端的情况，也就是系统由几个完全分离的cluster所构成。对于最为极端的情况，如果系统只能在cluster内部转移，而完全不会在cluster之间转移，那么转移矩阵A就会是分块矩阵的形式，比如下面的系统就可以分为两个完全不连通的cluster。

这样的一个矩阵存在着不止一个本征值为1的本征矢量，因为它的每个分块都可以看做一个转移矩阵，都会对应着一个稳态。比如对于上面的矩阵A，下面两个本征矢的本征值都为1。

如果我们得到的矢量都是这样的理想形式，那么聚类就很简单了，只要得到本征值为1的全部本征矢，把对应的元素大于0的数据点归为一类就可以。十分可惜的是，由于这两个本征矢量是简并的，它们线性叠加产生的矢量也是矩阵的本征值为1的本征矢量，比如这个矢量：

因此PCCA算法的思路，就是要从计算得到实际本征矢量，反推得到理想矢量，从而实现聚类。

如果将计算得到的k个本征值为1的本征列矢量并排合并，成为一个N*k的矩阵，那么矩阵的每一行可以看成对应与数据点的一个坐标。对于理想本征矢（对应下图蓝色基矢），数据点都是落在坐标轴上（因为除了所属的cluster所对应的那个本征值，其余的维度都是0），比如下图的红色和黄色的数据点。但是由于实际得到的本征矢量是理想本征矢的线性叠加，所以基矢就发生了旋转（对应黑色基矢）。

尽管如此，每个cluster的数据点落在一条直线上的性质并没有发生改变。现在的问题就变成了如何找到这些直线的方向。PCCA首先找到离原点最远的数据点，这个点相对于原点的矢量，就对应了一个理想本征矢。然后每一次都找与已知的理想矢量垂直（相对原点），而又离原点最远的数据点。只要找k次，就能找到所有的理想矢量。最后看数据点落在哪个方向上，就可以知道它们属于哪个cluster。实际情况下，矩阵并不会是完全的分块矩阵，所以除了第一个本征矢，其余本征矢量对应的本征值不会完全为1，而是接近于1。cluster之间的转移几率越小，这个算法的准确性自然越高。

聚类结果

PCCA算法的一个优点就是，它可以根据本征值来直接得到cluster的数目，有多少个接近于1的本征值就应该分多少个cluster。当然这也是一个限制，聚类数目不能随意给定，完全由数据性质决定，如果cluster的分界不明显，那么可能聚类就完全无效。

下面是PCCA对样品1的聚类结果。第一幅图是由算法自动判定聚类数目，正好为3，聚类十分成功；第二幅图是人为地要求分为4个cluster，结果算法就基本失效了。

PCCA是除了Chameleon算法外，对样品2聚类结果最好的一个算法。不过它并不能正确地判断出cluster的数目，总共划分出了7个cluster，这里显示了前5个。

PCCA可以自动判定cluster数目，而且也能得到非凸型的cluster，还可以适用于概率转移矩阵的聚类，看上去确实是一个性能比较好的聚类算法。不过，PCCA对数据的性质特征有比较强的预设，当数据性质偏离理想状况较远时，算法的稳定性有待考验。

4、SOM

之所以尝试SOM聚类，主要是因为这是基于神经网络的一种算法，而神经网络本身又是机器学习中的一个重要方法，所以就自己实践一下体会体会。

所谓SOM指的是Kohonen提出的竞争网络，全称是self-organizing map。这个算法有着非常多的改进和变种，在网络的拓扑结构、节点之间的反馈方式等方面各有不同。更一般地说，SOM应该是一个降维算法，它可以将高维的数据投影到节点平面上，实现高维数据可视化，然后也可以继续根据降维之后的数据再进行聚类，就像谱聚类一样。不过，因为这里仅仅是个人的算法尝试，所以我就使用最简单的方式，直接使用SOM进行聚类。

竞争网络，顾名思义就是网络节点相互竞争。对于每一个输入的数据点，网络节点都要进行竞争，最后只有一个节点获胜。获胜节点会根据赢得的数据点进行演化，变得与这个数据点更匹配。如果数据可以明显地分为多个cluster，节点变得跟某个cluster内的数据点更为匹配，一般而言就会跟其它cluster不太匹配，这样其它节点就会赢得了其它cluster的数据点。如此反复学习，每个节点就会变得只跟特定的一个cluster匹配，这样就完成了数据点的聚类。

SOM需要输入数据点的坐标矩阵，对应的，每个网络节点也有一个坐标，初始时刻随机赋值。每次输入一个数据点，与这个数据距离最近的节点获胜，获胜点的坐标向着这个数据点的方向偏移。聪明的看官们肯定发现了，这个简单化的SOM算法跟K-means算法思路基本一致，确实一些文章也提到，在节点数目偏少的情况下SOM的结果就类似于K-means。

聚类结果

SOM的聚类结果确实跟K-means比较类似，不过当聚类数目取为4时，经常也能正确的结果，而不会聚成4个cluster，这个跟学习时间以及节点的初始值有关。这应该是因为SOM的学习方式与K-means直接求平均不同。至于对样品2的聚类，SOM也跟K-means类似，就不贴出来了。

这个SOM聚类只是个人试水，并不能真正代表SOM聚类的最佳效果，仅作参考。

5、Affinity Propagation

Affinity Propagation（简称AP）是一个比较新的算法，在07年发表在Science上面，可见肯定是有一些独到之处的。

个人认为，AP算法的具体实现步骤没有很直观的物理意义，大致上就是一个网络自动演化的过程，实现起来也并不复杂。感兴趣的童鞋可以直接到算法作者的网页(http://www.psi.toronto.edu/index.php?q=affinity%20propagation)，上面也提供了各个版本的实现程序，可以直接拿来使用。

这里我就不详细地叙述算法的实现步骤了，只是介绍一下这个算法的一些特点。

1）AP只要求输入数据点之间相似性矩阵，而且还不需要是对称阵。从这个角度来说，适用范围非常大。

2）AP算法的核心是对于每个cluster找到一个代表数据点exemplar，使得cluster内其它数据点到这个点的距离平方和最小。这是AP算法的一个很大的优点，因为它不仅能完成聚类，还可以给出这个类别的代表。很多时候我们聚类也就是为了找出代表而已。

3）AP算法不能直接知道聚类的数目，它不仅不能判定合适的聚类数目，甚至在聚类完成前，用户都不知道这次会聚出多少个cluster来，只能自己慢慢调整参数，多次尝试……

AP算法的目标函数跟K-means类似，不过中心点不是平均值，而是真实的一个数据点。其实这个算法可以说是K-centers的一个高效实现，但归根到底得到的也就是K-centers最佳情况下的结果而已，跟K-means也类似，都是大小接近的凸型cluster，所以我就不贴结果了。可以说，当你想得到K-centers的结果时，AP算法是你的最佳选择，但如果你的目标不在于此，那就不要用这个方法了。(Via:中国统计网)