【经典机器学习算法】谱聚类算法及其实现（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∈ $R^m$ 。按照图论的思想，我们将这个 n 个数据向量当做 m 维空间中某一幅无向图上的一个个点，因为我们的目的是衡量这些点之间的相似性，所以本文把这幅图叫做相似图，记为 G=(V,E) ，其中 V={v1,v2,…,vn} 表示顶点， E 表示边的集合。连接两个顶点 vi 和 vj 的边的权重记为 $w_{ij}$ ，它们的相似性用 $s_{ij}$ 表示，该相似性的值越大，说明它们越相似，反之则越不相似。

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

这里 $w_{ij}$ 和 $s_{ij}$ 有直接关系

2.1 邻接矩阵

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

$W=\begin{bmatrix}w_{11}&w_{12}&\ldots&w_{1n}\\w_{21}&w_{22}&\ldots&w_{2n}\\\vdots&\vdots&\ldots&\vdots\\w_{n1}&w_{n2}&\ldots&w_{nn}\end{bmatrix}\quad(1)$

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

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

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

2.2 度与度矩阵

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

${d_i=\sum_{j=1}^nw_{ij}\quad(2)}$

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

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

${D=\begin{bmatrix}d_1&0&0&0&\ldots&0\\0&d_2&0&0&\ldots&0\\0&0&d_3&0&\ldots&0\\0&0&0&d_4&\ldots&0\\\vdots&\vdots&\vdots&\vdots&\ldots&\vdots\\0&0&0&0&\ldots&d_n\end{bmatrix}}\text{(3)}$

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

2.3 矩阵的相似

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

$W(A,B):=\sum_{\boldsymbol{v}_i\in A,\boldsymbol{v}_j\in B}w_{ij}\quad(4)$

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

子集的大小有两种定义：

子集内顶点的个数，记为|A|。
子集内所有顶点的度之和，记为： $\mathrm{vol}(A):=\sum_{\boldsymbol{v}_i\in A}d_i\text{}$ 。

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 $\epsilon$ 近邻法

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

$w_{ij}=\left\{\begin{array}{ll}0,&\text{if}s_{ij}>\epsilon\\\epsilon,&\text{if}s_{ij}\leq\epsilon\end{array}\right.\quad(5)$

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

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

3.2 k近邻法

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

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

$\left.w_{ij}=w_{ji}=\left\{\begin{array}{ll}0,&\boldsymbol{v}_i\not\in\mathrm{knn}(\boldsymbol{v}_j)\mathrm{~and~}\boldsymbol{v}_j\not\in\mathrm{knn}(\boldsymbol{v}_i)\\\frac{1}{s_{ij}},&\boldsymbol{v}_i\in\mathrm{knn}(\boldsymbol{v}_j)\mathrm{~or~}\boldsymbol{v}_j\in\mathrm{knn}(\boldsymbol{v}_i)\end{array}\right.\right.\quad(6)$

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

两个顶点vi与vj只同时在双方的k个近邻中，则令wij=wji，只要有一方不在另外一方的k个近邻中，则令wij=wji=0。综合： $\left.w_{ij}=w_{ji}=\left\{\begin{array}{ll}0,&\boldsymbol{v}_i\not\in\mathrm{knn}(\boldsymbol{v}_j)\mathrm{~or~}\boldsymbol{v}_j\not\in\mathrm{knn}(\boldsymbol{v}_i)\\\frac{1}{s_{ij}},&\boldsymbol{v}_i\in\mathrm{knn}(\boldsymbol{v}_j)\mathrm{~and~}\boldsymbol{v}_j\in\mathrm{knn}(\boldsymbol{v}_i)\end{array}\right.\right.\quad(7)$

3.3 高斯核函数

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

1、保证对称

2、和距离呈反函数

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

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

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

$w_{ij}=w_{ji}=e^{-\frac{1}{2}[(\boldsymbol{v}_{i}-\boldsymbol{v}_{j})^{T}\Sigma^{-1}(\boldsymbol{v}_{j}-\boldsymbol{v}_{j})]}\quad(8)$

注意，这里的 $\boldsymbol{v}_{i}^{T}\Sigma^{-1}\boldsymbol{v}_{j}$ 是一个标量，标量的转置仍然是它自身，所以公式(8)是一个对称的度量算子。为什么要求是对称的度量的算子，因为要保证租后得到的相似矩阵是相似的。

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

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

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

$L = D - W \quad(9)$

其中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，有：

$\begin{aligned} f^{T}L\boldsymbol{f}& =\boldsymbol{f}^T(D-W)\boldsymbol{f} \\ &=\boldsymbol{f}^TD\boldsymbol{f}-\boldsymbol{f}^TW\boldsymbol{f} \\ &=\sum_{i=1}^nd_if_i^2-\sum_{i,j=1}^nw_{ij}f_if_j \\ &=\frac12\left(\sum_{i=1}^nd_if_i^2-2\sum_{i,j=1}^nw_{ij}f_if_j+\sum_{j=1}^nd_jf_j^2\right) \\ &=\frac12\left[\sum_{i=1}^n(\sum_{j=1}^nw_{ij})f_i^2-2\sum_{i,j=1}^nw_{ij}f_if_j+\sum_{j=1}^n(\sum_{i=1}^nw_{ji})f_j^2\right] \\ &&\text{(10)} \\ &=\frac12\left[\sum_{i=1}^n(\sum_{j=1}^nw_{ij})f_i^2-2\sum_{i,j=1}^nw_{ij}f_if_j+\sum_{j=1}^n(\sum_{i=1}^nw_{ij})f_j^2\right] \\ &=\frac12\left(\sum_{i,j=1}^nw_{ij}f_i^2-2\sum_{i,j=1}^nw_{ij}f_if_j+\sum_{i,j=1}^nw_{ij}f_j^2\right) \\ &=\frac12\sum_{i,j=1}^nw_{ij}(f_i^2-2f_if_j+f_j^2) \\ &=\frac12\sum_{i,j=1}^nw_{ij}(f_i-f_j)^2 \end{aligned}$

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

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

$f^TLf\geq0\quad(11)$

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

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

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

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

黑塞矩阵（Hessian）：

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

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

$\begin{aligned} L_{1}& =(D-W)\mathbf{1} \\ &=D\mathbf{1}-W\mathbf{1} \\ &=\begin{bmatrix}d_1\\d_2\\\vdots\\d_n\end{bmatrix}-\begin{bmatrix}\sum_{j=1}w_{1j}\\\sum_{j=1}w_{2j}\\\vdots\\\sum_{j=1}w_{nj}\end{bmatrix}& (12) \\ &=\begin{bmatrix}d_1\\d_2\\\vdots\\d_n\end{bmatrix}-\begin{bmatrix}d_1\\d_2\\\vdots\\d_n\end{bmatrix} \\ &=\mathbf{0}=0*\mathbf{1} \end{aligned}$

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

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

$A=\begin{bmatrix}A_1&0&0&\ldots&0\\0&A_2&0&\ldots&0\\\vdots&\vdots&\vdots&\ldots&\vdots\\0&0&0&\ldots&A_n\end{bmatrix}\quad(13)$

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

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

一张图，如下：

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

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

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

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

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

$\mathbf{W(A,B):=\sum_{i\in A,j\in B}w_{ij}}$

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

$\mathrm{cut}(\mathbf{A}_1,\mathbf{A}_2,...,\mathbf{A}_\mathrm{k})=\frac12\sum_{\mathrm{i}=1}^\mathrm{k}\mathbf{W}(\mathbf{A}_\mathrm{i},\bar{\mathbf{A}}_\mathrm{i})$

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 )它还同时考虑最大化每个子图点的个数，即：

$\mathrm{RatiocCut}(\mathrm{A}_1,\mathrm{A}_2,...,\mathrm{A}_\mathrm{k})=\frac12\sum_{\mathrm{i}=1}^\mathrm{k}\frac{\mathrm{W}(\mathrm{A}_\mathrm{i},\bar{\mathrm{A}}_\mathrm{i})}{|\mathrm{A}_\mathrm{i}|}$

最小化这个函数即可。

5.2.2 Ncut切图

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

$\mathrm{NCut}(\mathrm{A}_1,\mathrm{A}_2,...,\mathrm{A}_\mathrm{k})=\frac12\sum_{\mathrm{i}=1}^\mathrm{k}\frac{\mathrm{W}(\mathrm{A}_\mathrm{i},\bar{\mathrm{A}}_\mathrm{i})}{\mathrm{vol}(\mathrm{A}_\mathrm{i})}$

5.2.3 总结

引入子图内连接强度：

$\text{intra connect}(A)=\sum_{u,v\in A}w(u,v)$

$vol(A_i)$ 的本质就可以用这个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)相反，谱聚类将数据从高维降到低维空间；尽管这两者都是为了使得投影后的数据线性可分，但是使用的方法却是相反的。

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

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