企业官网建设流程全解析

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 = Σ_{i=1}^{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-means++

K-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下降速率变缓的拐点，也就是所谓的“肘部”。

实际操作中，我一般会从K=1跑到K=8，记录每个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聚类结果（K=3）'); 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');

在合成数据上，这条曲线一般在K=3附近出现明显拐弯，后续K的WCSS下降变慢。如果你拿到的数据没有明显拐点，可以再结合轮廓系数曲线。操作方法是：对每个K值，保存对应的平均轮廓系数，画出K-轮廓系数的曲线，峰值所在位置就是推荐的K。我这里就不放完整循环了，思路就是把上面的脚本包一层循环，把sil存进数组。

4.3 在真实实验场景中的使用建议

实际使用这份代码时，我强烈建议把整个流程封装成一个函数，参数由上层脚本统一控制，方便批量跑实验。比如写一个script函数，输入数据集路径和K范围，输出每个K下的三个指标值。这样跑网格化实验的时候，不需要反复改主函数。

另外一点是，可视化部分我用的是gscatter，它是MATLAB专门用来按分组画散点图的函数，会自动分配颜色。如果你的数据维度超过2，可以先用PCA降到前两个主成分再画图，但注意此时坐标轴含义变成主成分分量，不要再把它当作原始特征解释。

5. 常见问题与避坑经验

5.1 典型问题速查表

下面这些是我在实际调试中遇到比较多的问题，整理成表格方便速查。

最后一个问题我在写代码时特别留意了，如果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存下来再取最优就行。这套代码我用了挺长时间，至今还没遇到比这个习惯更稳的整体方案。