K-Means 聚类算法

聚类（clustering）是一种寻找数据之间内在结构的技术。聚类把全体数据实例组织成一些相似组，而这些相似组被称作簇。处于相同簇中的数据实例彼此相同，处于不同簇中的实例彼此不同。K-Means、AgglomerativeClustering、DBSCAN、MeanShift、SpectralClustering 等是常用的聚类方法。 K-Means 算法又名 K 均值算法，K-Means 算法中的 K 表示的是聚类为 K 个簇，Means 代表取每一个聚类中数据值的均值作为该簇的中心，或者称为质心，即用每一个类的质心对该簇进行描述。

K-Means 算法原理

简介

K-Means 是一种广泛使用的无监督学习(没标准答案下去学习)算法，用于将数据点划分为 k 个簇（Cluster）。其目标是使同一簇内的数据点尽可能相似，而不同簇之间的数据点尽可能不同。

算法原理

K-Means 算法是一个迭代过程，主要包含以下步骤：

初始化 (Initialization): 随机选择 k 个数据点作为初始质心（Centroids）。
分配 (Assignment): 计算每个数据点到这 k 个质心的距离，并将每个数据点分配给最近的质心所代表的簇。
更新 (Update): 重新计算每个簇的质心。新的质心是该簇内所有数据点的平均值（均值）。
重复 (Repeat): 重复步骤 2 和 3，直到满足停止条件（例如：质心不再发生变化、达到最大迭代次数或误差平方和收敛）

数学公式

目标函数

K-Means 的目标是最小化簇内误差平方和(Within-Cluster Sum of Squares, WCSS)，也称为惯性（Inertia）。

$$ J =SEE= \sum_{i=1}^k\sum_{x\in C_i}||x-\mu_i||^2 $$

其中：

J 是目标函数值:整个数据集的误差平方和（Sum of Squared Error，SSE）。
k 是簇的数量。
$C_i$ 是第 i 个簇。
X 是簇 $C_i$中的数据点。
$\mu_i$是簇 $C_i$的质心。
$||x-\mu_i||^2$ 是数据点 x 到质心 $\mu_i$ 的欧几里得距离的平方

质心更新公式

在第 t+1 次迭代中，簇$C_i$ 的新质心$\mu_i^{(t+1)}$ 计算如下：

$$ \mu_i^{(t+1)} = \frac{1}{|C_i^{(t)}|} \sum_{x\in C_i^{(t)}} x $$

其中 $|C_i^{(t)}|$ 是第 t 次迭代中簇 $C_i$内数据点的数量。

简单样例说明

假设我们有以下一维数据点：[2, 4, 10, 12, 3, 20, 30, 11, 25]，我们想将其分为 k=2 个簇。

初始化: 随机选 2 和 30 为质心。
- μ1=2, μ2=30
分配:
- 靠近 2 的点: {2, 3, 4, 10, 11, 12} -> 簇 1
- 靠近 30 的点: {20, 25, 30} -> 簇 2
更新:
- 新 μ1=(2+3+4+10+11+12)/6=42/6=7
- 新 μ2=(20+25+30)/3=75/3=25
迭代: 使用新质心 7 和 25 重新分配，直到质心稳定

Python 代码样例

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101


import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from matplotlib.animation import FuncAnimation

# ==========================================
# 1. 数据准备
# ==========================================
print("Generating data...")
X, y_true = make_blobs(n_samples=300, centers=4, cluster_std=0.60, random_state=0)

# ==========================================
# 2. K-Means 算法手动实现 (为了记录每一步)
# ==========================================
class KMeansManual:
    def __init__(self, X, k=4):
        self.X = X
        self.k = k
        # 随机选择k个点作为初始质心
        indices = np.random.choice(X.shape[0], k, replace=False)
        self.centroids = X[indices]
        self.labels = None
        self.history = [] # 用于存储每一步的状态: (centroids, labels, title)

    def run(self, max_iters=10):
        # 记录初始状态
        self.history.append((self.centroids.copy(), None, "Initialization"))

        for i in range(max_iters):
            # --- 步骤 1: 分配 (Assignment) ---
            # 计算距离矩阵: (k, n_samples)
            # 利用广播机制计算每个点到每个质心的欧氏距离
            distances = np.linalg.norm(self.X[:, np.newaxis] - self.centroids, axis=2)

            # 获取最近质心的索引
            self.labels = np.argmin(distances, axis=1)

            # 记录分配后的状态
            self.history.append((self.centroids.copy(), self.labels.copy(), f"Iteration {i+1}: Assignment"))

            # --- 步骤 2: 更新 (Update) ---
            new_centroids = np.array([
                self.X[self.labels == j].mean(axis=0) if np.sum(self.labels == j) > 0 else self.centroids[j]
                for j in range(self.k)
            ])

            # 记录更新质心后的状态
            self.history.append((new_centroids.copy(), self.labels.copy(), f"Iteration {i+1}: Update Centroids"))

            # --- 检查收敛 ---
            if np.allclose(self.centroids, new_centroids):
                print(f"Converged at iteration {i+1}")
                break

            self.centroids = new_centroids

# 运行算法
k = 4
kmeans = KMeansManual(X, k=k)
kmeans.run()

# ==========================================
# 3. 动画可视化
# ==========================================
print("Starting animation...")
fig, ax = plt.subplots(figsize=(10, 7))

def update(frame_idx):
    ax.clear()
    centroids, labels, title = kmeans.history[frame_idx]

    # 绘制数据点
    if labels is None:
        # 初始状态，未分类，显示为灰色
        ax.scatter(X[:, 0], X[:, 1], c='gray', s=50, alpha=0.5)
    else:
        # 已分类，按类别着色
        ax.scatter(X[:, 0], X[:, 1], c=labels, s=50, cmap='viridis', alpha=0.6)

    # 绘制质心 (红色大叉)
    ax.scatter(centroids[:, 0], centroids[:, 1], c='red', s=200, marker='X', edgecolors='white', linewidths=2, label='Centroids')

    # 绘制质心移动轨迹 (可选，如果不是第一帧)
    if frame_idx > 0:
        prev_centroids, _, _ = kmeans.history[frame_idx-1]
        # 简单的绘制从上一个位置到当前位置的线
        for i in range(k):
            ax.plot([prev_centroids[i, 0], centroids[i, 0]], 
                    [prev_centroids[i, 1], centroids[i, 1]], 'k--', alpha=0.5)

    ax.set_title(title, fontsize=14)
    ax.set_xlabel("Feature 1")
    ax.set_ylabel("Feature 2")
    ax.legend()
    ax.grid(True, alpha=0.3)

# 创建动画
# interval=1000 表示每帧停留 1000ms (1秒)
ani = FuncAnimation(fig, update, frames=len(kmeans.history), interval=2000, repeat=True)

plt.show()

K-Means算法的问题

K如何确定

手肘法

算法所需要预设的参数至少有2个：簇个数K和初始簇质心

聚类的目标是使得每个样本点到距离其最近的聚类中心的总误差平方和(SSE)尽可能小。

理论上随着K的增加，SSE会单调递减，因为类数的增加意味着总有一部分样本点会因为归属到新的类簇而节约下一段距离，直到K=N时，情况将演变为每个样本自成一类，此时SSE值就会降为0,这也就没意义了。

根据学者们的长期实践经验，K值最大不应超过样本量的开平方根，即$K_{max}<= \sqrt N$ 。而确定了范围后，最优K值又应该怎么判断？一种简单的思路是：试图找到某一个K值，要求当K大于该值时，SSE的下降变化幅度（或速度）明显变小。换句话说，当K超过某一个数后，每个类簇的聚合程度不再获得显著提升，此时我们就可以认为已找到最佳K的取值。这也是手肘法通过画出不同K值与SSE值的折线图，若SSE值下降过程中存在“肘点”（下降速度骤减的拐点处），该点所对应的K值即合适的聚类数。不过遗憾的是，若SSE的下降是均匀的，传统的肘部图法也就失灵了。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48


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

# ==========================================
# 1. 数据准备
# ==========================================
print("Generating data...")
# 生成稍微复杂一点的数据，以便观察肘部
X, y_true = make_blobs(n_samples=500, centers=5, cluster_std=0.8, random_state=42)

# ==========================================
# 2. 计算不同 K 值的 WCSS (Inertia)
# ==========================================
wcss = [] # Within-Cluster Sum of Squares (簇内误差平方和)
k_range = range(1, 11) # 测试 k 从 1 到 10

print("Calculating WCSS for different k values...")
for k in k_range:
    kmeans = KMeans(n_clusters=k, random_state=42)
    kmeans.fit(X)
    wcss.append(kmeans.inertia_) # inertia_ 属性就是 WCSS
    print(f"k={k}, WCSS={kmeans.inertia_:.2f}")

# ==========================================
# 3. 绘制肘部图 (Elbow Plot)
# ==========================================
plt.figure(figsize=(10, 6))
plt.plot(k_range, wcss, marker='o', linestyle='-', color='b')

# 标注轴和标题
plt.title('Elbow Method For Optimal k (肘部法则)', fontsize=16)
plt.xlabel('Number of clusters (k)', fontsize=12)
plt.ylabel('WCSS (Inertia)', fontsize=12)
plt.xticks(k_range)
plt.grid(True)

# 标记可能的肘部点 (在这个例子中，我们知道真实中心是5，所以肘部应该在5附近)
# 这里为了演示，我们简单标注一下 k=5 的位置
optimal_k = 5
plt.annotate(f'Possible Elbow Point (k={optimal_k})', 
             xy=(optimal_k, wcss[optimal_k-1]), 
             xytext=(optimal_k+1, wcss[optimal_k-1]+1000),
             arrowprops=dict(facecolor='red', shrink=0.05),
             fontsize=12, color='red')

plt.show()

轮廓系数（Silhouette Coefficient）法

是一种用于评估聚类效果并帮助确定最优聚类数 K 的方法。它结合了凝聚度（cohesion）和分离度（separation）两个指标

轮廓系数的定义

对于数据集中的每一个样本点 $x_i$，其轮廓系数 $s(i)$ 定义为：

$$ s(i) = \frac{b(i) - a(i)}{\max\{a(i), b(i)\}} $$

其中： a(i)：样本 xi 到同簇内其他点的平均距离（衡量凝聚度）。 b(i)：样本 xi 到最近的其他簇中所有点的平均距离（衡量分离度）。

轮廓系数的取值范围为 [−1,1]： - 接近 1：表示样本聚类合理，簇内紧密、簇间分离良好。 - 接近 0：表示样本在两个簇边界上。 - 接近 -1：表示样本可能被分配到了错误的簇。

使用轮廓系数确定 K 值的步骤

对不同的 K 值（如 K=2,3,…,Kmax）运行 K-means 聚类。
对每次聚类结果，计算每个样本的轮廓系数，再求所有样本的平均轮廓系数。
选择使平均轮廓系数最大的 K 值作为最优聚类数.

注意：轮廓系数法假设簇是凸形且大小相近的，对于非球形或密度差异大的簇可能不适用。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57


import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.metrics import silhouette_score

# ==========================================
# 1. 数据准备
# ==========================================
print("Generating data...")
# 生成与之前类似的模拟数据，预设中心为 5
X, y_true = make_blobs(n_samples=500, centers=5, cluster_std=0.8, random_state=42)

# ==========================================
# 2. 计算不同 K 值的轮廓系数 (Silhouette Score)
# ==========================================
silhouette_scores = []
# 注意：轮廓系数至少需要 2 个簇才能计算，所以从 k=2 开始
k_range = range(2, 11) 

print("Calculating Silhouette Scores for different k values...")
for k in k_range:
    kmeans = KMeans(n_clusters=k, random_state=42)
    labels = kmeans.fit_predict(X)

    # 计算所有样本的平均轮廓系数
    # 轮廓系数范围是 [-1, 1]，越接近 1 表示聚类效果越好
    score = silhouette_score(X, labels)
    silhouette_scores.append(score)
    print(f"k={k}, Silhouette Score={score:.4f}")

# ==========================================
# 3. 绘制轮廓系数图
# ==========================================
plt.figure(figsize=(10, 6))
plt.plot(k_range, silhouette_scores, marker='o', linestyle='-', color='g', linewidth=2)

# 标注轴和标题
plt.title('Silhouette Method For Optimal k', fontsize=16)
plt.xlabel('Number of clusters (k)', fontsize=12)
plt.ylabel('Silhouette Score (Avg)', fontsize=12)
plt.xticks(k_range)
plt.grid(True, alpha=0.3)

# 标记最高分 (最佳 K 值)
best_k_idx = np.argmax(silhouette_scores)
best_k = k_range[best_k_idx]
best_score = silhouette_scores[best_k_idx]

plt.annotate(f'Best k={best_k}\nScore={best_score:.3f}', 
             xy=(best_k, best_score), 
             xytext=(best_k, best_score - 0.1), # 文本位置稍微往下一点
             arrowprops=dict(facecolor='red', shrink=0.05),
             fontsize=12, color='red', ha='center')

print(f"The optimal k according to Silhouette Score is {best_k}")
plt.show()

间隔统计法（gap statistic）

Gap统计量基于以下假设：如果聚类是有意义的，那么数据集中的样本点应该比随机数据更紧密地聚集在一起。因此，Gap统计量计算了实际数据集的WCSS与随机数据集WCSS的期望值之间的差异。

对于每一个K值，首先运行K-means算法，得到一个群内平方和。
然后，生成一组随机数据，并用相同的K值运行K-means算法。
比较真实数据的群内平方和和随机数据的结果，并计算他们之间的差距（称之为间隔值）。
对于多个K值，重复以上步骤，并选择拥有最大间隔值的K

初始的簇质心的选取

常用分析软件中的功能模块或函数包，基本上都已经代替使用者们自动预设了随机初始点，只需填入目标K值，就可以跑动算法。但实际上，K-Means对初始聚类中心的位置十分敏感，每次迭代，初始点的不同往往会导致不同的聚类结果。此外过于临近的初始中心点，有时还会导致模型的收敛时间变长（即Step4中迭代时间变长）。一种简单粗暴的解决方式是，选择不同的初始聚类中心，多次运行算法，挑出聚类效果更佳（SSE更小）、解释性更强的一组结果。

当然了，我们或许更想知道算法上的改进手段。一种常见的优化方法是采用最大距离法，如：首先选取数据集中距离最大的两个点作为初始聚类中心，将剩余数据对象依据到聚类中心点距离的远近分配到相应的簇中，并更新聚类中心，然后继续寻找与聚类中心距离最远的点作为下一个中心点……

与此类似地还有K-Means++ 算法，它是传统K-Means的改良版，同样是基于最大距离，这里结合加权概率的思想优化了对K个初始中心的选取，使得在选取第n+1(n+1<k)个聚类中心时，距离当前n个聚类中心越远的点会有更高的概率被选为第n+1个聚类中心。还有学者从点集密度的角度改进，又或者将优化搜索算法（如模拟退火、生物遗传算法等）运用到了聚类中心的选取中……

相似性与距离度量问题

特征量化后，不同个体的相似性反映在了向量之间的空间距离大小，常见的度量方法包括欧几里得距离、曼哈顿距离等等，有时我们还会用到余弦相似度等（如计算文档相似性）。而通常情况下，欧氏距离计算就可以满足我们对实现K-Means的需要。根据距离的度量方式容易发现，K-Means所划分出的类别是类球形的，换句话说，只有类球型分布的连续型样本数据，才能得到较好的聚类效果，而如果非数值型、样本类别极不平衡、非球形的分类，则聚类效果会受限。对于非理想情形的数据，有时我们就需要做一些灵活变通了。如，若数据为离散型，均值没有定义，我们可以采用K-众数(每个簇的质心不是均值，而是每个维度上出现频率最高的类别)的方法。如果样本类别极不平衡或者是非球型，则可能考虑更换聚类方法，如使用基于密度的聚类（经典的DBSCAN算法）、层次聚类法等等。

聚类时间问题

上面的流程中提到，当聚类中心不再改变时（数学上即要求SSE函数收敛），我们认为聚类过程结束，但是这并不是唯一的结束信号。为了节省计算时间，有时我们也会通过设置迭代次数、设置簇内平方和或SSE下降阈值，又或者替换为“直到仅有1%的点改变簇”这样的弱条件，来控制算法的进程。出于对问题复杂度和计算量的合理预判，若聚类中心的更新超过了迭代次数上限，或者代价函数SSE已经小于所设定的阈值，我们都有理由提前终止。

此外，为了提高收敛速度，还可以考虑采用二分K-Means法，将所有点作为一个簇，将该簇一分为二，然后选择能最大程度降低聚类代价函数的簇划分为两个簇，以此进行下去，直到簇的数目等于给定的个数K为止。值得一提的是，该法更突出的优点在于能够很好地解决K-Means收敛到局部最优的问题，帮助我们找到全局最优解

标准化问题

考虑到我们所研究的对象通常包含多列数据，这些数据代表不同方面的属性值，在单位和数量级上可能存在较大的差别，因此为了避免这些差异可能引发的计算精度下降等问题，对于连续属性，可以先对数据进行规范化处理（如零均值规范化、最大最小值规范化等），再进行距离的计算