K-Means Clustering —— 聚类算法
前言
在之前的文章中,我们讨论的算法都有一个共同点——它们都是有监督学习:数据既有特征 X,也有标签 y,模型的任务是学习从 X 到 y 的映射。
但现实世界中,大量的数据是没有标签的。比如电商平台有海量的用户行为数据,但没有事先标注好“这属于哪类用户”;新闻网站有无数篇文章,但没有事先分好“这是体育类还是政治类”。聚类(Clustering) 就是解决这类问题的最经典方法——它试图在无标签的数据中发现天然的分组结构。
K-Means 是聚类算法中最著名、最常用的算法之一。它于1955年由Stuart Lloyd提出,因其简单、直观、高效而经久不衰。尽管已经有半个多世纪的历史,K-Means至今仍然是数据科学中最常用的工具之一。
这篇文章,我们将从Lloyd迭代算法出发,一步步理解K-Means的数学原理,深入剖析K-Means与EM算法的深层联系,分析K-Means++的初始化优化,最后讨论肘部法则和轮廓系数的使用与局限。
一、K-Means的核心思想与Lloyd算法
1.1 K-Means在做什么?
假设我们有一堆数据点,想将它们分成 K 个组(簇)。K-Means的直觉极其简单:
每个簇应该有一个“中心”,每个数据点应该属于离它最近的那个中心。
用数学语言来表达,K-Means的目标是最小化所有数据点到其所属簇中心的距离平方和:
J=i=1∑N∥xi−μci∥2其中:
- N 是数据点的总数
- xi 是第 i 个数据点
- μci 是第 i 个数据点所属簇的中心
- ci∈{1,2,...,K} 是第 i 个数据点的簇分配
这个目标函数 J 也被称为畸变(Distortion) 或惯量(Inertia) ——它衡量了簇的“紧密程度”:值越小,簇内的点越聚集。
这里的优化变量有两个:簇中心 μ1,...,μK 和每个点的簇分配 c1,...,cN。如果同时优化这两个变量,问题是一个NP难问题——没有多项式时间的精确解法。
1.2 Lloyd算法:交替优化的经典策略
既然不能同时优化,那就交替优化。Lloyd算法(即标准K-Means算法)采用了一个极其简单的迭代策略:
步骤1:初始化。随机选择 K 个数据点作为初始簇中心。
步骤2:重复以下两步直到收敛:
(a) 分配(Assignment) —— E步的“硬”版本:将每个数据点分配给离它最近的簇中心。
ci=argjmin∥xi−μj∥2(b) 更新(Update) —— M步的“硬”版本:重新计算每个簇的中心(即簇内所有点的均值)。
μj=∣Cj∣1i∈Cj∑xi其中 Cj 是分配给第 j 个簇的所有数据点的集合。
收敛判断:当簇中心不再发生变化(或变化小于某个阈值)时,算法停止。
1import numpy as np2import matplotlib.pyplot as plt3from sklearn.datasets import make_blobs4
5class KMeans:6 """从零实现K-Means聚类算法(Lloyd算法)"""7
8 def __init__(self, n_clusters=3, max_iter=300, random_state=None):9 self.n_clusters = n_clusters10 self.max_iter = max_iter11 self.random_state = random_state12 self.centroids = None13 self.labels_ = None14 self.inertia_ = None15
16 def fit(self, X):17 np.random.seed(self.random_state)18 n_samples, n_features = X.shape19
20 # 1. 初始化:随机选择K个样本作为初始中心21 indices = np.random.choice(n_samples, self.n_clusters, replace=False)22 self.centroids = X[indices].copy()23
24 for _ in range(self.max_iter):25 # 2a. 分配(E步):计算每个点到所有中心的距离,分配到最近的中心26 distances = np.zeros((n_samples, self.n_clusters))27 for k in range(self.n_clusters):28 distances[:, k] = np.linalg.norm(X - self.centroids[k], axis=1)29 labels = np.argmin(distances, axis=1)30
31 # 2b. 更新(M步):重新计算每个簇的中心32 new_centroids = np.zeros_like(self.centroids)33 for k in range(self.n_clusters):34 if np.sum(labels == k) > 0:35 new_centroids[k] = X[labels == k].mean(axis=0)36 else:37 # 如果某个簇没有分配到点,保持原中心不变(或重新初始化)38 new_centroids[k] = self.centroids[k]39
40 # 检查是否收敛41 if np.allclose(self.centroids, new_centroids, rtol=1e-4):42 break43
44 self.centroids = new_centroids45
46 # 计算最终的簇分配和惯量47 distances = np.zeros((n_samples, self.n_clusters))48 for k in range(self.n_clusters):49 distances[:, k] = np.linalg.norm(X - self.centroids[k], axis=1)50 self.labels_ = np.argmin(distances, axis=1)51 self.inertia_ = np.sum(np.min(distances, axis=1) ** 2)52
53 return self54
55 def predict(self, X):56 distances = np.zeros((X.shape[0], self.n_clusters))57 for k in range(self.n_clusters):58 distances[:, k] = np.linalg.norm(X - self.centroids[k], axis=1)59 return np.argmin(distances, axis=1)60
61# 生成示例数据62X, y_true = make_blobs(n_samples=300, centers=3, cluster_std=0.8, random_state=42)63
64# 训练K-Means65kmeans = KMeans(n_clusters=3, random_state=42)66kmeans.fit(X)67
68print(f"收敛后的惯量: {kmeans.inertia_:.2f}")69print(f"簇标签: {np.unique(kmeans.labels_, return_counts=True)}")1.3 Lloyd算法的收敛性
Lloyd算法有一个重要的理论性质:它保证会收敛。
为什么?因为每一步都在单调地减小目标函数 J:
- 分配步骤:将每个点重新分配给最近的中心,J 不增
- 更新步骤:用均值更新中心,J 不增(因为均值是给定分配下使平方和最小的点)
由于 J≥0 有下界,且每一步都在减小,算法最终必然收敛。
然而,这个收敛是局部最优的。Lloyd算法只能保证收敛到某个局部极小值,而不是全局最优。不同的初始中心会导致不同的聚类结果。
二、EM视角:K-Means是EM算法的特例
2.1 EM算法回顾
EM(Expectation-Maximization)算法是处理含有隐变量的概率模型参数估计的通用框架。
EM算法的核心思想是:
当数据中存在隐变量(未被观测到的变量)时,我们无法直接最大化似然函数。EM算法通过交替进行“猜测隐变量的分布”(E步)和“基于猜测更新参数”(M步),逐步逼近最大似然估计。
EM算法的标准流程:
E步(Expectation) :基于当前参数 θ,计算隐变量的后验分布 P(Z∣X,θ)。
M步(Maximization) :基于E步得到的隐变量分布,用最大似然估计更新参数 θ。
2.2 K-Means作为EM的特例
那么,K-Means和EM是什么关系?
K-Means可以看作是EM算法的一个特例。
让我们把K-Means套进EM的框架中:
| EM框架 | K-Means的对应物 |
|---|---|
| 隐变量 Z | 每个数据点的簇分配 ci(属于哪个簇) |
| 参数 θ | 簇中心 μ1,...,μK |
| E步 | 将每个点硬分配到最近的簇中心 |
| M步 | 用簇内点的均值更新簇中心 |
关键区别在于 “硬”vs“软” :
- K-Means的E步:每个点以概率1分配给某一个簇(最近的),以概率0分配给其他簇——这是硬分配(Hard Assignment) 。
- GMM-EM的E步:每个点以概率分配给各个簇(后验概率)——这是软分配(Soft Assignment) 。
K-Means相当于假设高斯混合模型(GMM) 中每个高斯分量的协方差矩阵相同、各向同性(isotropic)、且方差趋于0。在这个极限情况下,后验概率退化为one-hot的硬分配,EM算法退化为K-Means。
从另一个角度来看,K-Means(Lloyd算法)可以看作是变分EM(Variational EM) 的一个特例,其中使用了截断后验(truncated posteriors) 作为变分分布。这种视角的一个重要优势是:不需要假设方差趋于0就能从理论上将K-Means与GMM联系起来。
1# 对比:硬分配(K-Means风格)vs 软分配(GMM-EM风格)2
3def hard_assignment(points, centroids):4 """硬分配:每个点完全属于一个簇"""5 distances = np.array([[np.linalg.norm(p - c) for c in centroids] for p in points])6 assignments = np.argmin(distances, axis=1)7 # 返回one-hot编码的硬分配8 hard = np.zeros((len(points), len(centroids)))9 hard[np.arange(len(points)), assignments] = 110 return hard11
12def soft_assignment(points, centroids, variances):13 """软分配:每个点以概率属于各个簇(GMM风格)"""14 n_points, n_clusters = len(points), len(centroids)15 probs = np.zeros((n_points, n_clusters))16 for k in range(n_clusters):17 # 计算每个点属于簇k的概率(假设高斯分布)18 diff = points - centroids[k]19 probs[:, k] = np.exp(-0.5 * np.sum(diff**2, axis=1) / variances[k])20 # 归一化为概率分布21 probs = probs / np.sum(probs, axis=1, keepdims=True)22 return probs23
24# 示例25points = np.array([[0, 0], [1, 1], [10, 10]])26centroids = np.array([[0, 0], [10, 10]])27variances = np.array([1.0, 1.0])28
29print("硬分配(K-Means风格):")30print(hard_assignment(points, centroids))31# 输出: 前两个点属于簇0,最后一个点属于簇132
33print("\n软分配(GMM-EM风格):")34print(soft_assignment(points, centroids, variances))35# 输出: 每个点有概率分布,如 [0.95, 0.05] 表示95%概率属于簇02.3 从硬分配到软分配:K-Means → GMM
K-Means和GMM-EM的对应关系揭示了一个重要的洞察:
K-Means是GMM-EM在“硬分配”极限下的特例。如果我们把“硬分配”放宽为“软分配”(即允许每个点以概率属于多个簇),就得到了高斯混合模型。
这个关系也解释了K-Means的一个核心假设:它假设所有簇都是球形且大小相近的。因为当所有高斯分量的协方差矩阵相同且各向同性时,决策边界是球形的——这正是K-Means用欧氏距离做最近邻分配所隐含的几何假设。
如果数据中的簇不是球形(比如拉长的椭圆),或者大小差异很大,K-Means的表现就会大打折扣。
三、K-Means++:聪明的初始化
3.1 为什么初始化如此重要?
Lloyd算法只能保证收敛到局部最优。不同的初始中心会导向完全不同的聚类结果。
一个糟糕的初始化可能导致:
- 收敛到很差的局部最优(簇分配不合理)
- 收敛速度极慢
- 某些簇没有分配到任何点(空簇问题)
传统的做法是:多次随机初始化,选择目标函数最小的结果。但这种方法计算量大,且不能保证找到好的初始点。
3.2 K-Means++的核心思想
K-Means++ 由David Arthur和Sergei Vassilvitskii于2007年提出,它用一种概率采样的方式选择初始中心,使得初始中心尽可能分散。
K-Means++的初始化步骤:
步骤1:从数据点中均匀随机选择第一个簇中心 μ1。
步骤2:对于每个数据点 x,计算它到最近已选中心的距离平方 D(x)2。
步骤3:以与 D(x)2 成正比的概率选择下一个簇中心(即距离已有中心越远的点,被选中的概率越大)——这称为 D2 采样(D²-sampling) 。
步骤4:重复步骤2-3,直到选出 K 个初始中心。
1def kmeans_plusplus_init(X, n_clusters, random_state=None):2 """K-Means++初始化算法的简化实现"""3 np.random.seed(random_state)4 n_samples, n_features = X.shape5
6 # 步骤1:随机选择第一个中心7 centroids = [X[np.random.choice(n_samples)]]8
9 # 步骤2-4:迭代选择剩余的中心10 for _ in range(1, n_clusters):11 # 计算每个点到最近已选中心的距离平方12 distances = np.array([13 min([np.linalg.norm(x - c) ** 2 for c in centroids])14 for x in X15 ])16 # 按概率(与距离平方成正比)选择下一个中心17 probs = distances / np.sum(distances)18 next_idx = np.random.choice(n_samples, p=probs)19 centroids.append(X[next_idx])20
21 return np.array(centroids)22
23# 对比:随机初始化 vs K-Means++24X, _ = make_blobs(n_samples=300, centers=3, cluster_std=0.8, random_state=42)25
26# 随机初始化27random_init = X[np.random.choice(len(X), 3, replace=False)]28print("随机初始化中心:\n", random_init)29
30# K-Means++初始化31kpp_init = kmeans_plusplus_init(X, 3, random_state=42)32print("\nK-Means++初始化中心:\n", kpp_init)33# K-Means++的中心通常更分散,覆盖数据的各个区域3.3 K-Means++的理论保证
K-Means++最吸引人的地方在于它的理论保证:
K-Means++能够以 O(logK) 的近似比逼近最优聚类代价。
通俗地说:K-Means++找到的初始中心,其对应的聚类代价(目标函数值)不会比全局最优值差太多——差距被控制在 O(logK) 倍以内。
这就是为什么scikit-learn的KMeans默认使用init='k-means++'——它用一个简单的概率采样策略,显著提高了找到高质量聚类的概率。
四、如何选择K?肘部法则与轮廓系数
K-Means要求用户预先指定簇的数量 K ——但在实际问题中,我们往往不知道数据应该分成几类。这就引出了聚类中最棘手的问题之一:如何选择 K?
4.1 肘部法则(Elbow Method)
肘部法则的思路很直观:
- 对 K=1,2,3,...,Kmax 分别运行K-Means
- 记录每个 K 对应的畸变(Distortion/Inertia) ——即所有点到其簇中心的距离平方和
- 绘制 K 与畸变的关系曲线
- 选择曲线**“肘部”**对应的 K 值
为什么畸变会随 K 增加而下降?因为簇越多,每个点离自己的中心就越近。但增加的簇带来的畸变下降会逐渐减小——那个“拐点”就是肘部。
1from sklearn.cluster import KMeans2
3def elbow_method(X, max_k=10):4 """肘部法则:计算不同K值下的惯量"""5 inertias = []6 K_range = range(1, max_k + 1)7
8 for k in K_range:9 kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)10 kmeans.fit(X)11 inertias.append(kmeans.inertia_)12
13 # 绘制肘部曲线14 plt.figure(figsize=(8, 5))15 plt.plot(K_range, inertias, 'bo-')16 plt.xlabel('簇的数量 K')17 plt.ylabel('惯量 (Inertia)')18 plt.title('肘部法则')19 plt.grid(True)20 plt.show()21
22 return inertias23
24# 生成数据并应用肘部法则25X, _ = make_blobs(n_samples=300, centers=4, cluster_std=0.6, random_state=42)26inertias = elbow_method(X, max_k=10)27# 观察:曲线在K=4处出现明显的“肘部”4.2 肘部法则的严重局限
肘部法则有一个致命的弱点:它严重缺乏理论支撑。有研究者甚至呼吁 “停止使用肘部法则” 。
肘部法则的主要问题:
问题一:主观性强。 什么是“肘部”?不同的人可能看到不同的拐点。在很多数据集上,畸变曲线是平滑下降的,根本没有明显的肘部。
问题二:对数据分布敏感。 当数据分布不均匀、有噪声、或簇的形状不是球形时,肘部法则的结果往往不可靠。
问题三:缺乏理论保证。 肘部法则只是一个启发式规则,没有任何统计理论保证它选择的 K 是最优的。
4.3 轮廓系数(Silhouette Coefficient)
轮廓系数由Rousseeuw于1987年提出,它衡量的是每个点与自身簇的紧密程度相对于与最近邻簇的分离程度。
对于第 i 个数据点,轮廓系数的计算分为三步:
步骤1:计算 a(i) —— 点 i 到同簇内其他点的平均距离(簇内凝聚度)。
步骤2:计算 b(i) —— 点 i 到最近的其他簇中所有点的平均距离(簇间分离度)。
步骤3:计算轮廓系数:
s(i)=max{a(i),b(i)}b(i)−a(i)轮廓系数的取值范围是 [−1,1] :
- s(i)≈1:点被很好地聚类(远亲不如近邻)
- s(i)≈0:点在两个簇的边界上
- s(i)≈−1:点被错误地聚类(应该属于其他簇)
通常,轮廓系数 > 0.7 表示聚类质量很好。
1from sklearn.metrics import silhouette_score2
3def silhouette_analysis(X, max_k=10):4 """轮廓系数分析:计算不同K值下的平均轮廓系数"""5 scores = []6 K_range = range(2, max_k + 1) # 轮廓系数至少需要2个簇7
8 for k in K_range:9 kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)10 labels = kmeans.fit_predict(X)11 score = silhouette_score(X, labels)12 scores.append(score)13
14 # 绘制轮廓系数曲线15 plt.figure(figsize=(8, 5))16 plt.plot(K_range, scores, 'ro-')17 plt.xlabel('簇的数量 K')18 plt.ylabel('平均轮廓系数')19 plt.title('轮廓系数分析')20 plt.axhline(y=0.7, color='green', linestyle='--', label='高质量聚类 (0.7)')21 plt.legend()22 plt.grid(True)23 plt.show()24
25 return scores26
27# 应用轮廓系数分析28scores = silhouette_analysis(X, max_k=10)29# 选择轮廓系数最高的K值30best_k = np.argmax(scores) + 231print(f"最佳K值: {best_k}, 轮廓系数: {scores[best_k-2]:.4f}")4.4 轮廓系数的局限
轮廓系数也有其局限:
计算开销大。 轮廓系数需要计算所有点对之间的距离,复杂度为 O(N2)。对于大规模数据集,计算可能非常耗时。
对簇形状敏感。 和K-Means一样,轮廓系数假设簇是凸的、球形的。对于非凸形状的簇,轮廓系数可能给出误导性的结果。
没有绝对标准。 虽然0.7被视为“好”的阈值,但这个阈值是经验性的,不是理论保证的。
4.5 如何科学地选择K?
基于以上的讨论,这里给出一些更可靠的选择K的方法:
| 方法 | 优点 | 缺点 |
|---|---|---|
| 领域知识 | 最可靠 | 并非总有领域知识 |
| 肘部法则 | 简单直观 | 主观、缺乏理论支撑 |
| 轮廓系数 | 有明确的统计解释 | 计算量大 |
| Gap Statistic | 有统计理论支撑 | 计算复杂 |
| 下游任务评估 | 最实用 | 需要定义具体任务 |
最实际的做法:结合多种方法,并用下游任务(如分类、异常检测)的绩效来验证聚类结果的有效性。毕竟,聚类的最终目的是服务于某个实际任务。
五、K-Means的假设与适用场景
5.1 K-Means的三个关键假设
理解K-Means的隐式假设,是正确使用它的前提:
假设一:簇是球形的(Spherical) 。K-Means使用欧氏距离,决策边界是超球面。如果数据中的簇是拉长的、椭圆的或任意形状的,K-Means可能表现不佳。
假设二:簇的大小相近(Equal Size) 。K-Means倾向于产生大小相近的簇。如果一个簇很大、另一个很小,大簇可能会“吞掉”小簇的部分点。
假设三:簇的密度相近(Equal Density) 。K-Means的分配基于距离,如果不同簇的密度差异很大,密度高的簇会被过度分割。
5.2 什么时候用K-Means?
K-Means最适合以下场景:
- 数据大致满足球形簇的假设
- 簇的大小和密度相近
- 数据维度不是特别高(维度灾难会影响距离度量的有效性)
- 需要快速、可扩展的聚类方案
- 数据量大,需要线性或近线性的时间复杂度
5.3 什么时候不用K-Means?
- 簇的形状是非凸的、拉长的、环形的——考虑DBSCAN或谱聚类
- 簇的大小差异巨大——考虑层次聚类或GMM
- 数据中存在大量噪声或离群点——K-Means对离群点敏感
- 数据维度很高(如 > 50维)——先做降维(如PCA)再聚类
六、完整代码示例
1import numpy as np2import matplotlib.pyplot as plt3from sklearn.datasets import make_blobs, make_circles, make_moons4from sklearn.cluster import KMeans5from sklearn.metrics import silhouette_score, adjusted_rand_score6from sklearn.preprocessing import StandardScaler7
8# 1. 生成不同形状的数据9np.random.seed(42)10
11# 数据集1:球形簇(K-Means擅长)12X1, y1 = make_blobs(n_samples=500, centers=4, cluster_std=0.6, random_state=42)13
14# 数据集2:环形数据(K-Means不擅长)15X2, y2 = make_circles(n_samples=500, factor=0.5, noise=0.05, random_state=42)16
17# 数据集3:月牙形数据(K-Means不擅长)18X3, y3 = make_moons(n_samples=500, noise=0.05, random_state=42)19
20# 2. 标准化(K-Means对尺度敏感)21scaler = StandardScaler()22X1_scaled = scaler.fit_transform(X1)23X2_scaled = scaler.fit_transform(X2)24X3_scaled = scaler.fit_transform(X3)25
26# 3. 对三个数据集运行K-Means27fig, axes = plt.subplots(2, 3, figsize=(15, 8))28datasets = [(X1_scaled, y1, '球形簇'),29 (X2_scaled, y2, '环形数据'),30 (X3_scaled, y3, '月牙形数据')]31
32for idx, (X, y_true, title) in enumerate(datasets):33 # 运行K-Means34 kmeans = KMeans(n_clusters=4 if idx == 0 else 2,35 random_state=42, n_init=10)36 labels = kmeans.fit_predict(X)37
38 # 可视化39 ax = axes[0, idx]40 ax.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis', alpha=0.6)41 ax.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1],42 c='red', marker='X', s=200, label='簇中心')43 ax.set_title(f'{title}\nK-Means聚类结果')44 ax.legend()45
46 # 评估47 ax2 = axes[1, idx]48 ax2.scatter(X[:, 0], X[:, 1], c=y_true, cmap='viridis', alpha=0.6)49 ax2.set_title(f'{title}\n真实标签')50
51 # 打印指标52 sil_score = silhouette_score(X, labels)53 ari_score = adjusted_rand_score(y_true, labels)54 print(f"{title}: 轮廓系数={sil_score:.4f}, ARI={ari_score:.4f}")55
56plt.tight_layout()57plt.show()58
59# 结论:K-Means在球形簇上表现优异,在环形和月牙形数据上表现较差七、总结
| 概念 | 核心内容 |
|---|---|
| Lloyd算法 | 交替执行“分配”(E步)和“更新”(M步),保证收敛到局部最优 |
| 目标函数 | J=∑∥xi−μci∥2,即畸变(Distortion)/惯量(Inertia) |
| K-Means与EM | K-Means是EM的特例——硬分配 + 各向同性高斯 |
| 硬分配 vs 软分配 | K-Means做硬分配(每个点100%属于一个簇),GMM-EM做软分配(概率分配) |
| K-Means++ | 用 D2 采样选择分散的初始中心,有 O(logK) 的近似保证 |
| 肘部法则 | 绘制 K-畸变曲线找“肘部”,缺乏理论支撑,不推荐作为唯一依据 |
| 轮廓系数 | s(i)=max{a(i),b(i)}b(i)−a(i),>0.7表示聚类质量好 |
| K-Means的假设 | 簇是球形的、大小相近、密度相近 |
核心要点回顾
-
Lloyd算法是K-Means的标准实现,通过交替优化簇分配和簇中心来最小化畸变。算法保证收敛,但只能收敛到局部最优。
-
K-Means是EM算法的特例——它相当于用硬分配(而非软分配)来处理高斯混合模型,且假设所有高斯分量的协方差矩阵相同且各向同性。理解这个联系,就能理解K-Means的局限性:它只适用于球形、大小相近的簇。
-
K-Means++通过概率采样(D2 采样)选择分散的初始中心,显著提高了找到高质量聚类的概率,且有 O(logK) 的理论保证。
-
肘部法则虽然有“肘部”这个直观概念,但严重缺乏理论支撑。轮廓系数提供了更具体的统计解释,但计算开销大。在实践中,应结合多种方法,并以下游任务的表现来验证聚类结果。
-
K-Means的适用场景是球形、大小相近、密度相近的簇。如果数据不满足这些假设(如环形、拉长、密度差异大),应考虑DBSCAN、谱聚类或GMM等其他算法。
Some information may be outdated