1、前言
上一篇谱聚类原理推导,我们已经对谱聚类的原理进行了推导,本文将利用谱聚类,基于Python,实现聚类。
2、实现流程
2.1、RatioCut
①求出邻接矩阵w,度矩阵D
②计算拉普拉斯矩阵L
③求出L的特征值,并取出前 k 1 k_1 k1个特征值所对应的特征向量组成Y。
④以Y作为样本,使用传统的聚类方法进行聚类,比如K_mean聚类,聚成 k 2 k_2 k2类
2.2、Ncut
①求出邻接矩阵w,度矩阵D
②计算拉普拉斯矩阵L
③归一化拉普拉斯矩阵( D − 1 2 L D − 1 2 D^{-\frac{1}{2}}LD^{-\frac{1}{2}} D−21LD−21)
③求出 D − 1 2 L D − 1 2 D^{-\frac{1}{2}}LD^{-\frac{1}{2}} D−21LD−21的特征值,并取出前 k 1 k_1 k1个特征值所对应的特征向量组成F。
④将F按行归一化。
④以F作为样本,使用传统的聚类方法进行聚类,比如K_mean聚类,聚成 k 2 k_2 k2类。
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) #导入模型训练集