1、从拉普拉斯矩阵说到谱聚类Spectral Hashing作者：July出处：结构之法算法之道blog目录2目录0引言31矩阵基础31.1理解矩阵的12点数学笔记. . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2一堆基础概念. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52拉普拉斯矩阵62.1 Laplacian matrix的定义. . . . . . . . . . . . . . . . . . . . . . . . . . . .
2、62.2拉普拉斯矩阵的性质. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83谱聚类93.1相关定义. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93.2目标函数. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103.3最小化RatioCut与最小化fLf等价. . . . . . . . . . . . . . . . . . .
3、 . 113.4谱聚类算法过程. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130引言30引言11月1日上午，机器学习班1第7次课，邹博讲聚类（PPT），其中的谱聚类引起了自己的兴趣，他从最基本的概念：单位向量、两个向量的正交、方阵的特征值和特征向量，讲到相似度图、拉普拉斯矩阵，最后讲谱聚类的目标函数和其算法流程。课后自己又琢磨了番谱聚类跟拉普拉斯矩阵，打算写篇博客记录学习心得，若有不足或建议，欢迎随时不吝指出，thanks。1矩阵基础在讲谱聚类之前，有必要了解一些矩阵方面的基础知识。1.1理解矩阵的12点数
4、学笔记如果对矩阵的概念已经模糊，推荐国内一人写的理解矩阵by孟岩系列，其中，抛出了很多有趣的观点，我之前在阅读的过程中做了些笔记，如下：1.简而言之：矩阵是线性空间里的变换的描述，相似矩阵则是对同一个线性变换的不同描述。那，何谓空间？本质而言，空间是容纳运动的一个对象集合，而变换则规定了对应空间的运动by孟岩在线性空间选定基后，向量刻画对象的运动，运动则通过矩阵与向量相乘来施加。然，到底什么是基？坐标系也。2.有了基，那么在(1)中所言的则应是：矩阵是线性空间里的变换的描述，相似矩阵则是对同一个线性变换在不同基（坐标系）下的不同描述。出来了两个问题，一者何谓变换，二者不同基（坐标系）如何理解？
5、事实上，所谓变换，即空间里从一个点（元素/对象）到另一个（元素对象）的跃迁，矩阵用来描述线性变换。基呢?通过前面已知，矩阵无非不过就是用来描述线性空间中的线性变换的一个东西而已，线性变换为名词，矩阵为描述它的形容词，正如描述同一个人长得好看可以用多个不同形容词“帅”“靓”描述，同一个线性变换也可以由多个不同的矩阵来描述，而由哪一个矩阵描述它，则由基（坐标系）确定。3.前面说了基，坐标系也，形象表述则为角度，看一个问题的角度不同，描述问题得到的结论也不同，但结论不代表问题本身，同理，对于一个线性变换，可以选1http:/ Ma = b，坐标点移动则是向量a经过矩阵M所描述的变换，变成了向量b；变
6、坐标系则是有一个向量，它在坐标系M的度量下结果为a，在坐标系I（I为单位矩阵，主对角为1，其它为0）的度量下结果为b，本质上点运动与变换坐标系两者等价。为何？如(5)所述，同一个变换，不同坐标系下表现不同矩阵，但本质相同。7. Ib，I在(6)中说为单位坐标系，其实就是我们常说的直角坐标系，如Ma = Ib，在M坐标系里是向量a，在I坐标系里是向量b，本质上就是同一个向量，故此谓矩阵乘法计算无异于身份识别。且慢，什么是向量？放在坐标系中度量，后把度量的结果（向量在各个坐标轴上投影值）按顺序排列在一起，即成向量。8. b在I坐标系中则是Ib，a在M坐标系中则是Ma，故而矩阵乘法M N，不过是N在
7、M坐标系中度量得到MN，而M本身在I坐标系中度量出。故Ma = Ib，M坐标系中的a转过来在I坐标系中一量，却成了b。如向量(x;y)在单位长度均为1的直角坐标系中一量，是(1;1)，而在X轴单位长度为2；Y轴单位长度为3的量则是(2;3)。9.何谓逆矩阵? Ma = Ib，之前已明了坐标点变换a ! b等价于坐标系变换M ! I，但具体M如何变为I呢，答曰：让M乘以M的逆矩阵。以坐标系:(2 00 3)1矩阵基础5为例，X轴单位度量长度变为原来的1/2，Y轴单位度量长度变为原来的1/3，即与矩阵(12 00 13)相乘，便成直角坐标系I。即对坐标系施加变换，即让其与变换矩阵相乘。1.2一堆基
8、础概念根据wikipedia的介绍2，在矩阵中，n阶单位矩阵，是一个的方形矩阵，其主对角线元素为1，其余元素为0。单位矩阵以In表示；如果阶数可忽略，或可由前后文确定的话，也可简记为I（或者E）。如下所示，便是一些单位矩阵：I1 = (1);I2 =(1 00 1);I3 =0BB1 0 00 1 00 0 11CCA;:;In =0BBBBB1 0 : 00 1 : 0. . . .0 0 : 11CCCCCA (1.1)单位矩阵中的第i列即为单位向量ei。单位向量同时也是单位矩阵的特征向量，特征值皆为1，因此这是唯一的特征值，且具有重数n。由此可见，单位矩阵的行列式为1，且迹数为n。单位向
9、量又是什么呢？数学上，赋范向量空间中的单位向量就是长度为1的向量。欧几里得空间中，两个单位向量的点积就是它们之间角度的余弦（因为它们的长度都是1）。一个非零向量u的正规化向量（即单位向量）u就是平行于的单位向量，记作：u = uu (1.2)这里u是u的范数（长度）。何谓点积？点积又称内积，两个向量a = a1;a2;:;an和b = b1;b2;:;bn的点积定义为：a b =ni=1aibi = a1b1 + a2b2 + : + anbn (1.3)这里的指示求和符号。例如，两个三维向量1;3; 5和4; 2; 1的点积是：1;3; 5 4; 2; 1 = (1)(4) + (3)( 2
10、) + ( 5)(1) = 3 (1.4)2http:/en.wikipedia.org/wiki/Matrix_(mathematics)2拉普拉斯矩阵6使用矩阵乘法并把（纵列）向量当作n 1矩阵，点积还可以写为：a b = aTb (1.5)这里aT的指示矩阵的转置。使用上面的例子，将一个1 3矩阵（就是行向量）乘以一个3 1向量得到结果(通过矩阵乘法的优势得到1 1矩阵也就是标量)：(1 3 5)0BB4211CCA= (3) (1.6)除了上面的代数定义外，点积还有另外一种定义：几何定义。在欧几里得空间中，点积可以直观地定义为：a b = jajjbjcos (1.7)这里jxj表示x
11、的模（长度）， 表示两个向量之间的角度。根据这个定义式可得：两个互相垂直的向量的点积总是零。若和都是单位向量（长度为1），它们的点积就是它们的夹角的余弦。正交是垂直这一直观概念的推广，若内积空间中两向量的内积（即点积）为0，则称它们是正交的，相当于这两向量垂直，换言之，如果能够定义向量间的夹角，则正交可以直观的理解为垂直。而正交矩阵（orthogonal matrix）是一个元素为实数，而且行与列皆为正交的单位向量的方块矩阵（方块矩阵，或简称方阵，是行数及列数皆相同的矩阵。）若数字 和非零向量v满足Av = v，则v为A的一个特征向量， 是其对应的特征值。换句话说，在v这个方向上，A做的事情无
12、非是把v沿其v的方向拉长/缩短了一点（而不是毫无规律的多维变换）， 则是表示沿着这个方向上拉伸了多少的比例。简言之，A对v做了手脚，使得向量v变长或变短了，但v本身的方向不变。矩阵的迹是矩阵的对角线元素之和，也是其个特征值之和。更多矩阵相关的概念可以查阅相关wikipedia，或矩阵分析与应用。2拉普拉斯矩阵2.1 Laplacian matrix的定义拉普拉斯矩阵（Laplacian matrix)，也称为基尔霍夫矩阵,是表示图的一种矩阵。给定一个有n个顶点的图G = (V;E)，其拉普拉斯矩阵被定义为:L = D W (2.1)2拉普拉斯矩阵7其中D为图的度矩阵，W为图的邻接矩阵。举个例子
13、。给定一个简单的图，如下：645321图1:一个简单的Graph把此“图”转换为邻接矩阵的形式，记为W：W =0BBBBBBBBBB0 1 0 0 1 01 0 1 0 1 00 1 0 1 0 00 0 1 0 1 11 1 0 1 0 00 0 0 1 0 01CCCCCCCCCCA(2.2)把W的每一列元素加起来得到N个数，然后把它们放在对角线上（其它地方都是零），组成一个N N的对角矩阵，记为度矩阵D，如下图所示：D =0BBBBBBBBBB2 0 0 0 0 00 3 0 0 0 00 0 2 0 0 00 0 0 3 0 00 0 0 0 3 00 0 0 0 0 11CCCCCC
14、CCCCA(2.3)根据拉普拉斯矩阵的定义L = D W，可得拉普拉斯矩阵L为：L =0BBBBBBBBBB2 1 0 0 1 01 3 1 0 1 00 1 2 1 0 00 0 1 3 1 11 1 0 1 3 00 0 0 1 0 11CCCCCCCCCCA(2.4)2拉普拉斯矩阵82.2拉普拉斯矩阵的性质介绍拉普拉斯矩阵的性质之前，首先定义两个概念，如下：1.对于邻接矩阵，定义图中A子图与B子图之间所有边的权值之和如下：W(A;B) =i2A;j2Bwij (2.5)其中，wij定义为节点i到节点j的权值，如果两个节点不是相连的，权值为零。2.与某结点邻接的所有边的权值和定义为该顶点的
15、度d，多个d形成一个度矩阵D（对角阵）di =nj=1wij (2.6)拉普拉斯矩阵L具有如下性质： L是对称半正定矩阵； L 1 = 0 1，即L的最小特征值是0，相应的特征向量是1。证明：L 1 = (D W) 1 = 0 = 0 1 (2.7)此外，别忘了，之前特征值和特征向量的定义：若数字 和非零向量v满足Av = v，则v为A的一个特征向量， 是其对应的特征值。 L有n个非负实特征值0 = 1 2 : n且对于任何一个属于实向量f 2 n，有以下式子成立：fLf = 12Ni;j=1wij(fi fj)2 (2.8)其中，L = D W，di =nj=1 wij，W(A;B) =i2
16、A;j2B wij。证明：fLf = fDf fWf =ni=1dif2i ni;j=1fifjwij= 12(ni=1dif2i 2ni;j=1fifjwij +nj=1djf2j ) = 12ni;j=1wij(fi fj)2(2.9)3谱聚类93谱聚类所谓聚类（Clustering），就是要把一堆样本合理地分成两份或者K份。从图论的角度来说，聚类的问题就相当于一个图的分割问题。即给定一个图G = (V;E)，顶点集V表示各个样本，带权的边表示各个样本之间的相似度，谱聚类的目的便是要找到一种合理的分割图的方法，使得分割后形成若干个子图，连接不同子图的边的权重（相似度）尽可能低，同子图内的边
17、的权重（相似度）尽可能高。物以类聚，人以群分，相似的在一块儿，不相似的彼此远离。至于如何把图的顶点集分割/切割为不相交的子图有多种办法，如：1. cut/Ratio Cut2. Normalized Cut3.不基于图，而是转换成SVD能解的问题目的是为了要让被割掉各边的权值和最小，因为被砍掉的边的权值和越小，代表被它们连接的子图之间的相似度越小，隔得越远，而相似度低的子图正好可以从中一刀切断。本文重点阐述上述的第一种方法，简单提一下第二种3，第三种本文不做解释，有兴趣的可以参考文末的参考文献条目13。3.1相关定义为了更好的把谱聚类问题转换为图论问题，定义如下概念（有些概念之前已定义，权当回
18、顾下）：无向图G = (V;E)，顶点集V表示各个样本，带权的边表示各个样本之间的相似度与某结点邻接的所有边的权值和定义为该顶点的度d，多个d形成一个度矩阵D（对角阵）di =nj=1wij (3.1)邻接矩阵，A子图与B子图之间所有边的权值之和定义如下：W(A;B) =i2A;j2Bwij (3.2)其中，定义wij为节点i到节点j的权值，如果两个节点不是相连的，权值为零。3整理者注：强烈建议读者阅读Jianbo Shi关于NCut的PAMI论文，相比RatioCut更有代表性。3谱聚类10相似度矩阵的定义。相似度矩阵由权值矩阵得到，实践中一般用高斯核函数（也称径向基函数核4）计算相似度，距
19、离越大，代表其相似度越小。s(xi;xj) = e xi xj2/(2 2) (3.3)子图A的指示向量如下：1A = (f1;:;fn)T 2 nfi = 1 if vi 2 Afi = 0 otherwise(3.4)3.2目标函数因此，如何切割图则成为问题的关键。换言之，如何切割才能得到最优的结果呢？举个例子，如果用一张图片中的所有像素来组成一个图，并把（比如，颜色和位置上）相似的节点连接起来，边上的权值表示相似程度，现在要把图片分割为几个区域（或若干个组），要求是分割所得的Cut值最小，相当于那些被切断的边的权值之和最小，而权重比较大的边没有被切断。因为只有这样，才能让比较相似的点被保
20、留在了同一个子图中，而彼此之间联系不大的点则被分割了开来。设A1;A2;:;Ak为图的几个子集（它们没有交集），为了让分割的Cut值最小，谱聚类便是要最小化下述目标函数：cut(A1;:;Ak) = 12ki=1W(Ai; Ai) (3.5)其中k表示分成k个组，Ai表示第i个组， Ai表示Ai的补集，W(Ai; Ai)表示第Ai组与第 Ai组之间的所有边的权重之和（换言之，如果要分成K个组，那么其代价就是进行分割时去掉的边的权值的总和）。为了让被切断边的权值之和最小，便是要让上述目标函数最小化。但很多时候，最小化cut通常会导致不好的分割。以分成2类为例，这个式子通常会将图分成了一个点和其余
21、的n 1个点。如图2所示，很明显，最小化的smallest cut不是最好的cut，反而把fA;B;C;Hg分为一边，fD;E;F;Gg分为一边很可能就是最好的cut。为了让每个类都有合理的大小，目标函数尽量让A1;A2;:;Ak足够大。改进后的目标函数为：RatioCut(A1;:;Ak) = 12ki=1W(Ai; Ai)jAij =ki=1cut(Ai; Ai)jAij (3.6)4整理者注：也就是我们常见的RBF Kernel3谱聚类11图2: Graph Cut示例其中jAj表示A组中包含的顶点数目，或：NCut(A1;:;Ak) = 12ki=1W(Ai; Ai)vol(Ai) =
22、ki=1cut(Ai; Ai)vol(Ai) (3.7)其中，vol(A) =i2A wij。3.3最小化RatioCut与最小化fLf等价下面，咱们来重点研究下RatioCut函数。目标函数：minA V RatioCut(A; A)定义向量f = (f1;:;fn)T 2 n，且：fi =8:jAj/jAj if vi 2 AjAj/ jAj if vi 2 A(3.8)根据之前得到的拉普拉斯矩阵矩阵的性质，已知：fLf = 12Ni;j=1wij(fi fj)2 (3.9)3谱聚类12现在把fi的定义式代入上式，我们将得到一个非常有趣的结论！推导过程如下：fLf = 12Ni;j=1wi
23、j(fi fj)2= 12i2A;j2 Awij0 jAjjAj +jAjjAj1A2+j2A;i2 Awij0 jAjjAj jAjjAj1A2= cut(A; A)( jAjjAj +jAjjAj + 2)= cut(A; A)( jAj + jAjjAj +jAj + jAjjAj)= jVj RatioCut(A; A)(3.10)是的，我们竟然从fLf推出了RatioCut，换句话说，拉普拉斯矩阵L和我们要优化的目标函数RatioCut有着密切的联系。更进一步说，因为jVj是一个常量，所以最小化RatioCut，等价于最小化fLf。同时，因单位向量1的各个元素全为1，所以直接展开可得
24、到约束条件：f2 =f2i = n且f1 =fi = 0，具体推导过程如下：ni=1fi =i2A jAjjAj i2 AjAjjAj = jAj jAjjAj jAjjAjjAj = 0f2 =ni=1f2i = jAjjAjjAj +jAjjAj jAj =jAj + jAj = n(3.11)最终我们新的目标函数可以由之前的minA V RatioCut(A; A)，写成：minf2nfLf subject tof ? 1;f = pn (3.12)其中：fi =8:jAj/jAj if vi 2 AjAj/ jAj if vi 2 A(3.13)且因f = pn，所以有：ff = n（
25、注：f是列向量的前提下，ff是一个值，实数值，ff是一个N N的矩阵）。继续推导前，再次提醒特征向量和特征值的定义：定义1 (特征值与特征向量)若数字 和非零向量v满足Av = v，则v为A的一个特征向量， 是其对应的特征值。假定Lf = f，此刻， 是特征值，f是L的特征向量。两边同时左乘f，得到fLf = ff，而ff = n，其中n为图中顶点的数量之和，因此fLf = n，因n是3谱聚类13个定值，所以要最小化fLf，相当于就是要最小化 。因此，接下来，我们只要找到L的最小特征值及其对应的特征向量即可。但到了这关键的最后一步，咱们却遇到了一个比较棘手的问题，即由之前得到的拉普拉斯矩阵的性
26、质“L最小的特征值为零，并且对应的特征向量正好为1”可知：其不满足的条件f ? 1，怎么办呢？根据论文“A Tutorial on Spectral Clustering”中所说的Rayleigh-Ritz理论，我们可以取第2小的特征值，以及对应的特征向量v。更进一步，由于实际中，特征向量v里的元素是连续的任意实数，所以可以根据v是大于0，还是小于0对应到离散情况下的f = (f1;:;fn)T 2 n，决定f是取jAj/jAj，还是取 jAj/ jAj。而如果能求取v的前K个特征向量，进行K-means聚类，得到K个簇，便从二聚类扩展到了K聚类的问题。而所要求的这前K个特征向量就是拉普拉斯矩
27、阵的特征向量（计算拉普拉斯矩阵的特征值，特征值按照从小到大顺序排序，特征值对应的特征向量也按照特征值递增的顺序排列，取前K个特征向量，便是我们所要求的前K个特征向量）！所以，问题就转换成了：求拉普拉斯矩阵的前K个特征值，再对前K个特征值对应的特征向量进行K-means聚类。而两类的问题也很容易推广到k类的问题，即求特征值并取前K个最小的，将对应的特征向量排列起来，再进行K-means聚类。两类分类和多类分类的问题，如出一辙。就这样，因为离散求解f = (f1;:;fn)T 2 n很困难，但RatioCut巧妙地把一个NP难度的问题转换成拉普拉斯矩阵特征值（向量）的问题，将离散的聚类问题松弛为连
28、续的特征向量，最小的系列特征向量对应着图最优的系列划分方法。剩下的仅是将松弛化的问题再离散化，即将特征向量再划分开，便可以得到相应的类别。不能不说妙哉！3.4谱聚类算法过程综上可得谱聚类的算法过程如下：1.根据数据构造一个Graph，Graph的每一个节点对应一个数据点，将各个点连接起来（随后将那些已经被连接起来但并不怎么相似的点，通过cut/RatioCut/NCut的方式剪开），并且边的权重用于表示数据之间的相似度。把这个Graph用邻接矩阵的形式表示出来，记为W。2.把W的每一列元素加起来得到N个数，把它们放在对角线上（其他地方都是零），组成一个N N的对角矩阵，记为度矩阵D，并把W D
29、的结果记为拉普拉斯矩阵L = D W。3谱聚类143.求出L的前k个特征值（前k个指按照特征值的大小从小到大排序得到）f gki=1，以及对应的特征向量fvgki=1。4.把这个特征（列）向量排列在一起组成一个N k的矩阵，将其中每一行看作k维空间中的一个向量，并使用K-means算法进行聚类。聚类的结果中每一行所属的类别就是原来Graph中的节点亦即最初的N个数据点分别所属的类别。或许你已经看出来，谱聚类的基本思想便是利用样本数据之间的相似矩阵（拉普拉斯矩阵）进行特征分解（通过Laplacian Eigenmap的降维方式降维），然后将得到的特征向量进行K-means聚类。此外，谱聚类和传统
30、的聚类方法（例如K-means）相比，谱聚类只需要数据之间的相似度矩阵就可以了，而不必像K-means那样要求数据必须是N维欧氏空间中的向量。3谱聚类15参考文献及推荐阅读1.孟岩之理解矩阵系列：http:/ wikipedia上关于拉普拉斯矩阵的介绍：http:/en.wikipedia.org/wiki/Laplacian_matrix；5.邹博之聚类PPT：http:/ Tutorial on Spectral Clustering”：http:/engr.case.edu/ray_soumya/mlrg/Luxburg07_tutorial_spectral_clustering.pd
31、f；7.知乎上关于矩阵和特征值的两个讨论：http:/ of Massive Datasets第10章：http:/infolab.stanford.edu/ullman/mmds/book.pdf；12. Tydsh: Spectral Clustering：http:/ H. Zha, C. Ding, M. Gu, X. He, and H.D. Simon. Spectral relaxation for K-meansclustering. Advances in Neural Information Processing Systems 14 (NIPS 2001).pp. 1057-1064, Vancouver, Canada. Dec. 2001；14.机器学习中谱聚类方法的研究：http:/
