kmeans聚类_KMeans与谱聚类

K-Means实现

import numpy as np
import matplotlib.pyplot as plt

def DataSet():
    data = np.loadtxt('data.txt', dtype=np.float32)
    return data

# 欧氏距离计算
def distEclud(x,y):
    return np.sqrt(np.sum((x-y)**2))  # 计算欧氏距离

# 为给定数据集构建一个包含K个随机质心的集合
def randCent(data,k):
    m,n = data.shape
    centroids = np.zeros((k,n))
    label = np.zeros(data.shape[0])
    for i in range(k):
        index = int(np.random.uniform(0,m))
        centroids[i,:] = data[index,:]
        label[index] = i+1
    return centroids,label

def KMeans(data,k):
    """
    :param data:数据
    :param k:类簇个数
    :return:
    """
    # 第1步 初始化centroids(聚类中心)
    centroids,label = randCent(data, k)
    clusterChange = True
    while clusterChange:
        for i in range(data.shape[0]):
            minDist = 100000.0
            minIndex = -1
            # 第2步 找出最近的质心
            for j in range(k):  # j代表第几个质心
                # 计算该样本到质心的欧式距离
                distance = distEclud(centroids[j, :], data[i, :])
                if distance < minDist:
                    # 把最小的距离和最小距离的质心的索引记录下来
                    minDist = distance
                    minIndex = j+1
            # 第3步:更新样本所属类别
            label[i] = minIndex
        min_attr = 0
        # 第4步：更新质心
        for j in range(1,k+1):
            pointsInCluster = data[np.nonzero(label == j)]
            # 获取簇类所有的点,np.nonzero得到数组array中非零元素的位置
            if distEclud(centroids[j-1, :],np.mean(pointsInCluster, axis=0))>min_attr:
                min_attr=distEclud(centroids[j-1, :],np.mean(pointsInCluster, axis=0))
            centroids[j-1, :] = np.mean(pointsInCluster, axis=0)  # 对矩阵的行求均值
        # 迭代退出条件
        if min_attr<0.1:
            clusterChange = False
    res = np.c_[data,label]
    return res,centroids

def showCluster(res, k, centroids):
    m, n = data.shape
    mark =['or','ob', 'og', 'ok', 'oy']
    if k > len(mark):
        print("k值太大了")
        return 1
    # 绘制所有的样本
    for i in range(m):
        markIndex = int(res[i,2]-1)
        plt.plot(data[i, 0], data[i, 1], mark[markIndex])
    mark = ['Dr', 'Db', 'Dg', 'Dk', 'Dy']
    # 绘制质心
    for i in range(k):
        plt.plot(centroids[i, 0], centroids[i, 1], mark[i])
    plt.show()

def calculate_error(centroids,res):
    print(centroids)
    for i in range(1,6):
        print(i,":",res[np.nonzero(res[:,2] == i)].shape[0])
    true_centroid = np.array([[1,4],[5.5,-4.5],[1,-1],[6,4.5],[9,0]])
    # 按照第一列的顺序排序
    true_centroid = true_centroid[np.argsort(true_centroid[:, 0])]

if __name__=="__main__":
    data = DataSet()
    res,centroids = KMeans(data,5)
    showCluster(res,5,centroids)
    calculate_error(centroids,res)

聚类结果可视化：

谱聚类（Normalized Spectral Clustering）实现

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans

def DataSet():
    data = np.loadtxt('data_2.txt', dtype=np.float32)
    return data

def calculate_w_ij(a,b,sigma=1):
    w_ab = np.exp(-np.sum((a-b)**2)/(2*sigma**2))
    #w_ab =  np.exp(-np.linalg.norm(a - b) ** 2 / (2 * sigma ** 2))
    return w_ab

def Construct_Matrix_W(data,k=5):
    # 计算临接矩阵
    rows = data.shape[0]
    W = np.zeros((rows,rows))
    for i in range(rows):
        for j in range(rows):
            if(i!=j):
                W[i][j] = calculate_w_ij(data[i,:],data[j,:])
        t = np.argsort(W[i,:])
        for ll in range(rows-k):
            W[i][t[ll]] = 0
    W = (W+W.T)/2
    return W

def Calculate_Matrix_L_sym(W):
    # 计算标准化的拉普拉斯矩阵
    degreeMatrix = np.sum(W, axis=1)
    # np.diag(v,k=0)
    # 以一维数组的形式返回方阵的对角线（或非对角线）元素，或将一维数组转换成方阵（非对角线元素为0）
    L = np.diag(degreeMatrix) - W
    # D^(-1/2) L D^(-1/2)
    sqrtDegreeMatrix = np.diag(1.0 / (degreeMatrix ** (0.5)))
    L_sym = np.dot(np.dot(sqrtDegreeMatrix, L), sqrtDegreeMatrix)
    return L_sym

def normalization(matrix):
    # 归一化
    sum = np.sqrt(np.sum(matrix**2,axis=1,keepdims=True))
    nor_matrix = matrix/sum
    return nor_matrix

def spKmeans(H,k=2):
    sp_kmeans = KMeans(n_clusters=k).fit(H)
    return sp_kmeans.labels_

def showCluster(res, k):
    m, n = data.shape
    mark =['or','ob', 'og', 'ok', 'oy']
    if k > len(mark):
        print("k值太大了")
        return 1
    # 绘制所有的样本
    for i in range(m):
        markIndex = int(res[i,2])
        plt.plot(data[i, 0], data[i, 1], mark[markIndex])
    plt.show()

if __name__=="__main__":
    data = DataSet()
    np.random.shuffle(data)

    W = Construct_Matrix_W(data)
    L_sym = Calculate_Matrix_L_sym(W)

    lam, H = np.linalg.eig(L_sym) # 矩阵特征值分解

    t = np.argsort(lam)
    H = np.c_[H[:,t[0]],H[:,t[1]]]
    H = normalization(H)

    labels = spKmeans(H,k=2)
    res = np.c_[data,labels]
    showCluster(res,k=2)

聚类结果可视化：

Reference

[1] On spectral clustering: analysis and an algorithm, NIPS, 2002.