谱聚类代码实现（Python）

1、前言

上一篇谱聚类原理推导，我们已经对谱聚类的原理进行了推导，本文将利用谱聚类，基于Python，实现聚类。

2、实现流程

2.1、RatioCut

①求出邻接矩阵w，度矩阵D

②计算拉普拉斯矩阵L

③求出L的特征值，并取出前$k_1$个特征值所对应的特征向量组成Y。

④以Y作为样本，使用传统的聚类方法进行聚类，比如K_mean聚类，聚成$k_2$类

2.2、Ncut

①求出邻接矩阵w，度矩阵D

②计算拉普拉斯矩阵L

③归一化拉普拉斯矩阵（$D^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}$）

③求出$D^{-\frac{1}{2}}L{D}^{-\frac{1}{2}}$的特征值，并取出前$k_1$个特征值所对应的特征向量组成F。

④将F按行归一化。

④以F作为样本，使用传统的聚类方法进行聚类，比如K_mean聚类，聚成$k_2$类。

3、代码实现

3.1、RatioCut

先看结果

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_circles
def kernel(x1,x2):
    '''
    高斯核函数
    :param x1: 输入x1
    :param x2: 输入x2
    :return:  返回相乘的值
    '''
    sigma=0.1 #高斯核函数参数
    L2=np.sum(x1**2,1).reshape(-1,1)+np.sum(x2**2,1)-2*x1@x2.T
    result=np.exp(-L2/(2*sigma**2))
    return result

class K_mean():
    '''
    k_mean 聚类
    '''
    def __init__(self):
        pass
    def train(self,Y,k):
        mean=np.mean(Y, axis=0) #求均值
        std = np.std(Y, axis=0) #标准差
        Y=(Y-mean)/std #标准化
        row=Y.shape[0] #一共有多少数据
        centers = Y[np.random.choice(row, size=k), :] #随机抽取k个点作为中心点
        result = np.zeros(shape=(row, k)) #初始化用于储存结果
        while True:
            old_center = centers.copy() #复制一份中心点
            for i in np.arange(k):  #迭代每一个中心点
                distance=Y-centers[i,:]  #每一个样本都减去中心点
                L2=np.linalg.norm(distance,ord=2,axis=1) #求L2范数
                result[:,i]=L2 #保存结果
            min_position=np.argmin(result,axis=1) #比较，查看该样本点属于哪一类

            #迭代更新每一个中心点
            for j in np.arange(k):
                #取出索引
                index=np.argwhere(min_position==j).reshape(-1)
                #按行取出对应的样本值并求出均值
                mean=Y[index,:].mean(axis=0)
                #更新样本点
                centers[j,:]=mean
            #如果中心点不再更新，则模型结束
            if (old_center==centers).all():
                break

        #绘图
        plot_figure(x,min_position)

class Spectral_cluster():
    '''谱聚类'''
    def __init__(self):
        pass
    def train(self,X,k):
        W=kernel(X,X) #邻接矩阵W
        D=np.diag(np.sum(W,axis=1)) #度矩阵D
        L=D-W #拉普拉斯矩阵

        eigenvalues,eigenvectors=np.linalg.eig(L) #计算特征值和特征向量

        min_k=np.argsort(eigenvalues)[:k] #取出前小k个特征值对应索引
        Y=eigenvectors[:,min_k] #取出前k列特征向量作为Y

        k_mean=K_mean() #初始化k_mean聚类方法
        k_mean.train(Y,2) #实现聚类，聚成两类
def plot_figure(x,label):
    #绘图函数
    color_map={0:"r",1:"g",2:"b"}
    color=[color_map[i]  for i in label]
    plt.scatter(x[:,0],x[:,1],c=color)
    plt.show()
if __name__ == '__main__':
    k=2  #保留前k个特征值
    x, y =make_circles(n_samples=1000, factor=0.5, noise=0.05) #生成数据集
    spectral_cluster=Spectral_cluster() #初始化模型
    spectral_cluster.train(x,k) #导入模型训练集

3.2、NCut

结果

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_circles
def kernel(x1,x2):
    '''
    高斯核函数
    :param x1: 输入x1
    :param x2: 输入x2
    :return:  返回相乘的值
    '''
    sigma=0.1 #高斯核函数参数
    L2=np.sum(x1**2,1).reshape(-1,1)+np.sum(x2**2,1)-2*x1@x2.T
    result=np.exp(-L2/(2*sigma**2))
    return result

class K_mean():
    '''
    k_mean 聚类
    '''
    def __init__(self):
        pass
    def train(self,Y,k):
        mean=np.mean(Y, axis=0) #求均值
        std = np.std(Y, axis=0) #标准差
        Y=(Y-mean)/std #标准化
        row=Y.shape[0] #一共有多少数据
        centers = Y[np.random.choice(row, size=k), :] #随机抽取k个点作为中心点
        result = np.zeros(shape=(row, k)) #初始化用于储存结果
        while True:
            old_center = centers.copy() #复制一份中心点
            for i in np.arange(k):  #迭代每一个中心点
                distance=Y-centers[i,:]  #每一个样本都减去中心点
                L2=np.linalg.norm(distance,ord=2,axis=1) #求L2范数
                result[:,i]=L2 #保存结果
            min_position=np.argmin(result,axis=1) #比较，查看该样本点属于哪一类

            #迭代更新每一个中心点
            for j in np.arange(k):
                #取出索引
                index=np.argwhere(min_position==j).reshape(-1)
                #按行取出对应的样本值并求出均值
                mean=Y[index,:].mean(axis=0)
                #更新样本点
                centers[j,:]=mean
            #如果中心点不再更新，则模型结束
            if (old_center==centers).all():
                break

        #绘图
        plot_figure(x,min_position)

class Spectral_cluster():
    '''谱聚类'''
    def __init__(self):
        pass
    def train(self,X,k):
        W=kernel(X,X) #邻接矩阵W
        D=np.diag(np.sum(W,axis=1)) #度矩阵D
        L=D-W #拉普拉斯矩阵

        #############
        #计算D^(-0.5)
        values,vector=np.linalg.eig(D)
        V=np.diag(values**(-0.5))
        D_two=vector*V*np.linalg.inv(vector)
        #############

        L=D_two@L@D_two #归一化拉普拉斯矩阵
        eigenvalues,eigenvectors=np.linalg.eig(L) #计算特征值和特征向量

        min_k=np.argsort(eigenvalues)[:k] #取出前小k个特征值对应索引
        Y=eigenvectors[:,min_k] #取出前k列特征向量作为Y

        Y=Y/(np.sum(Y,axis=1).reshape(-1,1)) #对所得Y按行归一化

        k_mean=K_mean() #初始化k_mean聚类方法
        k_mean.train(Y,2) #实现聚类，聚成两类
def plot_figure(x,label):
    #绘图函数
    color_map={0:"r",1:"g",2:"b"}
    color=[color_map[i]  for i in label]
    plt.scatter(x[:,0],x[:,1],c=color)
    plt.show()
if __name__ == '__main__':
    k=2  #保留前k个特征值
    x, y =make_circles(n_samples=1000, factor=0.5, noise=0.05) #生成数据集
    spectral_cluster=Spectral_cluster() #初始化模型
    spectral_cluster.train(x,k) #导入模型训练集