聚类--谱聚类

来源：互联网发布：android图片压缩算法编辑：程序博客网时间：2024/04/28 18:16

前言：关于谱聚类，已经有很多厉害的老师和大牛写过教程博客等，也有很不错的tutorial文章可供参考。此博文仅记述个人的一些总结、思考、疑问，算是对现有谱聚类学习资源的一个小补充。

1. 谱聚类简述

说到聚类，可能最先想到的就是经典的Kmeans算法。但是，Kmeans的应用是有前提条件的，它假设(目标式中的)误差服从标准正态分布，因此，Kmeans在处理非标准正态分布和非均匀样本集时，聚类效果会比较差。对于非专业滴同志们，这种统计学术语听起来总是让人觉得有点虚、不好懂，下面给两个图，直观感受一下。

图1.1 (a)非标准正态分布集 (b)非均匀分布集 (c) 处理two moon数据集kmeans可能的聚类结果

对于图1.1中给出的情形，用Kmeans聚类就显得有些力不从心了，比如对于two moon这个toy dataset，理想情况下的结果应该和(a)中所示一致，但kmeans聚类后的效果可能是(c)中的样子。

为了处理这种奇奇怪怪形状分布的数据，谱聚类(spectral clustering)该出场了。不同于Kmeans采用数据点和聚类中心的欧式距离作为相似性度量，谱聚类是基于数据的相似度矩阵，也即，只用给出数据的相似度矩阵而不用给出数据本身，谱聚类就可以得到最终的聚类结果。这里先提一句，谱聚类算法得到最终聚类结果，一般是通过Kmeans，就是说Kmeans成为了谱聚类算法的一部分。另外应该指出的是，谱聚类是一类基于谱方法的聚类算法，只不过平常提到谱聚类，一般都是指代最原始和经典的那两个(Ncut，Ratio Cut)，在很多文章和博客里，看到的算法伪代码就是它们的。还有，为什么取名叫谱聚类，个人觉着这个“谱”字，正表明了求解过程与特征分解有关。

2. 公式化表述

在这一节，我们来看一下图论中的标准切割(Ncut)和比例切割(RatioCut)，并由此得到标准化的和非标准化的谱聚类目标函数。下面咱们先简单了解一些图论中的概念（这方面基础内容不了解的话，可求助wikipedia/baidu，或见其它博文介绍，或者从公式推导角度直接往下看，也没多大问题）。

设G=(V,E)表示一个带权无向图，其中V={vi}ni=1是顶点，可以对应成n个数据点，E是用来连接顶点的边的集合，每条边上带有权值，vi和vj之间的权值大小用wij表示，由此可得到数据之间的相似度矩阵W，它的第i行第j列元素就是wij。定义顶点vi的度为di=∑nj=1wij，则度矩阵D就是由{di}ni=1作为对角元素组成的对角矩阵。传闻中的拉普拉斯矩阵(Laplacian matrix)就可以这样定义L=D−W，标准形式为L~=In−D−1/2WD−1/2，这里In表示大小为nxn的单位阵。

继续，来看看什么是“割–Cut”。给定一个顶点集A⊂V，其补集表示为A¯。评估集合A的大小有两种方式，一种是用A中的顶点个数|A|来衡量，一种是用A中的顶点的度之和vol(A)=∑i∈Adi来衡量。A和A¯之间的割可以这样来定义：cut(A,A¯)=∑i∈A,j∈A¯wij。

给了一堆概念和定义，一不小心就让人如坠云雾，为什么要扯出什么图论？这个“割”真的是用小刀割东西的那个“割”嘛？首先第一个问题，把数据点看作是一个图中的顶点，继而把对数据的处理，转化为对图的分析，这种隐去实际意义，抓住问题本质的过程，就是建模，模型建好了，在理论上求解完毕了，再反过来对应回一个个实际的问题。事实上，很多实际生活中的问题，就是通过图来分析的，比如实现手机地图软件里的实时导航功能，就是一个路径规划的问题—-两个目的地，即图中两个顶点，是否可达，即是否存在一系列边将这两个顶点连起来，最短路径又是哪条，即怎么走最近。第二个问题，就是你想的那个割，这里割的是图中的边。想一想，要把数据点聚类成若干个簇，对应就是把图中的顶点分成若干个组，我们很自然地希望组内顶点间的关联度大(wij大)，组间顶点的关联度小，也就是说，我们可以把权值小的那些边给割断，保留权值大的那些边，以达到分组的目的。

2.1 由Ratio Cut到非标准谱聚类

假设我们要把顶点集分成k个组—-聚类成k个簇，A1,...,Ak，满足Ai∩Aj=∅，A1∪,...,∪Ak=V。我们想让被割去的都是一些权值小的边，即保证了关联度大的顶点被分到同一个组内，那么我们可以最小化下面的比例切割目标式：

R a t i o C u t (A 1, . . ., A k) = \sum i = 1 k c u t ( A i , A i ¯ ) | A i | (2.1)

式(2.1)中，分母的作用是寻求一种balance，努力让各个簇的大小相近，避免有的簇包含很多数据点，而有的簇仅包含一个或几个数据点。这里提一句，还有一种最小割(mincut)，就是只将分子上的各项求和，没有考虑分母项，即对各簇的大小没有限制。式(2.1)这个目标式好提出，但是不好解啊—-NP hard问题，怎么办，那就绕点儿路，退而求它的松弛问题(relaxed problem)，凑合得个近似解好了。

定义k个簇指示向量(indicator vector) fj=(f1j,...,fnj)′，其中

f i j = ⎧ ⎩ ⎨ 1 / | A j | - - - \sqrt, 0, if v i \in A j otherwise (2.2)

上式中

i=1,…,n,j=1,…,k。将这k个指示向量按列排成矩阵，则得到指示矩阵(indicator matrix)

F∈Rn×k，根据前面给的定义，可以推出

F 满足：

F T F f' i L f i = I k = c u t ( A i , A i ¯ ) | A i | (2.3) (2.4)

这里式(2.4)的推导对于谱聚类的提出和理解还是比较重要的，它将一个图切割的问题，巧妙且简洁地用数据的拉普拉斯矩阵和簇指示向量来表示，继而聚类问题实际就变成了求解簇指示向量

f。为了不冲淡主线，证明部分作为附录补在文末，供有需要者参考。基于式(2.4)，最小化Ratio Cut就可以写成如下的问题：

min F \in R n \times k T r (F T L F) s . t . F defined as Eq.(2.2), F T F = I k . (2.5)

Tr()表示矩阵的迹，即对角元素之和。式(2.5)是带有正交约束，且变量为离散值的最小化问题，这里相比式(2.1)只是表述形式变了，依然是不好求解的，这时就需要发挥不要脸之精神，对于使问题变得困难的那个离散定义，我们假装看不见，假装看不见，然后问题就变成了：

min F \in R n \times k T r (F T L F) s . t . F T F = I k . (2.6)

这个就是Ratio Cut的松弛版问题(the relaxation of Ratio Cut)，原来指示向量的元素只能取固定的几个非负值，现在我们将它松弛成可以取任意实值。好了，我们给这个新的比较好解的问题取个名儿，对了，就是谱聚类！这里拉普拉斯矩阵L是非标准的形式，所以式(2.6)就是非标准谱聚类的目标函数。

2.2 由Ncut到标准谱聚类

标准切割(Normalized Cut, Ncut)和比例切割很类似，只是衡量簇的大小时是采用的另一种方式，Ncut的目标式是：

N C u t (A 1, . . ., A k) = \sum i = 1 k c u t ( A i , A i ¯ ) v o l ( A i ) (2.7)

同样，定义k个指示向量

fj=(f1j,...,fnj)′，其中

f i j = ⎧ ⎩ ⎨ 1 / v o l (A i) - - - - - - \sqrt, 0, if v i \in A j otherwise (2.8)

这时，我们可以推出

F 满足：

F T D F f' i L f i = I k = c u t ( A i , A i ¯ ) v o l ( A i ) (2.9) (2.10)

基于式(2.10)，再次假装看不见指示向量的离散定义，最小化松弛版的NCut就是：

min F \in R n \times k T r (F T L F) s . t . F T D F = I k . (2.11)

我们令

H=D1/2F，则上式可以写成：

min H \in R n \times k T r (H T L ~ H) s . t . H T H = I k . (2.12)

采用标准化的拉普拉斯矩阵

L~，我们得到了标准形式的谱聚类目标式。

3. 目标函数的求解

3.1 瑞利商(Rayleigh quotient)

设A∈Rn×n是一个对称阵，0≠x∈Rn，称

R (x) = x T A x x T x (3.1)

为矩阵

A的瑞利商。下面讨论对称阵的特征值与它的瑞利商的极值之间的关系。

定理：将A的特征值（都是实数）按从大到小的顺序排列为λ1≥λ2≥...≥λn，则有：

λ 1 = max x \neq 0 R (x), λ n = min x \neq 0 R (x)

证明：令

Λ=diag(λ1,...,λn)，对于对称阵

A，一定存在正交阵

Q使

A对角化，即

QTAQ=Λ，将

Q按列分块为

Q=(q1,...,qn)，则有

Aqj=λjqj，注意这里

q1,...,qn是

A的两两正交的单位特征向量。对于非零向量

x，存在一组数

c1,...cn，使得

x = c 1 q 1 + . . . + c n q n (| c 1 | 2 + . . . + | c n | 2 \neq 0) .

因此有

A x = c 1 λ 1 q 1 + . . . + c n λ n q n

于是

R (x) = x T A x x T x = λ 1 | c 1 | 2 + \dots + λ n | c n | 2 | c 1 | 2 + \dots + | c n | 2

由此可得

λn≤R(x)≤λ1，容易验证

R(q1)=λ1，

R(qn)=λn，证毕

■。

对于非标准谱聚类目标式(2.6)，因为fTifi=1为常数，所以

min f T i L f i \Leftrightarrow min R (f i)

即

fi的解就是拉普拉斯矩阵

L的最小特征值对应的特征向量，因此，

F就是由

L的前k个最小的特征值对应的特征向量按列组成（同样，

H就是

L~前k个最小的特征向量）。接下来，我们得从这个连续取值的指示矩阵

F出发，想办法得到最终的聚类结果。如果是原问题(Ratio Cut/Ncut)的离散解，那我们直接对照定义就可以得到数据所属的簇标号，现在这种连续的情况，我们可以把

F的每一行，看作是一个数据点在k维空间中的表达，然后采用Kmeans算法，对这n个k维的数据点进行聚类，得到最终的聚类结果。这里Kmeans操作，可以看作一个离散化的过程，所以谱聚类实际上做的，就是先将NP–hard的离散问题松弛为好求解的连续问题，再将连续问题的解离散化，得到最终的聚类结果。

在上一节中，我们假装看不见离散约束从而得到了谱聚类的目标式，大家心里一定会有困惑，哪能这样搞，那最后的结果肯定不是原问题的解了呀。的确，这么粗暴地简化问题，在理论上是不能保证我们最后得到的解，和原问题(Ratio Cut/Ncut)的解之间有什么必然联系的，只不过在实际问题中，谱聚类取得的效果的确还是不错的。

图3.1 一个反例：Ratio Cut和谱聚类得到截然不同的解，蓝色的线表示Ratio Cut切割方式，只用两刀，绿色线是非标准谱聚类的结果，需要切k刀，具体细节说明请参见参考文献[1]。

【注】：根据定义，可知拉普拉斯矩阵L是一个对称半正定阵，即所有特征值都为非负实数，且L的最小特征值为0。对于一个全连通图，L只有一个0特征值，且对应的特征向量为常1向量(constant one vector)，就是所有元素都取相同值的单位向量，这种指示向量其实是不具备区分信息的，或者说，它表示把原图划分为一个空集和图自身。所以有时候，我们直接把这个特征向量去掉，只对k-1个指示/特征向量进行分析。对于一个含c个连通分量的图，L有c个0特征值，对应的特征向量就是簇指示向量，可以直接区分出哪些数据点位于同一个连通分量内。另外，除去通过Kmeans得到最终的聚类结果，也有人用谱旋转(spectral rotation)得到，在此不做展开。

3.2 拉格朗日乘子法

在这一小节，我们直接从拉格朗日乘子法的角度来求解谱聚类。设Λ∈Rk×k为拉格朗日乘子矩阵，式(2.6)的拉格朗日函数为：

J (F) = T r (F T L F) - T r (Λ (F T F - I k)) (3.2)

求导并置零：

\partial J ( F ) \partial F = 2 L F - F (Λ + Λ T) = 0 (3.3)

记

Λ~=12(Λ+ΛT)，有

LF=FΛ~，所以

F就是

L的特征向量矩阵，

Λ~是相应的特征值组成的对角矩阵。

J(F)=Tr(FTFΛ~)=Tr(Λ~)，要最小化

J(F)，自然取

L的前k个最小的特征值，

F则由对应的特征向量组成。

4. 从图重构角度再理解

4.1 正交非负约束

谱聚类目标式中，指示向量可取任意实值，而原始Ratio Cut/Ncut中，指示向量元素是离散的非负值。如果我们在松弛过程中，带上非负约束，那么新的模型应该更接近于原始问题，公式化表述就是：

min F \in R n \times k T r (F T L F) s . t . F T F = I k ， F \geq 0. (4.1)

式(4.1)是一个带有正交非负约束的问题，它的解实际上已经是离散的了。由于非负性，

F中的元素大小直接反应了数据与簇之间的关系紧密程度，这个性质被广泛用于非负矩阵分解中。所以，如果我们能求出这个问题的解

F，那么最终的聚类结果可以直接根据

F每一行的最大元素的列索引得到，而不再需要额外的后处理步骤了，如前文提及的Kmeans。这是一个比较好的特点，因为这样做，算法的性能不再受制于其它算法的性能，要知道，Kmeans本身就对初始化十分敏感，不同的初始点，可能得到截然不同的聚类结果。

图4.1 正交非负约束下的F示意图。这两个约束使得F是一个离散取值的、稀疏的、具有物理意义的矩阵，物理意义就是上面说的元素值对应数据和簇的关系。示意图中，F是一个“9x3”的指示矩阵，表示有9个数据点，分别属于3个不同的簇，比如第一个点属于第一个簇，则F的第一行第一列元素非零（实心黑圆），第一行其它元素为0，指示向量（F的列）是单位向量。

4.1 图重构角度

在这一节，我们关注一种特殊情形，并得到谱聚类的图重构理解视角。如果我们构建了一个图，它的相似度矩阵W，恰好是一个双随机矩阵（各行和、列和均为1），则度矩阵D就简化为一个单位阵I，此时拉普拉斯矩阵L=I−W=L~，也就是说L已经是标准化的，此时Ratio Cut问题等价于Ncut问题。此时式(4.1)可写为如下形式：

min F T F = I, F \geq 0 ∥ W - F F T ∥ 2 F . (4.2)

证明：注意有

W是对称双随机阵，

Tr(WTW)和

Tr(FFTFFT)是常数项。

Eq.(4.1) \Leftrightarrow max F T F = I, F \geq 0 T r (F T W F) \Leftrightarrow min F T F = I, F \geq 0 - T r (W T F F T + F F T W) \Leftrightarrow min F T F = I, F \geq 0 T r [W T W - W T (F F T) - (F F T) W + (F F T) (F F T)] \Leftrightarrow Eq.(4.2) ■

式(4.2)就体现了图重构的思想。W是原本构建好的图（相似度矩阵确定后，图也确定了，故此处可用图代指），现在我们用簇指示矩阵来重构一个图，使得重构图和原始图W相近，这样使得重构图，不仅具有原图的信息，而且在优化过程中，获得了很好的结构—-分块对角阵。一个块就对应一个连通分量，往往一个连通分量内的数据属于同一个簇，或者一个连通分量直接就对应了一个簇，实际中，由于冗余特征和噪音的存在，构建的原图往往只有一个连通分量。所以在这种特殊情况下，谱聚类可以看作是在用一种具有结构信息的重构图，来近似原始构建的图，如图4.2所示。

图4.2 指示矩阵重构图示意

5. 与K-means的联系

Kmeans公式化表述的证明过程可参见另一篇谈PCA的博文，这里只贴结果，Kmeans的目标函数可以写成如下形式：

min F \in R n \times k T r (F T K F) s . t . F T F = I k ， F \geq 0. (5.1)

这里的

F也是指示矩阵，不过是这样定义的：如果数据

xi是第j个簇，则

Fij=1，否则为0。

X∈Rd×n表示数据矩阵，n为数据个数，

K=XTX可看作一种相似度矩阵，这里采用的是标准内积线性核。对于Kernel Kmeans，

K~=ϕ(X)Tϕ(X)，

ϕ()是一种高维映射。

与谱聚类对比，Kmeans的交替迭代过程，相当于是直接在求一个离散定义的带正交约束的问题，没有像谱聚类一样对原问题进行松弛。两者构建相似度矩阵、以及定义指示矩阵的方式不同，但从形式上来看，两者的确有一些共通之处。

6. 谱聚类的不足

因为要计算拉普拉斯矩阵的特征分解，谱聚类一个很大的问题是计算复杂度高，达到O(n3)，n为数据点个数。随着大数据时代的来临，如何改进谱聚类，使之能处理大规模数据，是一个很有价值的研究方向，比较典型的解决思路有通过Nystro¨m方法降低计算复杂度，有通过分步解决，先选取一部分代表点得到聚类标签，再根据代表点和其它点间的某种关系（常利于稀疏学习等技术构建这种关系），得到剩余的点的聚类标签。

另一个不足是out of sample的问题，即给定一个新的测试样例，无法方便地得到其聚类标签，必须重新对整个数据集（包括新加的样例）来一遍完整的谱聚类流程。这个问题的一种解决思路是，得到训练数据的指示矩阵F后，将F的行作为线性回归中的标签项（F−XTθ），求出对应的回归系数矩阵θ∈Rd×k，对于新来的样例，根据θ得到其k维空间中的表达（拟合标签），再直接对所有数据的k维表达进行聚类就行。

7. 一个小问题

本文写作过程中主要参考了A Tutorial on Spectral Clustering，在论文第九页，定义Ratio Cut和Ncut的目标函数时，我始终觉得有点问题，少了一个系数1/2，虽然这并不影响最终结果，因为对于最小最大问题，常系数是可以去掉的。不过单纯从加深理解的角度，还是可以看细一点，在此贴出我的问题，还希望正好知道这一点的人不吝赐教。

上面的Ratio Cut 和Ncut的目标式中，第二个等式我觉得有问题，举个例子，k=2时：

c u t (A 1, A 2) = 1 2 (W (A 1, A 2) + W (A 2, A 1)) = W (A 1, A 2)

所以两个目标式中的第三项，都比第二项要少一个系数

12，所以个人认为这里是个typo，或者是我的理解仍有误。

主要参考文献/资源

【1】Luxburg U V. A tutorial on spectral clustering[J]. Statistics and Computing, 2007, 17(4):395-416.
【2】http://blog.csdn.net/v_july_v/article/details/40738211
【3】http://blog.pluskid.org/?p=287
【4】http://www.cnblogs.com/xingshansi/p/6702188.html

附录—–Ratio Cut相关证明

先贴一个很重要的公式，这一项在一些地方被称作图正则项(graph regularization)，可以将有监督学习扩展为半监督学习，它可以从流形学习的观点来解释，这里不展开。

注意，这里

f是一个n维指示向量，即

fi表示一个实数，如果把

f换成指示矩阵

F，

fi表示一个n维向量，则只需将上式中实数的平方替换成向量的范数的平方，即

||fi−fj||22，等式依然是成立的。这个公式可以这样来证明：
这里写图片描述

在参考文献[1]中，讲Ratio Cut时，是用二聚类情形作为例子，进行了较为详细地引出和推证，对于更为一般的情形，原文中只是给出了结论，网络上现有的博文，大多也是直接贴上原文中在二聚类情形下的证明。在这里，我们给出一般情形下的证明（本文正文中也是只关注了一般情形），证明过程跟二聚类情形是很类似的。

我们要证明的结论是：f′iLfi=cut(Ai,Ai¯)|Ai|，i=1,...,k，k是簇的个数，fi是第i个子集Ai的指示向量，fi的第p个元素fip指示第p个数据点是否属于子集Ai。

f' i L f i = 1 2 \sum p, q = 1 n w p q (f i p - f i q) 2 = 1 2 ⎛ ⎝ \sum p \in A i, q \in A i ¯ w p q (1 | A i | - - - \sqrt - 0) 2 + \sum q \in A i, p \in A i ¯ w p q (0 - 1 | A i | - - - \sqrt) 2 ⎞ ⎠ = c u t ( A i , A i ¯ ) | A i |

注意到

f′iLfi=(FTLF)ii，所以

R a t i o C u t (A 1, . . ., A k) = \sum i = 1 k c u t ( A i , A i ¯ ) | A i | = \sum i = 1 k f' i L f i = T r (F T L F)

Ncut的相关结论也可按照上述过程推证，此处略。

0 0