【经典机器学习算法】谱聚类算法及其实现（python）

🌈 个人主页：十二月的猫-CSDN博客
🔥 系列专栏： 🏀深度学习_十二月的猫的博客-CSDN博客

💪🏻 十二月的寒冬阻挡不了春天的脚步，十二点的黑夜遮蔽不住黎明的曙光

目录

1. 前言

2. 前提知识

2.1  邻接矩阵

2.2 度与度矩阵

2.3 矩阵的相似

2.4 连通子图

3. 相似度的衡量方法

3.1 ​​近邻法

3.2 k近邻法

3.3 高斯核函数

4. 图拉普拉斯矩阵(Graph Laplacian Matrix)

4.1 非规范化的图拉普拉斯矩阵

 4.2 非规范化的图拉普拉斯矩阵的性质

5. 谱聚类（无向图切割）

5.1 谱聚类切割目标（优化目标-loss） 

5.2 谱聚类算法思想

5.2.1 RatioCut切图

 5.2.2 Ncut切图

5.2.3 总结

6. 谱聚类算法实现（基于python实现） 

7. 总结

---

1. 前言

在看一篇论文的过程中，遇到一个问题：

“已知数据集，要求将数据集分为几组，要求组间距离最大，组内距离最小”

这是一个无监督问题，在查阅资料后，认为聚类可以帮我解决这个问题

谱聚类的思想来源于图论，它把待聚类的数据集中的每一个样本看做是图中一个顶点，这些顶点连接在一起，连接的这些边上有权重，权重的大小表示这些样本之间的相似程度。同一类的顶点它们的相似程度很高，在图论中体现为同一类的顶点中连接它们的边的权重很大，不在同一类的顶点连接它们的边的权重很小。

于是谱聚类的最终目标就是找到一种切割图的方法，使得切割之后的各个子图内的权重很大，子图之间的权重很小。

​

可以看出，数据集总共分为2类左右，沿着图中蓝色线切割可以得到结果，这种切割的前提是这两个类之间的顶点，比如顶点i和j之间的权重最小，即Wij最小。

2. 前提知识

假设给定一个数据集X={x1,x2,…,xn}，其中每一个样本 xi∈ 。按照图论的思想，我们将这个 n 个数据向量当做 m 维空间中某一幅无向图上的一个个点，因为我们的目的是衡量这些点之间的相似性，所以本文把这幅图叫做相似图，记为 G=(V,E) ，其中 V={v1,v2,…,vn} 表示顶点， E 表示边的集合。连接两个顶点 vi 和 vj 的边的权重记为，它们的相似性用 表示，该相似性的值越大，说明它们越相似，反之则越不相似。

本文要求边的权重 wij≥0 ，权重等于0表示俩顶点无连接，则 n 个顶点的权重构成一个矩阵 W=(),i,j=1,2,…,n ，这个矩阵将在下文出现。

这里和有直接关系

2.1  邻接矩阵

对于一幅无向图G=(V,E)​​，学过图论或者数据结构的同学都知道，他有两个很重要的概念是图的邻接矩阵和顶点的度。所有顶点之间的权重构成一个n×n​​​矩阵，叫做邻接矩阵，也叫权重矩阵，即：

对于无向图，顶点vi​​与顶点vj​​之间的权重和顶点vj​​与顶点vi​​之间的权重是一样的，因而wij=wji​​，因此W​​是对称矩阵，即W=WT​​​。顶点自己到自己的权重是多少呢？这里先按下不表。这个邻接矩阵稍后将作为图的相似矩阵。注意这里的相似矩阵不是矩阵的相似。

相似矩阵：由点之间的相似值sij来组成的矩阵

矩阵的相似：两个矩阵，也就是两个图是否相似的定量衡量

2.2 度与度矩阵

在数据结构中，度定义为与该顶点直接连接的顶点的个数，或者是连接到该顶点的边的个数。不过不采用这个定义。对于某个顶点di,i=1,2,…,n​而是将度定义为：

从公式(2)可以看出，顶点vi的度其实就是邻接矩阵第i行的和(第i列的和也可以，因为W​是对称矩阵)。

度矩阵定义为n个度构成的对角矩阵，即：

相似矩阵对角线上的值：本行所有wij求和

2.3 矩阵的相似

给定顶点V的一个子集A⊂V，我们定义它的补为。再给定顶点V的一个子集B⊂V，我们定义它的补为，对于2个子集A和B，我们定义：

公式(4)表示两个子集中顶点之间的权重之和，注意这里不包含子集内顶点之间的权重。

子集的大小有两种定义：

• 子集内顶点的个数，记为|A|。
• 子集内所有顶点的度之和，记为：。

2.4 连通子图

对于一个非空子集A⊂V​，如果A中的任意两个顶点都至少存在一条路径将它们连接起来，并且A中的其它顶点也在这条路径上，则称A是连接的。如果子集A是连接的，并它与它的补A¯​不存在任何的连接。则称A​是一个连通子图。非空子集A1,A2,…,Ak构成图V的一个分割，用数学公式来写就是A1∪A2,…,∪Ak=V。

3. 相似度的衡量方法

wij：表示vi、vj两个点之间的权重

sij：表示vi、vj两个点之间的相似度

权重就是相似度，相似度越大权重越大

图中各个顶点的相似度衡量主要基于距离的度量，也就是说空间两个点的距离越近，则它们越相似，距离越远，则它们越不相似，即相似度与距离成反比，所以只要你使用的度量空间具有这种性质，都可以作为相似度的衡量方法。下面介绍三种相似度的衡量方法，同时也是相似矩阵的计算方法。

3.1 ​​近邻法

该方法采用欧式距离计算两个顶点的距离，然后设定一个阈值ϵ，使得：

从公式(5)可以看出，由此得到的相似矩阵其元素要么是0要么是ϵ，这种方法获得权重信息量太少了，一般很少使用。

缺陷：相似度不是一个连续的变量，且只有一个固定的值

3.2 k近邻法

该方法取与顶点最近的k​个顶点，该顶点与这k个顶点的权重都大于0，但这会导致最后所得的相似矩阵不一定是对称的，因为一个点vi在另外一个点vj的k个近邻中，并不能保证vj也在vi的k​​个近邻中。有两种可以保证所得的相似矩阵对称：

• 两个顶点vi​与vj​只要其中一个点在另外一个点的k​​个近邻中，则令wij=wji​，只有这两个顶点同时都不在任何一方的k​个近邻中，则令wij=wji=0​​​。综合可得：

方法本质：增加限制条件，保证其一定是对称的  

•  两个顶点vi与vj只同时在双方的k个近邻中，则令wij=wji，只要有一方不在另外一方的k个近邻中，则令wij=wji=0。综合：

3.3 高斯核函数

考虑到相似度计算的问题在于：

1、保证对称

2、和距离呈反函数

3、不论什么维度都要能够计算距离，从而计算相似度

到这里不难想到：高斯核函数

该方法将所有的顶点都连接起来。然后通过度量空间中某种对称度量算子来计算顶点之间的相似度。比如使用高斯核函数计算两个顶点之间的相似度：

注意，这里的是一个标量，标量的转置仍然是它自身，所以公式(8)​是一个对称的度量算子。为什么要求是对称的度量的算子，因为要保证租后得到的相似矩阵是相似的。

4. 图拉普拉斯矩阵(Graph Laplacian Matrix)

4.1 非规范化的图拉普拉斯矩阵

图拉普拉斯矩阵的定义比较简单，即：

 

其中D是公式(3)​的度矩阵，W​是公式(1)的权重矩阵(相似矩阵) 

举个例子，给定下面的图：

把此“图”转换为邻接矩阵的形式，记为：W

把的每一列元素加起来得到个数，然后把它们放在对角线上（其它地方都是零），组成一个N × N N \times NN×N对角矩阵，记为度矩阵D DD，如下图所示：

根据拉普拉斯矩阵的定义L = D − W L=D-WL=D−W，可得拉普拉斯矩阵 L LL为：

 4.2 非规范化的图拉普拉斯矩阵的性质

(1)对于任意的向量f∈Rn，有：

 

(2)L​是一个对称半正定矩阵。

因为经过相似矩阵W​的各种求法可知，其元素wij​是非负数，所以由公式(10)可知​：

恒成立。从而L是一个对称半正定矩阵。

补充一下正定矩阵的作用：

很多时候，我们在机器学习/深度学习/优化问题中需要计算最优解，要怎么判断我们所求的解就是最优解呢？

这里需要引入：黑塞矩阵（Hessian）

 黑塞矩阵（Hessian）：

1. 如果是正定矩阵，则临界点处是一个局部极小值
2. 如果是负定矩阵，则临界点处是一个局部极大值
3. 如果是不定矩阵，则临界点处不是极值

 (3)L的最小特征值为0，对应的特征向量为全1向量1。

 

所以，矩阵L的0特征值对应的特征向量为1。 

补充定理1：对于一个分块对角矩阵A​：

它的特征值等于各个分块矩阵Ai,i=1,2,…,n​的特征值。

5. 谱聚类（无向图切割）

一张图，如下：

将其分为几组，可以理解为：1、由单个点去聚合；2、由整张图去切割 

回收前面提到的“矩阵的相似”：

这里我们切割的目的就是：要让切割后的子图之间的相似程度最小，子图内的相似程度最大

切割子图之间的相似程度定义如下：

定义 A 和 B是图 G 中两个子图，则定义子图A和 B的切图权重为：

那么对于我们k个子图的集合：A 1 , A 2 , . . . , A k，我们定义切图 cut 为：

5.1 谱聚类切割目标（优化目标-loss） 

 那么如何切图可以让子图内的点权重和高，子图间的点权重和低呢？一个自然的想法就是最小化c u t ( A 1 , A 2 , . . . , A k )，但是可以发现，这种极小化的切图存在问题，如下图：

问题出现本质：没有考虑算法内聚性，没有让子图内的权重尽量高 

容易确保切割数量与cut函数的关系不是单调的，存在极值点：

1、当子图数量增加，则需要增加考虑子图间的cut值

2、当子图数量减少，需要增加考虑子图内部的连接强度

5.2 谱聚类算法思想

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

5.2.1 RatioCut切图

RatioCut切图为了避免上面出现的最小切图，对每个切图，不光考虑最小化cut( A 1,A 2 , ..,A k )它还同时考虑最大化每个子图点的个数，即：

最小化这个函数即可。

 5.2.2 Ncut切图

Ncut切图和RatioCut切图很类似，但是把Ratiocut的分母 ∣ A i ∣换成。由于子图样本的个数多并不一定权重就大，我们切图时基于权重也更合我们的目标，因此一般来说Ncut切图优于RatioCut切图。

5.2.3 总结

引入子图内连接强度：

的本质就可以用这个intra connect(A)来代替

本质上：除上intra connect(A)和|Ai|的目的都是考虑上子图内部的内聚性

6. 谱聚类算法实现（基于python实现） 

import numpy as np
from sklearn.cluster import KMeans
import matplotlib.pyplot as pltdef load_data(filename):"""载入数据:param filename: 文件名:return: numpy array 格式的数据"""data = np.loadtxt(filename, delimiter='\t')return datadef distance(x1, x2):"""获得两个样本点之间的欧几里得距离:param x1: 样本点1:param x2: 样本点2:return: 两个样本点之间的距离"""return np.linalg.norm(x1 - x2)def get_dist_matrix(data):"""获取距离矩阵:param data: 样本集合:return: 距离矩阵"""return np.linalg.norm(data[:, np.newaxis] - data[np.newaxis, :], axis=-1)def getW(data, k):"""获得邻接矩阵 W:param data: 样本集合:param k: KNN参数:return: 邻接矩阵 W"""n = len(data)dist_matrix = get_dist_matrix(data)W = np.zeros((n, n))for idx in range(n):# 获取最近k个邻居的索引idx_array = np.argsort(dist_matrix[idx])[1:k+1]  # 跳过自己W[idx, idx_array] = 1# 确保邻接矩阵是对称的return (W + W.T) / 2def getD(W):"""获得度矩阵:param W: 邻接矩阵:return: 度矩阵 D"""return np.diag(np.sum(W, axis=1))def getL(D, W):"""获得拉普拉斯矩阵:param D: 度矩阵:param W: 邻接矩阵:return: 拉普拉斯矩阵 L"""return D - Wdef getEigen(L, cluster_num):"""获得拉普拉斯矩阵的特征向量:param L: 拉普拉斯矩阵:param cluster_num: 聚类数目:return: 选定特征值对应的特征向量"""eigval, eigvec = np.linalg.eig(L)ix = np.argsort(eigval)[:cluster_num]  # 选择最小的cluster_num个特征值的索引return eigvec[:, ix]def plotRes(data, clusterResult, clusterNum):"""结果可视化:param data: 样本集:param clusterResult: 聚类结果:param clusterNum: 聚类个数"""scatterColors = ['black', 'blue', 'green', 'yellow', 'red', 'purple', 'orange']for i in range(clusterNum):color = scatterColors[i % len(scatterColors)]plt.scatter(data[clusterResult == i, 0], data[clusterResult == i, 1], c=color, marker='+')plt.title(f'Clustering Result with {clusterNum} clusters')plt.xlabel('Feature 1')plt.ylabel('Feature 2')plt.show()def cluster(data, cluster_num, k):"""聚类函数:param data: 输入数据:param cluster_num: 聚类数目:param k: KNN参数:return: 聚类标签"""W = getW(data, k)D = getD(W)L = getL(D, W)eigvec = getEigen(L, cluster_num)# 使用KMeans进行聚类clf = KMeans(n_clusters=cluster_num)label = clf.fit_predict(eigvec)  # 直接使用fit_predictreturn labelif __name__ == '__main__':cluster_num = 7knn_k = 5filename = '../data/Aggregation_cluster=7.txt'data = load_data(filename=filename)data = data[:, :-1]  # 去除最后一列（假设为标签列）label = cluster(data, cluster_num, knn_k)plotRes(data, label, cluster_num)

运行结果如下：

7. 总结

以上就是整个谱聚类的原理介绍、分析、实现和讨论。其本质呢还是从数据中构造某种相似矩阵(类比协方差矩阵)，然后对矩阵进行特征分解，为去掉冗余特征，再做投影(降维)，抓住主要成分，注意和PCA的区别，PCA的目的是用最少的特征尽可能地表示最多的信息(对应前几个最大的特征值)，而谱聚类是要求切图耗费的能量最少(对应前几个最小特征值)。

最后是谱聚类的一些问题：

(1)和k-means一样都要选择类别数/分组数k。

(2)选择相似性矩阵的度量方式，度量方式不同得到的图拉普拉斯矩阵不同，可能会导致不对称。

(3)可以看到，谱聚类在投影之后还是需要其他聚类方法介入，其实可以这么认为，谱聚类前面的这些工作可以看做是数据预处理的过程，而后再使用经典的聚类方法如k-means等。

(4)谱聚类对于非凸数据聚类很有用(请看前面的几个例子)。

(5)和支持向量机将数据投影到高维空间(kernel trick)相反，谱聚类将数据从高维降到低维空间；尽管这两者都是为了使得投影后的数据线性可分，但是使用的方法却是相反的。

---

 撰写文章不易，如果文章能帮助到大家，大家可以点点赞、收收藏呀~

十二月的猫在这里祝大家学业有成、事业顺利、情到财来