谱聚类 - 魏允臣

背景：谱聚类（Spectral Clustering）于 2006 年由 Ulrike von Luxburg 公布在论文 A Tutorial on Spectral Clustering 上。
本质：一种观点认为，谱聚类的本质就是通过拉普拉斯矩阵变换得到其特征向量组成的新矩阵的 K-Means 聚类，而其中的拉普拉斯矩阵变换可以被简单地看作是降维的过程。而谱聚类中的「谱」其实就是矩阵中全部特征向量的总称。

拉普拉斯矩阵

拉普拉斯矩阵（Laplacian Matrix），也称为基尔霍夫矩阵，是无向图的一种矩阵表示形式。对于一个有 $n$ 个顶点的图，其拉普拉斯矩阵定义为： $L_{n \times n} = D - A \tag{1}$ 其中， $D$ 为图的 度矩阵 ， $A$ 为图的 邻接矩阵 。

对于有边连接的两个点 $v_i$ 和 $v_j$ ， $w_{ij}>0$ ，对于没有边连接的两个点 $v_i$ 和 $v_j$ ， $w_{ij}=0$ 。对于图中的任意一个点 $v_i$ ，它的度 $d_i$ 定义为和它相连的所有边的权重之和，即 $d_i = \sum\limits_{j=1}^{n}w_{ij} \tag{2}$ 度矩阵是一个对角矩阵，主角线上的值由对应的顶点的度组成。

对于一张有 $n$ 个顶点的图 $L_{n \times n}$ ，其邻接矩阵 $A$ 为任意两点之间的权重值 $w_{ij}$ 组成的矩阵： $A=\begin{pmatrix} w_{11} & \cdots & w_{1n}\\ \vdots & & \vdots \\ w_{n1} & \cdots & w_{nn} \end{pmatrix} \tag{4}$ 对于构建邻接矩阵 $A$ ，一般有三种方法，分别是： $\epsilon$ -邻近法，全连接法和 K-近邻法。
第一种方法， $\epsilon$ -邻近法通过设置了一个阈值 $\epsilon$ ，再求解任意两点 $x_{i}$ 和 $x_{j}$ 间的欧式距离 $s_{ij}$ 来度量相似性。然后，根据 $s_{ij}$ 和 $\epsilon$ 的大小关系，定义 $w_{ij}$ 如下： $w_{ij}= \begin{cases} 0& {s_{ij} > \epsilon}\\ \epsilon& {{s_{ij} \leq \epsilon}} \end{cases} \tag{5}$ 第二种方法，全连接法通过选择不同的核函数来定义边权重，常用的有多项式核函数，高斯核函数和 Sigmoid 核函数。例如，当使用高斯核函数 RBF 时， $w_{ij}$ 定义如下： $w_{ij}=s_{ij}=exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2}) \tag{6}$ 除此之外，K-邻近法也可以被用于生成邻接矩阵。当我们令 $x_i$ 对应图中的一个节点 $v_i$ 时，如果 $v_j$ 属于 $v_i$ 的前 $k$ 个最近邻节点集合，则连接 $v_j$ 和 $v_i$ 。
但是，此种方法会产生不对称的有向图。因此，将有向图转换为无向图的方法有两种：
第一，忽略方向性，即如果 $v_j$ 属于 $v_i$ 的前 $k$ 个最近邻节点集，或 $v_i$ 属于 $v_j$ 的前 $k$ 个最近邻节点集，则连接 $v_j$ 和 $v_i$
$w_{ij}=w_{ji}= \begin{cases} 0& {x_i \notin KNN(x_j) \;and \;x_j \notin KNN(x_i)}\\ exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})& {x_i \in KNN(x_j)\; or\; x_j \in KNN(x_i}) \end{cases} \tag{7}$ 第二是当且仅当 $v_j$ 属于 $v_i$ 的前 $k$ 个最近邻节点集，且 $v_i$ 属于 $v_j$ 的前 $k$ 个最近邻节点集时连接 $v_j$ 和 $v_i$ 。第二种方法体现了相互性。 $w_{ij}=w_{ji}= \begin{cases} 0& {x_i \notin KNN(x_j) \;or\;x_j \notin KNN(x_i)}\\ exp(-\frac{||x_i-x_j||_2^2}{2\sigma^2})& {x_i \in KNN(x_j)\; and \; x_j \in KNN(x_i}) \end{cases} \tag{8}$ 那么，对于 1.1 中的无向图，它对应的邻接矩阵、度矩阵以及拉普拉斯矩阵如下。
那么，对于 $k$ 个子图点的集合： $A_1, A_2, \cdots, A_k$ 定义切图为：
$cut(A_1,A_2,\cdots, A_k) = \frac 1 2 \sum\limits_{i=1}^{k}W(A_i, \overline{A}_i ) \tag{11}$ 其中 $\overline{A}_i$ 为 $A_i$ 的补集，意为除 $A_i$ 子集外其他 $V$ 的子集的并集。

容易看出， $cut(A_1,A_2,\cdots, A_k)$ 描述了子图之间的相似性，其值越小则代表子图的差异性越大。但是，公式（11）在划分子图时并没有考虑到子图最少节点的数量，这就容易导致一个数据点或者很少的数据点被划为一个独立的子图，但这明显是不正确的。
为了接近这个问题，我们会引入一些正则化方法，其中最常用的就是 RatioCut 和 Ncut。
RatioCut 切图时不光考虑最小化 $cut(A_1,A_2,\cdots, A_k)$ ，它还同时考虑最大化每个子图点的个数，其定义如下：
$RatioCut(A_1,A_2,...A_k) = \sum\limits_{i=1}^{k}\frac{cut(A_i, \overline{A}_i )}{|A_i|} \tag{12}$ 其中， $|A_i|$ 表示子图 $A_i$ 中节点的个数。
Ncut 切图和 RatioCut 切图很类似，但是把 Ratiocut 的分母 $|A_i|$ 换成 ${assoc(A_i)}$ . 由于子图样本的个数多并不一定权重就大，我们切图时基于权重也更合我们的目标，因此一般来说 Ncut 切图优于 RatioCut 切图。
$NCut(A_1,A_2,...A_k) = \sum\limits_{i=1}^{k}\frac{cut(A_i, \overline{A}_i )}{assoc(A_i, V)} \tag{13}$ 其中， ${assoc(A_i, V)} = \sum_{V_j\in A_i}^{k} d_j$ ， $d_j = \sum_{i=1}^{n} w_{ji}$ 。

谱聚类流程及实现

根据数据构造无向图

G

，图中的每一个节点对应一个数据点，将相似的点连接起来，并且边的权重用于表示数据之间的相似度。

计算图的邻接矩阵

A

和度矩阵

D

，并求解对应的拉普拉斯矩阵

L

。

求出

L

的前

k

个由小到大排列的特征值

\{\lambda\}_{i=1}^k

以及对应的特征向量

\{v\}_{i=1}^k

。

把

k

个特征向量排列在一起组成一个

N\times k

的矩阵，将其中每一行看作

k

维空间中的一个向量，并使用 K-Means 算法进行聚类，并得到最终的聚类类别。

可以看到，谱聚类的最后还是会用到 K-Means，谱聚类的本质可以理解为通过拉普拉斯矩阵变换得到其特征向量组成的新矩阵的 K-Means 聚类，而其中的拉普拉斯矩阵变换可以被简单地看作是降维的过程。而谱聚类中的「谱」其实就是矩阵中全部特征向量的总称。

首先，我们按照谱聚类流程，先使用 Python 实现谱聚类。先生成一组示例数据。我们使用 make_circles() 生成环状数据。

import numpy as np
from sklearn import datasets
from matplotlib import pyplot as plt
%matplotlib inline
noisy_circles, _ = datasets.make_circles(n_samples=300, noise=0.1, factor=.3, random_state=10)
plt.scatter(noisy_circles[:,0], noisy_circles[:,1], color='b')
参考上面的谱聚类流程，我们知道首先需要计算由数据生成无向图的邻接矩阵  $A$ ，本次实验使用 K-近邻方法计算无向图中边对应的相似度权重。下面通过 knn_similarity_matrix(data, k) 函数实现：
def knn_similarity_matrix(data, k):
    zero_matrix = np.zeros((len(data), len(data)))
    w = np.zeros((len(data), len(data)))
    for i in range(len(data)):
        for j in range(i + 1, len(data)):
            zero_matrix[i][j] = zero_matrix[j][i] = np.linalg.norm(data[i] - data[j]) # 计算欧式距离
    for i, vector in enumerate(zero_matrix):
        vector_i  = np.argsort(vector)
        w[i][vector_i[1 : k + 1]] = 1
    w = (np.transpose(w) + w)/2
    return w
得到邻接矩阵  $A$ ，就可以计算得到度矩阵  $D$ ，以及对应的拉普拉斯矩阵  $L$ ，进而实现整个谱聚类方法：
def spectral_clustering(data, k, n):  
    # 计算近邻矩阵、度矩阵、拉普拉斯矩阵
    A_matrix = knn_similarity_matrix(data, k) 
    D_matrix = np.diag(np.power(np.sum(A_matrix, axis=1), -0.5))  
    L_matrix = np.eye(len(data)) - np.dot(np.dot(D_matrix, A_matrix), D_matrix)  
    # 计算特征值和特征向量
    eigvals, eigvecs = np.linalg.eig(L_matrix)
    # 选择前 n 个最小的特征向量
    indices = np.argsort(eigvals)[: n]  
    k_eigenvectors = eigvecs[:, indices]
    k_eigenvectors
    # 使用 K-Means 完成聚类
    clusters = KMeans(n_clusters=n).fit_predict(k_eigenvectors)
    return clusters
最后，我们定义参数 k=5, 并将数据集聚为 2 类。
sc_clusters = spectral_clustering(noisy_circles, k=5, n=2)
sc_clusters
Out[6]:
array([0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1,
              0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0,
              1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0,
              0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,
              1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1,
              0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1,
              0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1,
              1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1,
              1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1,
              1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1,
              1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1,
              1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1,
              1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0,
              0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1], dtype=int32)
根据聚类标签可视化数据集，可以看到谱聚类的效果是非常不错的。
plt.scatter(noisy_circles[:,0], noisy_circles[:,1], c=sc_clusters)
scikit-learn 中的谱聚类
scikit-learn 提供了谱聚类的实现方法，具体为 sklearn.cluster.SpectralClustering()：
sklearn.cluster.SpectralClustering(n_clusters=8, eigen_solver=None, random_state=None, n_init=10, gamma=1.0, affinity='rbf', n_neighbors=10, eigen_tol=0.0, assign_labels='kmeans', degree=3, coef0=1, kernel_params=None, n_jobs=1)
其主要参数有：
n_clusters：聚类簇数量。
eigen_solver：特征值求解器。
gamma：affinity 使用核函数时的核函数系数。
affinity：邻接矩阵计算方法，可选择核函数、k-近邻以及  $\epsilon$ -邻近法。
n_neighbors：邻接矩阵选择 k-近邻法时的 k 值。
assign_labels：最终聚类方法，默认为 K-Means。

谱聚类在很多时候会被用于图像分割，我们尝试使用谱聚类的 scikit-learn 实现来完成一个有趣的例子。该例子参考自 scikit-learn 官方文档。
首先，我们需要生成一幅示例图像，图像中有几个大小不等的圆，且存在噪点。
1. 生成 100px * 100px 的图像
2. 在图像中添加 4 个圆
3. 添加随机噪声点
l = 100
x, y = np.indices((l, l))
center1 = (28, 24)
center2 = (40, 50)
center3 = (67, 58)
center4 = (24, 70)
radius1, radius2, radius3, radius4 = 16, 14, 15, 14
circle1 = (x - center1[0]) ** 2 + (y - center1[1]) ** 2 < radius1 ** 2
circle2 = (x - center2[0]) ** 2 + (y - center2[1]) ** 2 < radius2 ** 2
circle3 = (x - center3[0]) ** 2 + (y - center3[1]) ** 2 < radius3 ** 2
circle4 = (x - center4[0]) ** 2 + (y - center4[1]) ** 2 < radius4 ** 2
img = circle1 + circle2 + circle3 + circle4
mask = img.astype(bool)
img = img.astype(float)
img += 1 + 0.2 * np.random.randn(*img.shape)
plt.figure(figsize=(5, 5))
plt.imshow(img)
from sklearn.feature_extraction import image
from sklearn.cluster import spectral_clustering
graph = image.img_to_graph(img, mask=mask) # 图像处理为梯度矩阵
graph.data = np.exp(-graph.data / graph.data.std()) # 正则化
labels = spectral_clustering(graph, n_clusters=4)
label_im = -np.ones(mask.shape)
label_im[mask] = labels
plt.figure(figsize=(5, 5))
plt.imshow(label_im)
可以看出，谱聚类的表现普遍优于 K-Means。

谱聚类也有一些缺点。例如，由于最后使用的 K-Means 聚类，所以要提前指定聚类数量。另外，当使用 K-近邻生成邻接矩阵时还需要指定最近邻样本数量，对参数是比较敏感的。
聚类方法选择