谱聚类

背景：谱聚类（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}$ 度矩阵是一个对角矩阵，主角线上的值由对应的顶点的度组成。
$D = \left( \begin{array}{ccc} d_1 & \ldots & \ldots \\ \ldots & d_2 & \ldots \\ \vdots & \vdots & \ddots \\ \ldots & \ldots & d_n \end{array} \right) \tag{3}$

邻接矩阵

对于一张有 $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 中的无向图，它对应的邻接矩阵、度矩阵以及拉普拉斯矩阵如下。

邻接矩阵为：

A= \left(\begin{array}{rrrrrr} 0 & 1 & 0 & 0 & 1 & 0\\ 1 & 0 & 1 & 0 & 1 & 0\\ 0 & 1 & 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 & 1 & 1\\ 1 & 1 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0\\ \end{array}\right)

度矩阵为：

D= \left(\begin{array}{rrrrrr} 2 & 0 & 0 & 0 & 0 & 0\\ 0 & 3 & 0 & 0 & 0 & 0\\ 0 & 0 & 2 & 0 & 0 & 0\\ 0 & 0 & 0 & 3 & 0 & 0\\ 0 & 0 & 0 & 0 & 3 & 0\\ 0 & 0 & 0 & 0 & 0 & 1\\ \end{array}\right)

拉普拉斯矩阵为：

L = D - A = \left(\begin{array}{rrrrrr} 2 & -1 & 0 & 0 & -1 & 0\\ -1 & 3 & -1 & 0 & -1 & 0\\ 0 & -1 & 2 & -1 & 0 & 0\\ 0 & 0 & -1 & 3 & -1 & -1\\ -1 & -1 & 0 & -1 & 3 & 0\\ 0 & 0 & 0 & -1 & 0 & 1\\ \end{array}\right) \tag{9}

无向图切图

聚类其实就是寻求一种把整个无向图切分成小块的过程，而切分出来的小块就是每一个类别（簇）。所以，谱聚类的目标是将无向图划分为两个或多个子图，使得子图内部节点相似而子图间节点相异，从而达到聚类的目的。

对于无向图 G 的切图，我们的目标是将图

G(V,E)

切成相互没有连接的

k

个子图，每个子图点的集合为：

A_1, A_2, \cdots, A_k

，它们满足

A_i \cap A_j = \varnothing

，且

A_1 \cup A_2 \cup \cdots \cup A_k = V

对于任意两个子图点的集合 $A, B \subset V$ , $A \cap B=\varnothing$ , 我们定义 $A$ 和 $B$ 之间的切图权重为：

$W(A, B) = \sum\limits_{i \in A, j \in B}w_{ij} \tag{10}$
那么，对于 $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')

from sklearn.cluster import KMeans

plt.scatter(noisy_circles[:,0], noisy_circles[:,1], c=KMeans(n_clusters=2).fit_predict(noisy_circles))

参考上面的谱聚类流程，我们知道首先需要计算由数据生成无向图的邻接矩阵 $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)

接下来，我们使用谱聚类方法完成图像边缘检测，并得到处理后的图像。

"""
1. 生成 100px * 100px 的图像
2. 在图像中添加 4 个圆
3. 添加随机噪声点
"""

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 会将少数离群点划为一类的现象。

D. Cai 等（DOI: 10.1109/TKDE.2005.198）曾经在 TDT2 和 Reuters-21578 基准聚类数据集上测试过谱聚类和 K-Means 的聚类准确率表现，结果如下：

	TDT2	TDT2	Reuters-21578	Reuters-21578
k 值	K-Means	谱聚类	K-Means	谱聚类
2	0.989	0.998	0.871	0.923
3	0.974	0.996	0.775	0.816
4	0.959	0.996	0.732	0.793
…	…	…	…	…
9	0.852	0.984	0.553	0.625
10	0.835	0.979	0.545	0.615

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

聚类方法选择

最后编辑于：2019.02.03 03:34:20

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 204,684评论 6赞 478
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 87,143评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 151,214评论 0赞 337
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,788评论 1赞 277
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,796评论 5赞 368
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,665评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 38,027评论 3赞 399
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,679评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 41,346评论 1赞 299
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,664评论 2赞 321
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,766评论 1赞 331
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,412评论 4赞 321
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 39,015评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,974评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,203评论 1赞 260
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 45,073评论 2赞 350
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,501评论 2赞 343

谱聚类