K-means大概是所有做数据分析、模式识别方向的人最早接触的聚类算法我自己从实验室那会儿就开始反复用它。它的核心逻辑确实很直白把数据点分到最近的簇中心再根据分到的点重新计算簇中心一直迭代到稳定。但真正拿MATLAB从头实现的时候你会发现很多细节文档里不会提初始中心怎么选才不容易陷入局部最优、K到底设多少、不同量纲的数据怎么处理、聚类结果用什么指标评价才算靠谱。这次我把一份基于MATLAB的K-means全套实现整理成了代码仓库不调用内置的kmeans函数从迭代逻辑到WCSS、轮廓系数、Calinski-Harabasz指数和Davies-Bouldin指数全部手写一遍附带完整可视化脚本和调试过程中踩过的坑。如果你正在用MATLAB做聚类实验或者想彻底弄懂K-means的每个齿轮是怎么转的这份代码可以直接拿去跑也可以当模板改成自己的版本。我会先解释整体设计思路然后逐段拆解实现细节再放完整可运行的代码和实验演示最后整理一份常见问题速查表。整个过程按我自己平时写实验代码的习惯来组织代码风格偏可读性不追求炫技工作用完全够。1. 项目概述与整体思路1.1 为什么非要自己手写K-means我见过不少同学直接调MATLAB内置的kmeans函数一行代码出结果方便是真方便。但一旦需要改距离度量、换初始化方式、加新的评价指标或者想搞清楚算法内部哪里出了问题内置函数就像一个黑盒调试起来很痛苦。更重要的是写论文或者答辩的时候评审经常会问“你的初始中心怎么选的”“收敛条件是什么”“空簇怎么处理的”这些问题如果没手写过一遍很难答得深入。自己从零实现K-means能让我们完全控制每条细节。比如默认的欧氏距离可以换成余弦距离随机初始化可以换成K-means收敛判断可以用中心变化量也可以用WCSS变化量。这些改动在内置函数里要么不支持要么需要额外的参数配置远不如自己代码里改一行来得直接。所以这个项目定了一个原则基础功能全部手写实验功能全部开放。这样既能当教学示例也能作为各类聚类实验的起点代码。1.2 这份代码要完成的核心功能输入非常简单一份特征矩阵X、一个聚类数K、若干可选参数。输出包括聚类标签、簇中心、迭代过程的WCSS记录以及三个主流评价指标。具体来说实现标准K-means迭代流程包含向量化距离计算支持随机初始化和K-means初始化两种方式内置空簇处理机制和双重收敛判断计算WCSS用于肘部法则选K计算轮廓系数、Calinski-Harabasz指数、Davies-Bouldin指数提供二维数据可视化脚本方便直接观察聚类效果整体代码分三个文件聚类主函数、初始化辅助函数、评价指标计算函数。后面我会把每个文件的完整代码和设计原因都写清楚。2. 算法原理与实现细节2.1 K-means核心原理速览K-means做的事情说白了就是一句话让每个样本到它所属簇中心的距离平方和最小。这个目标函数叫WCSS公式是WCSS Σ_{i1}^{n} ||x_i - μ_{cluster(i)}||^2其中μ是簇中心cluster(i)是样本i被分配到的簇编号。算法通过交替执行两个步骤来最小化这个目标E步把每个样本分给最近的簇中心M步重新计算每个簇的中心为该簇所有样本的均值。可以用一个生活化的类比来理解想象K个队长站在一片场地上所有队员选择自己最近的队长走过去然后每个队长重新站到自己队员队伍的中间位置重复这个过程。当队长的位置几乎不再移动或者队员归属不再改变时聚类就完成了。K-means虽然简单但有三个问题几乎是绕不开的初始中心选得不好会陷入局部最优不同特征量纲差异大会让距离计算失衡K值需要人为指定。这套代码里前两个问题都有对应处理方案K值选择会用肘部法则和评价指标一起判断。2.2 数据预处理标准化不是可选项做K-means最容易忽略的一步是标准化。如果特征A的取值范围是0到1特征B的取值范围是10000到50000那么计算欧氏距离时特征B的数值会完全主导结果特征A的信息几乎等于没参与。聚类结果自然也就被量纲大的特征带偏了。MATLAB里做标准化很简单调用zscore矩阵就行X zscore(X);标准化之后每个特征变成均值为0、标准差为1这样每个特征对距离的贡献才在同一个量级上。但这里有一个很多新手会踩的坑标准化用的均值方差必须在训练集上计算后面如果有新样本加入要用同一组参数来转换而不是重新计算zscore。否则新样本的转换方式和训练样本不一致聚类结果就不可比了。我在写这个项目时特意在测试脚本里先做了zscore然后才调用聚类函数。如果你把标准化环节丢给样本集自己去做每次跑出来的结果都可能因为样本组合不同而漂移。2.3 初始化策略随机选和K-meansK-means对初始中心非常敏感。随机初始化如果踩到某几个点都集中在同一个簇内部剩下的簇可能永远抢不到样本最后得到很差的局部最优。而且随机初始化本身有随机性同一份数据跑两次结果可能差异很大。一个常用的改进是K-means初始化思路很直接第一个中心随机选之后每个新中心都优先选择离已选中心更远的样本点。这样让初始中心在数据空间里尽量铺开从源头减少局部最优的概率。整个过程里的概率采样很关键。我们不是简单选“最远”的点而是按距离平方的占比作为概率分布做加权随机。这样既保证了远离已选中心的样本更容易被选中又保留了一定随机性不会每次都得到完全相同的结果。我在代码里实现了这个逻辑后面会放完整函数。2.4 从零实现的MATLAB聚类主函数下面这版函数是实现K-means的主体关键逻辑都写了注释。我先给完整代码再逐段解释。function [label, centroid, allWCSS] myKMeans(X, K, opts) % 从零实现K-means聚类算法 % 输入: % X : n x m 矩阵行为样本列为特征 % K : 聚类数量 % opts : 可选参数结构体包含 maxIter, tol, init % 输出: % label : n x 1 簇标签向量取值1到K % centroid : K x m 最终簇中心 % allWCSS : 每次迭代的WCSS值用于观察收敛过程 if nargin 3 opts struct(); end if ~isfield(opts, maxIter), opts.maxIter 200; end if ~isfield(opts, tol), opts.tol 1e-4; end if ~isfield(opts, init), opts.init kmeans; end n size(X, 1); % 初始化簇中心 centroid zeros(K, size(X, 2)); switch lower(opts.init) case random idx randperm(n, K); centroid X(idx, :); case kmeans centroid kmeansPlusPlusInit(X, K); otherwise error(不支持的初始化方式: random 或 kmeans); end allWCSS zeros(opts.maxIter, 1); for iter 1:opts.maxIter % 计算每个样本到所有簇中心的欧氏距离 D pdist2(X, centroid); % 分配样本到最近的簇 [~, label] min(D, [], 2); % 更新簇中心 for k 1:K memberIdx (label k); if any(memberIdx) centroid(k, :) mean(X(memberIdx, :), 1); else % 空簇处理把离该簇原中心最远的样本重新分配进来 [~, farIdx] max(min(D, [], 2)); centroid(k, :) X(farIdx, :); label(farIdx) k; end end % 重新计算WCSS用本轮分配结果 D2 pdist2(X, centroid); distToCenter min(D2, [], 2); allWCSS(iter) sum(distToCenter.^2); % 收敛判断WCSS变化小于阈值则提前终止 if iter 1 abs(allWCSS(iter) - allWCSS(iter-1)) opts.tol allWCSS allWCSS(1:iter); break; end end if iter opts.maxIter abs(allWCSS(end) - allWCSS(end-1)) opts.tol fprintf(达到最大迭代次数 %d可能未完全收敛。\n, opts.maxIter); end end代码里有两个关键设计。第一距离计算用了pdist2函数它会一次性返回n行K列的距离矩阵然后min操作直接找到每个样本最近的簇。这种写法比双重for循环快很多尤其样本量上万时差异非常明显。第二WCSS的计算在簇中心更新之后重新做了一次因为分配用的旧中心和新中心可能有细微差异用新中心计算的WCSS才是这一轮真正意义上的目标函数值。收敛判断这里我做了双重保险当WCSS变化量小于1e-4时提前退出如果达到maxIter还没收敛会打印提示这样至少你知道结果可能是未完全收敛的。在实际实验里maxIter设200通常足够但如果你数据量大或者中心初始得很差可能要多跑几次。2.5 K-means初始化函数实现/ 初始化部分单独抽出来写方便后期改成其他初始化策略。下面是完整实现function centroids kmeansPlusPlusInit(X, K) % K-means 初始化让初始中心之间尽量互相远离 % 输入输出与 myKMeans 中的说明一致 n size(X, 1); centroids zeros(K, size(X, 2)); % 随机选择第一个中心 firstIdx randi(n); centroids(1, :) X(firstIdx, :); % D2记录每个样本到最近已选中心的距离平方 D2 pdist2(X, centroids(1, :)).^2; for k 2:K sumD2 sum(D2); if sumD2 0 % 极端情况所有样本都距离已选中心为0直接随机选一个未选样本 chosenIdx randi(n); while any(all(centroids(1:k-1, :) repmat(X(chosenIdx, :), k-1, 1), 2)) chosenIdx randi(n); end centroids(k, :) X(chosenIdx, :); continue; end % 按距离平方占比作为概率分布采样 probs D2 / sumD2; cumProbs cumsum(probs); r rand(); chosenIdx find(cumProbs r, 1, first); centroids(k, :) X(chosenIdx, :); % 更新每个样本到已选中心的最短距离 D2 min(D2, pdist2(X, centroids(k, :)).^2); end end这段代码的精华在最后两行。每选完一个新中心我们都要重新计算每个样本到“所有已选中心”的最短距离平方。用min(D2, 新距离平方)可以原地更新表达的是“我之前到最近中心的距离和到新中心的距离哪个更小”。这个更新策略让K-means的计算代价保持在可接受范围。2.6 空簇处理与收敛判断的细节实际聚类过程中空簇并不少见尤其是K设得比较大、数据本身又分布不均的时候。空簇如果不管后面的mean计算会出现NaN整个结果就毁了。我这里的处理方式其实很粗暴从当前所有样本里找一个距离该空簇原中心最远的样本把它当成新中心重新开始。这样做的好处是不引入过多额外逻辑坏处是这个样本可能已经被其他簇分配了所以紧接着还专门把label修正了一下。更规范一些的做法是重新用K-means在这个簇的附近重新采样或者把小簇合并到邻近簇里。我在这里选择了最直接的方案如果空簇出现频率高我建议优先检查初始化时的K是不是设得太大而不是依赖空簇补救逻辑。收敛判断这里我也要做个说明我只卡了WCSS的变化量没有同时检查样本分配是否变化。其实两种方式都行但WCSS变化能更平滑地反映整体收敛状态样本分配可能在末尾阶段来回抖动WCSS却已经稳定了。实际项目中双条件判断也可以代码只是提供了一种实现基线。3. 聚类评价指标计算3.1 WCSS与肘部法则选K聚类做完不能光看图片觉得“好像挺对的”量化指标才有说服力。最基础也最常用的指标是WCSS也就是组内平方和。它衡量的是所有样本到自己所属簇中心的距离平方总和这个值越小说明簇内越紧凑。但WCSS有个特点K越大WCSS通常越小。当K等于样本数时每个点自己就是一个簇WCSS直接变成0。所以不能用WCSS最小作为选K的依据而是要找WCSS下降速率变缓的拐点也就是所谓的“肘部”。实际操作中我一般会从K1跑到K8记录每个K对应的WCSS然后画一条折线图。如果数据有明显簇结构曲线会在某个位置出现明显拐弯那个位置的K就是比较合理的选择。不过现实数据常常拐点不明显所以还需要结合其他指标做综合判断。3.2 轮廓系数直观衡量簇内外距离轮廓系数的思想很直观。对每个样本i先计算它到同簇其他样本的平均距离a(i)再计算它到最近其他簇所有样本的平均距离b(i)。样本的轮廓系数定义为s(i) (b(i) - a(i)) / max(a(i), b(i))结果落在-1到1之间。接近1表示样本被正确聚类接近-1表示样本更像属于隔壁簇接近0表示样本正好在两个簇的边界附近。整体轮廓系数把所有样本的s(i)取平均用来衡量整个聚类结果的质量值越大越好。轮廓系数的计算复杂度是O(n^2)因为需要遍历所有样本到所有样本的距离。我在代码里先用pdist2一次性算好了全距离矩阵再按标签批量取平均。如果样本量超过几万内存会开始紧张那就要改为只计算簇内和近邻簇的局部距离这部分性能优化我后面会在经验部分展开讲。3.3 Calinski-Harabasz指数组间分离与组内紧凑的比值Calinski-Harabasz指数简称CH指数计算的是组间离散度与组内离散度的比值。公式为CH [BSS / (K - 1)] / [WSS / (n - K)]其中BSS是组间平方和衡量各簇中心到全局中心的加权距离WSS是组内平方和和WCSS一样。CH值越大说明组间分离越明显、组内越紧凑聚类效果越好。注意它不像轮廓系数有上界所以一般用来比较不同K下的结果而不是绝对判断好坏。CH指数有个非常好的特性它不需要计算所有样本之间的两两距离只涉及每个簇中心和各簇内部方差所以计算速度非常快。我做实验时通常会把它和轮廓系数结合使用一个快一个信息量丰富互相印证。3.4 Davies-Bouldin指数簇内散布与簇间距离Davies-Bouldin指数也叫DB指数思路是对每一对簇计算一个“坏程度”比值R_ij (S_i S_j) / d(μ_i, μ_j)其中S_i是簇i的样本到其中心的平均距离相当于簇内半径d(μ_i, μ_j)是两个簇中心之间的距离。每一个簇i会选择所有j≠i中R_ij最大的那个作为它的最坏情况最终DB指数是所有簇最坏情况的平均。DB指数越小说明每个簇内部越紧凑、簇之间分开得越干净。DB指数和CH指数互补CH侧重组间方差占比DB直接考察簇内散布和中心距离的比值。三者一起看基本能覆盖聚类评价的大部分需求。极端情况比如某两个簇中心完全重合时DB会变成无穷大代码里会返回一个大数值提醒你聚类失效了。3.5 评价指标函数完整实现我把三个指标封装到一个函数里一次性返回结果这样在跑不同K时可以快速对比。完整代码如下function [sil, ch, db] clusterMetrics(X, label, centroid) % 计算三个聚类评价指标 % 输入: % X : n x m样本矩阵 % label : n x 1簇标签 % centroid: K x m簇中心 % 输出: % sil : 平均轮廓系数 % ch : Calinski-Harabasz指数 % db : Davies-Bouldin指数 n size(X, 1); K size(centroid, 1); if n 2 || K 2 error(样本数或聚类数过少无法计算评价指标); end % 全样本距离矩阵便于轮廓系数计算 D pdist2(X, X); % 1. 轮廓系数 a zeros(n, 1); % 到同簇样本平均距离 b zeros(n, 1); % 到最近其他簇样本平均距离 for i 1:n sameIdx find(label label(i)); if length(sameIdx) 1 a(i) 0; else a(i) sum(D(i, sameIdx)) / (length(sameIdx) - 1); end otherLabels unique(label); otherLabels otherLabels(otherLabels ~ label(i)); bestB inf; for k 1:length(otherLabels) otherIdx find(label otherLabels(k)); avgDist sum(D(i, otherIdx)) / length(otherIdx); if avgDist bestB bestB avgDist; end end b(i) bestB; end s (b - a) ./ max(a, b); s s(isfinite(s)); % 处理边界情况 sil mean(s); % 2. Calinski-Harabasz指数 globalMean mean(X, 1); BSS 0; WSS 0; for k 1:K idx find(label k); nk length(idx); if nk 0, continue; end clusterMean mean(X(idx, :), 1); BSS BSS nk * sum((clusterMean - globalMean).^2); WSS WSS sum(sum((X(idx, :) - clusterMean).^2)); end if WSS 0 ch inf; else ch (BSS / (K - 1)) / (WSS / (n - K)); end % 3. Davies-Bouldin指数 S zeros(K, 1); for k 1:K idx find(label k); if isempty(idx) S(k) 0; else distToCenter sqrt(sum((X(idx, :) - repmat(centroid(k, :), length(idx), 1)).^2, 2)); S(k) mean(distToCenter); end end M zeros(K, K); for i 1:K for j 1:K if i j M(i, j) 0; else M(i, j) sqrt(sum((centroid(i, :) - centroid(j, :)).^2)); end end end R zeros(K, 1); for i 1:K maxR 0; for j 1:K if j i, continue; end if M(i, j) 0 tmpR inf; else tmpR (S(i) S(j)) / M(i, j); end if tmpR maxR maxR tmpR; end end R(i) maxR; end db mean(R); end这段代码里轮廓系数的部分做了isinf过滤是因为极端情况下如果所有样本都被分到同一个簇b会变成inf。过滤掉非有限值再去平均至少不会让整个指标变成NaN。CH和DB两个指标内部用的都是向量化或简单循环计算效率足够应对一般实验规模。4. 完整代码与实操演示4.1 测试数据生成与标准流程脚本我准备了一段演示脚本它生成三个有明显分离趋势的高斯簇然后依次完成标准化、K-means聚类、指标计算和可视化。整段代码可以直接复制运行得到的结果和我在下文展示的基本一致。rng(42); % 固定随机种子保证可复现 % 生成三簇二维数据 X1 randn(150, 2) repmat([0, 0], 150, 1); X2 randn(150, 2) repmat([5, 2], 150, 1); X3 randn(150, 2) repmat([-3, 4], 150, 1); X [X1; X2; X3]; % 标准化 X zscore(X); % 设置参数并运行K-means K 3; opts struct(maxIter, 200, tol, 1e-4, init, kmeans); [label, centroid, allWCSS] myKMeans(X, K, opts); % 计算评价指标 [sil, ch, db] clusterMetrics(X, label, centroid); fprintf(最终WCSS: %.4f\n, allWCSS(end)); fprintf(平均轮廓系数: %.4f\n, sil); fprintf(Calinski-Harabasz指数: %.4f\n, ch); fprintf(Davies-Bouldin指数: %.4f\n, db); % 可视化聚类结果 gscatter(X(:,1), X(:,2), label); hold on; plot(centroid(:,1), centroid(:,2), kx, MarkerSize, 12, LineWidth, 2); legend off; title(K-means聚类结果K3); hold off; % 画WCSS收敛过程 figure; plot(1:length(allWCSS), allWCSS, o-); xlabel(迭代次数); ylabel(WCSS); title(WCSS收敛曲线);运行这段脚本后会看到聚类图三个簇基本被正确分开簇中心也落在数据密集区域。输出的指标值大概是WCSS在几百到一千的区间轮廓系数在0.6到0.75之间CH指数在几百到上千DB指数在0.5左右。每次运行因为K-means引入的随机性会有小幅波动但整体趋势稳定。4.2 肘部法则选K的完整流程很多场景下K不是预设好的需要数据告诉我们该选多少。这段代码跑K从1到8的数据绘制WCSS曲线。wcssValues zeros(1, 8); for K 1:8 [~, ~, allWCSS] myKMeans(X, K, opts); wcssValues(K) allWCSS(end); end figure; plot(1:8, wcssValues, o-); xlabel(K值); ylabel(WCSS); title(肘部法则选K);在合成数据上这条曲线一般在K3附近出现明显拐弯后续K的WCSS下降变慢。如果你拿到的数据没有明显拐点可以再结合轮廓系数曲线。操作方法是对每个K值保存对应的平均轮廓系数画出K-轮廓系数的曲线峰值所在位置就是推荐的K。我这里就不放完整循环了思路就是把上面的脚本包一层循环把sil存进数组。4.3 在真实实验场景中的使用建议实际使用这份代码时我强烈建议把整个流程封装成一个函数参数由上层脚本统一控制方便批量跑实验。比如写一个script函数输入数据集路径和K范围输出每个K下的三个指标值。这样跑网格化实验的时候不需要反复改主函数。另外一点是可视化部分我用的是gscatter它是MATLAB专门用来按分组画散点图的函数会自动分配颜色。如果你的数据维度超过2可以先用PCA降到前两个主成分再画图但注意此时坐标轴含义变成主成分分量不要再把它当作原始特征解释。5. 常见问题与避坑经验5.1 典型问题速查表下面这些是我在实际调试中遇到比较多的问题整理成表格方便速查。问题现象根本原因解决办法每次聚类结果不一样随机初始化带来的随机性固定随机种子rng或用K-means并多跑几次取最优结果明显很差簇混在一起初始中心选到了局部最优改用kmeans或增加重复次数出现空簇某些中心没有样本K设太大或初始化不好降低K或使用代码里的空簇重分配逻辑运行速度慢几万样本跑不动距离计算用了双层for循环改用pdist2对全部中心计算距离矩阵轮廓系数或DB指数出现NaN某个簇只有一个样本检查空簇和极端K值去掉单点簇WCSS不降反升收敛判断条件写错了确认中心更新后重新计算WCSS再比较K-means偶尔选中重复样本概率采样逻辑有边界问题检查D2更新是否用了min操作最后一个问题我在写代码时特别留意了如果D2更新时用的是赋值操作而不是min操作就会覆盖掉到其他新中心的更小距离导致采样偏向重复选点。5.2 亲手改代码比调用内置函数更有价值我有几个实操心得想单独拿出来说。第一刚拿到代码时不要急着跑大数据先拿二维合成数据验证逻辑是否正确画出来的图能直接告诉你聚类是否合理。第二多跑几次不同的K把WCSS、轮廓系数、CH指数、DB指数记录到一个表格里你会发现在不同K下这些指标波动的方式其实很有规律这对理解它们的意义非常有帮助。另外一个很实用的技巧是调试时可以临时去掉标准化步骤看看聚类结果飘成什么样。这个对比能直观地让你理解标准化为什么是必需步骤。我见过不少同学在真实数据集上跑了半天结果聚类出来的簇完全是按某个大数值特征切的而小数值特征基本没参与一查发现忘了zscore。5.3 扩展方向与可复用模块这套代码最值得留下的资产不是K-means本身而是那个名称为clusterMetrics的评价指标模块。换个聚类算法比如层次聚类或DBSCAN只要你能拿到标签和簇中心就能直接复用这套指标计算逻辑。把评价指标和算法解耦是我后来做大量对比实验时最受益的设计习惯。另外如果你做的是图像聚类或者文本聚类距离度量可能需要切换。代码里的pdist2只算欧氏距离但改成余弦距离也很容易先把每个样本归一化成单位向量再算欧氏距离就能得到余弦相似度。这种做法在向量检索里应用很广K-means的迭代流程完全不用动。最后说一个我常年保留的习惯跑正式实验时先用K-means初始化同一数据跑10次记录每次的WCSS选最小WCSS那次作为最终结果。这样做虽然多花了一点时间但结果稳定很多写报告、答辩的时候也不容易被追问复现性的问题。代码层面只需要在外层包一个循环把每次的WCSS存下来再取最优就行。这套代码我用了挺长时间至今还没遇到比这个习惯更稳的整体方案。