谱聚类原理总结
引入
聚类算法一般可以分为两类：
- Compactness。代表的算法有 K-means，GMM 等。但这类算法只能处理凸集，为了处理非凸的样本集，必须引⼊核技巧。
- Connectivity。这类以 spectral clustering 为代表。
举个例子，将下述数据采用聚类算法进行聚类，可以采用 GMM 或 K-Means 的方法：
然而对于下述数据却并不能使用上述两种算法：
此时可以考虑采用谱聚类（spectral clustering）的方法。
谱聚类算法(Spectral Clustering)
主要思想是把所有的数据看做空间中的点，这些点之间可以用边连接起来。距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高，通过对所有数据点组成的图进行切图，让切图后不同的子图间边权重和尽可能的低，而子图内的边权重和尽可能的高，从而达到聚类的目的。
1 基础知识
1.1 无向带权图
对于一个图 $G$ ，一般用 $V$ 代表点的集合 和用 $E$ 来描述边的集合。则 $G(V,E)$ 。其中， $V=(v_1, v_2,...v_n)$ 。对于 $V$ 中的任意两个点，可以有边连接，也可以没有边连接。定义权重 $w_{ij}$ 为点 $v_i$ 和点 $v_j$ 之间的权重。对于无向图，有 $w_{ij} = w_{ji}$。
权重：对于有边相连的两个点 $v_i$ 和 $v_j$，有 $w_{ij} > 0$ ， 对于没有边连接的两个点 $v_i$ 和 $v_j$ ，$w_{ij} = 0$。
度：对于图中的任意一个点 $v_i$ ，它的度 $d_i$ 定义为和它相连的所有边的权重之和，即 ：
$d_i = \sum\limits_{j=1}^{n}w_{ij}$
利用节点的权重值，可以得到图的邻接矩阵 $W$ ，是一个 $n \times n$ 的矩阵，第 $i$ 行的第 $j$ 个值对应的权重为 $w_{ij}$。
除此之外，对于点集 $V$ 的的一个子集 $A \subset V$ ，我们定义：
-
- $|A|: = 子集A中点的个数$
- $ vol(A): = \sum\limits_{i \in A}d_i$
利用每个节点的度，可以得到一个 $n \times n$ 的度矩阵 $D$，它是一个对角矩阵，只有主对角线有值，对应第 $i$ 行的第 $i个$ 点的度数，定义如下：
$\mathbf{D} =\left( \begin{array}{ccc}d_1 & \ldots & \ldots \\ \ldots & d_2 & \ldots \\ \vdots & \vdots & \ddots \\ \ldots & \ldots & d_n \end{array} \right)$
1.2 相似矩阵
一般情况下，邻接矩阵 $W$，通常可以自己输入权重，但在谱聚类中，我们只有数据点的定义，并没有直接给出这个邻接矩阵，那么怎么得到这个邻接矩阵呢？
基本思想：距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高。一般来说，可以通过样本点距离度量的相似矩阵 $S$ 来获得邻接矩阵 $W$。
构建邻接矩阵 $W$ 的方法有三类：
- $\epsilon$ - 邻近法
- $K$ 邻近法
- 全连接法
1.2.1 $\epsilon$-邻近法
首先，设置一个距离阈值 $\epsilon$ ，然后用欧式距离 $s_{ij}$ 度量任意两点 $x_i$ 和 $x_j$ 的距离。即相似矩阵的 $s_{ij} = ||x_i-x_j||_2^2$ ，然后根据 $s_{ij}$ 和 $\epsilon$ 的大小关系，来定义邻接矩阵 $W$ 如下：
$w_{ij}=\begin{cases}0& {s_{ij} > \epsilon}\\\epsilon& {{s_{ij} \leq \epsilon}}\end{cases}$
从上式可见，两点间的权重要不是 $\epsilon$，要不就是 $0$。显然这很难精确区分每个点之间的距离大小。所以很少使用 $\epsilon$- 邻近法。
1.2.2 $K$ 邻近法
基本思想：利用 KNN 算法遍历所有的样本点，取每个样本最近的 $k$ 个点作为近邻，只有和样本距离最近的 $k$ 个点之间的 $w_{ij} > 0$ 。
但是这种方法会造成重构之后的邻接矩阵 $W$ 非对称，且后面的算法需要对称邻接矩阵。为解决这种问题，一般采取下面两种方法之一：
Method1：第一种 $K$ 邻近法是只要一个点在另一个点的 $K$ 近邻中，则保留 $s_{ij}$：
$w_{ij}=w_{ji}=\begin{cases}0& {x_i \notin KNN(x_j) \;and \;x_j \notin KNN(x_i)}\\exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})& {x_i \in KNN(x_j)\; or\; x_j \in KNN(x_i})\end{cases}$
Method2：第二种 $K$ 邻近法是必须两个点互为 $K$ 近邻中，才保留 $s_{ij}$
$w_{ij}=w_{ji}=\begin{cases}0& {x_i \notin KNN(x_j) \;or\;x_j \notin KNN(x_i)}\\exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})& {x_i \in KNN(x_j)\; and \; x_j \in KNN(x_i})\end{cases}$
1.2.3 全连接法
相比前两种方法，第三种方法所有的点之间的权重值都大于 $0$，因此称之为全连接法。
基本思想：选择不同的核函数来定义边权重，常用的有多项式核函数，高斯核函数和 Sigmoid 核函数。
最常用的是高斯核函数RBF，此时相似矩阵和邻接矩阵相同：
$w_{ij}=s_{ij}=exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})$
在实际的应用中，使用第三种全连接法来建立邻接矩阵是最普遍的，而在全连接法中使用高斯径向核 RBF 是最普遍的。
2 拉普拉斯矩阵
拉普拉斯矩阵 $L= D-W$。$D$ 即为度矩阵，它是一个对角矩阵。而 $W$ 即为 邻接矩阵。（参考《图神经网络基础二：谱图理论》）
举例：
普拉斯矩阵的性质如下：
1、拉普拉斯矩阵是 对称矩阵;
2、$L$ 的行和为零 ；
3、$\mathbf{L} $ 有一个特征值为零 ；
4、$\mathrm{L} $ 是半正定矩阵 ；
5、对于任意向量 $f$，有：
$f^{T} L f=\frac{1}{2} \sum \limits _{i=1}^{N} \sum \limits _{j=1}^{N} w_{i j}\left(f_{i}-f_{j}\right)^{2}$
这一性质利用拉普拉斯矩阵的性质很容易可以得到：
$\begin{array}{l}f^{T} L f&=f^{T} D f-f^{T} W f \\&=\sum \limits _{i=1}^{N} d_{i} f_{i}^{2}-\sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N} w_{i j} f_{i} f_{j} \\&=\frac{1}{2}\left(\sum\limits_{i=1}^{N} d_{i} f_{i}^{2}-2 \sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N} w_{i j} f_{i} f_{j}+\sum\limits_{j=1}^{N} d_{j} f_{j}^{2}\right) \\&=\frac{1}{2}\left(\sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N} w_{i j} f_{i}^{2}-2 \sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N} w_{i j} f_{i} f_{j}+\sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N} w_{i j} f_{j}^{2}\right) \\&=\frac{1}{2} \sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N} w_{i j}\left(f_{i}-f_{j}\right)^{2}\end{array}$
对于上面的性质5，如果 $f$ 为网络中信号的值的向量，那么 $f^{T} L f$ 称为图信号的总变差 (Total Variation)，可以刻画图信号整体的平滑度
3 无向图切图
对于无向图 $G$ 的切图，目标是将图 $G(V,E)$ 切成相互没有连接的 $k$ 个子图，每个子图点集合为：$A_1,A_2,..A_k$，它们满足 $A_i \cap A_j = \emptyset$，且$A_1 \cup A_2 \cup ... \cup A_k = V$。
对于任意两个子图点集合 $A, B \subset V$ ，$A \cap B = \emptyset$，我们定义 $A$ 和 $B$ 之间的切图权重为：
$W(A, B) = \sum\limits_{i \in A, j \in B}w_{ij}$
那么对于 $k$ 个子图点的集合：$A_1,A_2,..A_k$，定义切图cut为：
$cut(A_1,A_2,...A_k) = \frac{1}{2}\sum\limits_{i=1}^{k}W(A_i, \overline{A}_i )$
其中 $\overline{A}_i $ 为 $A_i$ 的补集，意为除 $A_i$ 子集外其他 $V$ 的子集的并集。
每个子图就相当于聚类的一个类，找到子图内点的权重之和最高，子图间的点的权重之和最低的切图就相当于找到了最佳的聚类。实现这一点的一个很自然的想法是最小化 $cut$ 。然而这种方法存在问题，也就是最小化的 $cut$ 对应的切图不一定就是符合要求的最优的切图，如下图：
在上面的例子中，我们选择一个权重最小的边缘的点，比如 $C$ 和 $H$ 之间进行 $cut$，这样可以最小化 $cut(A_1,A_2,...A_k)$，但是却不是最优的切图，如何避免这种切图，并且找到类似图中 "Best Cut" 这样的最优切图呢？
接下介绍谱聚类使用的切图方法。
4 谱聚类之切图聚类
为避免最小切图导致的切图效果不佳，需要对每个子图的规模做出限定，一般来说，有两种切图方式，
- 第一种是 RatioCut
- 第二种是 Ncut
4.1 RatioCut 切图
RatioCut 切图为避免上述的最小切图，对每个切图，不光考虑最小化 $cut(A_1,A_2,...A_k)$ ，还考虑最大化每个子图点的个数，即：
$RatioCut(A_1,A_2,...A_k) = \frac{1}{2}\sum\limits_{i=1}^{k}\frac{W(A_i, \overline{A}_i )}{|A_i|}$
为最小化这个 RatioCut 函数，引入指示向量 $h_j \in \{h_1, h_2,..h_k\}\; j =1,2,...k$，对于任意一个向量 $h_j$， 它是一个 $n$ 维向量（$n$ 为样本数），定义 $h_{ij}$ 为：
$h_{ij}=\begin{cases}0& { v_i \notin A_j}\\\frac{1}{\sqrt{|A_j|}}& { v_i \in A_j}\end{cases}$
定义 $h_i^TLh_i$ 为：
$ \begin{array}{l} h_i^TLh_i & = \frac{1}{2}\sum\limits_{m=1}\sum\limits_{n=1}w_{mn}(h_{im}-h_{in})^2 \\& =\frac{1}{2}(\sum\limits_{m \in A_i, n \notin A_i}w_{mn}(\frac{1}{\sqrt{|A_i|}} - 0)^2 + \sum\limits_{m \notin A_i, n \in A_i}w_{mn}(0 - \frac{1}{\sqrt{|A_i|}} )^2\\& = \frac{1}{2}(\sum\limits_{m \in A_i, n \notin A_i}w_{mn}\frac{1}{|A_i|} + \sum\limits_{m \notin A_i, n \in A_i}w_{mn}\frac{1}{|A_i|}\\& = \frac{1}{2}(cut(A_i, \overline{A}_i) \frac{1}{|A_i|} + cut(\overline{A}_i, A_i) \frac{1}{|A_i|}) \\& = \frac{cut(A_i, \overline{A}_i)}{|A_i|} \end{array}$
可以看出，对于某一个子图 $i$，它的 RatioCut 对应于 $h_i^TLh_i$，那么 $k$个子图呢？对应的 RatioCut 函数表达式为：
$\begin{array}{l}RatioCut \left(A_{1}, A_{2}, \ldots A_{k}\right)&=\sum\limits _{i=1}^{k} h_{i}^{T} L h_{i}\\&=\sum\limits_{i=1}^{k}\left(H^{T} L H\right)_{i i} \\&=\operatorname{tr}\left(H^{T} L H\right)\end{array}$
上式中 $\operatorname{tr}\left(H^{T} L H\right) $ 为矩阵 $H^{T} L H$ 的迹， $H=\left(\begin{array}{llll}h_{1} & h_{2} & \cdots & h_{k}\end{array}\right) $ ，需要注意这里 的 $H$ 满足 $H^{T} H=I $ ，并且 $H$ 的元素只能取 $0$ 或者 $\frac{1}{\left|A_{i}\right|} $ 。
所以我们需要优化以下目标函数：
$\begin{array}{c} \underset{H}{\operatorname{argmin}} \quad \operatorname{tr}\left(H^{T} L H\right) \\\text { s.t. } H^{T} H=I\end{array}$
注意到 $H$ 矩阵里面的每一个指示向量都是 $n$ 维的，向量中每个变量的取值为 $0$ 或者 $\frac{1}{\sqrt{|A_j|}}$ ，就有 $2^n$ 种取值，有 $k$ 个子图的话就有 $k$ 个指示向量，共有 $k2^n$ 种 $H$，因此找到满足上面优化目标的 $H$ 是一个 $NP$ 难的问题。
观察到 $tr(H^TLH)$ 中每一个优化子目标 $h_i^TLh_i$，其中 $h$ 是单位正交基， $L$ 为对称矩阵，此时 $h_i^TLh_i$ 的最大值为 $L$ 的最大特征值，最小值是 $L$ 的最小特征值。类比于 PCA ，我们的目标是找到协方差矩阵（对应此处的拉普拉斯矩阵 $L$ ）的最大的特征值，而在我们的谱聚类中，我们的目标是找到目标的最小的特征值，得到对应的特征向量，此时对应二分切图效果最佳。也就是说，我们这里要用到维度规约的思想来近似去解决这个NP难的问题。
对于 $h_i^TLh_i$ ，我们的目标是找到最小的 $L$ 的特征值，而对于 $tr(H^TLH) = \sum\limits_{i=1}^{k}h_i^TLh_i$，则我们的目标就是找到 $k$ 个最小的特征值，一般来说，$k$ 远远小于 $n$，也就是说，此时我们进行了维度规约，将维度从 $n$ 降到了 $k$ ，从而近似可以解决这个NP难的问题。
通过找到 $L$ 的最小的 $k$ 个特征值，可以得到对应的 $k$ 个特征向量，这 $k$ 个特征向量组成一个 $nxk$ 维度的矩阵，即为我们的 $H$。一般需要对 $H$ 矩阵按行做标准化，即
$h_{ij}^{*}= \frac{h_{ij}}{(\sum\limits_{t=1}^kh_{it}^{2})^{1/2}}$
由于我们在使用维度规约的时候损失了少量信息，导致得到的优化后的指示向量 $h$ 对应的 $H$ 现在不能完全指示各样本的归属，因此一般在得到 $n\times k$ 维度的矩阵 $H$ 后还需要对每一行进行一次传统的聚类，比如使用 K-Means 聚类。
4.2 Ncut 切图
Ncut 切图和 RatioCut 切图很类似，但是把 Ratiocut 的分母 $|Ai|$ 换成 $vol(A_i)$ 。由于子图样本的个数多并不一定权重就大，我们切图时基于权重也更合我们的目标，因此一般来说 Ncut 切图优于 RatioCut 切图。 $NCut(A_1,A_2,...A_k) = \frac{1}{2}\sum\limits_{i=1}^{k}\frac{W(A_i, \overline{A}_i )}{vol(A_i)}$
对应的，Ncut 切图对指示向量 $h$ 做了改进。注意到 RatioCut 切图的指示向量使用的是 $\frac{1}{\sqrt{|A_j|}}$ 标示样本归属，而 Ncut 切图使用了子图权重 $\frac{1}{\sqrt{vol(A_j)}}$ 来标示指示向量 $h$，定义如下:
$h_{ij}=\begin{cases}0& { v_i \notin A_j}\\\frac{1}{\sqrt{vol(A_j)}}& { v_i \in A_j}\end{cases}$
那么我们对于 $h_i^TLh_i$ 有：
$ \begin{array}{l} h_i^TLh_i & = \frac{1}{2}\sum\limits_{m=1}\sum\limits_{n=1}w_{mn}(h_{im}-h_{in})^2 \\& =\frac{1}{2}(\sum\limits_{m \in A_i, n \notin A_i}w_{mn}(\frac{1}{\sqrt{vol(A_i)}} - 0)^2 + \sum\limits_{m \notin A_i, n \in A_i}w_{mn}(0 - \frac{1}{\sqrt{vol(A_i)}} )^2\\& = \frac{1}{2}(\sum\limits_{m \in A_i, n \notin A_i}w_{mn}\frac{1}{vol(A_i)} + \sum\limits_{m \notin A_i, n \in A_i}w_{mn}\frac{1}{vol(A_i)}\\& = \frac{1}{2}(cut(A_i, \overline{A}_i) \frac{1}{vol(A_i)} + cut(\overline{A}_i, A_i) \frac{1}{vol(A_i)}) \\& = \frac{cut(A_i, \overline{A}_i)}{vol(A_i)} \end{array}$
推导方式和 RatioCut 完全一致。也就是说，我们的优化目标仍然是
$\begin{array}{l}\operatorname{NCut}\left(A_{1}, A_{2}, \cdots, A_{k}\right)&=\sum\limits _{i=1}^{k} h_{i}^{T} L h_{i} \\&=\sum\limits_{i=1}^{k}\left(H^{T} L H\right)_{i i} \\&=\operatorname{tr}\left(H^{T} L H\right)\end{array}$
但是此时我们的 $H^TH \neq I$，而是 $H^TDH = I$。推导如下：
$\begin{array}{l}H^{T} D H&=\left(\begin{array}{c}h_{1}^{T} \\h_{2}^{T} \\\vdots \\h_{k}^{T}\end{array}\right)\left(\begin{array}{ccccc}d_{1} & & & \\& d_{2} & & \\& & \ddots & \\& & & d_{N}\end{array}\right)\left(\begin{array}{cccc}h_{1} & h_{2} & \cdots & h_{k}\end{array}\right) \\&=\left(\begin{array}{cccc}h_{11} d_{1} & h_{12} d_{2} & \cdots & h_{1 N} d_{N} \\h_{21} d_{1} & h_{22} d_{2} & \cdots & h_{2 N} d_{N} \\\vdots & \vdots & \ddots & \vdots \\h_{k 1} d_{1} & h_{k 2} d_{2} & \cdots & h_{k N} d_{N}\end{array}\right)\left(\begin{array}{cccc}h_{1} & h_{2} & \cdots & h_{k}\end{array}\right) \\&=\left(\begin{array}{ccccc}\sum\limits _{i=1}^{N} h_{1 i}^{2} d_{i} & \sum\limits_{i=1}^{N} h_{1 i} h_{2 i} d_{i} & \cdots & \sum\limits_{i=1}^{N} h_{1 i} h_{k i j} d_{i} \\\sum\limits_{i=1}^{N} h_{2 i} h_{1 i} d_{i} & \sum\limits_{i=1}^{N} h_{2 i}^{2} d_{i} & \cdots & \sum\limits_{i=1}^{N} h_{2 i} h_{k i} d_{i} \\\vdots & \vdots&\ddots & \vdots \\\sum\limits_{i=1}^{N} h_{k i} h_{1 i} d_{i} & \sum\limits_{i=1}^{N}h_{k i} h_{2 i} d_{i} & \cdots & \sum\limits_{i=1}^{N} h_{k i}^{2} d_{i}\end{array}\right)\end{array}$
对于对角线元素
$ h_i^TDh_i = \sum\limits_{j=1}^{n}h_{ij}^2d_j =\frac{1}{vol(A_i)}\sum\limits_{j \in A_i}d_j= \frac{1}{vol(A_i)}vol(A_i) =1$
由于 $h_{m i}$ 和 $h_{n i}$ 不可能同时非零（$v_i \in V_m \quad and \quad v_i \in V_n$ ），因此对于非对角线元素有:
$\sum \limits _{i=1}^{N} h_{m i} h_{n i} d_{i}=\sum \limits _{i=1}^{N} 0 \cdot d_{i}=0$
此时我们的优化目标最终为：
$\begin {array} {c} \underset{H}{arg\;min} \; tr(H^TLH) \\s.t.\;H^TDH=I \end {array} $
此时我们的 $H$ 中的指示向量 $h$ 并不是标准正交基，所以在 RatioCut里面的降维思想不能直接用。怎么办呢？其实只需要将指示向量矩阵 $H$ 做一个小小的转化即可。
我们令$H = D^{-1/2}F$, 则：
$H^TLH = F^TD^{-1/2}LD^{-1/2}F$
$H^TDH=F^TF = I$
此时优化目标为:
$\begin {array}{c}\underset{F}{\operatorname{arg\;\;min}}\;\; \operatorname{tr}\left(F^{T} D^{-1 / 2} L D^{-1 / 2} F\right) \\\text { s.t. } F^{T} F=I\end {array}$
可以发现这个式子和 RatioCut 基本一致，只是中间的 $L$ 变成了 $D^{-1/2}LD^{-1/2}$。这样我们就可以继续按照RatioCut的思想，求出 $D^{-1/2}LD^{-1/2}$ 的最小的前 $k$ 个特征值，然后求出对应的特征向量，并标准化，得到最后的特征矩阵 $F$，最后对 $F$ 进行一次传统的聚类（比如K-Means）即可。
一般来说， $D^{-1/2}LD^{-1/2}$ 相当于对拉普拉斯矩阵 $L$ 做了一次标准化，即
$\left(D^{-1 / 2} L D^{-1 / 2}\right)_{i j}=\frac{L_{i j}}{\sqrt{d_{i} * d_{j}}} $
5 谱聚类算法流程
5.1 Ncut 谱聚类算法流程
下面 Ncut 谱聚类算法流程。
输入：样本集D=$(x_1,x_2,...,x_n)$，相似矩阵的生成方式, 降维后的维度$k_1$, 聚类方法，聚类后的维度$k_2$
输出： 簇划分$C(c_1,c_2,...c_{k_2})$.
- 根据输入的相似矩阵的生成方式构建样本的相似矩阵 $S$ ；
- 根据相似矩阵 $S$ 构建邻接矩阵 $W$，构建度矩阵 $D$；
- 计算出拉普拉斯矩阵 $L$；
- 构建标准化后的拉普拉斯矩阵$D^{-1/2}LD^{-1/2}$；
- 计算$D^{-1/2}LD^{-1/2}$最小的 $k_1$ 个特征值所各自对应的特征向量 $f$；
- 将各自对应的特征向量$f$组成的矩阵按行标准化，最终组成 $n \times k_1$ 维的特征矩阵 $F$；
- 对 $F$ 中的每一行作为一个$k_1$维的样本，共 $n$ 个样本，用输入的聚类方法进行聚类，聚类维数为$k_2$ ；
- 得到簇划分$C(c_1,c_2,...c_{k_2})$.
5.2 谱聚类算法的优缺点
谱聚类算法的主要优点有：
- 谱聚类只需要数据之间的相似度矩阵，因此对于处理稀疏数据的聚类很有效。这点传统聚类算法比如K-Means很难做到
- 由于使用了降维，因此在处理高维数据聚类时的复杂度比传统聚类算法好。
谱聚类算法的主要缺点有：
- 如果最终聚类的维度非常高，则由于降维的幅度不够，谱聚类的运行速度和最后的聚类效果均不好。
- 聚类效果依赖于相似矩阵，不同的相似矩阵得到的最终聚类效果可能很不同。
因上求缘，果上努力~~~~ 作者：图神经网络，转载请注明原文链接：https://www.cnblogs.com/BlairGrowing/p/15814994.html
