谱聚类

1. 基本原理

它的主要思想：把所有数据看成空间中的点，这些点之间可以用变连接起来，距离较远的两个点之间的边权重较低，而距离较近的两个点之间的权重较高，通过对所有数据点组成的图进行切图，让切图后的不同的子图间边权重和尽可能小（即距离远），而子图内的边权重和尽可能高（即距离近）。

难点：

1. 如何构建图？
2. 如何切分图？

2. 谱聚类基础

2.1 无向权重图

对于一个图$G$，我们一般用点集合V={v1,v2,....,vn}V=\{v_1,v_2,....,v_n\}V={v1​,v2​,....,vn​}和边集合$E$来描述，即$G=(V,E)$。我们定义权重$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 _{j=1}^n{w}_{ij}$$

利用每个点度的定义，我们可以得到一个$n\times n$的度矩阵$D$，它是一个对角阵，只有主对角有值，对应第$i$行为第$i$个点的度；利用所有点之间的权重，我们可以得到图的邻接矩阵$W$，它也是一个$n\times n$矩阵，第$i$行的第$j$个值对应权重$w_{ij}$

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

2.2 拉普拉斯矩阵

拉普拉斯矩阵$L=D-W$，其性质如下：

1. 对称矩阵，由于$D和W$都为对称矩阵
2. 由于是对称矩阵，它的所有特征值都是实数
3. 对于任意向量$f$，有
   $$\begin{gathered} 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{gathered}$$
4. 由于拉普拉斯矩阵是半正定的，其对应的$n$个特征值都大于等于0。

3. 构建图——构建邻接矩阵

3.1 $ϵ$邻近法

通过设置一个阈值$ϵ$，然后利用欧氏距离$s_{ij}$度量任意两点$v_i和v_j$的距离，即$s_{ij}=||v_i-v_j|{|}_2^2$，然后根据$s_{ij}和ϵ$的大小关系，来定义邻接矩阵$W$：
$$w_{ij}=\{\begin{array}{ll}0,& s_{ij}>ϵ\\ ϵ,& s_{ij}\le ϵ\end{array}$$

从上式可知，两点间的权重要么$ϵ$，要么0，就没有其他信息了，距离远近度量很不明确，因此在实际应用中，很少采用。

3.2 $k$近邻法

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

1. 只要一个点在另一个点的K近邻中，就保留$w_{ij}$
   $$w_{ij}=w_{ji}=\{\begin{array}{ll}0& v_i\notin KNN(v_j)and{v}_j\notin KNN(v_i)\\ exp(-\frac{||v_i-x_j|{|}_2^2}{2{\sigma }^2})& v_i\in KNN(v_j)or{v}_j\in KNN(v_i)\end{array}$$
2. 必须两个点互为$K$近邻中，才能保留$w_{ij}$
   $$w_{ij}=w_{ji}=\{\begin{array}{ll}0& v_i\notin KNN(v_j)or{v}_j\notin KNN(v_i)\\ exp(-\frac{||v_i-x_j|{|}_2^2}{2{\sigma }^2})& v_i\in KNN(v_j)and{v}_j\in KNN(v_i)\end{array}$$

3.3 全连接法

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

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

4. 图的切分

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

对于任意两个子图点的集合$A,B\subset V,A\cap B=\varnothing$，我们定义$A和B$之间的切图权重为：
$$W(A,B)=\sum _{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 _{i=1}^{k}W(A_i,{\bar{A}}_i)$$
其中${\bar{A}}_i$为$A_i$的补集

那么如何切图可以让子图内的点权重和高，子图间的点权重和低呢？

一个自然的想法就是最小化$cut(A_1,A_2,...A_k)$, 但是可以发现，这种极小化的切图存在问题，如下图：

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

4.1 $RatioCut$切图

对于每个切图，不仅要考虑最小化$cut(A_1,A_2,...,A_k)$，还要考虑最大化每个子图样本的个数，即最小化$RatioCut$函数：
$$RatioCut(A_1,A_2,...,A_k)=\frac{1}{2}\sum _{i=1}^k\frac{W(A_i,{\bar{A}}_i)}{|A_i|}$$

我们引入指示向量hj∈{h1,h2,...,hk}h_j\in \{h_1,h_2,...,h_k\}hj​∈{h1​,h2​,...,hk​}，对于任意一个向量$h_j=(h_{1,j},h_{2,j},...,h_{n,j}{)}^T$，它是一个$n$维向量（$n$为样本数），我们定义$h_{ij}$为：
$$h_{ij}=\{\begin{array}{ll}0& v_i\notin A_j\\ \frac{1}{|A_j|}& v_i\in A_j\end{array}$$
对于$h_i^{T}L{h}_i$有：
$$\begin{gathered} h_i^{T}L{h}_i=\frac{1}{2}\sum _{m=1}\sum _{n=1}{w}_{mn}(h_{m,i}-h_{n,i}{)}^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,{\bar{A}}_i)\frac{1}{|A_i|}+cut(A_i,{\bar{A}}_i)\frac{1}{|A_i|})=\frac{cut(A_i,{\bar{A}}_i)}{|A_i|} \end{gathered}$$

由上式可知，$RatioCut$函数表达式可改写为：
$$RatioCut(A_1,A_2,...,A_k)=\sum _{i=1}^k{h}_i^{T}L{h}_i=\sum _{i=1}^k(H^{T}LH{)}_{ii}=tr(H^{T}LH)$$
其中$tr(H^{T}LH)$为矩阵的迹，即我们的$RatioCut$切图，实际上是最小化迹$tr(H^{T}LH)$。注意到$H^{T}H=I$，则我们的优化目标为：
$$\begin{gathered} arg\underset{H}{\min }tr(H^{T}LH) \\ s.t.H^{T}H=I \end{gathered}$$

注意观察$tr(H^{T}LH)$的每一个优化子目标$h_i^{T}L{h}_i$，其中$h_i$是单位正交基，$L$是对称矩阵，此时$h_i^{T}L{h}_i$是矩阵$L$的一个特征值。对于$h_i^{T}L{h}_i$，我们的目标是找到矩阵$L$的最小特征值，而对于$tr(H^{T}LH)={\sum }_{i=1}^k{h}_i^{T}L{h}_i$，我们的目标就是找到矩阵$L$的$k$个最小特征值。

4.2 $Ncut$切图

$Ncut$切图与$Ratio$切图类似，只是将$RatioCut$的分母$|A_i|$换成$vol(A_i)$。由于子图样本的个数多不一定权重就大，我们切图时基于权重也更符合我们的目标，因此一般来说$Ncut$优于$RatioCut$，定义如下：
$$Ncut(A_1,A_2,...,A_k)=\frac{1}{2}\sum _{i=1}^k\frac{W(A_i,{\bar{A}}_i)}{vol(A_i)}$$

对应的，$Ncut$切图对指示向量$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}$$

我们的优化目标依然是：（推导与$RatioCut$完全一致）
$$Ncut(A_1,A_2,...,A_k)=\sum _{i=1}^k{h}_i^{T}L{h}_i=\sum _{i=1}^k(H^{T}LH{)}_{ii}=tr(H^{T}LH)$$
但是此时我们的$H^{T}H\ne I$，而是$H^{T}DH=I$。推导如下：
$$h_i^{T}D{h}_i=\sum _{j=1}^n{h}_{ji}^2{d}_j=\sum _{j\in A_i}\frac{1}{vol(A_i)}{d}_j=\frac{1}{vol(A_i)}\sum _{j\in A_i}{d}_j=\frac{1}{vol(A_i)}vol(A_i)=1$$
也就是说，我们的优化目标最终为：
$$\begin{gathered} arg\underset{H}{\min }tr(H^{T}LH) \\ s.t.H^{T}DH=I \end{gathered}$$
此时我们的$H$中的指示向量$h$不是单位正交基，所以我们令$H=D^{-\frac{1}{2}}F$，则$H^{T}LH=F^T{D}^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}F,H^{T}DH=F^{T}F=I$，也就是优化目标变成了：
$$\begin{gathered} arg\underset{F}{\min }tr(F^T{D}^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}F) \\ s.t.F^{T}F=I \end{gathered}$$
可以发现这个式子和$RatioCut$基本一致，只是中间的$L$变成了$D^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}$。这样，我们可以按照$RatioCut$的思想，求出$D^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}$的$k$个最小特征值

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

5. 谱聚类算法流程

$输入：样本集D=(x_1,x_2,...,x_n)，邻接矩阵的生成方式,降维后的维度k_1,聚类方法，聚类后的维度k_2$

$输出：簇划分C(c_1,c_2,...c_{k_2})$

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

6. 实例演示

import numpy as np 
import matplotlib.pyplot as plt 

from sklearn import cluster, datasets
from sklearn.preprocessing import StandardScaler

np.random.seed(0)

# 构建数据
n_samples = 1500
noisy_circles = datasets.make_circles(n_samples=n_samples, factor=0.5, noise=0.05)
noisy_moons = datasets.make_moons(n_samples=n_samples, noise=0.05)
blobs = datasets.make_blobs(n_samples=n_samples, random_state=8)

data_sets = [
    (
        noisy_circles,
        {
            "n_clusters": 2
        }
    ),
    (
        noisy_moons,
        {
            "n_clusters": 2
        }
    ), 
    (
        blobs, 
        {
            "n_clusters": 3
        }
    )
]
colors = ["#377eb8", "#ff7f00", "#4daf4a"]
affinity_list = ['rbf', 'nearest_neighbors']

plt.figure(figsize=(17, 10))

for i_dataset, (dataset, algo_params) in enumerate(data_sets):
    # 模型参数
    params = algo_params

    # 数据
    X, y = dataset
    X = StandardScaler().fit_transform(X)

    for i_affinity, affinity_strategy in enumerate(affinity_list):
        # 创建SpectralCluster
        spectral = cluster.SpectralClustering(
            n_clusters=params['n_clusters'],
            eigen_solver='arpack', 
            affinity=affinity_strategy
        )

        # 训练
        spectral.fit(X)

        # 预测
        y_pred = spectral.labels_.astype(int)

        y_pred_colors = []

        for i in y_pred:
            y_pred_colors.append(colors[i])
        
        plt.subplot(3, 4, 4*i_dataset+i_affinity+1)
        plt.title(affinity_strategy)
        plt.scatter(X[:, 0], X[:, 1], color=y_pred_colors)

plt.show()

7. 谱聚类算法小结

优点：

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

缺点：

1. 如果最终聚类的维度非常高，则由于降维的幅度不够，谱聚类的运行速度和最后的聚类效果均不好
2. 聚类效果依赖于邻接矩阵，不同的邻接矩阵得到的最终聚类效果可能很不同