谱聚类初探

本学期初参加了我校的一个大学生科技活动。
导师甩给我的课题是：“谱聚类算法”，我自己选择的应用方向是“图像分割”。

在网上大量浏览相关帖子后，终于对该算法有了一个大概的掌握。
https://www.cnblogs.com/pinard/p/6221564.html
个人觉得上面这篇博客对算法总体的阐述比较好。
可以作为了解该算法的一篇入门文章。

所谓聚类，如名字一样，就是将杂乱无章的数据，按照彼此之间的相似度，分为几类。其中大量使用了矩阵变换的知识。

输入：

样本集D=(x1,x2,…,xn)(x1,x2,…,xn)。

我是直接从.mat文件读取数据，也可以用matplotlib.pyplot库的imread函数读取普通图片。

写作：

import matplotlib.pyplot as plt
pic = plt.imread('one.png')

输出：

簇划分C(c1,c2,…ck2)C(c1,c2,…ck2).
再使用

plt.show()

便可以将图像显示出来。

算法如下：

• 根据输入的相似矩阵的生成方式构建样本的相似矩阵S（网上好多帖子都提示使用“高斯相似”、但由于常量σ不好取值，故使用了“欧式距离公式”，来表示两点间相似度。）
• 计算出度矩阵D
• 计算出拉普拉斯矩阵L=D-W
• 计算出L前k个最小的特征向量v_1,…,v_k
• 将前k个特征向量组合成一个矩阵V，其中第i列对应v_i列向量。
• 该矩阵V的每一行对应代表x_i的低维度的表示y_i。
• 对所有y_i进行k-means聚类，聚成k类

代码如下：

该代码非个人创作，而是借鉴了各位前辈的思路和方法。

加之个人理解消化了一下，加入了详细的注释。

# coding=utf-8
import scipy.io as sio
from PIL import Image
from sklearn.cluster import KMeans
import numpy as np
import random
import copy
import matplotlib.pyplot as plt
import operator

def loaddata():
    #读取数据
    file = 'GaussianData.mat'
    #文件名
    file = 'ringData.mat'
    data = sio.loadmat(file)
    matrix = np.array(data['Dataset'])
    #matrix = np.array(L)
    return matrix
def distance(p1,p2):  
    # 欧式距离
    return  np.linalg.norm(p1-p2)
def getWbyKNN(data,k):  
    #data:点阵坐标信息
    points_num = len(data)
    #点数量：300
    dis_matrix = np.zeros((points_num,points_num))
    #创建矩阵300*300
    W = np.zeros((points_num,points_num))
    #创建矩阵300*300
    for i in range(points_num):
        for j in range(i+1,points_num):

            dis_matrix[i][j] = dis_matrix[j][i] = distance(data[i],data[j])
    #将每个点与其他点的相似度表示成矩阵
    for idx,each in enumerate(dis_matrix):
        index_array  = np.argsort(each)
        #排序
        W[idx][index_array[1:k+1]] = 1
        #找出距离第idx个点距离最近的k个点
        #距离最短的是自己，所以取1到k+1
    tmp_W = np.transpose(W)
    #转置
    W = (tmp_W+W)/2
    #转置相加除以2是为了让矩阵对称
    return W
def getD(W):
    #获得度矩阵
    points_num = len(W)
    D = np.diag(np.zeros(points_num))
    #转换成对角阵
    for i in range(points_num):
        D[i][i] = sum(W[i])
    return D
def getEigVec(L,cluster_num):
    #从拉普拉斯矩阵获得特征矩阵
    eigval,eigvec = np.linalg.eig(L)
    #特征值、特征向量
    lenght = len(eigval)
    dictEigval = dict(zip(eigval,range(0,lenght)))
    kEig = np.sort(eigval)[0:cluster_num]
    #排序
    ix = [dictEigval[k] for k in kEig]
    return eigval[ix],eigvec[:,ix]
    #取特征值和对应的列向量
def randRGB():
    #添加颜色，便于区分
    return (random.randint(0, 255)/255.0,
            random.randint(0, 255)/255.0,
            random.randint(0, 255)/255.0)
def plot(matrix,C,k):
    colors = []
    for i in range(k):
        colors.append(randRGB())
    for idx,value in enumerate(C):
        plt.plot(matrix[idx][0],matrix[idx][1],'o',color=colors[int(C[idx])])
    plt.show()

if __name__ == '__main__':
    cluster_num = 2
    KNN_k = 10
    data = loaddata()
    W = getWbyKNN(data,KNN_k)
    #相似度矩阵
    D = getD(W)
    #度矩阵
    L = D-W
    #拉普拉斯矩阵
    eigval,eigvec = getEigVec(L,cluster_num)
    clf = KMeans(n_clusters=cluster_num)
    #新建类对象
    s = clf.fit(eigvec.real)
    #开始聚类
    plot(data,s.labels_,cluster_num)