聚类算法 距离矩阵_谱聚类

作者：汤进​

本文详细阐述了聚类中的一种算法--谱聚类，并通过代码实现展示了谱聚类的所有细节，将实践和理论紧密的结合起来，值得细细把玩。

在聚类与K-Means一文中，ARGO就非监督学习的聚类算法做了比较系统的介绍，在其中提到了用图聚类的算法-谱聚类，这次ARGO将详细聊聊这个算法。本文将分以下几个部分进行介绍：

• 基础概念
• 怎么度量样本间的相似程度
• 谱聚类的优化目标
• 算法步骤总结
• 实例

一 基础概念

聚类的概念是指，将多维空间中的样本分为多个组，使得组内的样本尽可能相似，组间样本尽可能不同。有K-means，谱聚类，DBSCAN等多种聚类算法。

谱聚类是图论演变而来的聚类算法。将所有样本点

视为图的顶点，并使用邻接矩阵W表示图，矩阵中的每一个元素

表示

边的权重。当

到

的边不存在的时候，

的值取一个非常小，趋近或等于0的值；当

到

的邻接非常紧密，

则取一个相对比较大的值。

对于n维空间内，一个有m个顶点的图，它的顶点组成矩阵：

这样的图，它的邻接矩阵是一个m×m矩阵：

之后，谱聚类算法将这样的图切分成k个子图，使得子图间的权重尽可能小，子图内的权重尽可能大，而这k个子图就是聚类生成的k个类。可是如何才能建立上述的邻接矩阵图呢？这得从描述样本间的相似程度说起。

二 如何度量样本间的相似程度

K-means算法使用欧式距离度量样本间的相似度，距离越小相似程度越大。与此类似， 谱聚类使用边的权重度量样本的相似程度，权重越大相似度越大。所以谱聚类的关键是计算出样本间的边的权重所使用的方法。谱聚类的相似度量方式有很多，常见的方法有两种，一是高斯核函数，一是K近邻方式。

高斯核函数

通常使用高斯核函数RBF将欧式距离转换为权重，以度量相似程度。 举一个简单的例子，假设空间中有10个顶点如下图所示：

import matplotlib.pyplot as plt
import numpy as np

x = np.arange(0, 10, 1)
y = np.zeros(x.shape)

fig, ax = plt.subplots()
fig.set_size_inches(5, 5)

ax.scatter(x, y)
for i in range(10):
    ax.annotate('$x_{%d}$' % (i+1), xy=(x[i], y[i]), xytext=(x[i]+0.05, y[i]+0.001))

# Hide the right and top spines
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

plt.show()

对于

来说，要使距离它越近的点与它组成的边权重越大，可以使用RBF将其它点与它之间的欧式距离进行映射，得到

与相连的所有点的边的权重

结果如下图所示：

fig, rbf = plt.subplots()
fig.set_size_inches(6, 6)

x1 = [np.exp(-np.square(x[i]-x[5])/2) for i in range(10)]

xi = np.linspace(x.min(),x.max(),300)
from scipy.interpolate import spline
yi = spline(x,x1,xi)

rbf.plot( xi, yi )
for i in range(10):
    rbf.annotate('$w_{6%d}$' % (i+1), xy=(x[i], x1[i]), xytext=(x[i]+0.05, x1[i]+0.001))
rbf.spines['right'].set_visible(False)
rbf.spines['top'].set_visible(False)

rbf.set_xticks([i for i in range(10)])
rbf.set_xticklabels(['$x_{%d}$' % (i+1) for i in range(10)])

rbf.set_ylabel('$w_6$')
plt.show()

可以看出，与

越近，与之相连的边的权重越大。

K近邻法

除此之外，还可以使用K近邻法度量相似程度，对于某个顶点，找出距离它最近的k个顶点，这k个顶点的权重设为1，其它的权重设为0。

高斯核函数和K近邻法在scikit-learn的谱聚类算法都已经被支持，在本文最终的例子中，使用K近邻法作用相似程度度量方法计算邻接矩阵。

三 谱聚类的优化目标

谱聚类的目标是聚类之后，子图间权重尽可能小，子图内权重尽可能大。可以用

表示聚类之后得到k个组，即k个子图。对于某个子图

，它与其它子图

之间没有边连接。

对于两个子图A和B，定义A与B间的权重：

则

与其他子图的权重为

，对于一个样本

，可以定义

度：

进而得到样本集合的度矩阵

对于子图

，组内的权重可以表示为：

对于子图

，由于希望图内权重尽可能大，并且与其它图之间的权重尽可能小，可构建优化目标函数：

以子图

为例，

将

和

进行统一，定义一个m行的列向量

它的第 i 行表示为：

代入可得：

令

，L 就是有名的拉普拉斯矩阵，那么

对k个子图的优化目标进行求和，那么对于整个图，优化目标可以表示为

于是我们对于整张图，我们的优化目标变成了：

考虑到指示向量h与度d的定义，可以得到：

而这个式子很眼熟，与主成分分析PCA的优化目标仅仅的差异在于约束条件多了一个D，接下来把差异的部分去掉：

优化条件可转化为

将

看成一个矩阵B，那么优化目标就与PCA的优化目标完全一样了，要求最小值，只需要找到矩阵B的k个最小的特征值就完成了。

而这k个特征值对应的特征向量就组成了m×k的矩阵F，F就是谱聚类中的“谱”，F的m个行向量对应图中的m个顶点。最后， 使用传统的聚类算法如K-Means 将F的m个行向量进行聚类，就得到了最终的聚类结果。

四 算法步骤总结

根据前面的分析，可以总结出，算法的步骤如下：

1. 计算出邻接矩阵W
2. 根据W，计算出度矩阵D
3. 根据D与W，计算出拉普拉斯矩阵L
4. 根据L与D，计算出

1. 求出B最小的k个特征值对应的k个特征向量组成的矩阵F
2. 对F进行聚类（如k-means），得到最终的结果

五 实例

以scikit-learn众多聚类实例中的圆圈分类为例。有1500个样本，样本分布如图所示：

import time
import warnings

import numpy as np
import matplotlib.pyplot as plt

from sklearn import cluster, datasets, mixture
from sklearn.preprocessing import StandardScaler
from itertools import cycle, islice

np.random.seed(0)

n_samples = 1500
noisy_circles = datasets.make_circles(n_samples=n_samples, factor=.5,
                                      noise=.05)[0]

x1 = noisy_circles[:,0]
x2 = noisy_circles[:,1]

originX = np.array([x1, x2]).T

fig, ax = plt.subplots()
fig.set_size_inches(4, 4)

ax.scatter(x1, x2)

现在需要对数据集进行聚类。可以从分布图看出，样本聚集在2个同心圆上。所以，理想的分类结果应该是聚成2类，每个圆上的样本聚成一类。

按照谱聚类的步骤，先算出邻接矩阵W：

from sklearn.neighbors import kneighbors_graph

connectivity = kneighbors_graph(originX, n_neighbors=10, include_self=False, mode = 'distance')
W = 0.5 * (connectivity + connectivity.T)
W = W.toarray()

利用邻接矩阵W，计算出度矩阵D：

D = np.zeros(W.shape)
for i in range(W.shape[0]):
    for j in range(W.shape[0]):
        D[i, i] = D[i, i]+W[i,j]

使用W与D，计算出拉普拉斯矩阵L：

L = D-W

计算

tmp = np.zeros(D.shape)
for i in range(tmp.shape[0]):
    tmp[i,i] = 1.0/np.sqrt(D[i,i])
B = tmp@L@tmp

for i in range(B.shape[0]):
    for j in range(B.shape[1]):
        B[i,j] = '%.9f'%(B[i,j])

求出B的2（k=2）个最小的特征值，即这2个特征对应的特征向量组成的矩阵F：

from numpy.linalg import eig
eigenvalue,eigenvector=eig(B)

first = np.where(eigenvalue == np.min(eigenvalue))[0][0]
second = np.where(eigenvalue == np.max(eigenvalue))[0][0]
for i in range(eigenvalue.shape[0]):
    if i == first:
        continue;
    if eigenvalue[i]<eigenvalue[second]:
        second=i;

F = np.array([eigenvector[:,first], eigenvector[:,second]], np.float64).T

最后，用K-Means对F中1500个样本进行聚类，结果如图所示：

from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=2, random_state=0).fit(F)

colors = []
for i in range(kmeans.labels_.shape[0]):
    if kmeans.labels_[i] == 1:
        colors.append('red')
    else:
        colors.append('blue')

fig, ax = plt.subplots()
fig.set_size_inches(4, 4)

ax.scatter(x1, x2, c=colors, marker='o')

从图中可以看出，算法成功地将样本聚成2类，每一类分别聚集在2个圆上，与我们期待的结果一致。 实际上，scikit-learn已经封装好了算法的所有步骤，在实际使用的时候，只需要调用scikit-learn的SpectralClustering类就行了，以下是scikit-learn的计算方法：

from sklearn.cluster import SpectralClustering

originX = StandardScaler().fit_transform(originX)
clustering = SpectralClustering(n_clusters=2,
        eigen_solver='arpack',
        affinity="nearest_neighbors").fit(originX)


colors = []
for i in range(clustering.labels_.shape[0]):
    if clustering.labels_[i] == 1:
        colors.append('red')
    else:
        colors.append('blue')

fig, ax = plt.subplots()
fig.set_size_inches(4, 4)

ax.scatter(x1, x2, c=colors, marker='o')

最终的聚类结果如下：

今天我们的谱聚类就介绍完了。 ARGO，专注AI与金融，用Fintech改变世界。欢迎大家关注我们，与ARGO一同进步。