从零实现K-means聚类算法:核心原理、Python代码与实战应用

📅 2026/8/5 19:23:03
从零实现K-means聚类算法:核心原理、Python代码与实战应用
1. 从“物以类聚”到代码实现K-means聚类的核心思想我们每天都在做“分类”这件事。整理书架时你会把技术书、小说、工具书分开摆放整理照片时你可能会按时间、地点或人物来分组。这种“物以类聚”的直觉正是聚类算法的核心。在数据科学和机器学习领域K-means算法就是将这种直觉数学化、自动化的经典工具。它不需要你事先告诉它“这是小说那是工具书”而是让算法自己从一堆杂乱无章的数据点中找出内在的规律和分组。想象一下你有一片散布着星星的夜空肉眼看去杂乱无章。K-means算法就像一位天文学家它的任务是找出这些星星中哪些是聚在一起的星团。它不知道星团应该有多少个也不知道星团的边界在哪但它会通过计算星星之间的距离反复尝试最终将星空划分成几个清晰的区域每个区域内的星星彼此靠近而不同区域的星星则相对疏远。这个过程就是聚类。K-means之所以成为最流行、最易理解的聚类算法之一关键在于它的简洁和高效。它的核心思想可以用一句话概括通过迭代计算将N个数据点划分到K个簇中使得每个数据点到其所属簇的“中心点”的距离平方和最小。这个“中心点”在算法中被称为“质心”。整个算法的目标就是找到K个最优的质心位置以及每个数据点的归属从而让簇内的点尽可能相似距离质心近簇间的点尽可能不同距离其他质心远。今天我们就抛开复杂的数学公式直接从Python3的代码实现入手手把手带你从零构建一个可运行的K-means算法。我会在代码的每一行关键处解释其背后的数学原理和设计意图并分享在实际项目中调试参数、评估效果、避开常见陷阱的实战经验。无论你是刚入门机器学习的新手还是想巩固基础、了解底层实现的老手这篇从代码反推原理的实践指南都能让你对K-means有更“手感”的理解。2. 算法骨架与核心概念拆解距离、质心与迭代在动手写代码之前我们必须彻底理解K-means算法的几个核心构件。如果把算法比作一台机器那么输入的数据、距离度量方式、质心的初始化和更新规则就是这台机器的齿轮和轴承。理解它们你才能知道代码每一行在做什么以及为什么这么做。2.1 算法的输入与输出数据与标签K-means算法的输入非常简单一个形状为(n_samples, n_features)的二维数组X。n_samples代表你有多少个数据点n_features代表每个数据点有多少个特征维度。例如如果你想对顾客进行聚类每个顾客可能有“年龄”、“年收入”、“消费频率”三个特征那么n_features就是3。输出则是每个数据点所属的簇标签一个形状为(n_samples,)的一维数组标签通常是0到K-1的整数。这里有一个至关重要的前提K-means假设数据特征都是数值型的并且最好是连续值。如果你有类别型数据如“性别”、“城市”需要先进行独热编码等处理将其转化为数值。此外由于算法使用欧氏距离如果不同特征的数量级差异巨大例如“年龄”范围0-100“年收入”范围0-1,000,000直接计算距离会导致数量级大的特征主导结果。因此对数据进行标准化如Z-score标准化或归一化缩放到[0,1]区间是必不可少的预处理步骤。很多初学者忽略了这一步导致聚类结果完全失真这是第一个要避开的坑。2.2 距离的度量欧氏距离及其意义K-means默认使用欧氏距离来衡量数据点之间的相似度。对于两个点p和q其欧氏距离计算公式为distance sqrt((p1-q1)^2 (p2-q2)^2 ... (pn-qn)^2)这个公式的几何意义非常直观就是在多维空间中的直线距离。K-means的目标函数——最小化所有点到其所属质心的欧氏距离平方和——正是基于此。选择距离平方和而不是直接的距离和在数学上更便于求导和优化同时对大距离的点施加了更大的惩罚使得质心对异常值不那么敏感但依然敏感。注意虽然欧氏距离最常用但它并不是唯一选择。在处理文本数据如TF-IDF向量时余弦相似度可能更合适在某些特定领域曼哈顿距离也可能被使用。不过修改距离度量通常意味着目标函数和质心更新公式也需要调整这会衍生出K-medoids等变种算法。在我们的基础实现中我们坚守最经典的欧氏距离版本。2.3 质心簇的“引力中心”质心是每个簇的虚拟中心点它的坐标由属于该簇的所有数据点的坐标平均值计算得出。这就是“means”均值一词的由来。例如一个簇里有三个点(1,2), (2,3), (3,4)那么该簇的质心就是((123)/3, (234)/3) (2, 3)。质心不一定是一个真实存在的数据点它只是一个计算出来的代表位置。质心的初始化至关重要糟糕的初始质心可能导致算法收敛到局部最优解甚至收敛速度很慢。最常见的初始化方法是随机选择K个数据点作为初始质心。我们将在代码中实现这种方法并讨论其局限性。2.4 迭代的两步曲分配与更新K-means算法的迭代过程清晰得像一首二重奏不断重复两个步骤直到质心稳定分配步骤Assignment遍历每一个数据点计算它到当前K个质心的欧氏距离然后将该点分配给距离最近的那个质心所在的簇。这一步结束后每个数据点都有了一个新的、临时的簇标签。更新步骤Update对于每一个簇重新计算它的质心。新的质心坐标等于该簇内所有数据点各维度坐标的算术平均值。这两个步骤循环进行。如何判断算法已经“收敛”可以停止了呢通常有两种标准一是质心的位置在连续两次迭代中不再发生变化或变化小于一个极小的阈值二是数据点的簇分配不再发生变化。在代码中我们通常会设置一个最大迭代次数防止在无法收敛的情况下陷入无限循环。理解了这些核心概念我们的大脑里已经搭建起了算法的逻辑框架。接下来我们就用Python3代码为这个框架注入生命。3. 手把手实现从零构建K-means类我们不依赖sklearn完全从零开始构建一个KMeans类。这个过程会让你对算法的每一个细节都了如指掌。我们将这个类设计得与sklearn的接口类似包含fit和predict方法方便理解和使用。3.1 类的初始化与参数设计首先我们定义类的结构。一个健壮的K-means实现需要考虑哪些参数呢import numpy as np import matplotlib.pyplot as plt from sklearn.datasets import make_blobs # 用于生成演示数据 class MyKMeans: def __init__(self, n_clusters8, max_iter300, tol1e-4, random_stateNone): 初始化MyKMeans聚类器。 参数: n_clusters (int): 要形成的簇的数量以及要生成的质心数量。默认8。 max_iter (int): 单次运行的最大迭代次数。默认300。 tol (float): 关于两次连续迭代的质心差异的容差用于声明收敛。默认1e-4。 random_state (int): 确定质心初始化的随机数生成状态。 self.n_clusters n_clusters self.max_iter max_iter self.tol tol self.random_state random_state self.centroids None # 质心坐标形状 (n_clusters, n_features) self.labels_ None # 每个样本的簇标签形状 (n_samples,) self.inertia_ None # 样本到其最近质心的距离平方和用于评估聚类效果 def _init_centroids(self, X): 随机初始化质心。从数据点中随机选择n_clusters个点作为初始质心。 np.random.seed(self.random_state) # 随机选择不重复的索引 random_idx np.random.permutation(X.shape[0]) centroids X[random_idx[:self.n_clusters]] return centroids参数解读与经验谈n_clusters (K值)这是K-means算法最核心、也是最让人头疼的参数。算法本身无法知道数据应该分成几类K值需要你事先指定。选小了不同类别的数据会被强行合并选大了一个自然的类别又会被拆散。如何确定K值最常用的方法是“肘部法则”我们会在后面专门讨论。这里先将其作为一个必须由用户提供的参数。max_iter和tol这是控制算法停止的条件。max_iter是安全网防止在数据难以收敛时无限循环。tol定义了“质心稳定”的精度。通常1e-4是个不错的选择。在实际运行中如果数据量很大你可能需要适当增大tol如1e-3以提前停止用微小的精度损失换取显著的计算时间节省。random_state设置随机种子是为了让实验结果可复现。在调试和对比不同参数时固定随机种子至关重要。否则每次运行因初始化不同可能得到不同的结果会让你无法判断是参数的影响还是随机性的影响。inertia_这个属性非常重要它记录了所有样本到其所属质心距离的平方和也称为“簇内平方和”或“畸变程度”。这个值越小说明簇内样本越紧密。它是评估聚类效果和选择K值的关键指标。3.2 核心迭代循环分配与更新的代码实现接下来是算法的核心——fit方法。它接收数据X并通过迭代找到最优的质心和样本分配。def fit(self, X): 计算K-means聚类。 参数: X (array-like): 训练数据形状 (n_samples, n_features) X np.array(X) n_samples, n_features X.shape # 1. 初始化质心 self.centroids self._init_centroids(X) # 开始迭代 for i in range(self.max_iter): # 保存旧的质心用于收敛判断 old_centroids self.centroids.copy() # 2. 分配步骤计算每个样本到每个质心的距离并分配标签 distances self._calc_distances(X, self.centroids) # 找到每个样本距离最近的质心索引即簇标签 self.labels_ np.argmin(distances, axis1) # 3. 更新步骤根据新的样本分配重新计算质心 new_centroids np.zeros((self.n_clusters, n_features)) for k in range(self.n_clusters): # 获取属于第k簇的所有样本 cluster_k X[self.labels_ k] # 防止空簇如果某个簇没有样本则保留旧质心或重新初始化 if len(cluster_k) 0: # 策略重新随机初始化该质心 new_centroids[k] X[np.random.randint(0, n_samples)] else: # 计算新质心簇内所有样本的均值 new_centroids[k] cluster_k.mean(axis0) self.centroids new_centroids # 4. 检查收敛如果质心移动很小则停止迭代 centroid_shift np.linalg.norm(old_centroids - self.centroids) if centroid_shift self.tol: print(f迭代在第 {i1} 轮收敛。) break # 计算最终的 inertia_ final_distances self._calc_distances(X, self.centroids) # 取每个样本到其所属质心的距离即最小距离 min_distances np.min(final_distances, axis1) self.inertia_ np.sum(min_distances ** 2) return self def _calc_distances(self, X, centroids): 计算每个样本到每个质心的欧氏距离。使用向量化操作提高效率。 # 利用广播机制计算距离矩阵 # 形状: (n_samples, 1, n_features) 和 (1, n_clusters, n_features) 相减 # 结果形状: (n_samples, n_clusters, n_features) differences X[:, np.newaxis, :] - centroids[np.newaxis, :, :] squared_differences differences ** 2 # 沿特征轴求和并开方得到距离矩阵 (n_samples, n_clusters) distances np.sqrt(np.sum(squared_differences, axis2)) return distances代码细节与避坑指南向量化计算距离在_calc_distances方法中我们使用了NumPy的广播机制一次性计算所有样本到所有质心的距离避免了低效的Python层循环。这是实现性能的关键。对于大数据集这个操作可能会消耗大量内存距离矩阵大小为n_samples * n_clusters如果内存不足可能需要分块计算。空簇问题在更新质心的循环中我们加入了if len(cluster_k) 0的判断。这是一个非常重要的边界情况处理。如果某个质心在分配步骤后没有分配到任何样本它就成了“空簇”。如果不处理计算均值时会出错。我们的策略是随机选择一个数据点作为该簇的新质心。其他常见策略包括选择距离当前所有质心最远的点或者直接移除该簇减少K值。空簇的出现往往是K值设置过大或初始化太差的信号。收敛判断我们使用质心移动的欧氏范数np.linalg.norm来度量变化。当这个变化小于容差tol时认为算法已收敛。打印收敛信息有助于调试。3.3 预测与接口方法实现fit之后我们还需要predict方法用于对新数据或原有数据进行簇标签预测。同时为了完整性我们也实现一个计算inertia_的方法。def predict(self, X): 预测X中每个样本所属的簇。 参数: X (array-like): 形状 (n_samples, n_features) 的数据 返回: labels (array): 形状 (n_samples,) 的簇标签 X np.array(X) distances self._calc_distances(X, self.centroids) labels np.argmin(distances, axis1) return labels def fit_predict(self, X): 同时进行拟合和预测返回样本标签。 self.fit(X) return self.labels_现在我们的MyKMeans类已经具备了基本功能。让我们用一段简单的代码来测试它。4. 实战测试可视化与结果分析理论说得再好不如跑一遍代码看看。我们使用sklearn的make_blobs函数生成一个易于可视化的二维数据集它本身就会生成几个高斯分布的“团块”非常适合测试聚类算法。4.1 生成数据与运行算法# 1. 生成模拟数据 np.random.seed(42) # 生成500个样本4个中心点二维特征标准差为0.8 X, y_true make_blobs(n_samples500, centers4, n_features2, cluster_std0.8, random_state42) # 2. 创建并训练我们的K-means模型 kmeans MyKMeans(n_clusters4, max_iter300, tol1e-4, random_state42) kmeans.fit(X) labels kmeans.labels_ centroids kmeans.centroids print(f质心坐标:\n{centroids}) print(f簇内平方和 (inertia_): {kmeans.inertia_:.2f})4.2 可视化聚类结果可视化是理解聚类结果最直观的方式。我们将真实标签如果有的话、聚类结果和质心位置一起画出来。# 3. 可视化结果 fig, (ax1, ax2) plt.subplots(1, 2, figsize(14, 5)) # 子图1真实分布如果已知 ax1.scatter(X[:, 0], X[:, 1], cy_true, cmapviridis, s30, edgecolork, alpha0.7) ax1.scatter(centroids[:, 0], centroids[:, 1], cred, markerX, s200, labelPredicted Centroids) ax1.set_title(True Cluster Distribution) ax1.set_xlabel(Feature 1) ax1.set_ylabel(Feature 2) ax1.legend() ax1.grid(True, linestyle--, alpha0.5) # 子图2K-means聚类结果 scatter ax2.scatter(X[:, 0], X[:, 1], clabels, cmapviridis, s30, edgecolork, alpha0.7) ax2.scatter(centroids[:, 0], centroids[:, 1], cred, markerX, s200, labelCentroids) ax2.set_title(fK-means Clustering Result (K{kmeans.n_clusters})) ax2.set_xlabel(Feature 1) ax2.set_ylabel(Feature 2) ax2.legend() ax2.grid(True, linestyle--, alpha0.5) plt.tight_layout() plt.show()运行这段代码你应该能看到两个并排的散点图。左图显示了数据真实的四个分布中心用不同颜色表示右图显示了我们的K-means算法找到的四个簇和质心红色的“X”。在数据分离度较好的情况下我们的算法应该能非常准确地还原出数据的真实结构质心也会落在每个簇的密度中心附近。结果分析通过对比两个图你可以直观地评估聚类效果。如果右图中的颜色块聚类结果与左图基本一致且边界清晰说明聚类是成功的。inertia_的值给出了一个量化的评估这个值越小越好。但要注意inertia_会随着K值的增大而单调减小因为每个点离自己的质心更近所以不能单纯用它来比较不同K值下的模型。5. 超越基础K值选择、评估与高级话题一个能运行的K-means只是开始。在实际项目中更大的挑战在于如何确定K值如何评估聚类效果的好坏以及如何处理K-means的固有缺陷5.1 如何选择最佳的K值“肘部法则”实战K-means最大的痛点就是需要预先指定K值。肘部法则是最常用的启发式方法。其思想是随着K值增大簇内样本会更紧密inertia_会下降。下降的幅度会在某个点突然变缓这个拐点就像手肘的关节对应的K值可能就是最佳选择。def plot_elbow_method(X, max_k10): 绘制不同K值对应的inertia_曲线寻找肘部。 inertias [] K_range range(1, max_k1) for k in K_range: kmeans MyKMeans(n_clustersk, random_state42) kmeans.fit(X) inertias.append(kmeans.inertia_) plt.figure(figsize(8, 5)) plt.plot(K_range, inertias, bo-) plt.xlabel(Number of clusters (K)) plt.ylabel(Inertia (簇内平方和)) plt.title(Elbow Method For Optimal K) plt.xticks(K_range) plt.grid(True, linestyle--, alpha0.5) plt.show() # 使用之前生成的数据X plot_elbow_method(X, max_k10)运行后你会看到一条下降曲线。曲线开始陡峭下降然后逐渐平缓。那个“拐弯”的点比如从K3到K4下降幅度明显变小对应的K值可能是3或4就是肘部。这个方法很直观但有时拐点并不明显需要结合业务理解来判断。5.2 聚类效果评估当没有真实标签时在有真实标签的数据集上我们可以使用调整兰德指数ARI或标准化互信息NMI等外部指标来评估。但在无监督学习中我们通常没有真实标签。这时可以使用轮廓系数。轮廓系数结合了簇内凝聚度和簇间分离度。对于每个样本ia(i)样本i到同簇其他样本的平均距离凝聚度。b(i)样本i到其他某簇所有样本的平均距离的最小值分离度。样本i的轮廓系数s(i) (b(i) - a(i)) / max(a(i), b(i))。S(i)的取值范围在[-1, 1]之间。越接近1说明样本聚类越合理越接近-1说明样本可能被分错了簇接近0则说明样本在两个簇的边界上。所有样本的轮廓系数的平均值可以作为整个聚类结果的评价指标。from sklearn.metrics import silhouette_score # 计算我们之前K4时的轮廓系数 score silhouette_score(X, labels) print(fK4时轮廓系数为: {score:.3f}) # 我们可以计算不同K值下的轮廓系数选择最高的 best_k 0 best_score -1 for k in range(2, 11): # 轮廓系数要求至少2个簇 kmeans MyKMeans(n_clustersk, random_state42) labels_k kmeans.fit_predict(X) score_k silhouette_score(X, labels_k) print(fK{k}: 轮廓系数 {score_k:.3f}) if score_k best_score: best_score score_k best_k k print(f\n最佳K值基于轮廓系数: {best_k})轮廓系数越高越好。它和肘部法则结合使用能更可靠地确定K值。5.3 K-means的局限性非球形簇与噪声我们必须清醒地认识到K-means的假设和局限这是避免误用的关键假设簇是凸形的、各向同性的K-means使用欧氏距离天然地倾向于发现球状或超球状的簇。对于流形、环形或任意形状的簇K-means效果会很差。下图展示了K-means对环形数据的失败案例。对噪声和异常值敏感由于使用均值作为质心异常值会极大地拉偏质心的位置。需要指定K值如前所述这是一个需要先验知识或多次尝试的参数。初始化敏感性随机初始化可能导致不同的局部最优解。一个改进方案是运行多次如10次并选择inertia_最小的那次结果。# 生成一个环形数据展示K-means的局限 from sklearn.datasets import make_circles X_circle, _ make_circles(n_samples300, factor0.5, noise0.05, random_state42) kmeans_circle MyKMeans(n_clusters2, random_state42) labels_circle kmeans_circle.fit_predict(X_circle) plt.figure(figsize(6, 6)) plt.scatter(X_circle[:, 0], X_circle[:, 1], clabels_circle, cmapviridis, s30) plt.scatter(kmeans_circle.centroids[:, 0], kmeans_circle.centroids[:, 1], cred, markerX, s200, labelCentroids) plt.title(K-means Fails on Non-convex Clusters (Circles)) plt.xlabel(Feature 1) plt.ylabel(Feature 2) plt.legend() plt.grid(True, linestyle--, alpha0.5) plt.axis(equal) plt.show()运行这段代码你会看到K-means强行将两个环形簇按距离分成了内外两个“半圆”状的簇这完全违背了数据的真实结构。对于这类数据DBSCAN基于密度的聚类或谱聚类等算法是更好的选择。6. 性能优化与生产环境考量我们实现的MyKMeans是一个教学版本强调了可读性。但在处理真实世界的大规模数据时性能至关重要。以下是几个关键的优化方向距离计算的优化我们使用了向量化计算这已经比循环快很多。但对于超大数据集计算所有样本到所有质心的距离矩阵O(n_samples * n_clusters * n_features)依然昂贵。可以使用更高效的距离计算库如scipy.spatial.distance.cdist或者采用三角不等式进行加速的算法变种如Elkan K-means。sklearn的KMeans默认就使用了Elkan算法。初始化优化K-means我们使用的是随机初始化这可能导致收敛慢或效果差。K-means是一种智能初始化方案它选择彼此相距较远的点作为初始质心能显著提高收敛速度和最终结果的质量。其核心思想是第一个质心随机选后续每个质心被选中的概率与它到已选质心的最短距离的平方成正比。这保证了初始质心分散在数据空间中。Mini-Batch K-means对于海量数据如数百万样本即使算法是O(n)复杂度单次迭代也可能很慢。Mini-Batch K-means每次迭代只使用数据的一个随机子集mini-batch来更新质心极大地减少了计算量通常能以轻微的质量损失换取巨大的速度提升。并行化距离计算和样本分配是天然可并行的。可以利用多核CPU或GPU进行加速。在实际项目中我强烈建议直接使用sklearn.cluster.KMeans它已经集成了K-means初始化、Elkan/ Lloyd优化算法、并行计算等高级特性并且经过了高度优化和测试。我们自己实现的目的是为了深入理解而不是为了替代成熟的库。7. 从代码到应用K-means能做什么理解了原理和实现我们来看看K-means在现实世界中的典型应用场景这能帮你更好地将知识落地客户细分根据用户的购买历史、 demographics人口统计信息、行为数据等将客户分成不同的群组以便进行精准营销。例如发现“高价值低频次”用户和“低价值高频次”用户并采取不同的策略。图像压缩颜色量化一张彩色图片可能有数百万种颜色。使用K-means可以将所有像素的颜色聚类成K种比如64种然后用每个簇的质心颜色代替簇内所有像素的颜色。这样在视觉损失不大的情况下能大幅减少存储空间。这就是GIF图像常用的技术。文档聚类将文本文档转化为TF-IDF向量后可以使用K-means进行聚类自动发现讨论相似主题的文档集合。异常检测正常的数据点通常会形成紧密的簇而异常点则远离任何质心。通过计算每个点到最近质心的距离可以设定一个阈值来识别异常。推荐系统在协同过滤中可以先使用K-means对用户或物品进行聚类然后在簇内进行推荐这可以减少计算量称为“分群推荐”。在应用时牢记K-means的假设。如果你的数据不是球状簇或者有大量噪声不要强行使用K-means。先可视化你的数据或通过降维技术如PCA/t-SNE可视化对数据的结构有一个直观认识这是选择正确算法的第一步。最后分享一个我自己的经验永远不要完全相信自动聚类的结果。无论指标多好一定要结合业务知识去审视每一个簇。有时算法发现的“簇”只是数据分布的巧合有时有业务意义的细分可能因为特征选择不当而被算法忽略。机器学习模型是辅助决策的工具而不是替代人类判断的神谕。将聚类结果交给领域专家去解读和验证往往能碰撞出更有价值的洞见。