试用谱聚类算法进行聚类分 析

谱聚类

1. 基本原理

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

难点：

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

2. 谱聚类基础

2.1 无向权重图

对于一个图$G$，我们一般用点集合$V=\{v_1,v_2,....,v_n\}$和边集合$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}^nw_{ij}$

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

除此之外，对于点集$V$的一个子集$A \subset V$，我们定义：
$|A|:=子集A中点的个数 \\ vol(A):=\sum_{i \in A}d_i$

2.2 拉普拉斯矩阵

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

1. 对称矩阵，由于$D和W$都为对称矩阵
2. 由于是对称矩阵，它的所有特征值都是实数
3. 对于任意向量$f$，有
   $f^TLf=f^TDf-f^TWf=\sum_{i=1}^nd_if_i^2-\sum_{i,j=1}^nw_{ij}f_if_j\\=\frac {1} {2}(\sum_{i=1}^nd_if_i^2-2\sum_{i,j=1}^nw_{ij}f_if_j+\sum_{j=1}^nd_jf_j^2)=\frac {1} {2}\sum_{i,j=1}^nw_{ij}(f_i-f_j)^2$
4. 由于拉普拉斯矩阵是半正定的，其对应的$n$个特征值都大于等于0。

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

3.1 $\epsilon$邻近法

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

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

3.2 $k$近邻法

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

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

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}^kW(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|}$

我们引入指示向量$h_j\in \{h_1,h_2,...,h_k\}$，对于任意一个向量$h_j=(h_{1,j},h_{2,j},...,h_{n,j})^T$，它是一个$n$维向量（$n$为样本数），我们定义$h_{ij}$为：
$h_{ij}= \begin{cases} 0 & v_i \notin A_j \\ \frac {1} {|A_j|} & v_i \in A_j \end{cases}$
对于$h_i^TLh_i$有：
$h_i^TLh_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|}$

由上式可知，$RatioCut$函数表达式可改写为：
$RatioCut(A_1,A_2,...,A_k)=\sum_{i=1}^kh_i^TLh_i=\sum_{i=1}^k(H^TLH)_{ii}=tr(H^TLH)$
其中$tr(H^TLH)$为矩阵的迹，即我们的$RatioCut$切图，实际上是最小化迹$tr(H^TLH)$。注意到$H^TH=I$，则我们的优化目标为：
$arg\min_{H}tr(H^TLH)\\s.t. \quad H^TH=I$

注意观察$tr(H^TLH)$的每一个优化子目标$h_i^TLh_i$，其中$h_i$是单位正交基，$L$是对称矩阵，此时$h_i^TLh_i$是矩阵$L$的一个特征值。对于$h_i^TLh_i$，我们的目标是找到矩阵$L$的最小特征值，而对于$tr(H^TLH)=\sum_{i=1}^kh_i^TLh_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{cases} 0 & v_i \notin A_j \\ \frac {1} {\sqrt {vol(A_j)}} & v_i \in A_j \end{cases}$

我们的优化目标依然是：（推导与$RatioCut$完全一致）
$Ncut(A_1,A_2,...,A_k)=\sum_{i=1}^kh_i^TLh_i=\sum_{i=1}^k(H^TLH)_{ii}=tr(H^TLH)$
但是此时我们的$H^TH \not=I$，而是$H^TDH=I$。推导如下：
$h_i^TDh_i=\sum_{j=1}^nh_{ji}^2d_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$
也就是说，我们的优化目标最终为：
$arg\min_{H}tr(H^TLH) \\ s.t. \quad H^TDH=I$
此时我们的$H$中的指示向量$h$不是单位正交基，所以我们令$H=D^{-\frac {1} {2}}F$，则$H^TLH=F^TD^{-\frac {1}{2}}LD^{-\frac {1}{2}}F,H^TDH=F^TF=I$，也就是优化目标变成了：
$arg\min_{F}tr(F^TD^{-\frac {1}{2}}LD^{-\frac {1}{2}}F) \\ s.t. \quad F^TF=I$
可以发现这个式子和$RatioCut$基本一致，只是中间的$L$变成了$D^{-\frac {1}{2}}LD^{-\frac {1}{2}}$。这样，我们可以按照$RatioCut$的思想，求出$D^{-\frac {1}{2}}LD^{-\frac {1}{2}}$的$k$个最小特征值

一般来说，$D^{-\frac {1}{2}}LD^{-\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}}LD^{-\frac {1}{2}}$
4. 计算$D^{-\frac {1}{2}}LD^{-\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. 聚类效果依赖于邻接矩阵，不同的邻接矩阵得到的最终聚类效果可能很不同