谱聚类原理总结

原文链接

\qquad 谱聚类（spectral clustering）是广泛使用的聚类算法，比起传统的K-Means算法，谱聚类对数据分布的适应性更强，聚类效果也很优秀，同时聚类的计算量也小很多，更加难能可贵的是实现起来也不复杂。在处理实际的聚类问题时，个人认为谱聚类是应该首先考虑的几种算法之一。下面我们就对谱聚类的算法原理做一个总结。

1. 谱聚类概述

\qquad 谱聚类是从图论中演化出来的算法，后来在聚类中得到了广泛的应用。它的主要思想是把所有的数据看做空间中的点，这些点之间可以用边连接起来。距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高，通过对所有数据点组成的图进行切图，让切图后不同的子图间边权重和尽可能的低，而子图内的边权重和尽可能的高，从而达到聚类的目的。

2. 谱聚类基础之一：无向权重图

\qquad 对于一个图$G$,我们一般用点的集合$V$和边的集合$E$来描述。即$G(V,E)$,其中$V$即为我们数据集里面所有的点$(v_1,v_2,\dots ,v_n)$,对于$V$中的任意两个点，可以有边连接，也可以没有边连接。我们定义权重$w_{ij}$为点$v_i$和$v_j$之间的权重，由于我们是无向图，所以$w_{ij}=w_{ji}$。

\qquad 对于有边连接的两个点$v_i$和$v_j$，$w_{ij}>0$,对于没有边连接的两个点$v_i$和$v_j$,$w_{ij}=0$,对于图中的任意一个点$v_i$，它的度$d_i$定义为和它相连的所有边的权重之和，即
$$d_i=\sum _{j=1}^n{w}_{ij}$$
\qquad 利用每个点度的定义，我们可以得到一个$n\ast n$的度矩阵$D$,它是一个对角矩阵，只有主对角线有值，对应第$i$行的第$i$个点的度数，定义如下：
$$D=\left(\begin{array}{cccc}{d}_1& \dots & \dots & \dots \\ \dots & d_2& \dots & \dots \\ ⋮& ⋮& \ddots & ⋮\\ \dots & \dots & \dots & d_n\end{array}\right)$$

\qquad 利用所有点之间的权重值，我们可以得到图的邻接矩阵$W$，它也是一个$n\ast n$的矩阵，第$i$行的第$j$个值对应我们的权重$w_{ij}$。

\qquad 除此之外，对于点集$V$的的一个子集$A\subset V$，我们定义：
$$\begin{gathered} |A|:=子集A中点的个数 \\ vol(A)=\sum _{i\in A}{d}_i \end{gathered}$$

3. 谱聚类基础之二：相似矩阵

\qquad 邻接矩阵$W$，它是由任意两点之间的权重值$w_{ij}$组成的矩阵。通常我们可以自己输入权重，但是在谱聚类中，我们只有数据点的定义，并没有直接给出这个邻接矩阵，那么怎么得到这个邻接矩阵呢？

\qquad 基本思想是，距离较远的两个点之间的边权重值较低，而距离较近的两个点之间的边权重值较高，不过这仅仅是定性，我们需要定量的权重值。一般来说，我们可以通过样本点距离度量的相似矩阵$S$来获得邻接矩阵$W$。

\qquad 构建邻接矩阵$W$的方法有三类。$ϵ$-邻近法，K邻近法和全连接法。

\qquad 对于$ϵ$-邻近法，它设置了一个距离阈值$ϵ$，然后用欧式距离$s_{ij}$度量任意两点$x_i$和$x_j$的距离。即相似矩阵的$s_{ij}=||x_i-x_j|{|}^2$,然后根据$s_{ij}$和$ϵ$的大小关系，来定义邻接矩阵$W$如下：
$$W=\{\begin{array}{ll}0& s_{ij}>ϵ\\ ϵ& s_{ij}\le ϵ\end{array}$$

\qquad 从上式可见，两点间的权重要不就是$ϵ$,要不就是0，没有其他的信息了。距离远近度量很不精确，因此在实际应用中，我们很少使用$ϵ$-邻近法。

\qquad 第二种定义邻接矩阵$W$的方法是K邻近法，利用KNN算法遍历所有的样本点，取每个样本最近的k个点作为近邻，只有和样本距离最近的k个点之间的$w_{ij}>0$。但是这种方法会造成重构之后的邻接矩阵W非对称，我们后面的算法需要对称邻接矩阵。为了解决这种问题，一般采取下面两种方法之一：

\qquad 第一种K邻近法是只要一个点在另一个点的K近邻中，则保留$s_{ij}$
$$w_{ij}=w_{ji}=\{\begin{array}{ll}0& x_i\notin KNN(x_j)and{x}_j\notin KNN(x_i)\\ exp(-\frac{||x_i-x_j|{|}^2}{2{\sigma }^2})& x_i\in KNN(x_j)or{x}_j\in KNN(x_i)\end{array}$$

\qquad 第二种K邻近法是必须两个点互为K近邻中，才能保留$s_{ij}$
$$w_{ij}=w_{ji}=\{\begin{array}{ll}0& x_i\notin KNN(x_j)or{x}_j\notin KNN(x_i)\\ exp(-\frac{||x_i-x_j|{|}^2}{2{\sigma }^2})& x_i\in KNN(x_j)and{x}_j\in KNN(x_i)\end{array}$$

\qquad 第三种定义邻接矩阵$W$的方法是全连接法，相比前两种方法，第三种方法所有的点之间的权重值都大于0，因此称之为全连接法。可以选择不同的核函数来定义边权重，常用的有多项式核函数，高斯核函数和Sigmoid核函数。最常用的是高斯核函数RBF，此时相似矩阵和邻接矩阵相同：
$$w_{ij}=s_{ij}=exp(-\frac{||x_i-x_j|{|}^2}{2{\sigma }^2})$$

\qquad 在实际的应用中，使用第三种全连接法来建立邻接矩阵是最普遍的，而在全连接法中使用高斯径向核RBF是最普遍的。

4. 谱聚类基础之三：拉普拉斯矩阵

\qquad 单独把拉普拉斯矩阵(Graph Laplacians)拿出来介绍是因为后面的算法和这个矩阵的性质息息相关。它的定义很简单，拉普拉斯矩阵$L=D-W$。$D$即为我们第二节讲的度矩阵，它是一个对角矩阵。而$W$即为我们第二节讲的邻接矩阵，它可以由我们第三节的方法构建出。

\qquad 拉普拉斯矩阵有一些很好的性质如下：

\qquad 1）拉普拉斯矩阵是对称矩阵，这可以由$D$和$W$都是对称矩阵而得。

\qquad 2）由于拉普拉斯矩阵是对称矩阵，则它的所有的特征值都是实数。

\qquad 3）对于任意的向量$f$,我们有:
$$f^{T}Lf=\frac{1}{2}\sum _{i,j=1}^n{w}_{ij}(f_i-f_j{)}^2$$
这个利用拉普拉斯矩阵的定义很容易得到如下：

$D$矩阵只有对角线有元素，$W$矩阵对角线没有元素
$$\begin{array}{rl}{f}^{T}Lf& =f^{T}Df-f^{T}Wf=\sum _{i=1}^n{d}_i{f}_i^2-\sum _{i,j=1}^n{w}_{ij}{f}_i{f}_j\\ & =\frac{1}{2}(\sum _{i=1}^n{d}_i{f}_i^2-2\sum _{i,j=1}^n{w}_{ij}{f}_i{f}_j+\sum _{j=1}^n{d}_j{f}_j^2)\\ & =\frac{1}{2}\sum _{i,j=1}^n{w}_{ij}(f_i-f_j{)}^2\end{array}$$

\qquad 4) 拉普拉斯矩阵是半正定的，且对应的$n$个实数特征值都大于等于$0$，即$0\le {\lambda }_1\le {\lambda }_2\le \cdots \le {\lambda }_n$， 且最小的特征值为$0$，这个由性质3很容易得出。

5. 谱聚类基础之四：无向图切图

\qquad 对于无向图$G$的切图，我们的目标是将图$G(V,E)$切成相互没有连接的k个子图,每个子图点的集合为：$A_1,A_2,\dots ,A_k$,它们满足$A_i\bigcap A_j=\varnothing$,且$A_i\bigcup A_2\bigcup \cdots \bigcup A_k=V$。

\qquad 对于任意两个子图点的集合$A,B\subset V,A\bigcap B=\varnothing$,我们定义A和B之间的切图权重为：
$$W(A,B)=\sum _{i\in A,j\in B}{w}_{ij}$$
\qquad 那么对于我们k个子图点的集合：$A_1,A_2,\dots ,A_k$,我们定义切图cut为：
$$cut(A_1,A_2,\dots ,A_k)=\frac{1}{2}\sum _{i=1}^{k}W(A_i,\overline{A_i})$$
\qquad 其中$\overline{A_i}$为$A_i$的补集，意为除$A_i$集外其他$V$的子集的并集。

\qquad 那么如何切图可以让子图内的点权重和高，子图间的点权重和低呢？一个自然的想法就是最小化$cut(A_1,A_2,\dots ,A_k)$, 但是可以发现，这种极小化的切图存在问题，如下图：

\qquad 我们选择一个权重最小的边缘的点，比如C和H之间进行cut，这样可以最小化$cut(A_1,A_2,\dots ,A_k)$,但是却不是最优的切图，如何避免这种切图，并且找到类似图中"Best Cut"这样的最优切图呢？我们下一节就来看看谱聚类使用的切图方法。

6. 谱聚类之切图聚类

\qquad 为了避免最小切图导致的切图效果不佳，我们需要对每个子图的规模做出限定，一般来说，有两种切图方式，第一种是RatioCut，第二种是Ncut。下面我们分别加以介绍。

6.1 RatioCut切图

\qquad RatioCut切图为了避免第五节的最小切图，对每个切图，不光考虑最小化$cut(A_1,A_2,\dots ,A_k)$，它还同时考虑最大化每个子图点的个数，即：
$$RatioCut(A_1,A_2,\dots ,A_k)=\frac{1}{2}\sum _{i=1}^k\frac{W(A_i,\overline{A_i})}{|A_i|}$$

\qquad 那么怎么最小化这个RatioCut函数呢？牛人们发现，RatioCut函数可以通过如下方式表示。

\qquad 我们引入指示向量 h j ∈ { h 1 , h 2 , … , h k }   j = 1 , 2 , … , k h_j \in \{h_1,h_2,\dots,h_k\} \ j = 1,2,\dots,k hj​∈{h1​,h2​,…,hk​} j=1,2,…,k,对于任意一个向量$h_j$,它是一个n维向量（n为样本数），我们定义$h_{ij}$为：
$$h_{ij=}\{\begin{array}{ll}0& v_i\notin A_j\\ \frac{1}{\sqrt{|A_j|}}& v_i\in A_j\end{array}$$

\qquad 那么我们对于$h_i^{T}L{h}_i$有,：
$$\begin{array}{rl}{h}_i^{T}L{h}_i& =\frac{1}{2}\sum _{m=1}\sum _{n=1}{w}_{mn}(h_{mi}-h_{ni}{)}^2\\ & =\frac{1}{2}[\sum _{m\in A_i,n\notin A_i}{w}_{mn}(\frac{1}{\sqrt{|A_i|}}-0{)}^2+\sum _{m\notin A_i,n\in A_i}{w}_{mn}(0-\frac{1}{\sqrt{|A_i|}}{)}^2]\\ & =\frac{1}{2}(\sum _{m\in A_i,n\notin A_i}{w}_{mn}\frac{1}{|A_i|}+\sum _{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}$$

\qquad 上述第一式用了上面第四节的拉普拉斯矩阵的性质3.
$$f^{T}Lf=\frac{1}{2}\sum _{i,j=1}^n{w}_{ij}(f_i-f_j{)}^2$$
\qquad 第二式用到了指示向量的定义，另外第三式用到了上述第5节$cut$,
$$cut(A_1,A_2,\dots ,A_k)=\frac{1}{2}\sum _{i=1}^{k}W(A_i,\overline{A_i})$$
$$\begin{array}{rl}cut(A_i,\overline{A_i})& =\frac{1}{2}(W(A_i,\overline{A_i})+W(\overline{A_i},A_i))\\ & =\frac{1}{2}(\sum _{m\in A_i,n\notin A_i}{w}_{mn}+\sum _{m\notin A_i,n\in A_i}{w}_{mn})\\ & =W(A_i,\overline{A_i})\end{array}$$
$$cut(\overline{A_i},A_i)=\frac{1}{2}(\sum _{m\notin A_i,n\in A_i}{w}_{mn}+\sum _{m\in A_i,n\notin A_i}{w}_{mn})$$

\qquad 可以看出，对于某一个子图$i$,它的RatioCut对应于$h_i^{T}L{H}_i$。
$$RatioCut(A_1)=\frac{1}{2}\frac{W(A_1,\overline{A_1})}{|A_1|}\propto \frac{cut(A_i,\overline{A_i})}{|A_i|}=h_i^{T}L{H}_i$$
$$RatioCut(A_1,A_2,\dots ,A_k)=\sum _{i=1}^k{h}_i^{T}L{h}_i=\sum _{i=1}^k(H^{T}LH{)}_{ii}=tr(H^{T}LH)$$

\qquad 其中，$tr(H^{T}LH)$是矩阵的迹，即对角线的求和，$H$是归一化的正交基矩阵。也就是说，我们的RatioCut切图，实际上就是最小化我们的$tr(H^{T}LH)$，注意到$H^{T}H=I$，则我们的切图优化目标为：
$$\underset{H}{\underbrace{argmin}}tr(H^{T}LH)s.t.H^{T}H=I$$
\qquad 注意到我们H矩阵里面的每一个指示向量都是n维的，向量中每个变量的取值为0或者$\frac{1}{\sqrt{|A_i|}}$有$2^n$种取值，有k个子图的话就有k个指示向量，共有$k{2}^n$种H，因此找到满足上面优化目标的H是一个NP难的问题。

\qquad 注意观察$tr(H^{T}LH)$中每一个优化子目标$h_i^{T}L{h}_i$，其中$h$是单位正交基， L为对称矩阵，此时$h_i^{T}L{h}_i$的最大值为L的最大特征值，最小值是L的最小特征值。如果你对主成分分析PCA很熟悉的话，这里很好理解。在PCA中，我们的目标是找到协方差矩阵（对应此处的拉普拉斯矩阵L）的最大的特征值，而在我们的谱聚类中，我们的目标是找到目标的最小的特征值，得到对应的特征向量，此时对应二分切图效果最佳。也就是说，我们这里要用到维度规约的思想来近似去解决这个NP难的问题。

\qquad 对于$h_i^{T}L{h}_i$，我们的目标是找到最小的L的特征值，而对于$tr(H^{T}LH)=\sum _{i=1}^n{h}_i^{T}L{h}_i$，则我们的目标就是找到k个最小的特征值，一般来说，k远远小于n，也就是说，此时我们进行了维度规约，将维度从n降到了k，从而近似可以解决这个NP难的问题。

\qquad 通过找到L的最小的k个特征值，可以得到对应的k个特征向量，这k个特征向量组成一个nxk维度的矩阵，即为我们的H。一般需要对H矩阵按行做标准化，即:
$$h_{ij}^{\ast }=\frac{h_{ij}}{\sqrt{{\sum }_{t=1}^k{h}_{it}^2}}$$
\qquad 由于我们在使用维度规约的时候损失了少量信息，导致得到的优化后的指示向量h对应的H现在不能完全指示各样本的归属，因此一般在得到nxk维度的矩阵H后还需要对每一行进行一次传统的聚类，比如使用K-Means聚类.

6.2 Ncut切图

\qquad Ncut切图和RatioCut切图很类似，但是把Ratiocut的分母$|A_i|$换成
$vol(A_i)$。 由于子图样本的个数多并不一定权重就大，我们切图时基于权重也更合我们的目标，因此一般来说Ncut切图优于RatioCut切图。
$$Ncut(A_1,A_2,\dots ,A_k)=\frac{1}{2}\sum _{i=1}^k\frac{W(A,A_i)}{vol(A_i)}$$
\qquad 对应的，Ncut切图对指示向量$h$做了改进。注意到RatioCut切图的指示向量使用的是$\frac{1}{\sqrt{|A_j|}}$标示样本归属，而Ncut切图使用了子图权重$\frac{1}{\sqrt{vol(A_j)}}$来标示指示向量$h$，定义如下:
$$h_{ij}=\{\begin{array}{ll}0& v_i\notin A_j\\ \frac{1}{\sqrt{vol(A_j)}}& v_i\in A_j\end{array}$$
\qquad 那么我们对于$h_i^{T}L{h}_i$,有：
$$\begin{array}{rl}{h}_i^{T}L{h}_i& =\frac{1}{2}\sum _{m=1}\sum _{n=1}{w}_{mn}(h_{mi}-h_{ni}{)}^2\\ & =\frac{1}{2}[\sum _{m\in A_i,n\notin A_i}{w}_{mn}(\frac{1}{\sqrt{vol(A_i)}}-0{)}^2+\sum _{m\notin A_i,n\in A_i}{w}_{mn}(0-\frac{1}{\sqrt{vol(A_i)}}{)}^2]\\ & =\frac{1}{2}(\sum _{m\in A_i,n\notin A_i}{w}_{mn}\frac{1}{vol{A}_i}+\sum _{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}$$

\qquad 推导方式和RatioCut完全一致。也就是说，我们的优化目标仍然是
$$Ncut(A_1,A_2,\dots ,A_k)=\sum _{i=1}^k{h}_i^{T}L{h}_i=\sum _{i=1}^k(H^{T}LH{)}_{ii}=tr(H^{T}LH)$$
\qquad 但是此时我们的$H^{T}H̸=I$，而是$H^{T}DH=I$。推导如下：
$$h_i^{T}D{h}_i=\sum _{j=1}^n{h}_{ij}^2{d}_j=\frac{1}{vol(A_i)}\sum _{j\in A_i}{d}_j=\frac{1}{vol(A_i)}vol(A_i)=1$$
\qquad 也就是说，此时我们的优化目标最终为：
$$\underset{H}{\underbrace{argmin}}tr(H^{T}LH)s.t.H^{T}DH=I$$
\qquad 此时我们的H中的指示向量$h$并不是标准正交基，所以在RatioCut里面的降维思想不能直接用。怎么办呢？其实只需要将指示向量矩阵H做一个小小的转化即可。

\qquad 我们令$H=D^{-1/2}F$, 则：$H^{T}LH=F^T{D}^{-1/2}D{D}^{-1/2}{F}^T,H^{T}DH=F^{T}F=I$，也就是说优化目标变成了:
$$\underset{H}{\underbrace{argmin}}tr(F^T{D}^{-1/2}D{D}^{-1/2}{F}^T)s.t.F^{T}F=I$$
\qquad 可以发现这个式子和RatioCut基本一致，只是中间的$L$变成了$D^{-1/2}D{D}^{-1/2}$。这样我们就可以继续按照RatioCut的思想，求出$D^{-1/2}D{D}^{-1/2}$的最小的前k个特征值，然后求出对应的特征向量，并标准化，得到最后的特征矩阵$F$,最后对$F$进行一次传统的聚类（比如K-Means）即可。

\qquad 一般来说， $D^{-1/2}D{D}^{-1/2}$,相当于对拉普拉斯矩阵$L$做了一次标准化，即$\frac{L_{ij}}{\sqrt{d_i\ast d_j}}$

7. 谱聚类算法流程

\qquad 总结下谱聚类的基本流程了。一般来说，谱聚类主要的注意点为相似矩阵的生成方式（参见第二节），切图的方式（参见第六节）以及最后的聚类方法（参见第六节）。

\qquad 最常用的相似矩阵的生成方式是基于高斯核距离的全连接方式，最常用的切图方式是Ncut。而到最后常用的聚类方法为K-Means。下面以Ncut总结谱聚类算法流程。

\qquad 输入：样本集$D=(x_1,x_2,\dots ,x_n)$,相似矩阵的生成方式, 降维后的维度$k_1$, 聚类方法，聚类后的维度$k_2$
\qquad 输出： 簇划分$C(c_1,c_2,\dots ,c_{k_2})$

\qquad 1) 根据输入的相似矩阵的生成方式构建样本的相似矩阵S
\qquad 2）根据相似矩阵S构建邻接矩阵W，构建度矩阵D
\qquad 3）计算出拉普拉斯矩阵L
\qquad 4）构建标准化后的拉普拉斯矩阵$D^{-1/2}D{D}^{-1/2}$
\qquad 5）计算$D^{-1/2}D{D}^{-1/2}$最小的$k_1$个特征值所各自对应的特征向量$f$
\qquad 6) 将各自对应的特征向量$f$组成的矩阵按行标准化，最终组成$n\ast k_1$维的特征矩阵$F$
\qquad 7）对F中的每一行作为一个$k_1$维的样本，共n个样本，用输入的聚类方法进行聚类，聚类维数为$k_2$。
\qquad 8）得到簇划分$C(c_1,c_2,\dots ,c_{k_2})$

8. 谱聚类算法总结

\qquad 谱聚类算法是一个使用起来简单，如果你掌握了谱聚类，相信你会对矩阵分析，图论有更深入的理解。同时对降维里的主成分分析也会加深理解。

\qquad 下面总结下谱聚类算法的优缺点。

\qquad 谱聚类算法的主要优点有：

\qquad 1）谱聚类只需要数据之间的相似度矩阵，因此对于处理稀疏数据的聚类很有效。这点传统聚类算法比如K-Means很难做到

\qquad 2）由于使用了降维，因此在处理高维数据聚类时的复杂度比传统聚类算法好。

\qquad 谱聚类算法的主要缺点有：

\qquad 1）如果最终聚类的维度非常高，则由于降维的幅度不够，谱聚类的运行速度和最后的聚类效果均不好。

\qquad 2) 聚类效果依赖于相似矩阵，不同的相似矩阵得到的最终聚类效果可能很不同。