四种常用聚类及代码（二）：谱聚类（spectral clustering）

谱聚类在最近几年变得受欢迎起来，主要原因就是它实现简单，聚类效果经常优于传统的聚类算法（如K-Means算法）。刚开始学习谱聚类的时候，给人的感觉就是这个算法看上去很难，但是当真正的深入了解这个算法的时候，其实它的原理并不难，但是理解该算法还是需要一定的数学基础的。如果掌握了谱聚类算法，会对矩阵分析，图论和降维中的主成分分析等有更加深入的理解。

1、谱聚类概述

谱聚类(Spectral Clustering, SC)是一种基于图论的聚类方法——将带权无向图划分为两个或两个以上的最优子图，使子图内部尽量相似，而子图间距离尽量距离较远，以达到常见的聚类的目的。其中的最优是指最优目标函数不同，可以是割边最小分割，也可以是分割规模差不多且割边最小的分割。
谱聚类通过对样本数据的拉普拉斯矩阵的特征向量进行聚类，从而对样本数据聚类。谱聚类可以理解为将高维空间的数据映射到低维，然后在低维空间用其它聚类算法（如K-Means）进行聚类。

1.1、简单说明

谱聚类过程主要有两步，第一步是构图，将采样点数据构造成一张网图，表示为G(V,E)，V表示图中的点，E表示点与点之间的边，如下图：
四种常用聚类及代码（二）：谱聚类（spectral clustering）
第二步是切图，即将第一步构造出来的按照一定的切边准则，切分成不同的图，而不同的子图，即我们对应的聚类结果，举例如下：

2、数学准备

2.1、谱

方阵作为线性算子，它的所有特征值的全体统称为方阵的谱。方阵的谱半径为最大的特征值。矩阵A的谱半径是矩阵A^TA的最大特征值。

2.2、无向权重图

对于一个图G，我们一般用点的集合V和边的集合E来描述。即为G(V,E)。其中V即为我们数据集里面所有的点(v1,v2,…vn)。对于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 = $\sum$ _j w_ij

利用每个点度的定义，我们可以得到一个nxn的度矩阵D,它是一个对角矩阵，只有主对角线有值，对应第i行的第i个点的度数。
利用所有点之间的权重值，我们可以得到图的邻接矩阵W，它也是一个nxn的矩阵，第i行的第j个值对应我们的权重w_ij。

2.3、相似矩阵

相似矩阵就是样本点中的任意两个点之间的距离度量，在聚类算法中可以表示为距离近的点它们之间的相似度比较高，而距离较远的点它们的相似度比较低，甚至可以忽略。这里用三种方式表示相似度矩阵：一是ϵ-近邻法（ϵ-neighborhood graph），二是k近邻法（k-nearest nerghbor graph），三是全连接法（fully connected graph）。下面我们来介绍这三种方法。

2.3.1、ϵ-neighborhood graph：

样本点中任意两点之间的欧式距离：
四种常用聚类及代码（二）：谱聚类（spectral clustering）
用此方法构造的相似度矩阵表示如下：

该相似度矩阵由于距离近的点的距离表示为，距离远的点距离表示为0，矩阵种没有携带关于数据集的太多的信息，所以该方法一般很少使用，在sklearn中也没有使用该方法。

2.3.2、k-nearest nerghbor graph：

由于每个样本点的k个近邻可能不是完全相同的，所以用此方法构造的相似度矩阵并不是对称的。因此，这里使用两种方式表示对称的knn相似度矩阵，第一种方式是如果v_i在v_j的k个邻域中或者 v_j在v_i的k个邻域中，则w_ij=w_ji为v_i与v_j之间的距离，否则为0；第二种方式是如果v_i在v_j的k个邻域中并且 v_j在v_i的k个邻域中，则w_ij=w_ji为v_i与v_j之间的距离，否则为0。很显然第二种方式比第一种方式生成的相似度矩阵要稀疏。这两种方式用公式表达如下：

第一种方式：
四种常用聚类及代码（二）：谱聚类（spectral clustering）
第二种方式：

2.3.3、fully connected graph:

该方法就是在算法描述中的高斯相似度方法，公式如下：
四种常用聚类及代码（二）：谱聚类（spectral clustering）
该方法也是最常用的方法，在sklearn中默认的也是该方法，表示任意两个样本点都有相似度，但是距离较远的样本点之间相似度较低，甚至可以忽略。这里面的参数控制着样本点的邻域宽度，即越大表示样本点与距离较远的样本点的相似度越大，反之亦然。

2.4、拉普拉斯矩阵

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

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

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

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

3）对于任意的向量f,我们有
$f^TLf = \frac{1}{2}\sum\limits_{i,j=1}^{n}w_{ij}(f_i-f_j)^2$
证明：
$f^TLf = f^TDf - f^TWf = \sum\limits_{i=1}^{n}d_if_i^2 - \sum\limits_{i,j=1}^{n}w_{ij}f_if_j$
$=\frac{1}{2}( \sum\limits_{i=1}^{n}d_if_i^2 - 2 \sum\limits_{i,j=1}^{n}w_{ij}f_if_j + \sum\limits_{j=1}^{n}d_jf_j^2) = \frac{1}{2}\sum\limits_{i,j=1}^{n}w_{ij}(f_i-f_j)^2$
4) 拉普拉斯矩阵是半正定的，且对应的n个实数特征值都大于等于0，即0=λ1≤λ2≤…≤λn，且最小的特征值为0：
令 $\overline{1}$ 表示n*1的全1向量，则
四种常用聚类及代码（二）：谱聚类（spectral clustering）
由D和W的定义可以得出上式。

3、切图聚类

3.1、切图

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

对于任意两个子图点的集合A,B⊂V, A∩B=∅, 我们定义A和B之间的切图权重为：
$W(A, B) = \sum\limits_{i \in A, j \in B}w_{ij}$
那么对于我们k个子图点的集合：A₁,A₂,…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的补集，意为除Ai子集外其他V的子集的并集。

那么如何切图可以让子图内的点权重和高，子图间的点权重和低呢？一个自然的想法就是最小化cut(A₁,A₂,…A_k), 但是可以发现，这种极小化的切图存在问题，如下图：
四种常用聚类及代码（二）：谱聚类（spectral clustering）
我们选择一个权重最小的边缘的点，比如C和H之间进行cut，这样可以最小化cut(A₁,A₂,…A_k), 但是却不是最优的切图，如何避免这种切图，并且找到类似图中"Best Cut"这样的最优切图呢？

3.2、三种切图方法

3.2.1、最小切（mincut）

即最小化cut(A₁,A₂,…A_k)

3.2.2、RatioCut切图

RatioCut切图为了避免最小切图，对每个切图，不光考虑最小化cut(A₁,A₂,…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∈{h₁,h₂,…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}^T$ Lh_i,有：
$\begin{aligned} 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{aligned}$ hiTLhi=21m=1∑n=1∑wmn(him−hin)2=21(m∈Ai,n∈/Ai∑wmn(∣Ai∣1−0)2+m∈/Ai,n∈Ai∑wmn(0−∣Ai∣1)2=21(m∈Ai,n∈/Ai∑wmn∣Ai∣1+m∈/Ai,n∈Ai∑wmn∣Ai∣1=21(cut(Ai,Ai)∣Ai∣1+cut(Ai,Ai)∣Ai∣1)=∣Ai∣cut(Ai,Ai)
上述用了上面的拉普拉斯矩阵的性质以及指示向量的定义。可以看出，对于某一个子图i，它的RatioCut对应于 $h_{i}^T$ Lh_i,那么我们的k个子图呢？对应的RatioCut函数表达式为：
$RatioCut(A_1,A_2,...A_k) = \sum\limits_{i=1}^{k}h_i^TLh_i = \sum\limits_{i=1}^{k}(H^TLH)_{ii} = tr(H^TLH)$
其中tr(H^TLH)为矩阵的迹。也就是说，我们的RatioCut切图，实际上就是最小化我们的tr(H^TLH)。注意到H^TH=I,则我们的切图优化目标为:
$\underbrace{arg\;min}_H\; tr(H^TLH) \;\; s.t.\;H^TH=I$
注意到我们H矩阵里面的每一个指示向量都是n维的，向量中每个变量的取值为0或者 $\frac{1}{\sqrt{|A~j~|}}$ ，就有2ⁿ种取值，有k个子图的话就有k个指示向量，共有k2ⁿ种H，因此找到满足上面优化目标的H是一个NP难的问题。那么是不是就没有办法了呢？

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

对于 $h_{i}^T$ Lh_i，我们的目标是找到最小的L的特征值，而对于tr(H^TLH)= $\sum_{i=1}^k$ $h_{i}^T$ Lh_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现在不能完全指示各样本的归属，因此一般在得到nxk维度的矩阵H后还需要对每一行进行一次传统的聚类，比如使用K-Means聚类.

3.2.3、Ncut切图

Ncut切图和RatioCut切图很类似，但是把Ratiocut的分母|A_i|换成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}^T$ Lh_i,有：
$\begin{aligned} 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{aligned}$
推导方式和RatioCut完全一致。也就是说，我们的优化目标仍然是:
$NCut(A_1,A_2,...A_k) = \sum\limits_{i=1}^{k}h_i^TLh_i = \sum\limits_{i=1}^{k}(H^TLH)_{ii} = tr(H^TLH)$
但是此时我们的H^TH≠I,而是H^TDH=I。推导如下：
$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$
也就是说，此时我们的优化目标最终为：
$\underbrace{arg\;min}_H\; tr(H^TLH) \;\; s.t.\;H^TDH=I$
此时我们的H中的指示向量h并不是标准正交基，所以在RatioCut里面的降维思想不能直接用。怎么办呢？其实只需要将指示向量矩阵H做一个小小的转化即可。
我们令H=D^−1/2F, 则：H^TLH=F^TD^−1/2LD^−1/2F，H^TDH=F^TF=I,也就是说优化目标变成了:
$\underbrace{arg\;min}_F\; tr(F^TD^{-1/2}LD^{-1/2}F) \;\; s.t.\;F^TF=I$
可以发现这个式子和RatioCut基本一致，只是中间的L变成了D^−1/2LD^−1/2。这样我们就可以继续按照RatioCut的思想，求出D^−1/2LD^−1/2的最小的前k个特征值，然后求出对应的特征向量，并标准化，得到最后的特征矩阵F,最后对F进行一次传统的聚类（比如K-Means）即可。
一般来说， D^−1/2LD^−1/2相当于对拉普拉斯矩阵L做了一次标准化，即 $\frac{L_{ij}}{\sqrt{d_i*d_j}}$

4、谱聚类流程

第一步：数据准备，生成图的邻接矩阵；

第二步：归一化普拉斯矩阵；

第三步：生成最小的k个特征值和对应的特征向量；

第四步：将特征向量kmeans聚类(少量的特征向量)；

一般来说，谱聚类主要的注意点为相似矩阵的生成方式，切图的方式以及最后的聚类方法。
具体介绍一下：
最常用的相似矩阵的生成方式是基于高斯核距离的全连接方式，最常用的切图方式是Ncut。而到最后常用的聚类方法为K-Means。下面以Ncut总结谱聚类算法流程。

输入：样本集D=(x₁,x₂,…,x_n)，相似矩阵的生成方式, 降维后的维度k₁, 聚类方法，聚类后的维度k₂
输出：簇划分C(c₁,c₂,… $c_{k_{2}}$ )

1）根据输入的相似矩阵的生成方式构建样本的相似矩阵S

2）根据相似矩阵S构建邻接矩阵W，构建度矩阵D

3）计算出拉普拉斯矩阵L

4）构建标准化后的拉普拉斯矩阵D^−1/2LD^−1/2
　　　　
5）计算D^−1/2LD^−1/2最小的k₁个特征值所各自对应的特征向量f

6）将各自对应的特征向量f组成的矩阵按行标准化，最终组成n×k₁维的特征矩阵F

7）对F中的每一行作为一个k₁维的样本，共n个样本，用输入的聚类方法进行聚类，聚类维数为k₂。

8）得到簇划分C(c₁,c₂,… $c_{k_{2}}$ ).

5、谱聚类优缺点

优点

（1）当聚类的类别个数较小的时候，谱聚类的效果会很好，但是当聚类的类别个数较大的时候，则不建议使用谱聚类；

（2）谱聚类算法使用了降维的技术，所以更加适用于高维数据的聚类；

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

（4）谱聚类算法建立在谱图理论基础上，与传统的聚类算法相比，它具有能在任意形状的样本空间上聚类且收敛于全局最优解

缺点

（1）谱聚类对相似度图的改变和聚类参数的选择非常的敏感；

（2）谱聚类适用于均衡分类问题，即各簇之间点的个数相差不大，对于簇之间点个数相差悬殊的聚类问题，谱聚类则不适用；

6、python实现

import numpy
import scipy
from sklearn.cluster import KMeans


def laplacian(A):
    """Computes the symetric normalized laplacian.
    L = D^{-1/2} A D{-1/2}
    """
    D = numpy.zeros(A.shape)
    w = numpy.sum(A, axis=0)
    D.flat[::len(w) + 1] = w ** (-0.5)  # set the diag of D to w
    return D.dot(A).dot(D)


def k_means(X, n_clusters):
    kmeans = KMeans(n_clusters=n_clusters, random_state=1231)
    return kmeans.fit(X).labels_


def spectral_clustering(affinity, n_clusters, cluster_method=k_means):
    L = laplacian(affinity)
    eig_val, eig_vect = scipy.sparse.linalg.eigs(L, n_clusters)
    X = eig_vect.real
    rows_norm = numpy.linalg.norm(X, axis=1, ord=2)
    Y = (X.T / rows_norm).T
    labels = cluster_method(Y, n_clusters)
    return labels
    ```