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

最新推荐文章于 2025-03-21 10:14:07 发布

杨铖

最新推荐文章于 2025-03-21 10:14:07 发布

阅读量6.9w

点赞数 158

分类专栏：算法篇文章标签：谱聚类 SC clustering scala

本文链接：https://blog.youkuaiyun.com/yc_1993/article/details/52997074

版权

算法篇同时被 2 个专栏收录

3 篇文章

订阅专栏

scala

1 篇文章

订阅专栏

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

先把大纲列出来：

一、从狄多公主圈地传说说起
二、谱聚类的演算
  （一）、演算
      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)
       第二步是切图，即将第一步构造出来的按照一定的切边准则，切分成不同的图，而不同的子图，即我们对应的聚类结果，举例如下：
               切图4
                            图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 = {0, i f s i, j > ε ε, i f s i, j \leq ε

$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 $W$ 非对称，为克服这种问题，一般采取下面两种方法之一：
一是只要点

xi $x_{i}$ 在

xj $x_{j}$ 的K个近邻中或者

xj $x_{j}$ 在

xi $x_{i}$ 的K个近邻中，则保留

si,j $s_{i,j}$ ，并对其做进一步处理

W $W$ ，此时为：

W i, j = W j, i = ⎧ ⎩ ⎨ 0, i f x i \notin K N N (X j) a n d x j \notin K N N (X i) e - ∥ ∥ x i - x j ∥ ∥ 2 2 σ 2, i f x i \in K N N (X j) o r x j \in K N N (X i)

$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个近邻中，才会保留

si,j $s_{i,j}$ 并做进一步变换，此时

W $W$ 为：

W i, j = W j, i = ⎧ ⎩ ⎨ 0, i f x i \notin K N N (X j) o r x j \notin K N N (X i) e - ∥ ∥ x i - x j ∥ ∥ 2 2 σ 2, i f x i \in K N N (X j) a n d x j \in K N N (X i)

$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，一般使用高斯距离：

si,j=e−∥∥xi−xj∥∥22σ2 $s_{i,j}=e^{\frac{-\left \| x_{i}-x_{j} \right \|^{2}}{2\sigma ^{2}}}$ ，则重构之后的矩阵

W $W$ 与之前的相似矩阵

S $S$ 相同，为：

Wi,j=Si,j=[s]i,j $W_{i,j}= S_{i,j}=[s]_{i,j}$ 。
在了解三种构图方式后，还需要注意一些细节，对于第一二中构图，一般是重构基于欧氏距离的，而第三种构图方式，则是基于高斯距离的，注意到高斯距离的计算蕴含了这样一个情况：对于

∥∥xi−xj∥∥2 $\left \| x_{i}-x_{j} \right \|^{2}$ 比较大的样本点，其得到的高斯距离反而值是比较小的，而这也正是

S $S$ 可以直接作为

W $W$ 的原因，主要是为了将距离近的点的边赋予高权重。
得到邻接矩阵

W $W$ 后，需要做进一步的处理：
(1).计算阶矩(degree matrix)

D $D$ ：

D i, j = ⎧ ⎩ ⎨ ⎪ ⎪ 0, i f i \neq j \sum j w i, j, i f i = j

$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.$
其中其中

wi,j $w_{i,j}$ 为邻接矩阵

W $W$ 元素，

∑jwi,j $\sum_{j}w_{i,j}$ 表示将图中某点所连接的其他点的边的权重之和，可以看出，

D $D$ 为对角矩阵.
(2).计算拉普拉斯矩阵(Laplacians matrix)

L $L$ ：

L=D−W $L=D-W$
如此，在构图阶段基本就完成了，至于为什么要计算出拉普拉斯矩阵

L $L$ ，可以说

L=D−W $L=D-W$ 这种形式对于后面极大化问题是非常有利的(不知道前人是怎么想到的，不然

L $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$ ，则子集与子集之间连边的权重和为：

c u t (A 1, A 2, \dots, A k) = 1 2 \sum i k W (A i, A i ¯)

$cut(A_{1},A_{2},\cdots,A_{k})=\frac{1}{2}\sum_{i}^{k}W(A_{i},\bar{A_{i}})$
其中

Ai¯ $\bar{A_{i}}$ 为

Ai $A_{i}$ 的补集，意为除

Ai $A_{i}$ 子集外其他V的子集的并集，

W(Ai,Ai¯) $W(A_{i},\bar{A_{i}})$ 为

Ai $A_{i}$ 与其他子集的连边的和，为：

W (A i, A i ¯) = \sum m \in A i, n \in A i ¯ w m, n

$W(A_{i},\bar{A_{i}})=\sum_{m\in A_{i},n \in \bar{A_{i}}}w_{m,n}$
其中

wm,n $w_{m,n}$ 为邻接矩阵W中的元素。
由于我们切图的目的是使得每个子图内部结构相似，这个相似表现为连边的权重平均都较大，且互相连接，而每个子图间则尽量没有边相连，或者连边的权重很低，那么我们的目的可以表述为：

m i n c u t (A 1, A 2, \dots, A k)

$min\quad cut(A_{1},A_{2},\cdots,A_{k})$
但是稍微留意可以发现，这种极小化的切图存在问题，如下图：
切图5

图3 问题切图
如上图，如果用

mincut(A1,A2,⋯,Ak) $min\quad cut(A_{1},A_{2},\cdots,A_{k})$ 这种方法切图，会造成将V切成很多个单点离散的图，显然这样是最快且最能满足那个最小化操作的，但这明显不是想要的结果，所以有了RatioCut和Ncut，具体地：

R a t i o c u t (A 1, A 2, \dots, A k) = 1 2 \sum i k W ( A i , A i ¯ ) | A i |

$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 |}$

N c u t (A 1, A 2, \dots, A k) = 1 2 \sum i k W ( A i , A i ¯ ) v o l ( A i )

$Ncut(A_{1},A_{2},\cdots,A_{k})=\frac{1}{2}\sum_{i}^{k}\frac{W(A_{i},\bar{A_{i}})}{vol(A_{i})}$
其中

|Ai| $\left | A_{i} \right |$ 为

Ai $A_{i}$ 中点的个数，

vol(Ai) $vol(A_{i})$ 为

Ai $A_{i}$ 中所有边的权重和。
下面详细讲解这两种方法的计算：Ratiocut & Ncut

(1).Ratiocut

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

m i n R a t i o c u t (A 1, A 2, \dots, A k)

$min\quad Ratiocut(A_{1},A_{2},\cdots,A_{k})$
对上述极小化问题做一下转化，引入

{A1,A2,⋯,Ak} $\{A_{1},A_{2},\cdots ,A_{k}\}$ 的指示向量

hj={h1,h2,⋯,hi,⋯,hk},j=1,2,⋯,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 = ⎧ ⎩ ⎨ ⎪ ⎪ ⎪ ⎪ 1 ∣ ∣ A j ∣ ∣ - - - \sqrt, i f v i \in A j 0, i f v i \notin A 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.$
通俗理解就是，每个子集

Aj $A_{j}$ 对应一个指示向量

hj $h_{j}$ ，而每个

hj $h_{j}$ 里有N个元素，分别代表N个样本点的指示结果，如果在原始数据中第i个样本被分割到子集

Aj $A_{j}$ 里，则

hj $h_{j}$ 的第i个元素为

1∣∣Aj∣∣√ $\frac{1}{\sqrt{\left | A_{j} \right |}}$ ，否则为0。
进一步，计算，可以得到：

h T i L h i = h T i (D - W) h i = h T i D h i - h T i W h 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 2 i, m D m, m - \sum m = 1 \sum n = 1 h i, m h i, n w m, n = 1 2 (\sum m = 1 h 2 i, m D m, m - 2 \sum m = 1 \sum n = 1 h i, m h i, n w m, n + \sum n = 1 h 2 i, n D n, n)

$\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为对角矩阵，所以只有对角线上元素与向量

hi $h_{i}$ 相乘有意义；
第五步，代数变换。
进一步地，做如下：

h T i L h i = 1 2 (\sum m = 1 h 2 i, m D m, m - 2 \sum m = 1 \sum n = 1 h i, m h i, n w m, n + \sum n = 1 h 2 i, n D n, n) = 1 2 (\sum m = 1 \sum n = 1 h 2 i, m 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 2 i, n w n, m) = 1 2 \sum m = 1 \sum n = 1 w n, m (h i, m - h i, n) 2

$\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}$
同样分解如下：
第一步到第二步，由上文阶矩的定义可知，其对角线上元素

Dm,m=∑n=1wm,n $D_{m,m}=\sum_{n=1}w_{m,n}$
第二部到第三步，二次函数变换；
进一步，将上文

hj,i $h_{j,i}$ 代入上式，最终可以得到：

h T i L h i = 1 2 \sum m = 1 \sum n = 1 w m, n (h i, m - h i, n) 2 = 1 2 (\sum m \in A i, n \in A i ¯ w m, n (1 | A i | - - - \sqrt - 0) 2 + \sum m \in A i ¯, n \in A i w m, n (0 - 1 | A i | - - - \sqrt) 2) = 1 2 (\sum m \in A i, n \in A i ¯ w m, n 1 | A i | + \sum m \in A i ¯, n \in A i w m, n 1 | A i |) = 1 2 (c u t (A i, A i ¯) 1 | A i | + c u t (A i ¯, A i) 1 | A i |) = c u t ( A i , A i ¯ ) | A i | = R a t i o c u t (A i, A 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}$
分解如下：
第一步到第二步，将指示变量

hi,m $h_{i,m}$ ，

hi,n $h_{i,n}$ 对应

Ai $A_{i}$ 的元素分别代入；
第三，四，五步，分别做简单运算，以及将

∑m∈Ai,n∈Ai¯wm,n $\sum_{m\in A_{i},n \in \bar{A_{i}}}w_{m,n}$ 替换为

cut(Ai,Ai¯) $cut(A_{i},\bar{A_{i}})$ ；
可以看到，通过引入指示变量h，将L矩阵放缩到与Ratiocut等价，这无疑是在谱聚类中划时代的一举。
为了更进一步考虑进所有的指示向量，令

H={h1,h2,⋯,hk} $H=\{h_{1},h_{2},\cdots ,h_{k}\}$ ，其中

hi $h_{i}$ 按列排列，由

hi $h_{i}$ 的定义知每个

hi $h_{i}$ 之间都是相互正交，即

hi∗hj=0,i≠j $h_{i}*h_{j}=0,i \neq j$ ，且有

hi∗hi=1 $h_{i}*h_{i}=1$ ，可得到：

R a t i o c u t (A 1, A 2, \dots, A k) = \sum i = 1 k h T i L h i = \sum i = 1 k (H T L H) i i = T r (H T L H)

$\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表示对角线求和。
如此，便将极小化问题：

minRatiocut(A1,A2,⋯,Ak) $min\quad Ratiocut(A_{1},A_{2},\cdots,A_{k})$ ，转化为：

a r g m i n H T r (H T L H), s . t . H T H = I

$arg\mathop{min}\limits_{H} \quad Tr(H^{T}LH), \quad s.t. H^{T}H=I$
将之前切图的思路以一种严谨的数学表达表示出来，对于上面极小化问题，由于L矩阵是容易得到的，所以目标是求满足条件的H矩阵，以使得

Tr(HTLH) $Tr(H^{T}LH)$ 最小。
注意到一点，要求条件下的H，首先是求得组成H的每个

hi $h_{i}$ ，而每个

hi $h_{i}$ 都是Nx1的向量，且向量中每个值都是二值分布，即取值0或者

1/|Ai|−−−√ $1/\sqrt{\left | A_{i} \right |}$ ，那么对于每个元素，有2种选择，则当个

hi $h_{i}$ 就有

2N $2^{N}$ 种情况，对于整个H矩阵，则有

(C1k)N $(C_{k}^{1})^{N}$ 种情况，对于N很大的数据，这无疑是灾难性的NP-hard问题，显然我们不可能遍历所有的情况来求解，那么问题是不是没有解？答案自然是否定的，虽然我们没办法求出精确的，但是我们可以用另外一种方法近似代替

hi $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 = H i , j ( \sum k j = 1 H 2 i , j ) 1 / 2

$H_{i,j}^{*}=\frac{H_{i,j}}{(\sum_{j=1}^{k}H_{i,j}^{2})^{1/2}}$
现在H矩阵也完成了，剩下的就是对样本聚类了，要明白我们的目标不是求

Tr(HTLH) $Tr(H^{T}LH)$ 的最小值是多少，而是求能最小化

Tr(HTLH) $Tr(H^{T}LH)$ 的H，所以聚类的时候，分别对H中的行进行聚类即可，通常是kmeans，也可以是GMM，具体看效果而定。
       至于为什么是对H的行进行聚类，有两点原因：
       1.注意到H除了是能满足极小化条件的解，还是L的特征向量，也可以理解为W的特征向量，而W则是我们构造出的图，对该图的特征向量做聚类，一方面聚类时不会丢失原图太多信息，另一方面是降维加快计算速度，而且容易发现图背后的模式。
       2.由于之前定义的指示向量

hi $h_{i}$ 是二值分布，但是由于NP-hard问题的存在导致

hi $h_{i}$ 无法显式求解，只能利用特征向量进行近似逼近，但是特征向量是取任意值，结果是我们对

hi $h_{i}$ 的二值分布限制进行放松，但这样一来

hi $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的目标，也是极小化各子图连边的和，如下：

m i n N c u t (A 1, A 2, \dots, A k)

$min\quad Ncut(A_{1},A_{2},\cdots,A_{k})$
下面对该问题做一下转化，先引入

{A1,A2,⋯,Ak} $\{A_{1},A_{2},\cdots ,A_{k}\}$ 的指示变量

hj={h1,h2,⋯,hi,⋯,hk},j=1,2,⋯,k $h_{j}=\{h_{1},h_{2},\cdots,h_{i},\cdots ,h_{k}\},\quad j=1,2,\cdots,k$ ，同样，其中i表示样本下标，j表示子集下标，

hj,i $h_{j,i}$ 表示样本i对子集j的指示，不同的是，其具体为：

h j, i = ⎧ ⎩ ⎨ ⎪ ⎪ 1 v o l ( A i ) - - - - - - \sqrt, i f v i \in A j 0, i f v i \notin A 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个样本被分割到子集

Aj $A_{j}$ 里，则

hj $h_{j}$ 的第i个元素为

1/vol(Ai)−−−−−−√ $1/\sqrt{vol(A_{i})}$ ，否则为0。
进一步，计算

hTiLhi $h_{i}^{T}Lh_{i}$ ，可以得到：

h T i L h i = h T i D h i - h T i W h 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 2 i, m D m, m - \sum m = 1 \sum n = 1 h i, m h i, n w m, n = 1 2 (\sum m = 1 h 2 i, m D m, m - 2 \sum m = 1 \sum n = 1 h i, m h i, n w m, n + \sum n = 1 h 2 i, n D n, n) = 1 2 (\sum m = 1 \sum n = 1 h 2 i, m 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 2 i, n w n, m) = 1 2 \sum m = 1 \sum n = 1 w m, n (h i, m - h i, n) 2 = 1 2 (\sum m \in A i, n \in A i ¯ w m, n (1 v o l ( A i ) - - - - - - \sqrt - 0) 2 + \sum m \in A i ¯, n \in A i w m, n (0 - 1 v o l ( A i ) - - - - - - \sqrt) 2) = c u t ( A i , A i ¯ ) v o l ( A i ) = N c u t (A i, A 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}$
由于推导与上文类似，这里就综合在一起，主要是为了说明

hTiLhi $h_{i}^{T}Lh_{i}$ 如何得到

Ncut(Ai,Ai¯) $Ncut(A_{i},\bar{A_{i}})$ ，以便将

minNcut(A1,A2,⋯,Ak) $min\quad Ncut(A_{1},A_{2},\cdots,A_{k})$ 用严谨数学式子进行表达。
       这里对上述一大串式子分解一下：
       第一，二，三步，做一些简单的变换，如L=D-W代入，将向量内积转换成连加形式，同样，D是对角矩阵，可省去一些项；
       第四，五，六步，主要是做二次函数变换；
       第七步，将

hi,m $h_{i,m}$ ，

hi,n $h_{i,n}$ 对应

Ai $A_{i}$ 的元素分别代入；
第八，九步，得到

Ncut(Ai,Ai¯) $Ncut(A_{i},\bar{A_{i}})$ 形式；
此外，相比与Ratiocut，Ncut特有的一点性质是：

hTiDhi=1 $h_{i}^{T}Dh_{i}=1$ ，具体如下：

h T i D h i = \sum n = 1 h 2 i, n D n, n = \sum v n \in A i 1 v o l ( A i ) D n, n + \sum v n \notin A i 0 \cdot D n, n = 1 v o l ( A i ) \sum v n \in A i w v n + 0 = 1 v o l ( A i ) v o l (A 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={h1,h2,⋯,hk} $H=\{h_{1},h_{2},\cdots ,h_{k}\}$ ，则

HTDH=I $H^{T}DH=I$ ,其中I为单位对角矩阵。
同样，可以得到：

N c u t (A 1, A 2, \dots, A k) = \sum i = 1 k h T i L h i = \sum i = 1 k (H T L H) i i = T r (H T L H)

$\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}$
如此，便将极小化问题：

minNcut(A1,A2,⋯,Ak) $min\quad Ncut(A_{1},A_{2},\cdots,A_{k})$ ，转化为：

a r g m i n H T r (H T L H), s . t . H T D H = I

$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/2F $H=D^{-1/2}F$ ，

D−1/2 $D^{-1/2}$ 意为对D对角线上元素开方再求逆。将

H=D−1/2F $H=D^{-1/2}F$ 代入

argminHTr(HTLH),s.t.HTDH=I $arg\mathop{min}\limits_{H} \quad Tr(H^{T}LH), \quad s.t. H^{T}DH=I$ ，可以得到：

H T L H = (D - 1 / 2 F) T L D - 1 / 2 F = F T D - 1 / 2 L D - 1 / 2 F

$\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}$
以及：

H T D H = (D - 1 / 2 F) T D D - 1 / 2 F = F T D - 1 / 2 L D - 1 / 2 F = F T I F = I

$\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}$
至此，目标操作

argminHTr(HTLH),s.t.HTDH=I $arg\mathop{min}\limits_{H} \quad Tr(H^{T}LH), \quad s.t. H^{T}DH=I$ 可以修改为：

a r g m i n H T r (F T D - 1 / 2 L D - 1 / 2 F), s . t . F T F = 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/2LD−1/2 $D^{-1/2}LD^{-1/2}$ 这一步操作，就是对L矩阵进行normalize，而normalize可以理解为标准化，具体为：

Li,jvol(Ai)vol(Aj)√ $\frac{L_{i,j}}{\sqrt{vol( A_{i})vol(A_{j})}}$ ，这么做的一个好处就是对L中的元素进行标准化处理使得不同元素的量纲得到归一，具体理解就是，同个子集里，不同样本点之间的连边可能大小比较相似，这点没问题，但是对于不同子集，样本点之间的连边大小可能会差异很大，做这一步normalize操作，可以将L中的元素归一化在[-1,1]之间，这样量纲一致，对算法迭代速度，结果的精度都是有很大提升。
最后，F矩阵的求解可通过求

D−1/2LD−1/2 $D^{-1/2}LD^{-1/2}$ 的前k个特征向量组成，然后对F再做进一步的标准化：

F * i, j = F i , j ( \sum k j = 1 F 2 i , j ) 1 / 2

$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, i, j = 1, 2, \dots, N

$P_{i,j}=w_{i,j}/D_{i,j},\quad i,j=1,2,\cdots,N$
其中

wi,j $w_{i,j}$ ，

Di,j $D_{i,j}$ 分别为W，D矩阵中的元素。可以看出，其实

P=WD−1 $P=WD^{-1}$ 。为方便做进一步论述，假设一个初始分布概率

π=(π1,π2,⋯,πi,⋯,πN)T $\pi = (\pi_{1},\pi_{2},\cdots,\pi_{i},\cdots,\pi_{N})^{T}$ ，其中：

π i = D i , i v o l ( V )

$\pi_{i}=\frac{D_{i,i}}{vol(V)}$
对于最简单的双簇聚类，假设将原图V分割成

A1 $A_{1}$ 以及

A2 $A_{2}$ ，则目标是极小化下面概率：

m i n 1 2 (P (A 1 | A 2) + P (A 2 | A 1))

$min \frac{1}{2}(P(A_{1}|A_{2})+P(A_{2}|A_{1}))$
这种极小化概率的意义在于，样本在不同簇之间切换的可能性，应该尽量低。
那么，对于

P(A1|A2) $P(A_{1}|A_{2})$ ，可以写为：

P (A 1 | A 2) = P ( A 1 , A 2 P ( A 2 )

$P(A_{1}|A_{2})=\frac{P(A_{1},A_{2}}{P(A_{2})}$
对于

P(A2) $P(A_{2})$ ，可以理解为任意样本属于

A2 $A_{2}$ 的概率，则：

P (A 2) = v o l ( A 2 ) v o l ( A )

$P(A_{2})=\frac{vol(A_{2})}{vol(A)}$
可通过

A2 $A_{2}$ 簇里边和与整张图V边和做对比。
对于

P(A1,A2) $P(A_{1},A_{2})$ ，可以理解为，任取一点，其从

A2 $A_{2}$ 切换至

A1 $A_{1}$ 的概率，具体为：

P (A 1, A 2) = \sum i \in A 1, j \in A 2 π i P i, j = \sum i \in A 1, j \in A 2 D i , j v o l ( V ) w i , j D i , j = 1 v o l ( V ) \sum i \in A 1, j \in A 2 w i, j

$\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}$
则：

P (A 1 | A 2) = 1 v o l ( V ) \sum i \in A 1 , j \in A 2 w i , j v o l ( A 2 ) v o l ( V ) = \sum i \in A 1 , j \in A 2 w i , j v o l ( A 2 )

$\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}$
最终：

m i n 1 2 (P (A 1 | A 2) + P (A 2 | A 1)) = m i n 1 2 (\sum i \in A 1 , j \in A 2 w i , j v o l ( A 2 ) + \sum i \in A 1 , j \in A 2 w i , j v o l ( A 1 )) = m i n 1 2 \sum i = 1 2 W ( A i , A i ¯ ) v o l ( A i )

$\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,sigma,epsilon)
1. S = Euclidean(data) 
2. W = Gaussian( kNN(S,k,kNNType) , sigma ) 
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，具体看上文说明，sigma为构造W矩阵时的高斯函数参数，epsilon为kmeans或者GMMs中的更新步长。SC输出为聚出K个簇。

三、谱聚类的实现

实现过程涉及到的一些概念有：Similarity Matrix、kNN/mutual kNN、Laplacian Matrix、Normalized、Eigenvector(Jacobi methond)、kmeans/GMM；下面一一按序解析，使用语言为scala，这里需要说明一点，由于技术保密，这里有些东西只能介绍一些简单的版本，动手操作可以发现，上述计算在小样本还算过得去，样本一大简直无法入目，只能各位自己多去查看论文了，若日后有缘，会开篇做另外的优化阐述。

（一）、Similarity Matrix

// calculate the similarity matrix
def calculateSimilarityMatrix(SCInput: Array[Array[Double]]): Array[Array[Double]] = {
    // the Euclidean Distance
    def SCEuclideanDistance(SCEDInput: Array[Double]): Array[Double] = {
        SCInput.map(_.zip(SCEDInput))
               .map(a => a.map(b => math.pow(b._1 - b._2, 2)).sum)
     }
     SCInput.map(SCEuclideanDistance)
}

       这里不代入数据了，直接构造一个similarity Matrix了，如下：

                     图4 similarity Matrix

（二）、kNN/mutual kNN

下面是对上面similarity matrix一步构造的相似图做调整，将其转化为W（adjacency matrix）.

    // define the kNN function
    def kNN(k: Int, kNNType: String, sigma: Double = 1.0, SMatrix: Array[Array[Double]]): Array[Array[Double]] = {
        val len = SMatrix.length
        val AdjacencyMatrix = Array.ofDim[Double](len, len)

        // define the function of calculating the mutual adjacency matrix (for mutualkNN)
        @tailrec
        def calculateMutualAdjacencyMatrix(n: Int, S: Array[Array[Double]]): Array[Array[Double]] = {
            if (n < len) {
                // take out the smallest k values
                val kSmallestValue = S(n).zipWithIndex.sortWith((a, b) => a._1 < b._1).take(k + 1)
                val indexOfValue = kSmallestValue.map(_._2).distinct
                // calculate the Gaussian similarity value
                for (i <- indexOfValue) {
                    val GaussianSimilarity = math.exp(-S(n)(i) / 2 / sigma / sigma)
                    AdjacencyMatrix(n)(i) = GaussianSimilarity
                }
                calculateMutualAdjacencyMatrix(n + 1, S)
            } else {
                // judge mutual or not.
                for (i <- 0 until len) {
                    val notZero = AdjacencyMatrix(i).zipWithIndex.filter(_._1 != 0.0).map(_._2)
                    for (j <- notZero) if (AdjacencyMatrix(j)(i) == 0.0) AdjacencyMatrix(i)(j) = 0.0
                }
                AdjacencyMatrix
            }
        }

        // define the function of calculating the mutual adjacency matrix (for kNN)
        @tailrec
        def calculateAdjacencyMatrix(n: Int, S: Array[Array[Double]]): Array[Array[Double]] = {
            if (n < len) {
                // take out the smallest k values
                val kSmallestValue = S(n).zipWithIndex.sortWith((a, b) => a._1 < b._1).take(k + 1)
                val indexOfValue = kSmallestValue.map(_._2).distinct
                // calculate the Gaussian similarity value
                for (i <- indexOfValue) {
                    val GaussianSimilarity = math.exp(-S(n)(i) / 2 / sigma / sigma)
                    AdjacencyMatrix(n)(i) = GaussianSimilarity
                    AdjacencyMatrix(i)(n) = GaussianSimilarity
                }
                calculateAdjacencyMatrix(n + 1, S)
            } else {
                AdjacencyMatrix
            }
        }

        if (kNNType == "mutualkNN") {
            calculateMutualAdjacencyMatrix(0, SMatrix)
        } else {
            calculateAdjacencyMatrix(0, SMatrix)
        }
    }

       两种操作结果如下：kNN和mutualkNN。PS：由于是直接在IntelliJ直接截图，所以看起来可能没那么好看(ㄒoㄒ)
               kNN
                     图5 adjacency Matrix (kNN)

图6 adjacency Matrix (mutualkNN)

两个矩阵计算时k值均指定为2，sigma指定为1，可以观察出，上文矩阵 $\begin{bmatrix} 0.0 & 8.0 & 7.0 & 5.0 \\ 8.0 & 0.0 & 6.0 & 10.0 \\ 7.0 & 6.0 & 0.0 & 3.0 \\ 5.0 & 10.0 & 3.0 & 0.0 \end{bmatrix}$ ，其对应的kNN处理后的W矩阵，比较容易理解，即取任意样本最近的2个样本作为连点即可，而mutualkNN处理后的W矩阵，则需保证连点的两边均是该两点的最近的2个样本之一，举个例子，与第1个样本最近的两个点分别为第3，4个点，其距离分别为7.0和5.0，而第3个点最近的两个点分别为第2，4个点，其距离为6.0和3.0，所以这样下来，与第1个样本连边的点就只有第4个点，看图6可明白。

（三）、 Laplacian Matrix

下面利用上面W（adjacency matrix）矩阵，进一步计算D（degree matrix）矩阵和L（laplacian matrix）矩阵.

    // define the function of calculating the Laplacian Matrix
    def calculateLaplacianMatrix(adjacencyMatrix: Array[Array[Double]]): Array[Array[Double]] = {
        val len = adjacencyMatrix.length
        val laplacianMatrix = Array.ofDim[Double](len, len)

        // define the function of calculating the degree matrix
        def calculateDegreeMatrix(AM: Array[Array[Double]]): Array[Array[Double]] = {
            val degreeMatrix = Array.ofDim[Double](len, len)
            for (i <- 0 until len) {
                degreeMatrix(i)(i) = AM(i).sum
            }
            degreeMatrix
        }
        val degreeMatrix = calculateDegreeMatrix(adjacencyMatrix)

        // calculate the Laplacian matrix
        for (i <- 0 until len) {
            laplacianMatrix(i) = degreeMatrix(i).zip(adjacencyMatrix(i)).map(a => a._1 - a._2)
        }
        laplacianMatrix
    }

       这里我把D的计算和L的计算写在一起，一步计算到位，结果如下：
               degreeMatrix
                     图7 degree Matrix

laplacianMatrix
图8 laplacian Matrix

（四）、 Normalized

下面进一步对上面L（laplacian matrix）矩阵进行正则化处理.

    // define the function of calculating the normalized Laplacian Matrix
    def Normalized(laplacianMatrix: Array[Array[Double]], adjacencyMatrix: Array[Array[Double]]): Array[Array[Double]] = {
        val len = adjacencyMatrix.length

        // define the function of calculate the -1/2 power of the degree matrix
        def calculateAdjustDegreeMatrix(AM: Array[Array[Double]]): Array[Array[Double]] = {
            val adjustDegreeMatrix = Array.ofDim[Double](len, len)
            for (i <- 0 until len) {
                adjustDegreeMatrix(i)(i) = 1 / math.pow(AM(i).sum, 0.5)
            }
            adjustDegreeMatrix
        }
        val adjustDegreeMatrix = calculateAdjustDegreeMatrix(adjacencyMatrix)

        // calculate the normalized laplacian matrix
        def matrixProduct(left: Array[Array[Double]], right: Array[Array[Double]]): Array[Array[Double]] = {
            val len = left.length
            val output = Array.ofDim[Double](len,len)

            for(i <- 0 until len; j <- i until len){
                output(i)(j) = left(i).zip(right.map(_(j))).map(a => a._1 * a._2).sum
                output(j)(i) = output(i)(j)
            }
            output
        }
        val temp = matrixProduct(adjustDegreeMatrix,laplacianMatrix)
        val normalizedLaplacianMatrix = matrixProduct(temp,adjustDegreeMatrix)
        normalizedLaplacianMatrix
    }

       下面是上文L（laplacian matrix）矩阵的正则结果：

                     图9 Normalized laplacian Matrix

（五）、 Eigenvector(Jacobi methond)

下面进一步对上面正则化处理之后的L矩阵L’，取对应的特征向量，组成新的矩阵，特征向量的计算这里用的是串行的Jacobi旋转方法，内在逻辑就是对L’矩阵行列转换，得到L’的相似矩阵，且满足该相似矩阵为对角矩阵。值得一说的是，这种方法小样本可行，大样本效率很低，之后看情况，看是否再把优化的方法写出来，看缘分吧。。。

    // define the function of calculating the k smallest eigenvectors of normalized laplacian matrix with Jacobi method.
    def kSmallestEigenvectors(k: Int, normalizedLaplacian: Array[Array[Double]]): (Array[Double], Array[Array[Double]]) = {

        val len = normalizedLaplacian.length

        // initial the eigenvector matrix.
        val eigenvectorMatrix = Array.ofDim[Double](len, len)
        for (i <- 0 until len) eigenvectorMatrix(i)(i) = 1.0

        // initial the parameter of epsilon.
        val epsilon = (math.pow(10, -10), -1)

        // calculate the largest one Of normalized laplacian matrix(off-diagonal).
        def calculateLargestOfNormL(Input: Array[Array[Double]]): (Double, Int) = {
            val temp = Input.map(_.zipWithIndex)
            val largestOfRow = new Array[(Double, Int)](len - 1)
            for (i <- 0 until (len - 1)) {
                largestOfRow(i) = temp(i).filter(_._2 > i).sortWith((a, b) => a._1.abs > b._1.abs).head
            }
            val largestOfInput = largestOfRow.sortWith((a, b) => a._1.abs > b._1.abs).head
            if (epsilon._1 > largestOfInput._1.abs) epsilon else largestOfInput
        }
        var largestOfNormL = calculateLargestOfNormL(normalizedLaplacian)

        // Judge condition.
        var loop: Boolean = true

        // main loop.
        while (loop) {

            // the index of the largest value of normalized laplacian matrix.
            val mRow = largestOfNormL._2
            val mCol = normalizedLaplacian(largestOfNormL._2).indexOf(largestOfNormL._1)

            // calculate the new normalized Laplacian matrix: angle, sin, cos, new(row,col)
            // cache the temp value.
            val nL_ij = normalizedLaplacian(mRow)(mCol)
            val nL_ii = normalizedLaplacian(mRow)(mRow)
            val nL_jj = normalizedLaplacian(mCol)(mCol)
            val angle = if (nL_jj == nL_ii) math.Pi / 4.0 else 0.5 * math.atan2(2 * nL_ij, nL_jj - nL_ii)
            val sinAngle = math.sin(angle)
            val cosAngle = math.cos(angle)
            // update the normalized Laplacian matrix (ii, jj, ij, ji)
            normalizedLaplacian(mRow)(mRow) = nL_ii * cosAngle * cosAngle + nL_jj * sinAngle * sinAngle - 2 * nL_ij * cosAngle * sinAngle
            normalizedLaplacian(mCol)(mCol) = nL_ii * sinAngle * sinAngle + nL_jj * cosAngle * cosAngle + 2 * nL_ij * cosAngle * sinAngle
            normalizedLaplacian(mRow)(mCol) = (cosAngle * cosAngle - sinAngle * sinAngle) * nL_ij + sinAngle * cosAngle * (nL_ii - nL_jj)
            normalizedLaplacian(mCol)(mRow) = normalizedLaplacian(mRow)(mCol)
            for (i <- (0 until len).filter(a => a != mRow && a != mCol)) {
                val tempRi = normalizedLaplacian(mRow)(i)
                val tempCi = normalizedLaplacian(mCol)(i)
                normalizedLaplacian(mRow)(i) = cosAngle * tempRi - sinAngle * tempCi
                normalizedLaplacian(mCol)(i) = sinAngle * tempRi + cosAngle * tempCi
                normalizedLaplacian(i)(mRow) = normalizedLaplacian(mRow)(i)
                normalizedLaplacian(i)(mCol) = normalizedLaplacian(mCol)(i)
            }

            // update the eigenvector matrix.
            for (i <- 0 until len) {
                val eigenIR = eigenvectorMatrix(i)(mRow)
                val eigenIC = eigenvectorMatrix(i)(mCol)
                eigenvectorMatrix(i)(mRow) = eigenIR * cosAngle - eigenIC * sinAngle
                eigenvectorMatrix(i)(mCol) = eigenIC * cosAngle + eigenIR * sinAngle
            }

            // update the temp value again.
            largestOfNormL = calculateLargestOfNormL(normalizedLaplacian)

            // update the judge condition.
            if (largestOfNormL._2 == -1) loop = false

        }

        // return the k smallest eigenvalue and it's corresponding eigenvector matrix.
        val eigenvalue = new Array[Double](len)
        for (i <- 0 until len) eigenvalue(i) = normalizedLaplacian(i)(i)
        def kSmallest(eigenvalue: Array[Double], eigenvector: Array[Array[Double]]): (Array[Double], Array[Array[Double]]) = {
            val eigenvalueWithIndex = eigenvalue.zipWithIndex.sortWith((a, b) => a._1 < b._1).take(k)
            val ouput = Array.ofDim[Double](len, k)
            var n: Int = 0
            for ((v, i) <- eigenvalueWithIndex) {
                for (j <- 0 until len) {
                    ouput(j)(n) = eigenvector(j)(i)
                }
                n += 1
            }
            (eigenvalueWithIndex.map(_._1), ouput)
        }
        val kSmallestValueVector = kSmallest(eigenvalue, eigenvectorMatrix)

        // final ouput.
        (kSmallestValueVector._1, kSmallestValueVector._2)

    }

       为直观一点，下面拿矩阵 $\begin{bmatrix} 0.0 & 8.0 & 7.0 & 5.0 \\ 8.0 & 0.0 & 6.0 & 10.0 \\ 7.0 & 6.0 & 0.0 & 3.0 \\ 5.0 & 10.0 & 3.0 & 0.0 \end{bmatrix}$ 计算特征值和特征向量，结果如下：

                     图10 特征值（对角线）

图11 特征向量

结果可以和Matlab和R，python做比较，我只和R比对过，基本无差。

（六）、 kmeans/GMM

再做一步最终的聚类，便基本完成了算法。

    // define the kmeans function.
    def kmeans(K: Int, eigenvector: Array[Array[Double]]): Array[(Int, Array[Double])] = {

        val len = eigenvector.length
        val col = eigenvector(0).length

        // initial the random centers.
        var center: Map[Int, Array[Double]] = Map()
        var Karr: List[Int] = Nil
        while(Karr.length < K){
            val RandomK = Random.nextInt(len)
            if(!Karr.contains(RandomK)) Karr = Karr ::: List(RandomK)
        }
        for (i <- 0 until K) {
            center += (i -> eigenvector(Karr(i)))
        }

        // classify the points into K clusters with the present center.
        def classify(ct: Map[Int, Array[Double]], input: Array[Array[Double]]): Array[(Int, Array[Double])] = {

            // calculate the euclidean distance.
            val tempArr = input.map(a => {
                val euclidean = new Array[Double](K)
                for (i <- 0 until K) {
                    euclidean(i) = ct(i).zip(a).map(m => math.pow(m._1 - m._2, 2)).sum
                }
                euclidean.zipWithIndex
            })

            // tagging the points.
            val tagging = tempArr.map(a => {
                val tag = a.sortWith((x, y) => x._1 < y._1).head
                tag._2
            })

            // output.
            val output = tagging.zip(input)
            output
        }
        //val pointsWithTag = classify(center, eigenvector)

        // update the center.
        def updateCenter(oldCenter: Map[Int, Array[Double]], PWT: Array[(Int, Array[Double])]): Map[Int, Array[Double]] = {

            // groupby the result for computing.
            val clusters = PWT.groupBy(_._1)

            // update the newCenter
            var newCenter: Map[Int, Array[Double]] = Map()
            for (i <- 0 until K) {
                val clustersI = clusters.get(i) match {
                    case Some(s) => s.map(_._2)
                    //case None => new Array(col)
                }
                val n = clustersI.length
                val centerI = for (j <- 0 until col) yield clustersI.map(_ (j)).sum / n
                newCenter += (i -> centerI.toArray)
            }
            newCenter
        }
        //val newCenter = updateCenter(center,pointsWithTag)

        // initialize the blank variable.
        var pointsWithTag: Array[(Int, Array[Double])] = Array()
        var newCenter: Map[Int, Array[Double]] = Map()
        var movement: Seq[Double] = Seq()

        // loop.
        var loop: Boolean = true
        var j = 0
        while (loop) {
            j += 1

            // tagging the points and update the center.
            pointsWithTag = classify(center, eigenvector)
            newCenter = updateCenter(center, pointsWithTag)

            // the movement of the center.
            movement = for (i <- 0 until K) yield newCenter(i).zip(center(i)).map(a => (a._1 - a._2).abs).sum

            // judge the movement is small enough or not.  movement.exists(_ > math.pow(10, -10))
            if (movement.exists(_ > math.pow(10, -5))) {
                center = newCenter
            } else {
                loop = false
            }
        }
        pointsWithTag
        //newCenter
        //j
    }

按照上文pesudo-code的指引，一步一步操作，最终输出结果如下，spectral clustering和kmeans的对比（twocircles datase）：

图12 spectral clustering和kmeans的对比

四、一些参考资料

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