做数据分析的人十有八九用过K-Means但大多数人只是import一下就出图。真正麻烦的是那些import之后的事为什么运行两次结果不一样为什么某个簇总是空为什么归一化前后聚类天差地别这次我不打算糊弄过去直接从零手撕一遍K-Means聚类算法——不调sklearn不用现成的聚类函数只用Numpy把初始化、距离计算、样本分配、质心更新和收敛判定一步步写成代码。整个过程里踩过的坑、查过的资料、改过的bug全记录在下面。适合刚接触聚类算法的新手也适合把KMeans()当黑盒用了很久、实在受不了拍脑袋调参的老手。看完你会明白这个老朋友在每次迭代里到底偷偷做了什么。1. K-Means聚类算法到底在做什么思路先对齐代码才不会歪1.1 一句话直觉让“相似”变成“距离”K-Means处理的是没有标签的样本集合。聚类算法的目标是把这些样本自动分成K个组尽可能让同一个组里的样本相似不同组之间的样本差异明显。K-Means把“相似”这个词翻译成了数学语言“相似”就是“距离近”。更具体地说它只关心欧氏距离——两个点在多维空间里离得越近就越可能属于同一簇。你可以把它想象成数据版的“矿石分拣”。你有一堆混合矿石不知道里面含几种金属也不知道每块矿石具体属于哪一类。K-Means根据矿石的质地、颜色、重量这些特征把它们自动归成几堆。归完之后每一堆的中心就是这堆样本的代表也就是质心。整个过程没有老师告诉你正确答案纯靠数据本身的结构来分组所以这类算法也叫无监督学习。这里有个容易忽略的点K-Means不是一把梭直接算出来的。它先随便猜几个质心然后让每个样本“投靠”距离最近的质心再根据新组员的平均位置重新计算质心如此循环直到质心不再明显移动。这个“猜—分—算—再猜”的循环是整个算法的灵魂。搞懂这个循环后面所有代码都只是它的具体翻译。1.2 目标函数与硬分配K-Means凭什么能收敛只看直觉的话很容易把K-Means理解成“循环求平均”虽然没错但会漏掉最重要的东西——目标函数。K-Means的目标是让每个样本到它所属质心的距离平方和最小用公式写出来就是SSE Σ(k1..K) Σ(x_i ∈ C_k) ‖x_i - μ_k‖²其中C_k是第K个簇μ_k是该簇的质心也就是簇内样本的均值。这个指标通常叫簇内平方和SSEK-Means的每次迭代都在压低这个值。这个式子决定了K-Means的一个关键性质硬分配。每个样本在每一步里只属于离它最近的唯一簇不存在“有60%属于A簇、40%属于B簇”这种说法。argmin一锤定音样本要么是这簇的要么是那簇的没有中间地带。这一点听起来理所当然但跟后面要提到的高斯混合模型GMM形成了鲜明对比GMM是软聚类会给每个样本输出一个概率向量。目标函数存在还有另一层意义它解释了K-Means为什么一定收敛。经典的Lloyd算法本质上在交替优化两个变量——固定质心时给样本分配最近的簇固定簇时把质心更新成簇内均值。两步操作都在降低同一个SSE目标函数所以迭代过程是单调不增的最终必然收敛。只不过这个收敛不是“收敛到全局最优”而是收敛到某个局部极小值。SSE这个函数并不凸存在很多坑坑洼洼的谷底随机初始化不同跑到的谷底也不同。这是后面讲K-Means和n_init时反复要提到的前提。1.3 为什么偏偏用欧氏距离超球面偏好的由来很多刚上手的人会问能不能用曼哈顿距离能用余弦距离吗答案是能算但那个算法就不叫K-Means了数学性质也会彻底变掉。原因在于欧氏距离与“均值”天生是一对。给定一堆点能让“到所有点的欧氏距离平方和”最小的那个中心恰恰就是它们的算术平均。而K-Means的每一步更新都需要求均值两者刚好无缝衔接。如果把距离换成曼哈顿距离最优中心就变成了中位数算法演变成K-Medians或K-Medoids换了个物种。所以K-Means不是“非要用欧氏距离”而是“用均值做质心自然推导出了欧氏距离”。欧氏距离还有一个隐藏含义它对应的簇形状是各向同性的圆。用专业点的话说K-Means假设每个簇在特征空间里是一个超球面而且各个方向半径差不多。一旦数据里出现细长条状的簇、半包围的圆弧形簇K-Means就会生硬地把它拦腰切开或者跨越边界乱分。判断一个数据集适不适合用K-Means最简单的方法就是降维到2D看一眼如果散点图呈现几个圆团直接用如果形状妖娆趁早换DBSCAN或GMM。欧氏距离对量纲的敏感也是它在实际项目里最常被忽略的坑。想象两个特征“年龄”取值20到60“年消费额”取值5万到100万。计算欧氏距离时年龄那一位的贡献基本被消费金额淹没聚类结果跟年龄一点关系都没有。这不是算法错了是距离度量被量纲绑架了。所以在跑K-Means之前标准化几乎是最重要的一步预处理后面4.1会专门展开。2. 从零手撕一遍K-Means完整实现与收敛过程可视化2.1 数据准备先造一份能画出结果的2D数据集既然要手撕算法数据我也选择自己造。用Numpy直接生成三个高斯簇二维带真实中心只是用来验证效果算法本身完全不需要标签import numpy as np import matplotlib.pyplot as plt rng np.random.default_rng(42) centers_true np.array([[2.0, 2.0], [-2.0, 3.0], [3.0, -2.0]]) X np.vstack([ rng.normal(loccenter, scale0.6, size(100, 2)) for center in centers_true ]) rng.shuffle(X)为什么选2D数据因为可视化直观每一轮质心的移动都能画出来人眼直接看到算法“动起来”。实际业务里通常是几十维甚至上百维但算法逻辑完全一样只是没法用眼睛直接看罢了。scale0.6是给每个簇加适量噪声让三个团有轻微重叠但又不至于糊成一团这样K-Means既有挑战又能稳定收敛。shuffle这步容易被忽略如果不打乱样本顺序随机初始化质心时更容易连续抽到同一区域的样本增加空簇概率。虽然对这个干净数据集影响不大但这是防止“莫名其妙翻车”的好习惯。真实数据里样本顺序常常带有隐藏结构比如按时间排列、按来源排列先shuffle能避免初始化被顺序影响。2.2 核心代码逐段拆解初始化、分配、更新、收敛我把手撕版K-Means写成一个可直接运行的函数包含距离计算、空簇处理、质心位移记录。每个步骤都带注释方便对照原理def euclidean_dist(X, center): # X: (n, d), center: (d,) - 返回每个样本到该质心的距离 return np.sqrt(np.sum((X - center) ** 2, axis1)) def kmeans_custom(X, k, max_iter200, tol1e-4, seed42): rng np.random.default_rng(seed) # 1. 初始化从样本里随机挑k个样本当质心 initial_idx rng.choice(len(X), sizek, replaceFalse) centroids X[initial_idx].copy() for i in range(max_iter): # 2. 分配步骤计算每个样本到全部质心的距离 # distances 的形状是 (k, n)每行是“到某一个质心”的距离 distances np.array([euclidean_dist(X, c) for c in centroids]) labels np.argmin(distances, axis0) # 每个样本分给最近质心 # 3. 更新步骤对每个簇求均值作为新质心 new_centroids np.zeros_like(centroids) for c_idx in range(k): cluster_points X[labels c_idx] if len(cluster_points) 0: # 空簇处理随机补一个样本当质心避免实际簇数悄悄减一 new_centroids[c_idx] X[rng.choice(len(X))] else: new_centroids[c_idx] cluster_points.mean(axis0) # 4. 收敛判定质心位移小于阈值 shift np.linalg.norm(new_centroids - centroids) centroids new_centroids if shift tol: print(f迭代 {i 1} 轮后收敛总位移 {shift:.6f}) return centroids, labels print(f达到最大迭代次数 {max_iter}可能需要检查tol的合理性) return centroids, labels运行指定K3我这个测试里通常第3轮就收敛。三个质心大约落在(2.02, 1.98)、(-1.98, 2.94)、(2.95, -1.96)附近与真实中心(2,2)、(-2,3)、(3,-2)误差在0.1以内。这说明数据生成得干净K-Means在这种球状簇场景下简直是为它量身定做的。这里必须提醒一个新手最容易踩的代码级坑argmin的axis到底写几。在初版实现里distances的形状是(k, n)所以argmin(axis0)才能得到n个样本各自的簇编号。如果你手滑写成axis1取到的将是“对每个质心而言最近的样本索引”组成的列表形状直接变成(k,)结果完全不是簇标签而且大概率不报错。这类错误极其隐蔽我建议写完先打印distances.shape和labels.shape验证一下别直接往业务代码里塞。2.3 向量化提速用Numpy广播把循环压成矩阵运算教学版代码容易看懂但每次算距离都开一个循环数据量一上来就慢得可怕。实际生产需要的是同一逻辑的向量化版本。用Numpy广播一步算出所有样本到所有质心的距离矩阵def kmeans_vec(X, k, max_iter200, tol1e-4, seed42): rng np.random.default_rng(seed) idx rng.choice(len(X), sizek, replaceFalse) centroids X[idx].copy() for i in range(max_iter): # X: (n, d) - X[:, None, :] 形状(n, 1, d) # centroids: (k, d) - centroids[None, :, :] 形状(1, k, d) # 广播相减得到 (n, k, d)即每个样本对每个质心的分量差 diff X[:, None, :] - centroids[None, :, :] distances np.sqrt((diff ** 2).sum(axis2)) # 形状(n, k) labels np.argmin(distances, axis1) # 每个样本最近的质心 new_centroids np.zeros_like(centroids) for c_idx in range(k): cluster_points X[labels c_idx] if len(cluster_points) 0: new_centroids[c_idx] X[rng.choice(len(X))] else: new_centroids[c_idx] cluster_points.mean(axis0) if np.linalg.norm(new_centroids - centroids) tol: centroids new_centroids break centroids new_centroids return centroids, labels这个版本里distances变成(n, k)矩阵每个元素是第i个样本到第j个质心的距离argmin(axis1)就是按行找最小值得到n个样本的簇标签。对比教学版可以发现两个版本唯一的本质区别是距离矩阵的轴方向反了随之argmin的axis也反了。搞懂这一点你就同时看穿了大多数K-Means源码里让人迷惑的轴操作。数据量到几十万行、几十维特征时向量化版本比循环版快一到两个数量级这就是Numpy广播的威力。实际项目里如果还想再快可以用scipy.spatial.distance.cdist直接算距离矩阵它底层调了优化过的C实现比手写广播更快。不过在“手撕”阶段用广播更能看清形状变化的脉络因为每一步shape的变化就是算法逻辑的具象化。2.4 每轮迭代都在发生什么质心移动的可视化复盘代码跑通后最有意思的事情是看质心怎么一步步“归位”。为了画图我稍微改造一下函数把每轮迭代的质心都记录下来def kmeans_process(X, k, max_iter10, seed42): rng np.random.default_rng(seed) idx rng.choice(len(X), sizek, replaceFalse) centroids X[idx].copy() history [centroids.copy()] for i in range(max_iter): diff X[:, None, :] - centroids[None, :, :] distances np.sqrt((diff ** 2).sum(axis2)) labels np.argmin(distances, axis1) # 画图观察用的简化版省略了空簇处理逻辑 new_centroids np.array([X[labels c].mean(axis0) for c in range(k)]) history.append(new_centroids.copy()) if np.linalg.norm(new_centroids - centroids) 1e-4: centroids new_centroids break centroids new_centroids return np.array(history), labels, i 1 history, labels, iters kmeans_process(X, k3)然后用matplotlib画一个2x2的子图依次展示初始质心、第1轮后、第2轮后、最终结果。初始质心是从样本里随机挑的很可能三个点都挤在数据一侧第1轮分配完成后质心猛跳到三个簇的大致中心第2轮只做小修正第3轮基本要么触发收敛要么只需要再微调一步。我在实际跑图时观察到的典型现象是第一轮质心位移最大之后断崖式减小。这不是巧合而是坐标下降法的普遍特点——先快速找到大体区域再慢速精确逼近。理解这一点就知道为什么很多K-Means实现默认只迭代一两百轮就够了它真的不需要太多轮。如果发现收敛很慢根源几乎都不是“轮数不够”而是初始化差、数据没标准化或者K选得不合理。3. 手撕路上踩过的坑K-Means常见问题与排查实录3.1 空簇问题质心是怎么“消失”的第一个坑是空簇。现象是跑完后有一类标签没人用实际有效簇数从K变成了K-1。我第一次手撕时遇到这个还挺懵明明设置的K3为什么labels里只有0和1原因通常有三类。初始化质心位置太差比如三个初始质心全落在同一堆数据里另一个质心被孤立数据里有明显离群点离群点独自抢走一个质心或者K设得太大某些簇本来就只有零星几个样本随机初始化恰好没覆盖到。如果不处理空簇后果不只是少一个标签。求均值时用空簇样本取mean会直接报错或产生nan后面算SSE也全乱套。教科书版本往往回避这个细节但工程实现必须处理。我用的方案是检测到空簇就重新随机选一个样本作为该质心简单有效但不保证最优。更稳的做法有三种按我的推荐程度排列把新质心定为当前离其他质心最远的样本让空簇赶紧去覆盖空白区域。找到当前SSE最大的簇把它一分为二把新质心塞给空簇。直接用K-Means做初始化从源头上大幅降低空簇概率。我在小规模比赛和项目里用“最远样本填充”最多效果好、代码也就几行逻辑更接近生产版本。3.2 初始化玄学同数据两次结果为何天差地别第二个问题比空簇隐蔽得多同一份数据跑两遍K-Means结果可能完全对不上。这不是随机数种子的锅而是目标函数本身非凸。K-Means的SSE存在多个局部极小值。随机初始化相当于从某座山头出发往下走走到哪个谷底全看起点。有些谷底对应的聚类结果很离谱比如把两个本应分开的簇合并了或者把一个簇硬劈成两半。我第一次遇到时怀疑自己代码有bug查了一晚上才发现同样的代码换了个随机种子跑到另一个局部最优去了。工程界的解法早就标准化多跑几次取SSE最小的那次结果。sklearn里的n_init10就是干这个的。手撕版实现也很直接best_centroids None best_labels None best_sse np.inf for seed in range(10): centroids, labels kmeans_vec(X, k, seedseed) sse ((X - centroids[labels]) ** 2).sum() if sse best_sse: best_sse sse best_centroids centroids best_labels labels不过比“多次运行取最优”更优雅的方案是K-Means初始化。它的核心思路是让初始质心尽量离得远先随机选第一个质心之后每次以“距离平方占全场总距离平方的比例”作为概率挑选下一个质心。这个策略几乎把所有坏初始化压到了地板水平。手写K-Means初始化也不难def kmeans_plusplus_init(X, k, seed42): rng np.random.default_rng(seed) centroids [X[rng.choice(len(X))]] for _ in range(1, k): # 每个样本到最近已有质心的距离平方 diff X[:, None, :] - np.array(centroids)[None, :, :] d2 np.min((diff ** 2).sum(axis2), axis1) prob d2 / d2.sum() # 距离越远的样本越可能被选为质心 centroids.append(X[rng.choice(len(X), pprob)]) return np.array(centroids)这段代码的价值足够让一个手撕版K-Means的稳定性提升一个档次。我实际测试过原始随机初始化200个seed里大约有3到4次掉进明显差的局部解换K-Means之后同样200次试验一次都没翻车。这个差距落到业务里直接体现为“聚类结果跟业务case对得上”还是“结果看起来离谱”。3.3 收敛判定陷阱tol、max_iter与浮点精度的博弈还有一类坑不在算法逻辑而在边界条件。先说tol。tol控制的是“质心移动多少距离就认为收敛”。如果特征量级是几十万tol1e-4意味着质心几乎在原地踏步迭代会一直跑满max_iter才被迫停下如果特征已经标准化到零附近1e-4又相对合理。所以我建议跑K-Means前先标准化再根据质心位移曲线决定tol别直接照抄默认值。max_iter也有讲究。经典教材说“一般几十轮就收敛”但高维稀疏数据、初始化差时跑几百轮也不一定稳。我习惯设300到500。为什么不是一万因为每轮都要算一次全量距离矩阵数据大时多跑就是烧算力不如先把初始化做好。另一个容易被忽视的是浮点数精度。曾经我在判断收敛时把阈值写成1e-10结果跑满迭代次数也没触发收敛——因为质心在接近最优时每一步位移只有大约1e-7量级1e-10这个阈值过于苛刻。更麻烦的是当矩阵接近零时np.linalg.norm对极小位移会产生浮点抖动收敛判定不稳定。所以工程上我更倾向于用“簇分配标签不再变化”作为收敛信号。对K-Means来说标签不变意味着再怎么迭代结果也不会变这是最硬气的收敛判据if np.array_equal(labels, prev_labels): break当然常规情况下标签不变时质心也基本不动了两个判据等价。但标签判据对浮点噪声更鲁棒我后来的手撕代码里干脆同时允许两种退出条件代码量几乎没增加保险系数却高了不少。3.4 排查问题速查表现象、原因与处理方案把上面这些坑整理成一张表以后遇到问题直接对照现象典型原因我的处理方案聚类结果每次都不一样随机初始化陷入不同局部最优多次运行取最小SSE或用K-Means某个簇是空的初始化覆盖差、K太大、有离群点空簇填充最远样本或减小K或换K-Means迭代到max_iter才停tol相对量级太小、数据未标准化先标准化把tol调到1e-4量级聚类边界明显乱切数据簇形状非球、有重叠改用GMM、DBSCAN或先降维可视化结果对业务无解释性只看了统计最优没结合业务约束结合业务设定K人工复核每个簇如果遇到表里没有的怪问题一个笨但有效的排查方式把质心和簇内样本散点图画出来人眼扫一遍往往就有思路了。K-Means本身不复杂失败模式也不复杂大多数问题的根源不外乎“数据没洗干净”和“参数没对应上数据量级”这两处。4. 从手撕到实战标准化、K值选择与算法进阶路线4.1 跑之前先标准化欧氏距离的量纲陷阱K-Means基于欧氏距离而欧氏距离会被大数值特征主导。用客户分群举个例子特征里有“年龄”和“年消费额”年消费额动辄几万年龄只有几十距离计算几乎只看消费额年龄直接成了摆设。聚类结果看着像按消费额切了几刀实际上是量纲在帮忙做特征选择。标准做法通常是Z-score标准化把每个特征变成均值0、标准差1X_scaled (X - X.mean(axis0)) / X.std(axis0)还有一种Min-Max缩放把特征压到0到1之间对没有极端值的数据效果也不错。但Min-Max对离群点极敏感——一个极大的离群点会把其他所有数据压到0到0.01之间等于把正常样本全挤在一起。所以我个人对分布偏稳健的数据用Z-score更多。不过并不是所有场景都非要标准化。当所有特征本来就是同量纲时比如全是0-1的one-hot编码或全是百分比标准化反而可能破坏原始语义。判断标准很简单先问自己这个特征取值的绝对大小有没有业务意义如果只是量纲不同标准化如果量大本身代表某种强度保留原始尺度有时更合理。4.2 K值怎么定最靠谱肘部法、轮廓系数与实际业务约束K是K-Means里用户给定的唯一超参数也是最难拍板的一个。常见方法有两个半。第一个是肘部法。画一条“K值—SSE”曲线K越大SSE越小但下降速度会变缓。曲线上的那个拐点就是“肘”代表从那个K开始多分一簇带来的收益明显变小。手撕版实现很直接k_list range(2, 10) sse_list [] for k in k_list: centroids, labels kmeans_vec(X_scaled, k, seed42) sse ((X_scaled - centroids[labels]) ** 2).sum() sse_list.append(sse) plt.plot(k_list, sse_list, markero)第二个是轮廓系数。它同时考虑每个样本与自身簇内的紧密度以及到最近邻居簇的距离算出介于-1到1之间的分数全局平均轮廓系数越高说明簇结构越清晰。选K时一般取平均轮廓系数接近最大的那个值。不过轮廓系数在小数据集上非常吃计算量n超过几万就开始吃力那时更推荐用采样或者干脆靠业务经验。剩下半个方法是Gap Statistic它通过与随机数据对比来找K理论上更漂亮但计算量大实际项目里用得少。最后我还得说一句大实话数值指标永远只是参考真实业务里K经常由需求决定。你要做五个用户层级K就是5你要把客户分成能承载不同运营话术的几类K得看团队接不接受。最优K不等于最可用K这层判断只有人做得来。4.3 现成库到底多做了什么与sklearn KMeans的差异对照手撕完之后建议一定要做一次对照实验把结果跟sklearn的KMeans对比看看差多少。我测下来两者在常规数据上的SSE几乎一致差异主要来自几个工程细节sklearn默认初始化是k-means默认n_init10也就是跑10次完整聚类后取SSE最小的结果。它的默认tol1e-4、默认max_iter300处理空簇时用“最远点填充法”。底层是Cython实现还支持elkan变种利用三角不等式减少部分距离计算。我手撕版的“重新随机选点填充”方案跟它比粗糙一些但主流程完全对得上。如果你不想写代码Orange3这类可视化数据分析平台里也集成了聚类分析模块把K-Means封装成图形节点点几下鼠标就能跑。它的底层依然是“初始化—分配—更新—收敛”这套逻辑参数含义跟我上面写的完全一致。所以不管用什么工具理解这些参数都是通用的。数据规模上也要有数sklearn的KMeans在百万级样本、几十维特征上能跑得动但手撕版在同样规模下会因为全距离矩阵而内存爆炸。手撕版更像教学脚手架和问题排查工具真到生产级规模合理选择是Mini-Batch K-Means或者分布式实现。4.4 从K-Means走向哪K-Medoids、GMM高斯混合与更多聚类思路K-Means虽经典但有三个肉眼可见的短板均值对离群点敏感、只会生成球形簇、只做硬分配。对应这三个短板至少有三条进阶路线。第一条是K-Medoids。它不把质心定义为簇内均值而是从簇内样本里挑一个“最典型的代表”因此对离群点更稳。它适合特征不是连续值或均值没有意义的场景比如文本向量、评分数据。第二条是高斯混合模型GMM属于模型聚类算法。GMM给每个样本输出一个“属于各簇的概率向量”而不是一个硬标签同时允许每个簇有自己的协方差结构所以能拟合椭圆形的簇。从数学上看K-Means其实是GMM的特例——当GMM各组件协方差趋向于各向同性且共享时它的最大后验分配就退化为按欧氏距离最近分配也就是K-Means的行为。理解这一点你会发现K-Means处于整个概率聚类谱系的一端而不是孤立知识点。第三条是Mini-Batch K-Means每次只取一小批样本更新质心牺牲少量精度换速度千万级数据场景下我一般直接走这条线。还有一个方向值得提谱聚类、DBSCAN这些方法能处理K-Means完全做不了的形状比如嵌套的环形、细长的带状结构。选型逻辑很简单如果散点图画出来是几个圆团K-Means够了如果形状怪异就别硬扛先考虑密度聚类或图方法。这遍手撕最直接的收获是我终于敢在排查K-Means问题时说“我大概知道它哪里出了问题”而不是只能拍脑袋调参。另一个小技巧分享给刚开始学的人跑任何聚类之前先画一张降维后的散点图人眼扫一遍判断簇数目和形状。这一步花五分钟却能省掉后面无数瞎调参的时间。K-Means看起来简单但把它从“可运行的库”变成“自己写的代码”你才能真正感觉到它在做什么以及那些参数在什么情况下说了算。先把基础版跑稳后续我会接着写K-Means的细节和K值自动选择的实操对比感兴趣的话可以继续跟。