by pluskid, on 2009-02-25, in 
本文是“漫谈 Clustering 系列”中的第 6 篇，参见本系列的其他文章。
如果说 K-means 和 GMM 这些聚类的方法是古代流行的算法的话，那么这次要讲的 Spectral Clustering 就可以算是现代流行的算法了，中文通常称为“谱聚类”。由于使用的矩阵的细微差别，谱聚类实际上可以说是一“类”算法。
Spectral Clustering 和传统的聚类方法（例如 K-means）比起来有不少优点：
- 和 K-medoids 类似，Spectral Clustering 只需要数据之间的相似度矩阵就可以了，而不必像 K-means 那样要求数据必须是 N 维欧氏空间中的向量。
- 由于抓住了主要矛盾，忽略了次要的东西，因此比传统的聚类算法更加健壮一些，对于不规则的误差数据不是那么敏感，而且 performance 也要好一些。许多实验都证明了这一点。事实上，在各种现代聚类算法的比较中，K-means 通常都是作为 baseline 而存在的。 :p
- 计算复杂度比 K-means 要小，特别是在像文本数据或者平凡的图像数据这样维度非常高的数据上运行的时候。
突然冒出这么一个要求比 K-means 要少，结果比 K-means 要好，算得还比 K-means 快的东西，实在是让人不得不怀疑是不是江湖骗子啊。所以，是骡子是马，先拉出来溜溜再说。不过，在 K-medoids 那篇文章中曾经实际跑过 K-medoids 算法，最后的结果也就是一个 accuracy ，一个数字又不能画成图表之类的，看起来实在是没意思，而且 K-means 跑起来实在是太慢了，所以这里我还是稍微偷懒一下，直接引用一下一篇论文里的结果吧。
结果来自论文 Document clustering using locality preserving indexing 这篇论文，这篇论文实际上是提的另一种聚类方法（下次如果有机会也会讲到），不过在它的实验中也有 K-means 和 Spectral Clustering 这两组数据，抽取出来如下所示：
| k | TDT2 | Reuters-21578 | ||
| K-means | SC | K-means | SC | |
| 2 | 0.989 | 0.998 | 0.871 | 0.923 |
| 3 | 0.974 | 0.996 | 0.775 | 0.816 |
| 4 | 0.959 | 0.996 | 0.732 | 0.793 |
| … | ||||
| 9 | 0.852 | 0.984 | 0.553 | 0.625 |
| 10 | 0.835 | 0.979 | 0.545 | 0.615 |
其中 TDT2 和 Reuters-21578 分别是两个被广泛使用的标准文本数据集，虽然在不同的数据集上得出来的结果并不能直接拿来比较，但是在同一数据集上 K-means 和 SC (Spectral Clustering) 的结果对比就一目了然了。实验中分别抽取了这两个数据集中若干个类别（从 2 类到 10 类）的数据进行聚类，得出的 accuracy 分别在上表中列出（我偷懒没有全部列出来）。可以看到，Spectral Clustering 这里完胜 K-means 。
这么强大的算法，再戴上“谱聚类”这么一个高深莫测的名号，若不是模型无比复杂、包罗宇宙，就肯定是某镇山之宝、不传秘籍吧？其实不是这样，Spectral Clustering 不管从模型上还是从实现上都并不复杂，只需要能求矩阵的特征值和特征向量即可──而这是一个非常基本的运算，任何一个号称提供线性代数运算支持的库都理应有这样的功能。而关于 Spectral Clustering 的秘籍更是满街都是，随便从地摊上找来一本，翻开便可以看到 Spectral Clustering 算法的全貌：
- 根据数据构造一个 Graph ，Graph 的每一个节点对应一个数据点，将相似的点连接起来，并且边的权重用于表示数据之间的相似度。把这个 Graph 用邻接矩阵的形式表示出来，记为 $W$ 。一个最偷懒的办法就是：直接用我们前面在 K-medoids 中用的相似度矩阵作为 $W$ 。
- 把 $W$ 的每一列元素加起来得到 $N$ 个数，把它们放在对角线上（其他地方都是零），组成一个 $N\times N$ 的矩阵，记为 $D$ 。并令 $L = D-W$ 。
- 求出 $L$ 的前 $k$ 个特征值（在本文中，除非特殊说明，否则“前 $k$ 个”指按照特征值的大小从小到大的顺序）$\{\lambda\}_{i=1}^k$ 以及对应的特征向量 $\{\mathbf{v}\}_{i=1}^k$ 。
- 把这 $k$ 个特征（列）向量排列在一起组成一个 $N\times k$ 的矩阵，将其中每一行看作 $k$ 维空间中的一个向量，并使用 K-means 算法进行聚类。聚类的结果中每一行所属的类别就是原来 Graph 中的节点亦即最初的 $N$ 个数据点分别所属的类别。
就是这么几步，把数据做了一些诡异的变换，然后还在背后偷偷地调用了 K-means 。到此为止，你已经可以拿着它上街去招摇撞骗了。不过，如果你还是觉得不太靠谱的话，不妨再接着往下看，我们再来聊一聊 Spectral Clustering 那几个“诡异变换”背后的道理何在。
其实，如果你熟悉 Dimensional Reduction （降维）的话，大概已经看出来了，Spectral Clustering 其实就是通过 Laplacian Eigenmap 的降维方式降维之后再做 K-means 的一个过程──听起来土多了。不过，为什么要刚好降到 $k$ 维呢？其实整个模型还可以从另一个角度导出来，所以，让我们不妨先跑一下题。
在 Image Processing （我好像之前有听说我对这个领域深恶痛绝？）里有一个问题就是对图像进行 Segmentation （区域分割），也就是让相似的像素组成一个区域，比如，我们一般希望一张照片里面的人（前景）和背景被分割到不同的区域中。在 Image Processing 领域里已经有许多自动或半自动的算法来解决这个问题，并且有不少方法和 Clustering 有密切联系。比如我们在谈 Vector Quantization 的时候就曾经用 K-means 来把颜色相似的像素聚类到一起，不过那还不是真正的 Segmentation ，因为如果仅仅是考虑颜色相似的话，图片上位置离得很远的像素也有可能被聚到同一类中，我们通常并不会把这样一些“游离”的像素构成的东西称为一个“区域”，但这个问题其实也很好解决：只要在聚类用的 feature 中加入位置信息（例如，原来是使用 R、G、B 三个值来表示一个像素，现在加入 x、y 两个新的值）即可。
另一方面，还有一个经常被研究的问题就是 Graph Cut ，简单地说就是把一个 Graph 的一些边切断，让他被打散成一些独立联通的 sub-Graph ，而这些被切断的边的权值的总和就被称为 Cut 值。如果用一张图片中的所有像素来组成一个 Graph ，并把（比如，颜色和位置上）相似的节点连接起来，边上的权值表示相似程度，那么把图片分割为几个区域的问题实际上等价于把 Graph 分割为几个 sub-Graph 的问题，并且我们可以要求分割所得的 Cut 值最小，亦即：那些被切断的边的权值之和最小，直观上我们可以知道，权重比较大的边没有被切断，表示比较相似的点被保留在了同一个 sub-Graph 中，而彼此之间联系不大的点则被分割开来。我们可以认为这样一种分割方式是比较好的。
实际上，抛开图像分割的问题不谈，在 Graph Cut 相关的一系列问题中，Minimum cut （最小割）本身就是一个被广泛研究的问题，并且有成熟的算法来求解。只是单纯的最小割在这里通常并不是特别适用，很多时候只是简单地把和其他像素联系最弱的那一个像素给分割出去了，相反，我们通常更希望分割出来的区域（的大小）要相对均匀一些，而不是一些很大的区块和一些几乎是孤立的点。为此，又有许多替代的算法提出来，如 Ratio Cut 、Normalized Cut 等。
不过，在继续讨论之前，我们还是先来定义一下符号，因为仅凭文字还是很难表述清楚。将 Graph 表示为邻接矩阵的形式，记为 $W$，其中 $w_{ij}$ 是节点 $i$ 到节点 $j$ 的权值，如果两个节点不是相连的，权值为零。设 $A$ 和 $B$ 为 Graph 的两个子集（没有交集），那么两者之间的 cut 可以正式定义为：
\[ \displaystyle \text{cut}(A, B) = \sum_{i\in A, j\in B} w_{ij} \]先考虑最简单的情况，如果把一个 Graph 分割为两个部分的话，那么 Minimum cut 就是要最小化 $\text{cut}(A, \bar{A})$ （其中 $\bar{A}$ 表示 $A$ 的补集）。但是由于这样经常会出现孤立节点被分割出来的情况，因此又出现了 RatioCut ：
\[ \displaystyle \text{RatioCut}(A, B) = \frac{\text{cut}(A, \bar{A})}{|A|} + \frac{\text{cut}(A, \bar{A})}{|\bar{A}|} \]以及 NormalizedCut ：
\[ \displaystyle \text{NCut}(A, B) = \frac{\text{cut}(A, \bar{A})}{\text{vol}(A)} + \frac{\text{cut}(A, \bar{A})}{\text{vol}(\bar{A})} \]其中 $|A|$ 表示 $A$ 中的节点数目，而 $\text{vol}(A) = \sum_{i\in A}w_{ij}$ 。两者都可以算作 $A$ 的“大小”的一种度量，通过在分母上放置这样的项，就可以有效地防止孤立点的情况出现，达到相对平均一些的分割。事实上，Jianbo Shi 的这篇 PAMI paper：Normalized Cuts and Image Segmentation 正是把 NormalizedCut 用在图像分割上了。
搬出 RatioCut 和 NormalizedCut 是因为它们和这里的 Spectral Clustering 实际上有非常紧密的联系。看看 RatioCut ，式子虽然简单，但是要最小化它却是一个 NP 难问题，不方便求解，为了找到解决办法，让我们先来做做变形。
令 $V$ 表示 Graph 的所有节点的集合，首先定义一个 $N$ 维向量 $f$：
\[ \displaystyle f_i = \left\{\begin{array}{ll}\sqrt{|\bar{A}|/|A|} & \text{if } v_i \in A \\ -\sqrt{|A|/|\bar{A}|} & \text{if } v_i \in \bar{A} \end{array}\right. \]再回忆一下我们最开始定义的矩阵 $L=D-W$ ，其实它有一个名字，叫做 Graph Laplacian ，不过，我们后面可以看到，其实有好几个类似的矩阵都叫做这个名字：
Usually, every author just calls “his” matrix the graph Laplacian.
其实也可以理解，就好象现在所有的厂家都说自己的技术是“云计算”一样。这个 $L$ 有一个性质就是：
\[ \displaystyle f'Lf = \frac{1}{2}\sum_{i,j=1}^N w_{ij}(f_i-f_j)^2 \]这个是对任意向量 $f$ 都成立的，很好证明，只要按照定义展开就可以得到了。把我们刚才定义的那个 $f$ 带进去，就可以得到
\[ \displaystyle \begin{aligned} f'Lf &= \frac{1}{2}\sum_{i,j=1}^N w_{ij}(f_i-f_j)^2 \\ &= \frac{1}{2}\sum_{i\in A, j\in\bar{A}} w_{ij}\left(\sqrt{\frac{|\bar{A}|}{|A|}}+\sqrt{\frac{|A|}{|\bar{A}|}}\right)^2 + \sum_{i\in \bar{A}, j\in A} w_{ij}\left(-\sqrt{\frac{|\bar{A}|}{|A|}}-\sqrt{\frac{|A|}{|\bar{A}|}}\right)^2 \\ &= \text{cut}(A,\bar{A})\left(\frac{|\bar{A}|}{|A|}+\frac{|A|}{|\bar{A}|}+2\right) \\ &=\text{cut}(A,\bar{A})\left(\frac{|A|+|\bar{A}|}{|A|} + \frac{|A|+|\bar{A}|}{|\bar{A}|}\right)\\ &=|V|\cdot\text{RatioCut}(A,\bar{A}) \end{aligned} \]另外，如果令 $\mathbf{1}$ 为各个元素全为 1 的向量的话，直接展开可以很容易得到 $f'\mathbf{1} = \sum f_i = 0$ 和 $\|f\|^2 = \sum f_i^2 = n$ 。由于 $|V|$ 是一个常量，因此最小化 RatioCut 就等价于最小化 $f'Lf$ ，当然，要记得加上附加条件 $f \bot \mathbf{1}$ 以及 $\|f\| = \sqrt{n}$ 。
问题转化到这个样子就好求了，因为有一个叫做 Rayleigh quotient 的东西：
\[ \displaystyle R(A, x) = \frac{x'Ax}{x'x} \]他的最大值和最小值分别等于矩阵 $A$ 的最大的那个特征值和最小的那个特征值，并且极值在 $x$ 等于对应的特征向量时取到。由于 $f'f = \sqrt{n}$ 是常数，因此最小化 $f'Lf$ 实际上也就等价于最小化 $R(L, f)$ ，不过由于 $L$ 的最小的特征值为零，并且对应的特征向量正好为 $\mathbf{1}$ （我们这里仅考虑 Graph 是联通的情况），不满足 $f \bot \mathbf{1}$ 的条件，因此我们取第二个小的特征值，以及对应的特征向量 $v$ 。
到这一步，我们看起来好像是很容易地解决了前面那个 NP 难问题，实际上是我们耍了一个把戏：之前的问题之所以 NP 难是因为向量 $f$ 的元素只能取两个值 $\sqrt{|\bar{A}|/|A|}$ 和 $-\sqrt{|A|/|\bar{A}|}$ 中的一个，是一个离散的问题，而我们求的的特征向量 $v$ 其中的元素可以是任意实数，就是说我们将原来的问题限制放宽了。那如何得到原来的解呢？一个最简单的办法就是看 $v$ 的每个元素是大于零还是小于零，将他们分别对应到离散情况的 $\sqrt{|\bar{A}|/|A|}$ 和 $-\sqrt{|A|/|\bar{A}|}$ ，不过我们也可以采取稍微复杂一点的办法，用 $k=2$ 的 K-means 来将 $v$ 的元素聚为两类。
到此为止，已经有 Spectral Clustering 的影子了：求特征值，再对特征向量进行 K-means 聚类。实际上，从两类的问题推广到 k 类的问题（数学推导我就不再详细写了），我们就得到了和之前的 Spectral Clustering 一模一样的步骤：求特征值并取前 k 个最小的，将对应的特征向量排列起来，再按行进行 K-means 聚类。分毫不差！
用类似的办法，NormalizedCut 也可以等价到 Spectral Clustering 不过这次我就不再讲那么多了，感兴趣的话（还包括其他一些形式的 Graph Laplacian 以及 Spectral Clustering 和 Random walk 的关系），可以去看这篇 Tutorial ：A Tutorial on Spectral Clustering 。
为了缓和一下气氛，我决定贴一下 Spectral Clustering 的一个简单的 Matlab 实现：
function idx = spectral_clustering(W, k) D = diag(sum(W)); L = D-W; opt = struct('issym', true, 'isreal', true); [V dummy] = eigs(L, D, k, 'SM', opt); idx = kmeans(V, k); end |
最后，我们再来看一下本文一开始说的 Spectral Clustering 的几个优点：
- 只需要数据的相似度矩阵就可以了。这个是显然的，因为 Spectral Clustering 所需要的所有信息都包含在 $W$ 中。不过一般 $W$ 并不总是等于最初的相似度矩阵——回忆一下，$W$ 是我们构造出来的 Graph 的邻接矩阵表示，通常我们在构造 Graph 的时候为了方便进行聚类，更加强到“局部”的连通性，亦即主要考虑把相似的点连接在一起，比如，我们设置一个阈值，如果两个点的相似度小于这个阈值，就把他们看作是不连接的。另一种构造 Graph 的方法是将 n 个与节点最相似的点与其连接起来。
- 抓住了主要矛盾，忽略了次要的东西，Performance 比传统的 K-means 要好。实际上 Spectral Clustering 是在用特征向量的元素来表示原来的数据，并在这种“更好的表示形式”上进行 K-means 。实际上这种“更好的表示形式”是用 Laplacian Eigenmap 进行降维的后的结果，如果有机会，下次谈 Dimensionality Reduction 的时候会详细讲到。而降维的目的正是“抓住主要矛盾，忽略次要的东西”。
- 计算复杂度比 K-means 要小。这个在高维数据上表现尤为明显。例如文本数据，通常排列起来是维度非常高（比如，几千或者几万）的稀疏矩阵，对稀疏矩阵求特征值和特征向量有很高效的办法，得到的结果是一些 k 维的向量（通常 k 不会很大），在这些低维的数据上做 K-means 运算量非常小。但是对于原始数据直接做 K-means 的话，虽然最初的数据是稀疏矩阵，但是 K-means 中有一个求 Centroid 的运算，就是求一个平均值：许多稀疏的向量的平均值求出来并不一定还是稀疏向量，事实上，在文本数据里，很多情况下求出来的 Centroid 向量是非常稠密，这时再计算向量之间的距离的时候，运算量就变得非常大，直接导致普通的 K-means 巨慢无比，而 Spectral Clustering 等工序更多的算法则迅速得多的结果。
说了这么多，好像有些乱，不过也只能到此打住了。最后再多嘴一句，Spectral Clustering 名字来源于 Spectral theory ，也就是用特征分解来分析问题的理论了。
UPDATE 2011.11.23: 有不少同学问我关于代码的问题，这里更新两点主要的问题：
- 关于 LE 降维的维度和 Kmeans 聚类的类别数的关系：我上面的代码里，取成了一样的，但是这两者并不是要求一定要一样的。最初 Spectral Cluster 是分析聚两类的情况，就降到 1 维，然后用 thresholding 的方法来分成两类。对于K 类的情况，自然的类比就是降到 K-1 维，这也是和 LDA 保持一致。因为 Laplacian 矩阵的特征向量有一个全一的向量（对应于特征值 0 ），所以可以求 K 个特征向量然后把特征值 0 对应的那个特征向量去掉。但是在实际中并不一定要保持这两者一致的，也就是说，这个降维的维度可以作为一个参数进行调节的，选择效果好的参数。
- 关于示例代码：注意除非我在这里注明了是发布某个 software 或者 package 什么的，否则这里贴的代码主要都是为了示例作用，为了只显示出算法的主要部分，因此通常会省略掉一些实现细节，可以看成是“可执行的伪代码”，不推荐直接用这里的代码去跑实验之类的（包括其他 post 里的相关代码）。除非你想自己试验一下具体实现和改进，否则可以直接找网上现成的专用的 package ，比如 Spectral Clustering 可以考虑这个包。
kmeans都能比sc慢，那这kmeans实现的也太烂了
@sth4nth
呵呵，看来你没有仔细看我说的原因吧？不知道你说这话是有什么理论依据还是自己做过实验比较的呢？
实验就不用说了，写过code的都知道，从理论上kmeans就不可能比sc慢。
最简单的,对于kmeans我们也可以首先计算出所有pair wise inner product,然后用kernel kmeans进行迭代，这肯定比计算kernel矩阵的特征向量快多了。这相当是naive kmeans的dual。而在原空间上kmeans跟本不需要计算所有pairwise distance。更不用提kmeans还可以通过三角不等式还有基于spatial indexing的方法来加速了
@sth4nth
恩，多谢你给的信息。确实你说的方法我都没有去了解过。看来在文章里面下的结论确实有些武断了。
您好，
有两点疑问想向您请教
1：spectral clustering的话，我用matlab实现过，测试的相似矩阵是400*400的（400张图片的相似度），但是matlab运行时间很不稳定，刚才测试的结果是104秒左右，但是第2次就跑了半个小时。
2：我换了小点的数据12*12的，这会速度还可以接受，但是结果非常不稳定（12个图片的相似矩阵是没问题的），有时结果出错1个（这个是最好的结果），有时会多达5个。
我明白 这个是因为随机性。但是如果这样的话，这么不稳定，我们如何检测聚类结果呢？
再跟别的方法，比如分层聚类比较的时候，假如 偶尔一次结果比分层好 偶尔一次比分层差，这个问题如何解决呢？
写了很长，可能是废话，本人刚刚硕士一年，接触数据挖掘的时间不长，如果说问了白痴的问题，也请您花点时间为我解释下好吗？
@yao
你好！谱聚类的效果还不错吧，虽然有随机性，但是总体效果还挺好的，不过 12×12 的数据这么少，做出来稳定性不好也是自然的吧。400×400 感觉数据量也不大的，运行时间一个是解特征值，一个是做 k-means ，你可以看看瓶颈在哪里，运行时间和你的矩阵的稀疏程度有关的。不过话说回来，做这方面的东西，除非是要花心思优化算法的实现，否者一般的实验跑几个小时是很正常的事吧。 🙂
谢谢回复
如果可以的话，加你MSN 聊可以吗？
总体看效果是不错，但是能否在统计学基础上，判断spectral和kmeans还有分层的结果谁更好，有没有什么方法呢？
请问你比对过 矩阵不做稀疏处理，还有矩阵做了稀疏处理，表现如何吗？
现在运行时间不是问题了 400*400的 2秒左右（当然是经过稀疏处理，取最近邻居那种）。
只是 每次的结果不稳定，想证明那个方法更好，找不到路而已，因为spectral也有表现非常差的时候。
您说的 “12×12 的数据这么少，做出来稳定性不好也是自然的吧。400×400 感觉数据量也不大的”我非常感兴趣，能具体解释下马? 基于什么原理使得矩阵越大结果越稳定呢？非常急切的期待您的回答。
我的msn yaojiannan1985@hotmail.com 我做图像聚类的，能交个朋友吗？
@yao
你好，我的 MSN 在 About 页面里应该有，不过我不怎么上的，而且最近几天也比较忙。
其实这些东西我也是正在学，不过我可以就我自己了解的给你大致解释一下。K-means 和 Spectral Clustering 我想并没有能保证哪一个一定在所有的时候都比另一个要强，总的来说，Spectral Clustering 通常比 K-means 要好的，况且 K-means 本身也是有一定的随机性的，同样的参数前一次跑出来的结果和后一次跑出来的结果可能都不一样。但是 Spectral Clustering 既然是用 Laplacian Eigenmap 来做降维，那么它首先基于了 Laplacian Eigenmap 的一个假设，就是数据本身是分布在一个低维流形上的，如果数据本身不符合这个假设或者甚至完全抵触的话，那就不能那个保证一定能得到满意的结果了。
对于数据量的问题，像这里讨论的机器学习算法，实际上可以算作是基于统计的模型的，就好比说一个人的行为是不可预测的，但是一个群体的行为却是可以预测的，既然要想从数据里学习到合适的模型，那数据太少的话总是会出现偏差的，这就会导致不稳定了，就是说，结果要碰运气了。
@sth4nth
你用LDA试试看……kmeans有很多优化，sc也有很多优化。
@cerror
没有看过 LDA ，确实听说两者都有好多优化，不过好像在跑实验的时候都是直接用的最原始的解法。 :/
收藏了，慢慢吸收！我转载到我QQ空间里不算侵权吧？
@heshizhu
CC 协议的，不用于商业目的的话，注明出处就好了。 🙂
您好，我想请教个问题：如果只知道任意两个对象的距离(而不知道其他的对象标示），可以用那些聚类算法，不可以用那些聚类算法，最好用哪个聚类算法呢？K-Mean算法肯定不能用吧！
@heshizhu
你好，像谱聚类这样的算法就可以满足你的条件，主要是看聚类算法所需要的输入参数的结构就可以了。通常用谱聚类就可以了，不过也还是要看具体的应用场景的。
哦，非常感谢！
前段时间看Normalized Cut时，思想基本还能看懂，但是理论部分太抽象了，看了你的解释之后，基本上有点明白了，非常感谢！
作为最后k-means clustering的对象，为什么解出来的Eigenvectors要以Eigenvalue的大小来排列。有什么依据来着??
@Chrysalis
因为目标函数的值就等于特征值的和，所以取特征值最小的那些对应的特征向量。
@pluskid
答得不是我问的，我问的是选好特征向量之后，为什么要按照特征值大小来排列进而做最后的clustering。现在想应该是没有这个规定，因为最后做clustering是K-means 如果是用欧式距离的话 选定的那些特征向量如何排列都不会对距离产生影响。
当然在如何选特征向量上，是因该按照特征值大小来选取的，我理解的对不对??
spectral clustering最tricky的地方，就是选择了适当的vector “f”和对f的relaxation, 从而使得”对graph partition的objective function求极值问题”转化到了”关于Rayleigh Quotient求极值”的问题。
文中的Rayleigh Quotient拼写有错误..
@Chrysalis
因为选的时候已经排好序了，没有必要再把他打乱一下…… -.-///
好的，fixed typo……而且发现我一直都读错了，原来源自这里啊
你好，请问能说下在对高维文本数据的处理中，你所提到的高效的对稀疏矩阵求特征值和特征向量的方法具体有哪些么？谢谢
@dreaminyale
你好，我说的是 matlab 具有处理这样情况的能力，至于具体是什么算法，请去查阅数值计算相关的内容 （比如 Power Method 之类的）。
你好，请问你能说下Graph Cut和Spectral Clustering到底是什么关系么？谢谢！
@dreaminyale
你好，我的理解是 Graph Cut 和 Spectral Clustering 在数学推导上是一样的，也就是说在本质上是一样的，应该可以看作同样的数学模型在不同领域的应用或者说是从不同的理解角度来看待同一个问题。
我的理解，Spectral Clustering跟clustering并没有多少关系，重点还是Laplacian Eigenmap的降维。
降维的作用是把（高维）欧式空间的距离转换为流形上的距离。如果数据确实分布在某个低维流形上，那么这样做对于kmeans这类基于基于距离的算法是有意义的。
降维方法，或者说流形学习的方法，似乎目前并不完善。比方说，需要把数据降到多少维，很多时候都会难以把握。更夸张的是，JLMR上有一篇文章 “Manifold Learning: The Price of Normalization”，专门讲Laplacian Eigenmap, LLE, Hessian LLE, Diffusion Map等基于”spectral”的方法是如何的不靠谱。
同意前半部分观点。如果数据是欧氏的，那么K-means应该是最优的。(不知道有没这样的结论)但问题往往是数据不是欧氏，所以大家的做法是利用数据本身的结构把数据映射到欧氏空间中，再做K-means。
“Manifold Learning: The Price of Normalization”是从isometry的角度来看降维的，而spectral methods更多的考虑是映射的光滑性，因此它的观点是spectral methods往往都不是isometry。
可是聚类本身就是基于metric的……
如果高维到低维不能isometry，那降维后的聚类就没有太大意义吧
降维和聚类各有各的准则，而且这些准则往往是不一样的。比如单纯降维考虑的是isometry。聚类的话有很多准则，比如spectral clustering考虑的是smoothness(or locality), diffusion map考虑的是preserve diffusion distance。虽然说聚类本身是基于metric，但是这个metric不是standard的，需要去定义。
@C6H5NO2
恩，我也有这样的感觉，谢谢你提供的这个信息，我去看一看这篇文章，似乎很有趣的样子。 🙂
要用到 收藏了
感谢，请LZ继续加油~！
我看了一些关于sc的论文吧，但是没找到一篇确实提出可行的image segmentation的方法，我想找一篇是确切地提出怎么算识别率的方法，请楼主帮帮忙吧，解决一下这个问题。
我的mail是 mengqiuwang7@gmail.com,希望知道sc的image segmentation的识别率计算方法的高手跟我联系哈，主要就是你聚出了那几类，但是如何确定那就是你要多那一类呢，又没有中心点。
@mengqiuwang
你好，不是太明白你的问题，image segmentation 的识别率是什么意思？中心点又是什么？
可能我没有表述清楚吧，我的意思是说，假如有很多幅图片吧，每幅图片里面有河流、花草、山川，那么你对这幅图片进行了sc，最后形成了五类吧，你如何确定这五类里面哪一类是河流，哪一类是山川，哪一类又是花草呢？
Laplace Eigenmap的降维结果是一一对应的，可以聚了类找回点的原像，然后用原像的总体特征分类，比如总体色值，山偏绿，海偏蓝，或者其他一些如纹理。至于什么特征是山，什么特征是海，cluster是不可能知道的，因为根本没概念，只能有事先的参数或样本决定
@mengqiuwang
这个是没有办法做的吧，除非你有事先标记好的训练数据用于训练分类模型的话，但是这个属于分类的问题了，和 clustering 或者 segmentation 应该相对独立的，可以看作后续的处理过程。
这个大概要参考 Image Classification 方面的论文了，这个我也没有做过，估计可能要提取 SIFT 特征然后做 visual vocabulary 之类的方法来做，总之都是要有训练数据才行的。
你好，我关注你的日志好久了，请问你有谱聚类算法的matlab实现代码吗？我正在做这个毕业设计，我做的是改进NJW的算法，希望你可以发给我谱聚类的matlab源代码好吗？呵呵非常感谢你，我是真的愁的不知道怎么办好了。
邮箱：sunyxrizhao@yahoo.com.cn
@syx
你好，我这里没有现成的可以直接用 code ，不过你 Google spectral clustering matlab 能找到很多啦。
我是过来ym pluskid的，跟着何晓飞在做吗？博士还是硕士？
@gaohaidong
寒……硕士啊……
硕士也搞得这么理论。。。是明年上毕也吗？
@gaohaidong
我不是很清楚哪年，反正现在研一嘛…… =.=
不知楼主是否看过Jb Shi 和Malik 写的大名鼎鼎的Normailized Cut 算法的源代码，我看后感觉代码和论文有些地方有出入，代码反映出来的算法是基于边缘的图像分割算法，跟区域关系不太大。不知楼主认同这个观点不？
这个倒是没有看过，所以不知道具体是什么情况。
同学你好，仔细拜读了你的clustering系列，让我也萌生了写些自己总结的想法。
想问的是，spectral clustering中，最后的一步是k means吧，但此时每个数据点的维度是否只是1呢（即就是在第二小的特征向量上作聚类）？ 总觉得到这里有些怪怪的。
谢谢指教！
你好，维度不一定是 1 的，比如可以取 k 个最小特征值对应的特征向量来聚类，这样维度是 k 。
如果是分两类(k=2)，并且是用Laplacian求second smallest eigenvalue，那么其实降低后的维度就是1，推广到k类维度就是k-1。这是因为Laplacian的smallest eigenvector是值全为1的向量。
同学你好 请问下 我在有些资料上看的是前k个最大特征值对应的特征向量
再用kmeans聚类
或者是只取第2小的特征值对应的特征向量
是什么原因啊？难道是解的方程不同么？
是的，解的特征方程不同，便有最大最小的变种。
[…] 漫谈 Clustering (4): Spectral Clustering by pluskid, on 2009-02-25, in Machine Learning 47 comments […]
有个笔误：
其中 \\bar{A} 表示 A 的子集
不是 子集 应该是 补集
啊，确实，多谢指出，已经修正了。
请问关于谱聚类的阶次选取有何切实可行的办法
毕竟实际问题很难判断降到几维比较合适
请能谈谈吗？谢谢
随便说说：还是应该从问题出发，分析一下数据是什么样的结构。
要不就用cross validation选取合适的参数（比如到底是几维）
聚类的cross validation何解？
这个问题很困难啊，我是不知道有什么在实际中能有什么特别灵的办法。通常实际问题有它特定的结构可以分析。有时候通常会采用的一个启发式的办法就是把维度设成聚类类别数或者类别数减一。不过就聚类来说，选取合适的类别数本身就是很困难的问题啦！
我运行了博主的matlab程序，频繁的在特征值分解处出warning
opt = struct(‘issym’, true, ‘isreal’, true);
[eV eD] = eigs(L, D, k, ‘SM’, opt);
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 9.628061e-017.
一些背景：我的similarity矩阵比较稀疏。
请问，一，这会不会影响结果的准确性？是否有办法消除warnings？
二，有没有其他的spectral clustering程序可供参考？
多谢多谢
你好，网上应该有很多其他 spectral clustering 的程序的，我这个为了示例作用，写得比较简单，不推荐直接拿去实际中使用。不过对于你说的那个 warning ，是数据的问题，也和你的 similarity 矩阵构造有关系，如果确定前面的步骤都是合理的话，仅在求特征值一步，去除那个 warning 可以参考我最近的一篇 blog 文章：求最小的几个特征值。
我以前专门讨论过谱聚类和Kmeans聚类，后来发现其实这两个算法是一回事，只不过k取不同的值而已。具体分析在这里：Image Clustering using Local Discriminant Models and Global Integration
哦？居然是一样的？好，我去下载下来看看！thx
我还提供了code
嗯，赞！找到了！不过你主页上的照片失效了
拜读了，一直以为只有按照 similarity 矩阵计算出来的那几个才是 Laplacian 矩阵呢，原来只要满足半正定并且向量 1 是对应特征值 0 的特征向量就可以叫做 Laplacian matrix 啊？
顺便，你那个 HumanEva 数据集里的小人的图是用什么工具画的啊？
我让同学画的。任何一个拉普拉斯矩阵都对应一个图，这个图的权可以认为是similarity，比如在appendix里我就写出了那个图的边的权，虽然不是很严谨。更严谨的也可以写出来，但是太麻烦了，我就没放上去
唔，对呀！了解了，thx!
好ym啊啊。。。
我去……只要8行= =
看了学长的两篇博客，感觉都写的很好，解释的比老师讲的还清楚些。
请问spectral_clustering（）是matlab自带的嘛？
请 google spectral clustering matlab ，有很多工具包的。
请问为什么要对L=D-W取特征值特征向量呢?
你好，有几个问题没弄明白。第四步，求出k个特征向量后，把他们按列排成NxK的矩阵U，然后原始的数据点i就对应矩阵的第i行，然后对这N个K维向量进行Kmeans聚类。 这是为什么呢？ 按你说的用降维的思想理解的话，应该求出K个特征向量后，然后把原始数据点映射到这K个特征向量张成的空间，再进行kmeans聚类啊。
原始数据点为什么能对应于特征向量组成的矩阵的行向量？谢谢。
你好，其实降维的目标就是计算降维之后的点的坐标，而不是求低维空间的一组张成基。这里每一个特征向量都包含了所以数据点在那一个维度上的坐标。
我还是不太明白为什么这么映射，你可以给出一些证明吗。或者能推荐一些阐述这个问题的资料。
你好，因为问题的目标就是求数据的低维坐标，而特征值刚好就是这个问题的解，所以自然就是坐标了。如果你还不明白，可以参考这篇关于降维的文章，或者直接看 Laplacian Eigenmaps 的原始论文也是可以的。
你可以推算一下
U的每一行应该就是原始数据点在新空间中的坐标
你好，请问”3.求出L的前k个特征值（在本文中，除非特殊说明，否则“前k个”指按照特征值的大小从小到大的顺序） 以及对应的特征向量 。 “这一步中的k值如何确定？是由人为指定还是由计算机来处理，这个k需要和k_means中的类别数k一样吗？？谢谢！
K 是人工确定的，不一定要和 KMeans 中的 K 一样，有的时候会取成 KMeans 中的 K 再减去 1 ，不过通常也会作为一个参数来调。
博主，求问下，我现在在用c++实现这个算法，可是当找到k个特征值的时候，如果有复数了我们应该怎么处理，复数与实数怎么比较大小？
你好，如果你保证了你的矩阵是实对称的话，不应该有复数的特征值才对。如果是数值计算上的问题的话，可能需要你调整你使用的数值计算库了。
谢谢您这么迅速的回答我的提问。
仁兄伟才！ 受益匪浅，茅塞顿开
pluskid,为什么f’1=0，还有||f||^2=n，这两条等式都想不明白~~~能指导下否？
把 f_i 的定义式子直接带进去就得到了啊。
聚类技术和运筹学中的指派问题有没有区别呢？一直没想明白这一点。
请教～！
呃，我不知道运筹学中的指派问题是怎么定义的呢。如果是相似的话，也许是类似的问题，只是可能考察的角度和侧重点不一样。
如何找相似度函数呢？
可以看一下这篇文章第二部分，不好的一点是英文的，自己翻译下吧~~A Tutorial on Spectral Clustering
请问楼主是否知道用c++编程时，如何控制k？我从网上找了一些代码，多时计算了全部的特征值和特征向量，感觉这样就无法控制了k了。求解释呀，谢谢
这个不知道啊，你可以找一下有没有 power methods 之类的实现……
学长，我想问下，对于一个二分图（Bipartite Graph)的切割，如何将本文的方法套上去，大概怎么样的思路？
思路上的话，大致就是你要确定好你的目标（想要得到什么样的分割），然后得到一个目标函数，然后尝试把这个目标函数（一般是离散优化没法解）relax 为特征值问题阿之类的来进行近似。
博主你好，我的邻接矩阵用的是exp（-（距离平方）），但是最后算出来的结果不对，请问除了相似矩阵之外你用的是什么公式吗？
ym 学长，写的很不错
博主，您好！我有两个问题：1.照你这么说，特征向量构成的U矩阵中的第i行就对应着原始数据xi？第i行属于哪类，xi就属于哪个类？2.这里说的是取L的前k个最小特征值对应的特征向量，但为什么很对论文中和百度百科里面说的都是去前k个最大特征值对应的特征向量呢？
多谢pluskid的文章。请问一下有没有什么好的方法能自动统计cluster数k的?也就是说如果在聚类前并不知道最为合理的类别数为多少，有没有什么好的办法？
总的来说没有什么 out-of-box 直接能 work 的方法，这个问题本身也不是 well defined 的。:-/
那有没有什么好的办法能估计cluster数目的呢，或者需要iteratively的不断试验？
也可以试试不需要指定类别数目的算法，比如 Mean Shift，Affinity Propagation 之类的。
这其实可以从物理上的简正模式来理解，特征矢量表示的是每个点的运动方向，前k个矢量对应振动最大的k个模式，这样就比从图来理解直观多了
如果按照简正运动模式的想法来分析，做K-means之前还应该对各个特征矢量按其本征值进行加权，乘以eigenvalue^(-0.5)，因为前面的维度的区分性要高于后面的维度
我想请教一下，为什么求L的前k个本征矢要用[V dummy] = eigs(L, D, k, ‘SM’, opt)，而不是[V dummy]=eigs(L,k,’sm’)？好像也不是为了解决奇异阵的问题啊……
@pluskid 博主，你能回答下我的问题么？因为我到现在还是不明白到底Spectral Clustering中是取前k个最小特征值还是最大特征值的特征向量！谢谢博主！
如果用 Laplacian Matrix 就是取小特征值，如果用 similarity matrix 就是取最大的。
@pluskid 非常谢谢您的回复！像A为Gaussian similarity function构造的亲合矩阵，那L=D-1/2AD1/2 (不好意思，1/2是上标) 是similarity matrix ？因为它取的是L的前k个最大特征值！
这样的情况还是需要取小的特征值的。
哦，知道了！谢谢！那像上面的A应该是similarity matrix，L是Laplacian Matrix，对吧？
是的。
请问楼主，有没有什么办法能比较好的解决data normalization的？也就是说如果feature space的不同dimension上的scale差别很大，如果直接用Euclidean distance的话，实际上对scale最大的那个dimension的bias会非常大，从而降低cluster的效果。有没有什么通用的办法能比较好的解决这个问题的呢？非常感谢。
直接 normalize 呀，比如说全都缩放成均值 0 标准差 1 的情况，有各种 normalization 的标准。如果你想尝试下 fancy 一些的方法的话，可以搜索下 distance metric learning 相关的东西。
多谢，怕的是fancy的办法效果不好。不过用min-max，或者z-score的话不知道效果如何。如果是classification就简单多了。
还有就是请问一下楼主，用欧式距离和cosine similarity哪个有没有什么指导原则一类大的东西啊？
分析得很到位，具体关于laplacian eigenmap方法，感觉我还还需要再看看论文，
请问楼主有什么办法可以处理非对称相似度矩阵的吗？也就是说如果计算出的相似度矩阵为非对称矩阵，或者说是有向图的normal cut，有什么比较好的办法吗？或者有没有其他比较好的clustering的办法？
想问一下博主，如何去评价聚类效果的好坏
没有 groundtruth 的话，比较难弄。
那如果有ground truth的话，就是知道哪些样本会是同一类，比如社交网络的聚类，知道哪些节点会在同一个圈子，那我实验聚类出来，有哪些方法可以评价聚类效果。有一种是用正确率来评价，就是把所有聚类结果中分到同一类的最多数样本加起来，然后除于总样本数来得到。好像有一种是用边的正确率来评价，把聚类结果中同一类的节点彼此用边相连，与ground truth同一类节点彼此相连的网络来比较，统计哪些边正确，我不知道是不是可以这样的。
有 ground truth 的话，可以用 accuracy 或者 mutual information 之类的来衡量，需要先把聚类的类别做一个映射到正确的类别上。
博主好，我用下面的代码遇到一些问题：
function idx = spectral_cluster(W, k)
D = diag(sum(W));
L = D – W;
%opt = struct(‘issym’, true, ‘isreal’, true);
[V, dummy] = eigs(L, D, k);
[idx] = kmeans(V, k, ’emptyaction’, ‘singleton’,’replicates’, 20);
Unable to compute the specified eigenvalues because infinite eigenvalue(s) exist
想不明白是怎么回事，请问你知不知道呀
Laplacian Eigenmap 的降维方式降维和PCA降维有什么区别和联系吗？
感觉
Laplacian Eigenmap就是把协方差矩阵中方差小的特征拎了出来。组成子空间，然后向子空间进行投影？
PCA则是把协方差矩阵中方差最大的特征拎了出来，组成子空间，然后向子空间进行投影？
差别很大，LE 用的不是协方差矩阵。
博主，在计算谱聚类的时候应该是选取拉普拉斯矩阵矩阵从次小特征值开始的k个特征值对应的特征向量。应为(1,1,1,…)’一定是L的特征向量，刚好特征值是0所以L的最小特征值是0，选取的特征向量应该剔除0所对应的特征向量。
你说得没错，所以一般 k 类聚类实际上是取了 k-1 个特征向量，因为丢掉了 constant 的那个。但是也有些描述是没有专门去丢掉那个的，就是留在那里也没有什么影响，反正是 constant 的。
kid你好, 是不是可以理解为component个数等于特征值为0的个数等于k, 而最后计算的时候取的是除0以外的k个次小特征值? 此外, 不知道博主对于similarity graph的构造方法有没有深入研究, 如果选择构造全连通图的话, 除了Gaussian function外还有没有一些通用的函数呢?
Hi, 唔，我对 similarity graph 构造也没有什么特别的研究，我见过的一般主要就是三种：1. 全连通 Gaussian decay；2. k 近邻图；3. ε-近邻图。
其实LE方法包括了pca的情况，这里主要体现在拉普拉斯矩阵表达上面。问题最后归纳为广义特征值的求解。具体见你导师05年pami上的一篇文章
请问一下，在介绍Spectral Clustering的步骤里面，取L的前k个最小特征值对应的特征向量，这里的k是指要聚类的类数量吗？还是自定义一个值？比如要聚三类的话，只要只考虑前三的特征值？
可以这么取，看成是由 2 类问题到多类问题的直接推广。不过实际中一般会把它当做一个参数来调，选择效果最好的一个维度。
Vol(A) 的定义是不是不够严谨？没有写j的范围。还是可以简写成那样？
学长，L 的最小特征值为0，且对应的特征向量为1 这个能否说明下怎么得出来的 -.-
根据 L=D-W
你好 最近在看在这方面的东西，想问下那个凸球形的样本空间是个什么意思，谢谢
什么凸球形的样本空间？
博主你好，看了你的博客，感觉写的很棒。觉得下面这些地方写的有点问题：
1. f’Lf=2|V|RatioCut(A,\bar{A})这个地方，前面的常数系数中，没有2呀。推倒的时候，公式中少了一个1/2。所以最终出来的关系式错了。
2.“到此为止，已经有 Spectral Clustering 的影子了：求特征值，再对特征向量进行 K-means 聚类。”感觉这句写的也有点问题，可能会造成误解。当然如果理解了文中的算法的话，也不是问题了。O(∩_∩)O哈哈~
阅读了这篇文章，对降维也有了新认识，像LE这样的降维，不仅考虑原始空间和降维后的空间的关系，而且它在降维后的空间中保持原始数据空间中的局部性（通过权值调整）。像LE的目标函数，在基于图的算法中也是经常用（比如基于图的半监督学习）。
你好，非常感谢指出错误，那个地方确实系数少了 1/2，现在已经修正过来了！
博主你好，请教一下关于LE降维的问题。
LE降维中，给定数据矩阵X和相似度矩阵W，就能得到降维后的数据Y。
但是LE和PCA等线性降维不同，PCA能得到一个映射矩阵，把X中的每一个样本映射到低维空间组成Y。
对于PCA来说，再来一个新数据，只需要直接乘上映射矩阵就可以降维了；而在LE中，如果再来一个新的数据应该怎么处理呢？
可以参考 out of sample extension for lle isomap… by Y. Bengio 和里面的参考文献
文章中有提到Spectral Clustering和kmeans的比较
有一个地方说的是Spectral Clustering的计算复杂度，我理解的就是时间复杂度要比kmeans小，我从一个地方看到kmeans的时间复杂度是nkp，那么Spectral Clustering的时间复杂度是多少？我用graphlab跑了10000个数据点，实际情况却是kmeans完胜Spectral Clustering（快了将近两个数量级！），请问楼主能解答一下吗
主要的时间花在了用svd降维处理上！！
恩，分解拉普拉斯矩阵的复杂度是随着数据点的数目增长的。
师兄，您好！我搞不清楚similarity matrix W是怎么得到的呢？我今天看一篇文章‘Recursive discriminant regression analysis to find homogeneous groups’，但是我不知道理解的对不对，W是 Gram matrix么？
谢谢！
你好，这里的 W 一般是通过距离计算 K 近邻构造出来的，并不是 gram matrix。
请问大侠，Spectral Clustering之后，此方法能否得到同一类的特征谱。。具体比如得到A类在x坐标的哪几个点都会有震动
我是初学者，最近在准备搞这个，冒昧再问一句，原始文件是很多xy坐标，构成N张波普图，对应到分类的算法里面，是怎样的数据结构？
xyxyxyxyxy第一个波谱
xyxyxyxyxy第二个波谱
？是这样吗
学长，请问下，对于只有两类的Spectrum CLuster, 最小化目标函数 f’Lf也就是最小化RatioCut的值，可是为什么当f是L的第二个特征向量时(第一个特征值是0)，该目标函数就是最小的？
实际上特征值零所对应的 constant 特征向量才是最优的，只是我们不想去特征值 0，因为所有的东西都被 map 到一起了就没有太大意义了，所以取了第二个……
恩恩 懂了~谢谢学长
学长好，我在用谱聚类算法进行文本聚类的时候，实验中用最小的前k个特征值对应的特征向量进行聚类的效果非常差，查到几乎没有分类能力了，而用最大的前k个效果非常好，这是怎么回事呢？请指点！万分感谢~
如果你用的是 similarity 矩阵而不是 Laplacian 矩阵的话这是正常的行为。
可是，我确实是用了Laplacian变换(L = D^(-1/2) .* similarity .* D^(-1/2))对similarity矩阵进行了处理，还有什么细节会影响结果吗？请学长指点！
我把laplacian变换公式，改成(L = D^(-1/2) .* (D-similarity) .* D^(-1/2))和L = D-similarity，结果还是一样的。
我有个建议就是将你这个漫谈系列，做成PDF发布，这样方便转载
[…] 谱聚类的基本概念可以参考pluskid的文章： 漫谈 Clustering (4): Spectral Clustering […]
[…] Laplacian Eigenmaps看问题的角度和LLE十分相似。它们都用图的角度去构建数据之间的关系。图中的每个顶点代表一个数据，每一条边权重代表数据之间的相似程度，越相似则权值越大。并且它们还都假设数据具有局部结构性质。LE假设每一点只与它距离最近的一些点相似，再远一些的数据相似程度为0，降维后相近的点尽可能保持相近。而LLE假设每一个点都能通过周围邻域数据的线性组合来描述，并且降维后这一线性关系尽可能保持不变。不过这里我不主要介绍LLE，主要介绍LE，因为它在spectral clustering中被使用。 […]
hi~觉得这篇文章非常有用，谢谢！
但是我有一个问题，看了下前面有人提过了，但是学长没有解答。
为什么在MATLAB code里，要计算L和D的generalised eigenvector呢，我按照一开始给的步骤，计算laplacian的特征值和特征向量，与你的计算结果不符…不知道特征值与generalised特征值有什么不一样，已经晕了…可否解答一下？
想弱弱的问一下大神，开头说SC跟kmeans相比的优点，“
和 K-medoids 类似，Spectral Clustering 只需要数据之间的相似度矩阵就可以了，而不必像 K-means 那样要求数据必须是 N 维欧氏空间中的向量。”
可是相似度矩阵一般不还是由各个observation的n维向量互相求欧式距离得到的吗(有的还要加上阈值判定)，这样子看来SC还是不能避免数据要求有n维向量呀。数据直接能给出相似度的矩阵，不多吧？
另外还有一个问题想请教的，PCA+Kmeans 跟SC比起来，实现好像差不太多呀，您觉得呢？
你好世界
很多相似度不是靠欧式距离计算来的，比如时间序列(time series,ts)，你可以用DTW计算两个ts的相似度，但很难计算两个ts的平均值，因为它们极可能是不等长的。
这篇文章于认识谱聚类的原理特别有用，思路很清晰，感谢！
很有帮助，谢谢作者！
文章很精彩，感谢！
有一个地方没太懂，想再请教下：在证明f垂直于1那里，是因为A和A补的点数相同么，所以刚好加和后结果为0？
这个在 A Tutorial on Spectral Clustering 里有更详细的推倒，你可以看看