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

十二月的猫

于 2024-10-01 16:46:59 发布

阅读量594

点赞数 10

分类专栏：机器学习深度学习文章标签：机器学习算法聚类

本文链接：https://blog.csdn.net/m0_67656158/article/details/142668953

版权

深度学习同时被 2 个专栏收录

15 篇文章 1 订阅

订阅专栏

机器学习

4 篇文章 0 订阅

订阅专栏

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

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

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

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

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

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

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

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 plt

def load_data(filename):
    """
    载入数据
    :param filename: 文件名
    :return: numpy array 格式的数据
    """
    data = np.loadtxt(filename, delimiter='\t')
    return data

def 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) / 2

def 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 - W

def 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_predict
    return label

if __name__ == '__main__':
    cluster_num = 7
    knn_k = 5
    filename = '../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)相反，谱聚类将数据从高维降到低维空间；尽管这两者都是为了使得投影后的数据线性可分，但是使用的方法却是相反的。

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

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