谱聚类（spectral clustering)及其实现详解

Preface
开了很多题，手稿都是写好一直思考如何放到CSDN上来，一方面由于公司技术隐私，一方面由于面向对象不同，要大改，所以一直没贴出完整，希望日后可以把开的题都补充全。

---

先把大纲列出来：

一、从狄多公主圈地传说说起
二、谱聚类的演算
  （一）、演算
      1、谱聚类的概览
      2、谱聚类构图
      3、谱聚类切图
        （1）、RatioCut
        （2）、Ncut
        （3）、一点题外话
  （二）、pseudo-code
三、谱聚类的实现(scala)
  （一）、Similarity Matrix
  （二）、kNN/mutual kNN
  （三）、Laplacian Matrix
  （四）、Normalized
  （五）、Eigenvector(Jacobi methond)
  （六）、kmeans/GMM
四、一些参考文献

---

一、从狄多公主圈地传说说起

       谱聚类（spectral clustering）的思想最早可以追溯到一个古老的希腊传说，话说当时有一个公主，由于其父王去世后，长兄上位，想独揽大权，便杀害了她的丈夫，而为逃命，公主来到了一个部落，想与当地的酋长买一块地，于是将身上的金银财宝与酋长换了一块牛皮，且与酋长约定只要这块牛皮所占之地即可。聪明的酋长觉得这买卖可行，于是乎便答应了。殊不知，公主把牛皮撕成一条条，沿着海岸线，足足围出了一个城市。
       故事到这里就结束了，但是我们要说的才刚刚开始，狄多公主圈地传说，是目前知道的最早涉及Isoperimetric problem（等周长问题）的，具体为如何在给定长度的线条下围出一个最大的面积，也可理解为，在给定面积下如何使用更短的线条，而这，也正是谱图聚类想法的端倪，如何在给定一张图，拿出“更短”的边来将其“更好”地切分。而这个“更短”的边，正是对应了spectral clustering中的极小化问题，“更好”地切分，则是对应了spectral clustering中的簇聚类效果。
       谱聚类最早于1973年被提出，当时Donath 和 Hoffman第一次提出利用特征向量来解决谱聚类中的f向量选取问题，而同年，Fieder发现利用倒数第二小的特征向量，显然更加符合f向量的选取，同比之下，Fieder当时发表的东西更受大家认可，因为其很好地解决了谱聚类极小化问题里的NP-hard问题，这是不可估量的成就，虽然后来有研究发现，这种方法带来的误差，也是无法估量的，下图是Fielder老爷子，于去年15年离世，缅怀。

二、谱聚类的演算

（一）、演算

1、谱聚类概览

       谱聚类演化于图论，后由于其表现出优秀的性能被广泛应用于聚类中，对比其他无监督聚类（如kmeans），spectral clustering的优点主要有以下：

1.过程对数据结构并没有太多的假设要求，如kmeans则要求数据为凸集。
2.可以通过构造稀疏similarity graph，使得对于更大的数据集表现出明显优于其他算法的计算速度。
3.由于spectral clustering是对图切割处理，不会存在像kmesns聚类时将离散的小簇聚合在一起的情况。
4.无需像GMM一样对数据的概率分布做假设。

       同样，spectral clustering也有自己的缺点，主要存在于构图步骤，有如下：

1.对于选择不同的similarity graph比较敏感（如 epsilon-neighborhood， k-nearest neighborhood，fully connected等）。
2.对于参数的选择也比较敏感（如 epsilon-neighborhood的epsilon，k-nearest neighborhood的k，fully connected的 ）。

       谱聚类过程主要有两步，第一步是构图，将采样点数据构造成一张网图，表示为G(V,E)，V表示图中的点，E表示点与点之间的边，如下图：
              
                            图1 谱聚类构图(来源wiki)
       第二步是切图，即将第一步构造出来的按照一定的切边准则，切分成不同的图，而不同的子图，即我们对应的聚类结果，举例如下：
              
                            图2 谱聚类切图
       初看似乎并不难，但是…，下面详细说明推导。

2、谱聚类构图

       在构图中，一般有三种构图方式：
       1. $\varepsilon$-neighborhood
       2. k-nearest neighborhood
       3. fully connected
       前两种可以构造出稀疏矩阵，适合大样本的项目，第三种则相反，在大样本中其迭代速度会受到影响制约，在讲解三种构图方式前，需要引入similarity function，即计算两个样本点的距离，一般用欧氏距离：$s_{i,j} = \left \| x_{i}-x_{j} \right \| ^{2}$, $s_{i,j}$表示样本点$x_{i}$与$x_{j}$的距离，或者使用高斯距离$s_{i,j}=e^{\frac{-\left \| x_{i}-x_{j} \right \|^{2}}{2\sigma ^{2}}}$，其中$\sigma$ 的选取也是对结果有一定影响，其表示为数据分布的分散程度，通过上述两种方式之一即可初步构造矩阵$S:S_{i,j}=[s]_{i,j}$，一般称 为Similarity matrix(相似矩阵)。
       对于第一种构图$\varepsilon$ -neighborhood，顾名思义是取$s_{i,j}\leq \varepsilon$的点，则相似矩阵$S$可以进一步重构为邻接矩阵(adjacency matrix)$W$ :

$$W_{i,j}=\left\{\begin{matrix} 0,\quad if \quad s_{i,j}>\varepsilon\\ \varepsilon,\quad if \quad s_{i,j}\leq \varepsilon \end{matrix}\right.$$
       可以看出，在$\varepsilon$-neighborhood重构下，样本点之间的权重没有包含更多的信息了。
       对于第二种构图k-nearest neighborhood，其利用KNN算法，遍历所有的样本点，取每个样本最近的k个点作为近邻，但是这种方法会造成重构之后的邻接矩阵$W$非对称，为克服这种问题，一般采取下面两种方法之一：
       一是只要点$x_{i}$在$x_{j}$的K个近邻中或者$x_{j}$在$x_{i}$的K个近邻中，则保留$s_{i,j}$，并对其做进一步处理$W$，此时 为：
$$W_{i,j}=W_{j,i}=\left\{\begin{matrix} 0,\quad if \quad x_{i}\notin KNN(X_{j}) \quad and \quad x_{j}\notin KNN(X_{i})\\ e^{\frac{-\left \| x_{i}-x_{j} \right \|^{2}}{2\sigma ^{2}}}, \quad if \quad x_{i}\in KNN(X_{j}) \quad or \quad x_{j}\in KNN(X_{i}) \end{matrix}\right.$$
       二是必须满足点 在 的K个近邻中且 在 的K个近邻中，才会保留$s_{i,j}$并做进一步变换，此时$W$为：
$$W_{i,j}=W_{j,i}=\left\{\begin{matrix} 0,\quad if \quad x_{i}\notin KNN(X_{j}) \quad or \quad x_{j}\notin KNN(X_{i})\\ e^{\frac{-\left \| x_{i}-x_{j} \right \|^{2}}{2\sigma ^{2}}}, \quad if \quad x_{i}\in KNN(X_{j}) \quad and \quad x_{j}\in KNN(X_{i}) \end{matrix}\right.$$
       对于第三种构图fully connected，一般使用高斯距离： $s_{i,j}=e^{\frac{-\left \| x_{i}-x_{j} \right \|^{2}}{2\sigma ^{2}}}$，则重构之后的矩阵$W$与之前的相似矩阵$S$相同，为：$W_{i,j}= S_{i,j}=[s]_{i,j}$。
       在了解三种构图方式后，还需要注意一些细节，对于第一二中构图，一般是重构基于欧氏距离的 ，而第三种构图方式，则是基于高斯距离的 ，注意到高斯距离的计算蕴含了这样一个情况：对于$\left \| x_{i}-x_{j} \right \|^{2}$比较大的样本点，其得到的高斯距离反而值是比较小的，而这也正是$S$可以直接作为$W$的原因，主要是为了将距离近的点的边赋予高权重。
       得到邻接矩阵$W$后，需要做进一步的处理：
       (1).计算阶矩(degree matrix)$D$：
$$D_{i,j}=\left\{\begin{matrix} 0,\quad if \quad i \neq j\\ \sum_{j}w_{i,j},\quad if \quad i = j \end{matrix}\right.$$
       其中其中$w_{i,j}$为邻接矩阵$W$元素，$\sum_{j}w_{i,j}$表示将图中某点所连接的其他点的边的权重之和，可以看出，$D$为对角矩阵.
       (2).计算拉普拉斯矩阵(Laplacians matrix)$L$：$L=D-W$
       如此，在构图阶段基本就完成了，至于为什么要计算出拉普拉斯矩阵$L$，可以说$L=D-W$这种形式对于后面极大化问题是非常有利的(不知道前人是怎么想到的，不然$L$也不会被命名为拉普拉斯了吧)。

3、谱聚类切图

       谱聚类切图存在两种主流的方式：RatioCut和Ncut，目的是找到一条权重最小，又能平衡切出子图大小的边，下面详细说明这两种切法。
       在讲解RatioCut和Ncut之前，有必要说明一下问题背景和一些概念，假设V为所有样本点的集合，$\{A_{1},A_{2},\cdots ,A_{k}\}$表示V的子集集合，其中$A_{1} \cup A_{2} \cup \cdots \cup A_{k}=V 且A_{1} \cap A_{2}\cap \cdots\cap A_{k}= \varnothing$，则子集与子集之间连边的权重和为：

$$cut(A_{1},A_{2},\cdots,A_{k})=\frac{1}{2}\sum_{i}^{k}W(A_{i},\bar{A_{i}})$$
       其中$\bar{A_{i}}$为$A_{i}$的补集，意为除$A_{i}$子集外其他V的子集的并集，$W(A_{i},\bar{A_{i}})$为$A_{i}$与其他子集的连边的和，为：
$$W(A_{i},\bar{A_{i}})=\sum_{m\in A_{i},n \in \bar{A_{i}}}w_{m,n}$$
       其中$w_{m,n}$为邻接矩阵W中的元素。
       由于我们切图的目的是使得每个子图内部结构相似，这个相似表现为连边的权重平均都较大，且互相连接，而每个子图间则尽量没有边相连，或者连边的权重很低，那么我们的目的可以表述为：
$$min\quad cut(A_{1},A_{2},\cdots,A_{k})$$
       但是稍微留意可以发现，这种极小化的切图存在问题，如下图：
              
                            图3 问题切图
       如上图，如果用$min\quad cut(A_{1},A_{2},\cdots,A_{k})$这种方法切图，会造成将V切成很多个单点离散的图，显然这样是最快且最能满足那个最小化操作的，但这明显不是想要的结果，所以有了RatioCut和Ncut，具体地：
$$Ratiocut(A_{1},A_{2},\cdots,A_{k})=\frac{1}{2}\sum_{i}^{k}\frac{W(A_{i},\bar{A_{i}})}{\left | A_{i} \right |}$$
$$Ncut(A_{1},A_{2},\cdots,A_{k})=\frac{1}{2}\sum_{i}^{k}\frac{W(A_{i},\bar{A_{i}})}{vol(A_{i})}$$
       其中$\left | A_{i} \right |$ 为$A_{i}$中点的个数，$vol(A_{i})$为$A_{i}$中所有边的权重和。
       下面详细讲解这两种方法的计算：Ratiocut & Ncut

       (1).Ratiocut

       Ratiocut切图考虑了目标子图的大小，避免了单个样本点作为一个簇的情况发生，平衡了各个子图的大小。Ratiocut的目标同样是极小化各子图连边和，如下：

$$min\quad Ratiocut(A_{1},A_{2},\cdots,A_{k})$$
       对上述极小化问题做一下转化，引入$\{A_{1},A_{2},\cdots ,A_{k}\}$的指示向量$h_{j}=\{h_{1},h_{2},\cdots,h_{i},\cdots ,h_{k}\},\quad j=1,2,\cdots,k$. 其中i表示样本下标，j表示子集下标， 表示样本i对子集j的指示，具体为：
$$h_{j,i}=\left\{\begin{matrix} \frac{1}{\sqrt{\left | A_{j} \right |}},\quad if \quad v_{i}\in A_{j}\\ 0, \quad if \quad v_{i} \notin A_{j} \end{matrix}\right.$$
       通俗理解就是，每个子集$A_{j}$对应一个指示向量$h_{j}$，而每个$h_{j}$里有N个元素，分别代表N个样本点的指示结果，如果在原始数据中第i个样本被分割到子集$A_{j}$里，则$h_{j}$的第i个元素为$\frac{1}{\sqrt{\left | A_{j} \right |}}$，否则为0。
       进一步，计算 ，可以得到：
$$\begin{align} h_{i}^{T}Lh_{i} & = h_{i}^{T}(D-W)h_{i} \\ 　　& =h_{i}^{T}Dh_{i} - h_{i}^{T}Wh_{i}\　　& = \sum_{m=1}\sum_{n=1}h_{i,m}h_{i,n}D_{m,n} - \sum_{m=1}\sum_{n=1}h_{i,m}h_{i,n}w_{m,n}\　　& = \sum_{m=1}h_{i,m}^{2}D_{m,m} - \sum_{m=1}\sum_{n=1}h_{i,m}h_{i,n}w_{m,n}\　　& = \frac{1}{2}(\sum_{m=1}h_{i,m}^{2}D_{m,m} -2\sum_{m=1}\sum_{n=1}h_{i,m}h_{i,n}w_{m,n} +\sum_{n=1}h_{i,n}^{2}D_{n,n} ) 　　\end{align}$$
       上述过程稍微复杂，分解如下：
       第一步，第二步转化，利用上文拉普拉斯矩阵定义L=D-W；
       第三步，m,n分别表示第m,n个样本，其最大值均为N，D/W均为对称矩阵，此步主要将矩阵内积离散化为元素连加求和形式；
       第四步，由于D为对角矩阵，所以只有对角线上元素与向量$h_{i}$相乘有意义；
       第五步，代数变换。
       进一步地，做如下：
$$\begin{align} h_{i}^{T}Lh_{i} & = \frac{1}{2}(\sum_{m=1}h_{i,m}^{2}D_{m,m} -2\sum_{m=1}\sum_{n=1}h_{i,m}h_{i,n}w_{m,n} +\sum_{n=1}h_{i,n}^{2}D_{n,n} ) \\ 　　& = \frac{1}{2}(\sum_{m=1}\sum_{n=1}h_{i,m}^{2}w_{m,n} -2\sum_{m=1}\sum_{n=1}h_{i,m}h_{i,n}w_{m,n} +\sum_{n=1}\sum_{m=1}h_{i,n}^{2}w_{n,m} )\　　& = \frac{1}{2}\sum_{m=1}\sum_{n=1}w_{n,m}(h_{i,m}-h_{i,n})^{2} 　　\end{align}$$
       同样分解如下：
       第一步到第二步，由上文阶矩的定义可知，其对角线上元素$D_{m,m}=\sum_{n=1}w_{m,n}$
       第二部到第三步，二次函数变换；
       进一步，将上文$h_{j,i}$代入上式，最终可以得到：
$$\begin{align} h_{i}^{T}Lh_{i} & = \frac{1}{2}\sum_{m=1}\sum_{n=1}w_{m,n}(h_{i,m}-h_{i,n})^{2} \\ 　　& = \frac{1}{2}(\sum_{m\in A_{i},n \in \bar{A_{i}}}w_{m,n}(\frac{1}{\sqrt{\left | A_{i} \right |}}-0)^{2}+\sum_{m\in \bar{A_{i}},n \in A_{i}}w_{m,n}(0-\frac{1}{\sqrt{\left | A_{i} \right |}})^{2})\　　& = \frac{1}{2}(\sum_{m\in A_{i},n \in \bar{A_{i}}}w_{m,n}\frac{1}{\left | A_{i} \right |}+\sum_{m\in \bar{A_{i}},n \in A_{i}}w_{m,n}\frac{1}{\left | A_{i} \right |})\　　& = \frac{1}{2}(cut(A_{i},\bar{A_{i}})\frac{1}{\left | A_{i} \right |} + cut(\bar{A_{i}},A_{i})\frac{1}{\left | A_{i} \right |} )\　　& = \frac{cut(A_{i},\bar{A_{i}})}{\left | A_{i} \right |}\　　& = Ratiocut(A_{i},\bar{A_{i}}) 　　\end{align}$$
       分解如下：
       第一步到第二步，将指示变量$h_{i,m}$，$h_{i,n}$对应$A_{i}$的元素分别代入；
       第三，四，五步，分别做简单运算，以及将$\sum_{m\in A_{i},n \in \bar{A_{i}}}w_{m,n}$替换为$cut(A_{i},\bar{A_{i}})$；
       可以看到，通过引入指示变量h，将L矩阵放缩到与Ratiocut等价，这无疑是在谱聚类中划时代的一举。
       为了更进一步考虑进所有的指示向量，令$H=\{h_{1},h_{2},\cdots ,h_{k}\}$，其中$h_{i}$按列排列，由$h_{i}$的定义知每个$h_{i}$之间都是相互正交，即 $h_{i}*h_{j}=0,i \neq j$，且有$h_{i}*h_{i}=1$，可得到：
$$\begin{align} Ratiocut(A_{1},A_{2},\cdots,A_{k}) & =\sum_{i=1}^{k} h_{i}^{T}Lh_{i} \\ 　　& =\sum_{i=1}^{k} (H^{T}LH)_{ii}\　　& = Tr(H^{T}LH) 　　\end{align}$$
       Tr表示对角线求和。
       如此，便将极小化问题：$min\quad Ratiocut(A_{1},A_{2},\cdots,A_{k})$，转化为：
$$arg\mathop{min}\limits_{H} \quad Tr(H^{T}LH), \quad s.t. H^{T}H=I$$
       将之前切图的思路以一种严谨的数学表达表示出来，对于上面极小化问题，由于L矩阵是容易得到的，所以目标是求满足条件的H矩阵，以使得$Tr(H^{T}LH)$最小。
       注意到一点，要求条件下的H，首先是求得组成H的每个$h_{i}$，而每个$h_{i}$都是Nx1的向量，且向量中每个值都是二值分布，即取值0或者$1/\sqrt{\left | A_{i} \right |}$，那么对于每个元素，有2种选择，则当个$h_{i}$就有$2^{N}$种情况，对于整个H矩阵，则有$(C_{k}^{1})^{N}$种情况，对于N很大的数据，这无疑是灾难性的NP-hard问题，显然我们不可能遍历所有的情况来求解，那么问题是不是没有解？答案自然是否定的，虽然我们没办法求出精确的 ，但是我们可以用另外一种方法近似代替$h_{i}$，而这就是文章开头Fielder提出来的，当k=2的时候，可以用L的倒数第二小的特征向量(因为最小的特征向量为1，其对应的特征值为0，不适合用来求解)来代替 ，关于具体为什么可以这么做，可以参考Rayleigh-Ritz method(对应论文为：A short theory of the Rayleigh-Ritz method)，因为涉及泛函变分，这里就不多说了。
       那么当k取任意数字时，只需取L矩阵对应的最小那k个（毕竟倒数第二小只有1个），即可组成目标H，最后，对H做标准化处理，如下：
$$H_{i,j}^{*}=\frac{H_{i,j}}{(\sum_{j=1}^{k}H_{i,j}^{2})^{1/2}}$$
       现在H矩阵也完成了，剩下的就是对样本聚类了，要明白我们的目标不是求$Tr(H^{T}LH)$的最小值是多少，而是求能最小化$Tr(H^{T}LH)$的H，所以聚类的时候，分别对H中的行进行聚类即可，通常是kmeans，也可以是GMM，具体看效果而定。
       至于为什么是对H的行进行聚类，有两点原因：
       1.注意到H除了是能满足极小化条件的解，还是L的特征向量，也可以理解为W的特征向量，而W则是我们构造出的图，对该图的特征向量做聚类，一方面聚类时不会丢失原图太多信息，另一方面是降维加快计算速度，而且容易发现图背后的模式。
       2.由于之前定义的指示向量$h_{i}$是二值分布，但是由于NP-hard问题的存在导致$h_{i}$无法显式求解，只能利用特征向量进行近似逼近，但是特征向量是取任意值，结果是我们对$h_{i}$的二值分布限制进行放松，但这样一来$h_{i}$如何指示各样本的所属情况？所以kmeans就登场了，利用kmeans对该向量进行聚类，如果是k=2的情况，那么kmeans结果就与之前二值分布的想法相同了，所以kmeans的意义在此，k等于任意数值的情况做进一步类推即可。
       以上是Ratiocut的内容。

       (2).Ncut

       Ncut切法实际上与Ratiocut相似，但Ncut把Ratiocut的分母$|A_{i}|$换成$vol(A_{i})$ ，这种改变与之而来的，是L的normalized，这种特殊称谓会在下文说明，而且这种normalized，使得Ncut对于spectral clustering来说，其实更好，下文会说明。
       同样，Ncut的目标，也是极小化各子图连边的和，如下：

$$min\quad Ncut(A_{1},A_{2},\cdots,A_{k})$$
       下面对该问题做一下转化，先引入$\{A_{1},A_{2},\cdots ,A_{k}\}$的指示变量$h_{j}=\{h_{1},h_{2},\cdots,h_{i},\cdots ,h_{k}\},\quad j=1,2,\cdots,k$，同样，其中i表示样本下标，j表示子集下标，$h_{j,i}$表示样本i对子集j的指示，不同的是，其具体为：
$$h_{j,i}=\left\{\begin{matrix} \frac{1}{\sqrt{vol(A_{i})}},\quad if \quad v_{i}\in A_{j}\\ 0, \quad if \quad v_{i} \notin A_{j} \end{matrix}\right.$$
       如果在原始数据中第i个样本被分割到子集$A_{j}$里，则$h_{j}$的第i个元素为$1/\sqrt{vol(A_{i})}$，否则为0。
       进一步，计算$h_{i}^{T}Lh_{i}$，可以得到：
$$\begin{align} h_{i}^{T}Lh_{i} & =h_{i}^{T}Dh_{i} - h_{i}^{T}Wh_{i}\　　& = \sum_{m=1}\sum_{n=1}h_{i,m}h_{i,n}D_{m,n} - \sum_{m=1}\sum_{n=1}h_{i,m}h_{i,n}w_{m,n}\　　& = \sum_{m=1}h_{i,m}^{2}D_{m,m} - \sum_{m=1}\sum_{n=1}h_{i,m}h_{i,n}w_{m,n}\　　& = \frac{1}{2}(\sum_{m=1}h_{i,m}^{2}D_{m,m} -2\sum_{m=1}\sum_{n=1}h_{i,m}h_{i,n}w_{m,n} +\sum_{n=1}h_{i,n}^{2}D_{n,n} )\　　& = \frac{1}{2}(\sum_{m=1}\sum_{n=1}h_{i,m}^{2}w_{m,n} -2\sum_{m=1}\sum_{n=1}h_{i,m}h_{i,n}w_{m,n} +\sum_{n=1}\sum_{m=1}h_{i,n}^{2}w_{n,m} )\　　& = \frac{1}{2}\sum_{m=1}\sum_{n=1}w_{m,n}(h_{i,m}-h_{i,n})^{2}\　　& = \frac{1}{2}(\sum_{m\in A_{i},n \in \bar{A_{i}}}w_{m,n}(\frac{1}{\sqrt{vol( A_{i})}}-0)^{2}+\sum_{m\in \bar{A_{i}},n \in A_{i}}w_{m,n}(0-\frac{1}{\sqrt{vol( A_{i})}})^{2})\　　& = \frac{cut(A_{i},\bar{A_{i}})}{vol(A_{i})} \　　& = Ncut(A_{i},\bar{A_{i}}) 　　\end{align}$$
       由于推导与上文类似，这里就综合在一起，主要是为了说明$h_{i}^{T}Lh_{i}$如何得到$Ncut(A_{i},\bar{A_{i}})$，以便将$min\quad Ncut(A_{1},A_{2},\cdots,A_{k})$用严谨数学式子进行表达。
       这里对上述一大串式子分解一下：
       第一，二，三步，做一些简单的变换，如L=D-W代入，将向量内积转换成连加形式，同样，D是对角矩阵，可省去一些项；
       第四，五，六步，主要是做二次函数变换；
       第七步，将$h_{i,m}$， $h_{i,n}$对应$A_{i}$的元素分别代入；
       第八，九步，得到$Ncut(A_{i},\bar{A_{i}})$ 形式；
       此外，相比与Ratiocut，Ncut特有的一点性质是：$h_{i}^{T}Dh_{i}=1$，具体如下：
$$\begin{align} h_{i}^{T}Dh_{i} & =\sum_{n=1}h_{i,n}^{2}D_{n,n} \\ 　　& =\sum_{v_{n} \in A_{i}}\frac{1}{vol(A_{i})}D_{n,n} + \sum_{v_{n} \notin A_{i}}0\cdot D_{n,n} \　　& = \frac{1}{vol(A_{i})}\sum_{v_{n} \in A_{i}}w_{v_{n} } + 0 \　　& = \frac{1}{vol(A_{i})}vol(A_{i}) \　　& = 1 　　\end{align}$$
       进一步，令$H=\{h_{1},h_{2},\cdots ,h_{k}\}$，则$H^{T}DH=I$ ,其中I为单位对角矩阵。
       同样，可以得到：
$$\begin{align} Ncut(A_{1},A_{2},\cdots,A_{k}) & =\sum_{i=1}^{k} h_{i}^{T}Lh_{i} \\ 　　& =\sum_{i=1}^{k} (H^{T}LH)_{ii}\　　& = Tr(H^{T}LH) 　　\end{align}$$
       如此，便将极小化问题：$min\quad Ncut(A_{1},A_{2},\cdots,A_{k})$，转化为：
$$arg\mathop{min}\limits_{H} \quad Tr(H^{T}LH), \quad s.t. H^{T}DH=I$$
       除了H的限制条件，Ncut的极小化与Ratiocut的极小化基本无差异，但是就是这微小的条件，使得Ncut实质上是对L做了normalize（不得不叹数学真是十分奇妙）。
       下面进一步说明，令$H=D^{-1/2}F$， $D^{-1/2}$意为对D对角线上元素开方再求逆。将 $H=D^{-1/2}F$代入$arg\mathop{min}\limits_{H} \quad Tr(H^{T}LH), \quad s.t. H^{T}DH=I$，可以得到：
$$\begin{align} H^{T}LH & =(D^{-1/2}F)^{T}LD^{-1/2}F \\ 　　& = F^{T}D^{-1/2}LD^{-1/2}F 　　\end{align}$$
        以及：
$$\begin{align} H^{T}DH & =(D^{-1/2}F)^{T}DD^{-1/2}F \\ 　　& = F^{T}D^{-1/2}LD^{-1/2}F \　　& = F^{T}IF \　　& = I 　　\end{align}$$
        至此，目标操作$arg\mathop{min}\limits_{H} \quad Tr(H^{T}LH), \quad s.t. H^{T}DH=I$可以修改为：
$$arg\mathop{min}\limits_{H} \quad Tr(F^{T}D^{-1/2}LD^{-1/2}F), \quad s.t. F^{T}F=I$$
        其中，$D^{-1/2}LD^{-1/2}$这一步操作，就是对L矩阵进行normalize，而normalize可以理解为标准化，具体为：$\frac{L_{i,j}}{\sqrt{vol( A_{i})vol(A_{j})}}$ ，这么做的一个好处就是对L中的元素进行标准化处理使得不同元素的量纲得到归一，具体理解就是，同个子集里，不同样本点之间的连边可能大小比较相似，这点没问题，但是对于不同子集，样本点之间的连边大小可能会差异很大，做 这一步normalize操作，可以将L中的元素归一化在[-1,1]之间，这样量纲一致，对算法迭代速度，结果的精度都是有很大提升。
       最后，F矩阵的求解可通过求$D^{-1/2}LD^{-1/2}$的前k个特征向量组成，然后对F再做进一步的标准化：
$$F_{i,j}^{*}=\frac{F_{i,j}}{(\sum_{j=1}^{k}F_{i,j}^{2})^{1/2}}$$
       之后再对F的行进行kmeans或GMM聚类即可。
       以上是Ncut的内容。

       (3).一点题外话

       写到这里，如果只是应用spectral clustering，则此部分可以忽略，直接看下文pseudo-code部分即可，但是对于喜欢深入探究，不妨看一看。
       值得一提的是，从概率的视角出发，与上文推导也是不谋而合，而且得到的结论，与Ncut更是异曲同工。这种概率视角在多数论文里，称之为随机游走(Random walks)，在随机数学里，常见于马尔可夫模型。这部分的详细出处可以参考Lovaszl(1993: Random Walks on Graphs: A Survey)，以及Meila & Shi (2001:A Random Walks Views of Spectral Segmentation)。
       在随机游走框架下，通常都会构建一个转移概率矩阵(transition matrix)，同样，利用上文的邻接矩阵，可以得到该转移概率矩阵P为：

$$P_{i,j}=w_{i,j}/D_{i,j},\quad i,j=1,2,\cdots,N$$
       其中$w_{i,j}$，$D_{i,j}$ 分别为W，D矩阵中的元素。可以看出，其实$P=WD^{-1}$。为方便做进一步论述，假设一个初始分布概率$\pi = (\pi_{1},\pi_{2},\cdots,\pi_{i},\cdots,\pi_{N})^{T}$，其中：
$$\pi_{i}=\frac{D_{i,i}}{vol(V)}$$
       对于最简单的双簇聚类，假设将原图V分割成$A_{1}$以及 $A_{2}$，则目标是极小化下面概率：
$$min \frac{1}{2}(P(A_{1}|A_{2})+P(A_{2}|A_{1}))$$
       这种极小化概率的意义在于，样本在不同簇之间切换的可能性，应该尽量低。
       那么，对于$P(A_{1}|A_{2})$，可以写为：
$$P(A_{1}|A_{2})=\frac{P(A_{1},A_{2}}{P(A_{2})}$$
       对于$P(A_{2})$，可以理解为任意样本属于$A_{2}$的概率，则：
$$P(A_{2})=\frac{vol(A_{2})}{vol(A)}$$
可通过$A_{2}$簇里边和与整张图V边和做对比。
对于$P(A_{1},A_{2})$，可以理解为，任取一点，其从$A_{2}$切换至 $A_{1}$的概率，具体为：
$$\begin{align} P(A_{1},A_{2}) & =\sum_{i\in A_{1},j \in A_{2}}\pi_{i}P_{i,j}\　　& = \sum_{i\in A_{1},j \in A_{2}}\frac{D_{i,j}}{vol(V)}\frac{w_{i,j}}{D_{i,j}} \　　& = \frac{1}{vol(V)} \sum_{i\in A_{1},j \in A_{2}}w_{i,j} 　　\end{align}$$
       则：
$$\begin{align} P(A_{1}|A_{2}) & =\frac{\frac{1}{vol(V)} \sum_{i\in A_{1},j \in A_{2}}w_{i,j}}{\frac{vol(A_{2})}{vol(V)}} \& = \frac{ \sum_{i\in A_{1},j \in A_{2}}w_{i,j}}{vol(A_{2})} 　　\end{align}$$
       最终：
$$\begin{align} min \frac{1}{2}(P(A_{1}|A_{2})+P(A_{2}|A_{1})) & =min \frac{1}{2}(\frac{ \sum_{i\in A_{1},j \in A_{2}}w_{i,j}}{vol(A_{2})} + \frac{ \sum_{i\in A_{1},j \in A_{2}}w_{i,j}}{vol(A_{1})}) \& = min \frac{1}{2}\sum_{i=1}^{2}\frac{W(A_{i},\bar{A_{i}})}{vol(A_{i})} 　　\end{align}$$
       至此，可以看到Ncut在概率的视角上，也有了相应的立足之地,Ratiocut的推导也同理。

（二）、pesudo-code

       Spectral clustering的pseudo-code有很多种，这里只讲最常用的normalized版本，也就是Ncut作为切图法的版本。
       具体为：
       1.构造S矩阵(similarity matrix)，同时指定要聚类的簇数k；
       2.利用S矩阵构造W矩阵(adjacent matrix)；
       3.计算拉普拉斯矩阵L，其中L=D-W；
       4.对L矩阵标准化，即令 ；
       5.计算normalize后的L矩阵的前k个特征向量，按特征值升序排列；
       6.对k个向量组成的矩阵的行进行看kmeans/GMM聚类；
       7.将聚类结果的各个簇分别打上标记，对应上原数据，输出结果。

input: (data,K,k,kNNType,GaussianDelta,epsilon)
1. S = Euclidean(data) 
2. W = Gaussian( kNN(S,k,kNNType) , GaussianDelta ) 
3. L = D - W
4. L‘ = normalized(L)
5. EV = eigenvector(L‘,K)
6. while( (newCenter - oldCenter) > epsilon){
      newCenter = kmeans(EV,K)
   }
output：K clusters of SC

       其中，data为样本，对应code为二维数组，K为要聚类的簇数，k为kNN的邻接个数，kNNType为SC中的kNN函数类型，一般为kNN/mutalkNN，具体看上文说明，GaussianDelta为构造W矩阵时的高斯函数参数，epsilon为kmeans或者GMMs中的更新步长。SC输出为聚出K个簇。

四、一些参考文献

1. Meila,shi: A Random Walks View of Spectral Segmentation
2. Harry Yserentant: A short theory of the Rayleigh-Ritz method
3. Ulrike von Luxburg: A Tutorial on Spectral Clustering
4. Andrew Y.Ng, Michael I.Jordan, Yair Weiss: On Spectral Clustering Analysis and an algorithm
5. L. LOVASZ: Random Walks on Graphs: A Survey
6. Fan R.K. Chung: Spectral Graph Theory
7. Xiao-Dong Zhang: The Laplacian eigenvalues of graphs: a survey
8. Bojan Mohar: THE LAPLACIAN SPECTRUM OF GRAPHS
9. pluskid: http://blog.pluskid.org/?p=287
10. Wiki: https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
11. G. E. FORSYTHE AND P. HENRICI：THE CYCLIC JACOBI METHOD FOR COMPUTING THE PRINCIPAL VALUES OF A COMPLEX MATRIX
12. Jacobi Transformations of a Symmetric Matrix