谱聚类（spectral clustering）原理总结

文章目录

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

一、谱聚类概述

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

乍一看，这个算法原理的确简单，但是要完全理解这个算法的话，需要对图论中的无向图，线性代数和矩阵分析都有一定的了解。下面我们就从这些需要的基础知识开始，一步步学习谱聚类。

二、谱聚类基础

2.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$。对于图中的任意一个点vi，它的度di定义为和它相连的所有边的权重之和，即
$d_i=∑_{j=1}^nw_{ij}$
利用每个点度的定义，我们可以得到一个nxn的度矩阵D,它是一个对角矩阵，只有主对角线有值，对应第i行的第i个点的度数，定义如下：
$D = \begin{pmatrix} d_1 & \cdots & \cdots \\ \cdots & d_2 & \cdots \\ \vdots & \vdots & \ddots\\ \cdots & \cdots & d_n \end{pmatrix}$
利用所有点之间的权重值，我们可以得到图的邻接矩阵W，它也是一个nxn的矩阵，第i行的第j个值对应我们的权重 $w_{ij}$。
除此之外，对于点集$V$的的一个子集$A⊂V$，我们定义：
$\begin{aligned} |A|: =&子集A中点的个数 \\ vol(A): =& ∑_{i∈A}d_i \end{aligned}$

2.2 相似矩阵

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

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

$W$的方法有三类:

$ϵ$-邻近法，

K邻近法和

全连接法。

1. 对于$ϵ$-邻近法:
   它设置了一个距离阈值$ϵ$，然后 用欧式距离$s_{ij}$度量任意两点$x_i$和$x_j$的距离。即相似矩阵的$s_{ij}=||x_i−x_j||_2^2$, 然后根据$s_{ij}$ 和 $ϵ$的大小关系，来定义邻接矩阵 $W$如下：
   $w_{ij}=\left\{\begin{matrix} 0 & s_{ij}>ϵ\\ ϵ & s_{ij}≤ϵ \end{matrix}\right.$
   从上式可见，两点间的权重要不就是$ϵ$,要不就是0，没有其他的信息了。距离远近度量很不精确，因此在实际应用中，我们很少使用$ϵ$-邻近法。
2. K邻近法:
   利用KNN算法遍历所有的样本点，取每个样本最近的k个点作为近邻，只有和样本距离最近的k个点之间的 $w_{ij}$>0。但是这种方法会造成重构之后的邻接矩阵 $W$
   第一种K邻近法是只要一个点在另一个点的K近邻中，则保留$s_{ij}$
   $w_{ij}=w_{ji}=\begin{cases} 0 & x_i∉KNN(x_j) ~ \text{and} ~x_j∉KNN(x_i) \\ \exp(−\frac{||x_i−x_j||_2^2}{2σ^2}) & x_i∈KNN(x_j)~\text{or}~x_j∈KNN(x_i) \end{cases}$
   第二种K邻近法是必须两个点互为K近邻中，才能保留$s_{ij}$
   $w_{ij}=w_{ji}=\begin{cases} 0 & x_i∉KNN(x_j) ~ \text{or} ~x_j∉KNN(x_i) \\ \exp(−\frac{||x_i−x_j||_2^2}{2σ^2}) & x_i∈KNN(x_j)~\text{and}~x_j∈KNN(x_i) \end{cases}$
3. 全连接法
   相比前两种方法，第三种方法所有的点之间的权重值都大于0，因此称之为全连接法。可以选择不同的核函数来定义边权重，常用的有多项式核函数，高斯核函数和Sigmoid核函数。最常用的是高斯核函数RBF，此时相似矩阵和邻接矩阵相同：
   $w_{ij}=s_{ij}= \exp(−\frac{||x_i−x_j||_2^2}{2σ^2})$
   在实际的应用中，使用第三种全连接法来建立邻接矩阵是最普遍的，而在全连接法中使用高斯径向核RBF是最普遍的。

2.3 拉普拉斯矩阵

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

1）拉普拉斯矩阵是 对称矩阵，这可以由

$D$和

$W$都是对称矩阵而得。

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

$f$,我们有

$f^TLf=\frac{1}{2}∑_{i,j=1}^nw_{ij}(f_i−f_j)^2$

$\begin{aligned} f^TLf&=f^TDf−f^TWf=∑_{i=1}^nd_if_i^2−∑_{i,j=1}^nw_{ij}f_if_j \\ &=\frac{1}{2}\left ( ∑_{i=1}^nd_if_i^2 - 2 ∑_{i,j=1}^nw_{ij}f_if_j + ∑_{j=1}^nd_jf_j^2 \right ) \\ &=\frac{1}{2}∑_{i,j=1}^nw_{ij}(f_i−f_j)^2 \end{aligned}$

$0=λ_1≤λ_2≤...≤λ_n$， 且最小的特征值为 0，这个由性质3很容易得出。

2.4 无向图切图

对于无向图$G$的切图，我们的目标是将图$G(V,E)$切成相互没有连接的$k$个子图，每个子图点的集合为：$A_1,A_2,..A_k$，它们满足$A_i∩A_j=∅,且A_1∪A_2∪...∪A_k=V.$

对于任意两个子图点的集合$A,B⊂V, A∩B=∅,$我们定义$A$和$B$之间的切图权重为：$W(A,B)=∑_{i∈A,j∈B}w_{ij}$
那么对于我们k个子图点的集合：$A_1,A_2,..A_k$，我们定义切图cut为：$cut(A_1,A_2,..A_k)=\frac{1}{2}∑_{i=1}^kW(A_i,\bar{A}_i)$
其中$\bar{A}_i$为$A_i$的补集，意为除$A_i$子集外其他$V$的子集的并集。

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

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

三、谱聚类切图聚类

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

3.1 RatioCut切图

RatioCut切图为了避免第五节的最小切图，对每个切图，不光考虑最小化 $cut(A_1,A_2,..A_k)$，它还同时考虑最大化每个子图点的个数，即：
$RatioCut(A_1,A_2,..A_k)=\frac{1}{2}∑_{i=1}^k\frac{W(A_i,\bar{A}_i)}{|A_i|}$
那么怎么最小化这个RatioCut函数呢？RatioCut函数可以通过如下方式表示。
我们引入指示向量 $h_j∈\{h_1,h_2,..h_k\} j=1,2,...k,$对于任意一个向量$h_j$, 它是一个n维向量（n为样本数），我们定义$h_{ij}$为：
$h_{ij}=\begin{cases} 0 & v_i∉A_j\\ \frac{1}{\sqrt{|A_j|}} &v_i∈A_j \end{cases}$
那么我们对于$h_i^TLh_i$，有：
$\begin{aligned} h_i^TLh_i & =\frac{1}{2}∑_{m=1}∑_{n=1}w_{mn}(h_{im}−h_{in})^2 \\ &=\frac{1}{2}\left ( ∑_{m∈A_i,n∉A_i}w_{mn}(\frac{1}{\sqrt{|A_i|}}-0)^2 + ∑_{m∉A_i,n∈A_i}w_{mn}(0- \frac{1}{\sqrt{|A_i|}})^2 \right ) \\ &=\frac{1}{2} \left ( ∑_{m∈A_i,n∉A_i}w_{mn}\frac{1}{|A_i|} + ∑_{m∉A_i,n∈A_i}w_{mn}\frac{1}{|A_i|} \right ) \\ &=\frac{1}{2} \left (\text{cut}(A_i,\bar{A}_i)\frac{1}{|A_i|} + \text{cut}(\bar{A}_i,A_i)\frac{1}{|A_i|} \right ) \\ &=\frac{\text{cut}(A_i,\bar{A}_i)}{|A_i|} \end{aligned}\tag{1}$

上述公式第一行用了上面第四节的拉普拉斯矩阵的性质3. 第二行用到了指示向量的定义。可以看出，对于某一个子图$i$，它的RatioCut 对应于$h_i^TLh_i$,那么我们的k个子图呢？对应的RatioCut函数表达式为：
$RatioCut(A_1,A_2,..A_k)=∑_{i=1}^kh_i^TLh_i=∑_{i=1}^k(H^TLH)_{ii}=tr(H^TLH)$
其中 $tr(H^TLH)$ 为矩阵的迹。也就是说，我们的RatioCut切图，实际上就是最小化我们的$tr(H^TLH)$。注意到 $H^TH=I$,则我们的切图优化目标为:
$\underbrace{argmin}_{\text{H}}tr(H^TLH)\Rightarrow s.t.H^TH=I$
注意到我们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很熟悉的话，这里很好理解。在PCA中，我们的目标是找到协方差矩阵（对应此处的拉普拉斯矩阵L）的最大的特征值，而在我们的谱聚类中，我们的目标是找到目标的最小的特征值，得到对应的特征向量，此时对应二分切图效果最佳。也就是说，我们这里要用到维度规约的思想来近似去解决这个NP难的问题。

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

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

3.2 Ncut切图

Ncut切图和RatioCut切图很类似，但是把Ratiocut的分母$|A_i|$换成$vol(A_i)$. 由于子图样本的个数多并不一定权重就大，我们切图时基于权重也更合我们的目标，因此一般来说Ncut切图优于RatioCut切图。
$NCut(A_1,A_2,..A_k)=\frac{1}{2}∑_{i=1}^k\frac{W(A_i,\bar{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∉A_j\\ \frac{1}{\sqrt{vol(A_j)}} &v_i∈A_j \end{cases}$
那么我们对于$h_i^TLh_i$，有：
$\begin{aligned} h_i^TLh_i & =\frac{1}{2}∑_{m=1}∑_{n=1}w_{mn}(h_{im}−h_{in})^2 \\ &=\frac{1}{2}\left ( ∑_{m∈A_i,n∉A_i}w_{mn}(\frac{1}{\sqrt{vol(A_i)}}-0)^2 + ∑_{m∉A_i,n∈A_i}w_{mn}(0- \frac{1}{\sqrt{vol(A_i)}})^2 \right ) \\ &=\frac{1}{2} \left ( ∑_{m∈A_i,n∉A_i}w_{mn}\frac{1}{vol(A_i)} + ∑_{m∉A_i,n∈A_i}w_{mn}\frac{1}{vol(A_i)} \right ) \\ &=\frac{1}{2} \left (\text{cut}(A_i,\bar{A}_i)\frac{1}{vol(A_i)} + \text{cut}(\bar{A}_i,A_i)\frac{1}{vol(A_i)} \right ) \\ &=\frac{\text{cut}(A_i,\bar{A}_i)}{vol(A_i)} \end{aligned}\tag{1}$
推导方式和RatioCut完全一致。也就是说，我们的优化目标仍然是：
$NCut(A_1,A_2,..A_k)=∑_{i=1}^kh_i^TLh_i=∑_{i=1}^k(H^TLH)_{ii}=tr(H^TLH)$
但是此时我们的$H^TH \neq I$，而是 $H^TDH=I$。推导如下：
$h_i^TDh_i=∑_{j=1}^nh_{ij}^2d_j=\frac{1}{vol(A_i)}∑_{j∈A_i}d_j=\frac{1}{vol(A_i)}vol(A_i)=1$

也就是说，此时我们的优化目标最终为：$\underbrace{arg \min}_{\text{F}}tr(H^TLH)\Rightarrow s.t.H^TH=I$

此时我们的H中的指示向量h并不是标准正交基，所以在RatioCut里面的降维思想不能直接用。怎么办呢？其实只需要将指示向量矩阵H做一个小小的转化即可。

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

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

四、谱聚类算法流程

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

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

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

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

谱聚类算法总结
谱聚类算法是一个使用起来简单，但是讲清楚却不是那么容易的算法，它需要你有一定的数学基础。如果你掌握了谱聚类，相信你会对矩阵分析，图论有更深入的理解。同时对降维里的主成分分析也会加深理解。

谱聚类算法的主要优点有：
1）谱聚类只需要数据之间的相似度矩阵，因此对于处理稀疏数据的聚类很有效。这点传统聚类算法比如K-Means很难做到
2）由于使用了降维，因此在处理高维数据聚类时的复杂度比传统聚类算法好。

谱聚类算法的主要缺点有：
1）如果最终聚类的维度非常高，则由于降维的幅度不够，谱聚类的运行速度和最后的聚类效果均不好。
2）聚类效果依赖于相似矩阵，不同的相似矩阵得到的最终聚类效果可能很不同。

代码实现：

import numpy as np
from sklearn.cluster import KMeans


# (1) 距离矩阵
def euclidDistance(x1, x2, sqrt_flag=False):
    res = np.sum((x1-x2)**2)
    if sqrt_flag:
        res = np.sqrt(res)
    return res

def calEuclidDistanceMatrix(X):
    X = np.array(X)
    S = np.zeros((len(X), len(X)))
    for i in range(len(X)):
        for j in range(i+1, len(X)):
            S[i][j] = 1.0 * euclidDistance(X[i], X[j])
            S[j][i] = S[i][j]
    return S

# (2) 邻接矩阵
def myKNN(S, k, sigma=1.0):
    N = len(S)
    A = np.zeros((N,N))

    for i in range(N):
        dist_with_index = zip(S[i], range(N))
        dist_with_index = sorted(dist_with_index, key=lambda x:x[0])
        neighbours_id = [dist_with_index[m][1] for m in range(k+1)] # xi's k nearest neighbours

        for j in neighbours_id: # xj is xi's neighbour
            A[i][j] = np.exp(-S[i][j]/2/sigma/sigma)
            A[j][i] = A[i][j] # mutually
    return A

#（3）标准化的拉普拉斯矩阵
def calLaplacianMatrix(adjacentMatrix):

    # compute the Degree Matrix: D=sum(A)
    degreeMatrix = np.sum(adjacentMatrix, axis=1)

    # compute the Laplacian Matrix: L=D-A
    laplacianMatrix = np.diag(degreeMatrix) - adjacentMatrix

    # normailze
    # D^(-1/2) L D^(-1/2)
    sqrtDegreeMatrix = np.diag(1.0 / (degreeMatrix ** (0.5)))
    return np.dot(np.dot(sqrtDegreeMatrix, laplacianMatrix), sqrtDegreeMatrix)

# （4）特征值分解
lam, H = np.linalg.eig(Laplacian) # H'shape is n*n

# (5) Kmeans
def spKmeans(H):
    sp_kmeans = KMeans(n_clusters=2).fit(H)
    return sp_kmeans.labels_

​​https://github.com/SongDark/SpectralClustering​​

鸣谢

https://code.google.com/p/pspectralclustering/