聚类——谱聚类算法以及Python实现

谱聚类(spectral cluster)可以视为一种改进的Kmeans的聚类算法。常用来进行图像分割。缺点是需要指定簇的个数，难以构建合适的相似度矩阵。优点是简单易实现。相比Kmeans而言，处理高维数据更合适。

核心思想

构建样本点的相似度矩阵(图)，将图切割成K个子图，使得各个子图内相似度最大，子图间相似度最弱

算法简介

构建相似度矩阵的拉普拉斯矩阵。对拉普拉斯矩阵进行特征值分解，选取前K(也是簇的个数)个特征向量(按特征值从小到大的顺序)构成K维特征空间，在特征空间内进行Kmeans聚类。概括地讲，就是将原始数据映射到特征空间进行Kmeans聚类。因此，谱聚类适合于簇的个数比较小的情况下。（个人的理解是当维度比较高，簇的个数比较小时，可以将其视为一种降维的方法）
拉普拉斯矩阵可以分为规范化的($L_{norm}$)和未规范化的拉普拉斯矩阵($L$)。

$$L = D-W\\L_{norm}=D^{-1/2}LD^{-1/2}$$ 其中D为图的度矩阵(对角矩阵，节点边权重之和)，W为相似度矩阵。

算法流程

• Input: 训练数据集data，簇的个数， 阈值epsilon, 最大迭代次数maxstep, 相似度计算方法及参数
• Output: 标签数组
• Step1:构建相似度矩阵，再构建拉普拉斯矩阵，对拉普拉斯矩阵进行特征值分解，将样本数据点映射到特征空间。
• Step2: 再特征空间内进行Kmeans聚类。

代码

"""
谱聚类算法
核心思想：构建样本点的图，切分图，使得子图内权重最大，子图间权重最小
"""
import numpy as np
from kmeans import KMEANS


class Spectrum:
    def __init__(self, n_cluster, epsilon=1e-3, maxstep=1000, method='unnormalized',
                 criterion='gaussian', gamma=2.0, dis_epsilon=70, k=5):
        self.n_cluster = n_cluster
        self.epsilon = epsilon
        self.maxstep = maxstep
        self.method = method  # 本程序提供规范化以及非规范化的谱聚类算法
        self.criterion = criterion  # 相似性矩阵的构建方法
        self.gamma = gamma  # 高斯方法中的sigma参数
        self.dis_epsilon = dis_epsilon  # epsilon-近邻方法的参数
        self.k = k  # k近邻方法的参数

        self.W = None  # 图的相似性矩阵
        self.L = None  # 图的拉普拉斯矩阵
        self.L_norm = None  # 规范化后的拉普拉斯矩阵
        self.D = None  # 图的度矩阵
        self.cluster = None

        self.N = None

    def init_param(self, data):
        # 初始化参数
        self.N = data.shape[0]
        dis_mat = self.cal_dis_mat(data)
        self.cal_weight_mat(dis_mat)
        self.D = np.diag(self.W.sum(axis=1))
        self.L = self.D - self.W
        return

    def cal_dis_mat(self, data):
        # 计算距离平方的矩阵
        dis_mat = np.zeros((self.N, self.N))
        for i in range(self.N):
            for j in range(i + 1, self.N):
                dis_mat[i, j] = (data[i] - data[j]) @ (data[i] - data[j])
                dis_mat[j, i] = dis_mat[i, j]
        return dis_mat

    def cal_weight_mat(self, dis_mat):
        # 计算相似性矩阵
        if self.criterion == 'gaussian':  # 适合于较小样本集
            if self.gamma is None:
                raise ValueError('gamma is not set')
            self.W = np.exp(-self.gamma * dis_mat)
        elif self.criterion == 'k_nearest':  # 适合于较大样本集
            if self.k is None or self.gamma is None:
                raise ValueError('k or gamma is not set')
            self.W = np.zeros((self.N, self.N))
            for i in range(self.N):
                inds = np.argpartition(dis_mat[i], self.k + 1)[:self.k + 1]  # 由于包括自身，所以+1
                tmp_w = np.exp(-self.gamma * dis_mat[i][inds])
                self.W[i][inds] = tmp_w
        elif self.criterion == 'eps_nearest':  # 适合于较大样本集
            if self.dis_epsilon is None:
                raise ValueError('epsilon is not set')
            self.W = np.zeros((self.N, self.N))
            for i in range(self.N):
                inds = np.where(dis_mat[i] < self.dis_epsilon)
                self.W[i][inds] = 1.0 / len(inds)
        else:
            raise ValueError('the criterion is not supported')
        return

    def fit(self, data):
        # 训练主函数
        self.init_param(data)
        if self.method == 'unnormalized':
            w, v = np.linalg.eig(self.L)
            inds = np.argsort(w)[:self.n_cluster]
            Vectors = v[:, inds]
        elif self.method == 'normalized':
            D = np.linalg.inv(np.sqrt(self.D))
            L = D @ self.L @ D
            w, v = np.linalg.eig(L)
            inds = np.argsort(w)[:self.n_cluster]
            Vectors = v[:, inds]
            normalizer = np.linalg.norm(Vectors, axis=1)
            normalizer = np.repeat(np.transpose([normalizer]), self.n_cluster, axis=1)
            Vectors = Vectors / normalizer
        else:
            raise ValueError('the method is not supported')
        km = KMEANS(self.n_cluster, self.epsilon, self.maxstep)
        km.fit(Vectors)
        self.cluster = km.cluster
        return


if __name__ == '__main__':
    from sklearn.datasets import make_blobs
    from itertools import cycle
    import matplotlib.pyplot as plt

    data, label = make_blobs(centers=3, n_features=10, cluster_std=1.2, n_samples=500, random_state=1)
    sp = Spectrum(n_cluster=3, method='unnormalized', criterion='gaussian', gamma=0.1)
    sp.fit(data)
    cluster = sp.cluster

    # km = KMEANS(4)
    # km.fit(data)
    # cluster_km = km.cluster

    # def visualize(data, cluster):
    #     color = 'bgrym'
    #     for col, inds in zip(cycle(color), cluster.values()):
    #         partial_data = data[inds]
    #         plt.scatter(partial_data[:, 0], partial_data[:, 1], color=col)
    #     plt.show()
    #     return

    # visualize(data, cluster)

    def cal_err(data, cluster):
        # 计算MSE
        mse = 0
        for label, inds in cluster.items():
            partial_data = data[inds]
            center = partial_data.mean(axis=0)
            for p in partial_data:
                mse += (center - p) @ (center - p)
        return mse / data.shape[0]

    print(cal_err(data, cluster))
    # print(cal_err(data, cluster_km))

我的GitHub
注：代码尚未进行严格测试，如有不当之处，请指正。