K均值（K-means）聚类算法（Python3实现代码）

helloWorld

77031人浏览 · 2023-05-10 18:07:42

helloWorld · 2023-05-10 18:07:42 发布

本文的代码与数据集已上传至github：https://github.com/bitacron/MachineLearning
数据集地址：https://pan.baidu.com/s/1OYoD06nxvG5kcH35C25cXQ?pwd=p599

一、K-means算法基本思想

1、K-means简介

K-means标准算法是1957 年史都华·劳埃德（Stuart Lloyd）作为一种脉冲码调制的技术所提出，但直到1982年才被贝尔实验室公开出版在《IEEE Transactions on Information Theory》中，原始文献为《Least squares quantization in PCM》（DOI: 10.1109/TIT.1982.1056489）。
术语“k-均值”于1967年才被詹姆斯·麦昆(James MacQueen) 在文献《Some Methods for Classification and Analysis of Multiva》中提出，描述 K-means算法的完整理论并进行了详细的研究。作为最经典的划分聚类算法，K-means算法的实现并不复杂，具有较高的可伸缩性，同时K-means算法具有良好的可靠性和高效性，是一种广泛应用的聚类算法。

2、K-means算法流程

K-means算法接受一个参数K用以决定结果中簇的数目。算法开始时，要在数据集中随机选择K个数据对象用来当做k个簇的初始中心，而将剩下的各个数据对象就根据他们和每个聚类簇心的距离选择簇心最近的簇分配到其中。然后重新计算各个聚类簇中的所有数据对象的平均值，并将得到的结果作为新的簇心；逐步重复上述的过程直至目标函数收敛为止。

下面介绍该算法的具体步骤：

对于给定的一组数据，随机初始化K个聚类中心（簇中心）
计算每个数据到簇中心的距离（一般采用欧氏距离），并把该数据归为离它最近的簇。
根据得到的簇，重新计算簇中心。
对步骤2、步骤3进行迭代直至簇中心不再改变或者小于指定阈值。
K-means算法的流程图

3、K-means伪代码

输入：聚类数目k，数据集 $X=\left\{\right.x_1,x_2,…x_i,…,x_n\left\}\right.$ ，停止阈值ε，模糊因子m，最大迭代次数T；
输出：聚类中心 $C=[c_1,c _2,…,c_k]$ ；

begin
01.C ← Ø, U ← Ø, t ← 0 // 初始化聚类中心C、隶属度矩阵U和迭代次数t
02.while t < T do 
03.    for i ← 1 to n do
04.        for j ← 1 to k do
05.            calculate dij  //根据公式计算第i个数据和第j个聚类中心的距离
06.            if dij = Min{D(Xi ,Zj)} then // 根据距离归类
07.                Xi ∈ Cj  // 将数据xi归到第j类
08.        end for
09.    end for
10.    for j ← 1 to k do 
11.        calculate cj  //根据公式重新计算聚类中心点
12.    end for
13.    calculate J  //根据公式重新计算目标函数J
14.    t += 1
15.    if ||J(t) - J(t-1)|| ≤ ε then
16.        return C
17.    end if
18.end while
19.return C
end

4、K-means时间复杂度分析

在有n个样本的d维数据集X中，数据集被划分成k个簇，如果本次聚类经过了t次迭代，每次迭代需要计算一次目标函数，计算过程中需要分别求出d维数据中k个簇的聚类中心与n个数据的距离，算法在此过程的时间复杂度为O(n∙k∙d)。综上所述K-means算法的时间复杂度为O(n∙k∙t∙d)。但是由于在一般情况下类簇数目k、数据维度d以及算法的迭代次数t均可认为是常量，所以模糊C均值算法的时间复杂度可以简化为O(n)。由此可见，K-means算法的时间复杂度随着数据点的数量增加呈线性增长，同时也受迭代次数、聚类簇数目以及数据维度的影响。

5、K-means优缺点

K-means优点：

K-means算法收敛速度快，时间复杂度低，计算效率较高；
K-means算法算法可解释性好，原理简单易于理解；
K-means算法调参简单，需要调的参数只有聚类簇数K。

K-means缺点：

K-means算法的参数簇数目(K值)需要手动指定，K值的确定直接影响聚类的结果，通常情况下，K值需要经验和对数据集的理解指定，因此指定的数值未必理想，聚类的结果也就无从保证。
K-means算法对噪声和异常值敏感，聚类结果容易受到噪声点的影响；
K-means算法欧氏距离进行相似性度量，适用于簇是凸形或球形的数据集，在非凸形数据集中难以达到良好的聚类效果。；
K-means算法的初始中心点选取上采用的是随机的方法。K-means算法极为依赖初始中心点的选取：一旦错误地选取了初始中心点，对于后续的聚类过程影响极大，很可能得不到最理想的聚类结果，同时聚类迭代的次数也可能会增加。而随机选取的初始中心点具有很大的不确定性，也直接影响着聚类的效果；
K-means算法结果不一定是全局最优，只能保证局部最优。

二、K-means实现（Python3）

本文使用的数据集为UCI数据集，分别使用鸢尾花数据集Iris、葡萄酒数据集Wine、小麦种子数据集seeds进行测试，本文从UCI官网上将这三个数据集下载下来，并放入和python文件同一个文件夹内即可。同时由于程序需要，将数据集的列的位置做出了略微改动。数据集具体信息如下表：

数据集	样本数	属性维度	类别个数
Iris	150	4	3
Wine	178	3	3
Seeds	210	7	3

数据集在我主页资源里有，免积分下载，如果无法下载，可以私信我。

1、Python3代码实现

"""
KMeans.py

K-means聚类算法实现

"""

import time
from collections import namedtuple

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.optimize import linear_sum_assignment
from sklearn.decomposition import PCA
from sklearn.metrics import (
    accuracy_score,
    adjusted_rand_score,
    confusion_matrix,
    f1_score,
    normalized_mutual_info_score,
    rand_score,
)
from sklearn.preprocessing import LabelEncoder

# ==========================================================
# 数据集配置
# ==========================================================

# path        : 数据集路径
# clusters    : 聚类簇数
# label_col   : 标签列索引(None表示无标签)
# id_col      : ID列索引(None表示无ID列)
# header      : 表头行(0表示第一行为表头，None表示无表头)
# sep         : 分隔符
# skiprows    : 跳过前几行
# comment     : 注释符(None表示无注释符)

Dataset = namedtuple(
    "Dataset",
    [
        "path",
        "clusters",
        "label_col",
        "id_col",
        "header",
        "sep",
        "skip_rows",
        "comment"
    ]
)

DATASETS = {
    # ==========================
    # UCI数据集
    # ==========================

    # 最后一列是标签，英文逗号作为分割符，无特征名表头
    "iris": Dataset("dataset/uci/iris.data", 3, -1, None, None, ",", 0, None),
    # 最后一列是标签，空格作为分割符，无特征名表头
    "seeds": Dataset("dataset/uci/seeds_dataset.txt", 3, -1, None, None, r"\s+", 0, None),
    # 最后一列是标签，英文逗号作为分割符，无特征名表头
    "glass": Dataset("dataset/uci/glass.data", 6, -1, 0, None, ",", 0, None),
    # 第一列是标签，英文逗号作为分割符，无特征名表头
    "wine": Dataset("dataset/uci/wine.data", 3, 0, None, None, ",", 0, None),
    # 第一列是ID，第二列是标签，英文逗号作为分割符，无特征名表头
    "wdbc": Dataset("dataset/uci/wdbc.data", 2, 1, 0, None, ",", 0, None),
    # 第一列是ID，第二列是标签，英文逗号作为分割符，无特征名表头
    "ecoli": Dataset("dataset/uci/ecoli.data", 8, -1, 0, None, r"\s+", 0, None),

    # ==========================
    # 人工合成数据集
    # ==========================

    # 空格分隔，最后一列是标签
    "flame": Dataset("dataset/synthetic/Flame.txt", 2, -1, None, None, ",", 0, None),
    "jain": Dataset("dataset/synthetic/Jain.txt", 2, -1, None, None, ",", 0, None),
    "spiral": Dataset("dataset/synthetic/Spiral.txt", 3, -1, None, None, ",", 0, None),
    "panelB": Dataset("dataset/synthetic/panelB.txt", 3, None, None, None, ",", 0, None),
    "panelC": Dataset("dataset/synthetic/panelC.txt", 3, None, None, None, ",", 0, None),
    "r15": Dataset("dataset/synthetic/R15.txt", 15, -1, None, None, ",", 0, None),
    "d31": Dataset("dataset/synthetic/D31.txt", 31, -1, None, None, ",", 0, None),
    "aggregation": Dataset("dataset/synthetic/Aggregation.txt", 7, -1, None, None, ",", 0, None),
    "compound": Dataset("dataset/synthetic/Compound.txt", 6, -1, None, None, ",", 0, None),
    "pathbased": Dataset("dataset/synthetic/Pathbased.txt", 3, -1, None, None, ",", 0, None),
}


# ==========================================================
# 数据加载
# ==========================================================

def load_dataset(dataset_name):
    """
    加载数据集

    Parameters
    ----------
    dataset_name : str
        数据集名称

    Returns
    -------
    features : ndarray
        特征矩阵

    labels_true : list or None
        真实标签

    num_clusters : int
        聚类簇数
    """

    config = DATASETS[dataset_name]
    print("加载数据集:   "f"数据集={dataset_name} ")
    df = pd.read_csv(config.path, sep=config.sep, header=config.header, skiprows=config.skip_rows,comment=config.comment)

    # 提取真实标签
    if config.label_col is not None:
        labels_true = df.iloc[:, config.label_col].to_numpy()
    else:
        labels_true = None

    # 构造需要删除的列(包括标签列、id列)
    drop_cols = []

    if config.id_col is not None:
        drop_cols.append(df.columns[config.id_col])

    if config.label_col is not None:
        drop_cols.append(df.columns[config.label_col])

    # 得到特征矩阵
    if drop_cols:
        # 一次性删除，防止错乱
        features = df.drop(columns=drop_cols).values
    else:
        features = df.values

    # 输出样本数n_samples和属性数n_attributes(维度dimensions)
    print(f"样本数量={features.shape[0]}  " f"属性数量(维度)={features.shape[1]}")
    return features, labels_true, config.clusters


# ==========================================================
# K-means核心算法
# ==========================================================

def initialize_centroids(features, num_clusters, random_state=42):
    """
    初始化聚类中心（从数据集中随机选择k个点作为初始质心）

    Parameters
    ----------
    features : ndarray, shape=(n_samples, n_features)
        数据矩阵 X

    num_clusters : int
        聚类簇数 k

    random_state : int, default=42
        随机种子

    Returns
    -------
    centroids : ndarray, shape=(num_clusters, n_features)
        初始聚类中心
    """
    rng = np.random.default_rng(random_state)
    indices = rng.choice(features.shape[0], num_clusters, replace=False)
    return features[indices].copy()


def assign_clusters(features, centroids):
    """
    将数据点分配给最近的聚类中心

    Parameters
    ----------
    features : ndarray, shape=(n_samples, n_features)
        数据矩阵 X

    centroids : ndarray, shape=(num_clusters, n_features)
        聚类中心矩阵 V

    Returns
    -------
    cluster_labels : ndarray, shape=(n_samples,)
        每个样本的聚类标签
    """
    # 计算每个样本到每个聚类中心的欧氏距离
    distances = np.linalg.norm(features[:, None, :] - centroids[None, :, :], axis=2)
    # 返回距离最近的聚类中心索引
    return np.argmin(distances, axis=1)


def calculate_centroids(features, cluster_labels, num_clusters):
    """
    计算每个簇的新质心（簇内数据点的均值）

    Parameters
    ----------
    features : ndarray, shape=(n_samples, n_features)
        数据矩阵 X

    cluster_labels : ndarray, shape=(n_samples,)
        聚类标签

    num_clusters : int
        聚类簇数 k

    Returns
    -------
    centroids : ndarray, shape=(num_clusters, n_features)
        更新后的聚类中心
    """
    centroids = np.zeros((num_clusters, features.shape[1]))
    for i in range(num_clusters):
        # 获取属于第i簇的所有样本
        cluster_points = features[cluster_labels == i]
        if len(cluster_points) > 0:
            centroids[i] = cluster_points.mean(axis=0)
        else:
            # 如果某个簇没有样本，重新随机初始化（避免空簇）
            rng = np.random.default_rng()
            centroids[i] = features[rng.choice(features.shape[0])]
    return centroids


def kmeans_clustering(features, num_clusters, tol=1e-5, max_iter=100, random_state=42):
    """
    K-means聚类主函数

    Parameters
    ----------
    features : ndarray
        特征矩阵(符号：X)

    num_clusters : int
        聚类簇数(符号：k)

    tol : float, default=1e-5
        终止阈值容差tolerance(符号：epsilon)

    max_iter : int, default=100
        最大迭代数(符号：T)

    random_state : int, default=42
        随机种子

    Returns
    -------
    cluster_labels : ndarray
        聚类标签

    cluster_centers : ndarray
        聚类中心
    """
    start = time.time()
    cluster_labels = np.zeros((num_clusters, features.shape[1]))  # 提前初始化，避免警告
    # Step1: 初始化聚类中心
    centroids = initialize_centroids(features, num_clusters, random_state)
    iteration = 0

    while iteration < max_iter:
        # Step2: 分配样本到最近的聚类中心
        cluster_labels = assign_clusters(features, centroids)

        # Step3: 保存旧聚类中心
        old_centroids = centroids.copy()

        # Step4: 更新聚类中心
        centroids = calculate_centroids(features, cluster_labels, num_clusters)

        # Step5: 判断收敛
        if np.linalg.norm(centroids - old_centroids) < tol:
            break

        iteration += 1

    elapsed_time = time.time() - start
    print(f"总计迭代次数: {iteration};  耗时: {elapsed_time:.4f} s")

    return cluster_labels, centroids


# ==========================================================
# 标签对齐
# ==========================================================

def align_labels(labels_true, labels_pred):
    """
    使用匈牙利算法进行标签匹配
    """
    cm = confusion_matrix(labels_true, labels_pred)
    row_ind, col_ind = linear_sum_assignment(-cm)
    mapping = {}
    for true_label, pred_label in zip(row_ind, col_ind):
        mapping[pred_label] = true_label

    labels_pred_aligned = np.array([
        mapping[label]
        for label in labels_pred
    ])

    return labels_pred_aligned


# ==========================================================
# 聚类评价指标
# ==========================================================

def clustering_indicators(labels_true, labels_pred):
    """
    计算聚类评价指标
    """

    if isinstance(labels_true[0], str):
        # 真实标签转化为数字标签
        labels_true = LabelEncoder().fit_transform(labels_true)

    labels_pred_aligned = align_labels(labels_true, labels_pred)

    f1 = f1_score(labels_true, labels_pred_aligned, average="macro")
    accuracy = accuracy_score(labels_true, labels_pred_aligned)
    nmi = normalized_mutual_info_score(labels_true, labels_pred)
    ri = rand_score(labels_true, labels_pred)
    ari = adjusted_rand_score(labels_true, labels_pred)

    return f1, accuracy, nmi, ri, ari


# ==========================================================
# 聚类结果可视化
# ==========================================================

def draw_cluster(features, cluster_centers, labels_pred):
    """
    绘制聚类结果散点图
    """

    features = np.asarray(features)
    cluster_centers = np.asarray(cluster_centers)

    if features.shape[1] > 2:
        pca = PCA(n_components=2)
        features_2d = pca.fit_transform(features)
        cluster_centers_2d = pca.transform(cluster_centers)
    else:
        features_2d = features
        cluster_centers_2d = cluster_centers

    plt.figure(figsize=(8, 6))
    # 做散点图
    plt.scatter(features_2d[:, 0], features_2d[:, 1], marker='o', c='black', s=7)  # 原图
    plt.show()
    plt.scatter(features_2d[:, 0], features_2d[:, 1], c=labels_pred, cmap="nipy_spectral", s=7, marker="o")

    plt.scatter(cluster_centers_2d[:, 0], cluster_centers_2d[:, 1], marker="x", color="m", s=30)
    # 设置x和y坐标轴刻度的标签字体和字号
    # plt.xticks(fontproperties='Times New Roman', fontsize=10.5)
    # plt.yticks(fontproperties='Times New Roman', fontsize=10.5)
    # plt.xlabel("x - label", fontdict={'family': 'Times New Roman', 'size': 10.5}, loc="right")
    # plt.ylabel("y - label", fontdict={'family': 'Times New Roman', 'size': 10.5}, loc="top")
    plt.title("K-means Clustering Result")

    plt.show()


# ==========================================================
# 主程序
# ==========================================================

if __name__ == "__main__":
    # 初始化本次数据集名称
    _dataset_name = "seeds"
    # 加载数据集
    _features, _labels_true, _num_clusters = load_dataset(_dataset_name)

    _tol = 1e-5
    _max_iter = 100
    _random_state = 42
    print("开始K-means:   " f"k={_num_clusters}  " f"n={len(_features)}   " f"epsilon={_tol}  " f"T={_max_iter}  " f"random_state={_random_state}")

    # 进行K-means算法
    _labels_pred, _cluster_centers = kmeans_clustering(_features, _num_clusters, _tol, _max_iter, _random_state)

    # 开始计算聚类指标
    if _labels_true is not None:
        # 有标签：计算聚类指标
        F1, ACC, NMI, RI, ARI = clustering_indicators(_labels_true, _labels_pred)
        print("聚类指标: " f"F1={F1:.6f}  " f"ACC={ACC:.6f}  " f"NMI={NMI:.6f}  " f"RI={RI:.6f}  " f"ARI={ARI:.6f}")
    else:
        # 数据集无标签：无法计算聚类指标
        print("Dataset without labels，unable to calculate clustering index")

    # 输出散点图
    draw_cluster(_features, _cluster_centers, _labels_pred)

2、聚类结果分析

本文选择了F值（F-measure，FM）、准确率（Accuracy，ACC）、标准互信息（Normalized Mutual Information，NMI）和兰德指数（Rand Index，RI）作为评估指标，其值域为[0,1]，取值越大说明聚类结果越符合预期。

F值结合了精度（Precision）与召回率（Recall）两种指标，它的值为精度与召回率的调和平均，其计算公式见公式：

$Precision=\frac{TP}{TP+FP}$

$Recall=\frac{TP}{TP+FN}$

$F-measure=\frac{2Recall \times Precision}{Recall+Precision}$

ACC是被正确分类的样本数与数据集总样本数的比值，计算公式如下：

$ACC=\frac{TP+TN}{TP+TN+FP+FN}$

其中，TP（True Positive）表示将正类预测为正类数的样本个数，TN （True Negative）表示将负类预测为负类数的样本个数，FP（False Positive）表示将负类预测为正类数误报的样本个数，FN（False Negative）表示将正类预测为负类数的样本个数。

NMI用于量化聚类结果和已知类别标签的匹配程度，相比于ACC，NMI的值不会受到族类标签排列的影响。计算公式如下：

$NMI=\frac{I\left(U,V\right)}{\sqrt{H\left(U\right)H\left(V\right)}}$

其中H(U)代表正确分类的熵，H(V)分别代表通过算法得到的结果的熵。

其具体实现代吗如下：
由于数据集中给定的正确标签可能为文本类型而不是数字标签，所以在计算前先判断数据集的标签是否为数字类型，如果不是，则转化为数字类型

聚类结果存在标签置换（Label Permutation）问题。聚类算法生成的簇标签是无序的，计算 ACC 和 F 值时若直接与真实标签比较，标签未对齐会导致评价结果失真。而 NMI、ARI、RI 不会又此问题，所以一直稳定。所以需要采用匈牙利算法建立最优标签映射关系，从而进行标签匹配（Label Matching）

# 调用confusion_matrix函数，使用匈牙利算法，进行标签匹配（Label Matching）
def align_labels(labels_true, labels_pred):
    cm = confusion_matrix(labels_true, labels_pred)

    row_ind, col_ind = linear_sum_assignment(-cm)

    mapping = {}
    for true_label, pred_label in zip(row_ind, col_ind):
        mapping[pred_label] = true_label

    labels_pred_aligned = np.array(
        [mapping[label] for label in labels_pred]
    )

    return labels_pred_aligned


# 计算聚类指标
def clustering_indicators(labels_true, labels_pred):
    if type(labels_true[0]) != int:
        labels_true = LabelEncoder().fit_transform(df[columns[len(columns)-1]])  # 如果标签为文本类型，把文本标签转换为数字标签
    # 标签对齐
    labels_pred_aligned = align_labels(labels_true, labels_pred)
    f_measure = f1_score(labels_true, labels_pred_aligned, average='macro')  # F值
    accuracy = accuracy_score(labels_true, labels_pred_aligned)  # ACC
    normalized_mutual_information = normalized_mutual_info_score(labels_true, labels_pred)  # NMI
    rand_index = rand_score(labels_true, labels_pred)  # RI
    adjusted_rand_index = adjusted_rand_score(labels_true, labels_pred)  # ARI
    return f_measure, accuracy, normalized_mutual_information, rand_index, adjusted_rand_index