企业官网建设流程全解析

做聚类分析的朋友一定有过这种经历：课上老师讲了K-means的原理，作业要求自己写一遍，但翻开任何一本教材，公式都能看懂，落到MATLAB代码上却总觉得隔了一层。等你好不容易调通了聚类，要写实验报告或者论文时，又卡在评价指标上——轮廓系数怎么算？CH指数到底在衡量什么？为什么这块的结果和别人论文里的对不上？

这篇博文就围绕“在MATLAB里从零实现K-means算法，并把常用的聚类评价指标一并计算出来”这件事展开，给出可以直接运行的代码、关键参数的取舍逻辑、以及我在实际项目中踩过的一些坑。适合正在做课程设计、准备算法竞赛、或者需要在科研里快速验证聚类效果的同学；哪怕你只是想把散落各处的代码整合成一套能复用的分析工具，这篇文章也值得花几分钟读完。

1. 为什么我坚持在MATLAB里手写K-means，而不是只调内置函数

MATLAB自带的kmeans函数非常成熟，参数丰富，甚至支持多种距离度量、重复采样和并行计算。但我在带项目、审阅学生作业的过程中发现一个问题：大量使用者只会填一个矩阵和一个K值，得到标签之后直接画图，一旦问“你的初始中心怎么选的？收敛条件是什么？空簇怎么处理？”就答不上来。这在课程设计里或许能蒙混过关，但到了实际业务场景，数据一脏、分布一复杂，内置函数返回的结果就可能悄悄“出错”，而你完全没有意识。

手写K-means最大的价值不是重复造轮子，而是把算法的假设、边界和代价全部暴露在你能看到的地方。举例来说，K-means隐含的核心假设是“簇呈凸形、密度相近”，它优化的是每个样本到所属簇中心的欧氏距离平方和。用生活化的类比：一群人要分成几个小组，每个小组选一个“临时队长”，规则是“每个人加入离自己最近的队长那一队，然后队长改成自己队伍所有人的平均位置，再重新分队”，如此反复直到队长位置不再变化。这个过程里有三个决策点极易出错：初始队长怎么选、队员到队长的“距离”用什么定义、“不再变化”的判定阈值是多少。

自己实现一遍之后，这些决策点就变成了你代码里的参数和分支。而且MATLAB的矩阵运算天然适合K-means的向量化实现，写出来并不比调库慢多少。更关键的是，你需要计算轮廓系数、CH指数、DB指数、Dunn指数这些评价指标时，MathWorks虽然在某些工具箱里有零散函数，但它们往往要求特定的输入格式，不同版本命名还不一样——与其依赖它们，不如基于距离矩阵自己写一套，所有中间量都能打印出来检查，出了问题也容易定位。

出于这些考虑，我在自己的代码仓库里保留了一套自写的K-means实现，并用它作为评价指标计算的基础。下面从主程序开始逐部分拆解。

2. 从零搭建K-means主程序：输入、迭代与收敛控制

我的myKmeans函数设计目标是“足够透明但足够高效”，输入一个n×d的特征矩阵和聚类数K，输出标签向量、最终簇中心、每次迭代的目标函数值（SSE），方便后续绘制收敛曲线。

function [labels, centroids, sseHistory] = myKmeans(X, K, opts) % myKmeans - 自实现K均值聚类 % 输入: % X : n×d 数值矩阵，每一行是一个样本 % K : 聚类数量 % opts : 可选结构体，包含： % MaxIter 最大迭代次数（默认200） % Tol 中心点移动距离阈值（默认1e-4） % InitMethod 初始化方法 'random' 或 'kmeanspp'（默认'random'） % Distance 'euclidean' 或 'cosine'（默认'euclidean'） % 输出: % labels : n×1 标签向量，取值1..K % centroids : K×d 最终簇中心 % sseHistory : 每次迭代的各类内误差平方和之和 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, 'InitMethod'), opts.InitMethod = 'random'; end if ~isfield(opts, 'Distance'), opts.Distance = 'euclidean'; end [n, d] = size(X); if K > n error('K不能大于样本数量'); end % 初始化簇中心 if strcmp(opts.InitMethod, 'kmeanspp') centroids = initKmeansPP(X, K); else % 随机选K个不同样本作为初始中心 rng('shuffle'); % 每次运行结果不同；复现实验时可固定种子 idx = randperm(n, K); centroids = X(idx, :); end sseHistory = zeros(1, opts.MaxIter); labels = zeros(n, 1); for iter = 1:opts.MaxIter % 分配步骤：计算每个样本到所有中心的距离 distMatrix = pdist2(X, centroids, opts.Distance); % n×K [minDist, labels] = min(distMatrix, [], 2); % 更新步骤：计算每个簇的均值作为新中心 newCentroids = zeros(K, d); for k = 1:K clusterPoints = X(labels == k, :); if ~isempty(clusterPoints) newCentroids(k, :) = mean(clusterPoints, 1); else % 空簇处理：常见做法是保留原中心，或者重新随机初始化 % 我这里保留原中心，并在命令行给出提示 fprintf('警告: 第 %d 轮迭代中簇 %d 为空，保留上一轮中心\n', iter, k); newCentroids(k, :) = centroids(k, :); end end % 计算本轮SSE sseHistory(iter) = sum(minDist.^2); % 收敛判定：所有中心点的最大位移小于阈值 shift = max(sqrt(sum((newCentroids - centroids).^2, 2))); centroids = newCentroids; if shift < opts.Tol sseHistory = sseHistory(1:iter); break; end end end

这段代码有几个值得展开的设计细节。

第一，用pdist2一次性算所有样本到所有中心的距离，在MATLAB里比双重循环快一个数量级。如果数据量超过十万行，再考虑分块计算避免内存暴涨；这个优化后面在“大样本处理”一节再说。第二，收敛判定我用的是中心点最大位移小于阈值，比“标签不再变化”更灵敏。实际操作中，经常出现中心还在缓慢漂移但标签已经不变的情况，此时提前终止会让SSE还没压到最低；用位移阈值更接近理论收敛点。第三，空簇处理是初学者最容易忽略的。某些中心初始化得太差，或者数据分布特殊，迭代中会出现某个簇没有分配到任何样本，此时mean一个空集全是NaN，NaN会像病毒一样传染到后面的距离矩阵。我选择保留上一轮中心，同时在命令窗口打警告，保证算法不会崩溃，并且让你知道自己数据里存在这样的不稳定现象。

调用示例我习惯用鸢尾花数据集，因为它是MATLAB内置的、带真实类别标签的经典数据，非常适合验证：

% 加载数据并做z-score标准化 load fisheriris X = meas; % 150×4 X = zscore(X); % 标准化，防止量纲影响距离 K = 3; opts.InitMethod = 'kmeanspp'; opts.MaxIter = 100; [labels, centroids, sse] = myKmeans(X, K, opts); % 输出聚类结果与真实类别的交叉矩阵 trueLabel = species; tabulate(labels)

标准化这一步很重要却常被忽略。鸢尾花数据里PetalLength的标准差远大于SepalLength，如果不标准化，距离计算会被花瓣长度主导，萼片的信息几乎不起作用。类似地，任何量纲不一致的特征，在做欧氏距离聚类前都应该z-score，这是K-means使用中性价比最高的一条经验。

3. 四种聚类评价指标的MATLAB实现与选择逻辑

聚类评价指标有两个方向：有真实标签时用外部指标（比如纯度、兰德指数、NMI），没有真实标签时用内部指标。实际业务中大部分情况是没有标签的，所以内部指标价值更大。内部指标的共同思路是“绕开真实答案，只看聚类结果的紧密程度和分离程度”。我挑了四个最常用、也最能互相补充的指标，逐一给出MATLAB实现。

3.1 轮廓系数：单个样本视角的“聚得紧不紧、分得开不开”

轮廓系数对每个样本计算一个值，取值在-1到1之间。记样本(i)所属簇为(A)，该样本到簇内所有其他样本的平均距离记为(a(i))；到最近的其他簇(C)内所有样本的平均距离记为(b(i))。则：

[ s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))} ]

直观理解：如果(s(i))接近1，说明它离自己簇内其他成员很近、离邻居簇很远，聚类效果好；接近-1说明它更可能被分在了错误的簇；接近0说明它处于两个簇的边界。所有样本的(s(i))取平均，就是全局轮廓系数。

MATLAB实现时注意，最朴素的双重三重循环在样本上万时会非常慢，我建议先算一个完整的n×n距离矩阵，然后利用labels构造同簇/异簇的掩膜矩阵，用矩阵运算把系数批量求出来：

function S = silhouetteCoefficient(X, labels) n = size(X, 1); D = pdist2(X, X); % 对称距离矩阵 D(1:n+1:end) = NaN; % 对角线置NaN，忽略自身距离 S = zeros(n, 1); for i = 1:n same = (labels == labels(i)); same(i) = false; a_i = mean(D(i, same), 'omitnan'); otherClusters = unique(labels(labels ~= labels(i))); b_i = inf; for c = otherClusters' mask = (labels == c); b_i = min(b_i, mean(D(i, mask), 'omitnan')); end S(i) = (b_i - a_i) / max(a_i, b_i); end end

这里用了'omitnan'参数确保对角线上的NaN不干扰均值。外层循环走一遍所有样本，内层对每个样本遍历其他簇，整体复杂度(O(n^2 + nK))，在几千样本内跑起来很顺畅。如果你的样本量到了几十万，这个指标就不太现实了，可以考虑随机抽样子集计算，或者换用下面的CH指数。

3.2 CH指数（Calinski-Harabasz）：类内紧、类间散的方差比

CH指数的计算公式为：

[ CH = \frac{\text{trace}(B)/(K-1)}{\text{trace}(W)/(n-K)} ]

其中(B)是簇间离差矩阵，(W)是簇内离差矩阵，trace是矩阵的迹。说得直白点，分子衡量“各簇中心到全局中心的距离平方和”，值越大说明簇间分得越开；分母衡量“每个样本到本簇中心的距离平方和”，值越小说明簇内越紧凑。两者一除，再加样本量与K的自由度修正，就得到CH指数。这个指标计算量很轻，只需要簇中心和每簇内部距离，不需要所有样本间的两两距离，因此特别适合大样本。

function ch = calinskiHarabaszIndex(X, labels, centroids) n = size(X, 1); K = size(centroids, 1); globalMean = mean(X, 1); % 簇间离差矩阵和 B = 0; for k = 1:K nk = sum(labels == k); centroidDiff = centroids(k, :) - globalMean; B = B + nk * (centroidDiff * centroidDiff'); end % 簇内离差平方和 W = 0; for k = 1:K points = X(labels == k, :); diff = points - centroids(k, :); W = W + sum(sum(diff .* diff, 2)); end ch = (trace(B) / (K - 1)) / (trace(W) / (n - K)); end

注意trace(B)这里因为centroidDiff是行向量，centroidDiff * centroidDiff'是一个(d×d)矩阵，取迹就得到该簇中心到全局中心的欧氏距离平方再乘以样本数。这个指标的优点是不依赖全距离矩阵，而且随着K增大会自然给出惩罚——K太多时每个簇越来越小，分子虽然也可能增大，但分母的自由度惩罚会让CH增速放缓或下降，适合用来选K。

3.3 DB指数（Davies-Bouldin）：簇内分散度与簇间距离的最坏比值

DB指数的思路非常工程化：对任意两个簇(i)和(j)，定义

[ R_{ij} = \frac{S_i + S_j}{M_{ij}} ]

其中(S_i)是簇(i)内样本到中心(i)的平均距离（簇内分散度），(M_{ij})是两个簇中心的距离（簇间分离度）。DB指数就是对每个簇找最坏的（最大的）(R_{ij})，再对所有簇求平均：

[ DB = \frac{1}{K} \sum_{i=1}^{K} \max_{j \neq i} R_{ij} ]

DB指数越小越好。直观理解就是每个簇都得保证“自己别太松散，同时离最近的那个邻居也不要太近”。它比CH指数更强调局部最坏情况，不会因为某个簇分得特别干净就掩盖另一个簇的交叉纠缠。

function db = daviesBouldinIndex(X, labels, centroids) K = size(centroids, 1); S = zeros(K, 1); for k = 1:K points = X(labels == k, :); diff = points - centroids(k, :); S(k) = mean(sqrt(sum(diff .* diff, 2))); end R = zeros(K); for i = 1:K for j = 1:K if i ~= j R(i, j) = (S(i) + S(j)) / norm(centroids(i, :) - centroids(j, :)); end end end db = mean(max(R, [], 2)); end

如果数据本身存在严重离群点，某个簇的(S_i)会非常大，导致所有跟它有关的(R_{ij})都偏大，DB指数整体被抬高。这其实是设计意图——DB指数对噪声比较敏感，能帮你发现聚类结果里的离群样本。

3.4 Dunn指数：最坏情况下也能保持分离

Dunn指数定义稍微复杂一些，它同时关注“簇之间的最小距离”和“簇内部的直径”：

[ D = \frac{\min_{1 \le i < j \le K} d(C_i, C_j)}{\max_{1 \le l \le K} \text{diam}(C_l)} ]

分子是所有不同簇间最小的“簇间距离”，分母是所有簇中最大的簇直径（簇内最远的两个样本之间的距离）。Dunn指数越大越好，因为它要求最坏情况下的簇间距离也要大于最大簇直径。这个指标的理论意义很强——它对应聚类有效性分析里的“分离且紧凑”理想，但是计算复杂度很高，需要每个簇内所有样本间的距离，样本量一大就不太现实。

function dunn = dunnIndex(X, labels) K = max(labels); n = size(X, 1); D = pdist2(X, X); D(1:n+1:end) = NaN; minInter = inf; for i = 1:K for j = i+1:K maskI = (labels == i); maskJ = (labels == j); subDist = D(maskI, maskJ); minInter = min(minInter, min(subDist(:))); end end maxDiam = 0; for k = 1:K mask = (labels == k); subDist = D(mask, mask); subDist(1:size(subDist,1)+1:end) = NaN; maxDiam = max(maxDiam, max(subDist(:), [], 'omitnan')); end dunn = minInter / maxDiam; end

实际使用中Dunn指数波动大，经常出现几个K值下数值都很接近的情况，所以我一般不单独拿它选K，而是把它当作辅助判断：如果某个聚类方案的Dunn显著高于其他方案，那说明这个方案在最坏情况下（最紧凑的簇距离最近的邻居）仍然分离良好，可信度很高。

3.5 指标之间的互补关系

下表是我在实际工作中总结的四指标适用场景：

单看任何一个指标都有可能被骗。举个常见现象：轮廓系数在K=5时最大，但CH指数在K=3时最大。这时候不要急着相信某个指标，而是回到业务本身——如果任务本来只是想把数据分成大方向上的几类，SME（领域专家）能解释K=3的每个簇，K=5的某些簇是同一个簇硬拆出来的，那就该选K=3。指标是参考坐标系，不是真理。

4. 从评价指标倒推最优K值：肘部法则与指标曲线的实操细节

聚类分析里最痛苦的问题不是怎么写K-means，而是K到底取几。我见过很多人直接“猜一个”，画个散点图自圆其说，这在课程报告里也许能混过去，在论文审稿人面前就站不住脚了。我平时选K有一套固定流程，跑一遍所有候选K对应的指标，再综合判断。

4.1 肘部法则与SSE曲线

肘部法则看的是“SSE随K增大的下降曲线”。K越大，每个簇越可以变得更紧密，SSE必然下降；但下降到某个拐点后，再增加K带来的SSE收益会显著变小，这个拐点就是“肘部”。比如K从1到2降幅很大，2到3降幅也大，3到4开始趋缓，那K=3就是合理的。问题是很多真实数据的曲线根本没有明显拐点，这时候需要结合其他指标来定。

% 对K=1..10分别运行K-means，记录各指标 KRange = 1:10; sseAll = zeros(1, length(KRange)); silAll = zeros(1, length(KRange)); chAll = zeros(1, length(KRange)); dbAll = zeros(1, length(KRange)); for i = 1:length(KRange) K = KRange(i); opts.InitMethod = 'kmeanspp'; [labels, centroids, sse] = myKmeans(X, K, opts); sseAll(i) = sse(end); silAll(i) = mean(silhouetteCoefficient(X, labels)); chAll(i) = calinskiHarabaszIndex(X, labels, centroids); dbAll(i) = daviesBouldinIndex(X, labels, centroids); end figure; subplot(2,2,1); plot(KRange, sseAll, '-o'); title('SSE肘部法则'); xlabel('K'); ylabel('SSE'); subplot(2,2,2); plot(KRange, silAll, '-o'); title('平均轮廓系数'); xlabel('K'); ylabel('Silhouette'); subplot(2,2,3); plot(KRange, chAll, '-o'); title('CH指数'); xlabel('K'); ylabel('CH'); subplot(2,2,4); plot(KRange, dbAll, '-o'); title('DB指数'); xlabel('K'); ylabel('DB');

注意这段代码里每次运行K-means都应该设置InitMethod='kmeanspp'并固定随机种子，否则因为初始化不同，同一K值每次跑出来的SSE可能差很多，曲线全是锯齿，拐点根本看不出来。我做对比实验时会额外加一行rng(42)，保证整个候选K范围内的初始化条件公平一致。

4.2 图形判断之外的经验：换模型不如换数据，换数据不如换预处理

有一类数据，无论怎么选K，所有指标都很差，SSE、轮廓系数都提示“聚类结构不明显”。这时候不要急着在K-means框架内较劲，先检查输入数据究竟适不适合聚类。我总结的经验顺序是：

• 第一步检查标准化。量纲差异是否已经被消除？有些变量标准差比其他变量大几百倍，等于其他变量被直接忽略了。
• 第二步做降维可视化。数据如果超过三维，先做PCA，看前两个主成分的散点图。PCA的方差占比高，说明数据确实有结构；K-means结果和PCA结果联合起来看，才能真正判断聚类是否合理。
• 第三步检查是否存在离群点。K-means对离群点极其敏感，几个极端样本就能把某个簇的中心拉飞。先尝试用isoutlier或去除标准化后距离超过3σ的样本，再跑聚类，效果经常天壤之别。

我自己经历过一个案例：某个业务数据集共12个特征，第一次跑聚类，轮廓系数只有0.2左右，怎么调K都不行；后来画PCA图发现有两个特征的极端值把整个空间拉伸了，去掉那3%的离群样本再聚类，轮廓系数直接跳到0.45。这说明什么？预处理质量决定了聚类上限，评价指标只是帮你找到当前上限下的最优K。

4.3 K值不确定性下的稳健性判断

如果多个指标给出不同K，我还会加一步“稳定性验证”。做法是：对同一个K值，用不同的随机种子（或者不同的初始中心选择）重复跑10次，看标签结果是否稳定。计算两两运行之间的调整兰德指数（ARI）或归一化互信息（NMI），如果稳定地在0.9以上，说明这个K对应的聚类结构是数据固有的；如果每次跑出来都大不一样，说明这个K更像是初始化运气撑起来的。

MATLAB里可以用rng(seed)循环跑，把每次的labels存成矩阵列，然后计算列与列之间的相似度。具体代码不长，但需要在NMI函数上花一点功夫，这里就不展开了。由易到难的替代方案是看多次运行里SSE的方差，方差小说明结果稳定；方差大说明很可能落入了不同的局部最优解，这时候换个初始值策略比纠结指标更有效。

5. 工程化改进与踩坑记录：初始化、距离度量与数据预处理

这一节记录我实际使用K-means时踩过的坑和做过的改进，很多都和菜单栏里随手点的“默认配置”有关。

5.1 K-means++初始化：为什么它几乎成了默认选项

随机初始化可能让人连续几次跑进很差的局部最优解，SSE只是最优解的1.5倍甚至更高。K-means++的思路很简单：初始化第一个中心时随机选；之后依次选择新的中心，使得“距离已有中心越远的样本越容易被选为中心”。这一步近似保证了初始中心能铺满整个数据空间，大幅降低后续陷入局部最优的概率。

function centroids = initKmeansPP(X, K) n = size(X, 1); centroids = zeros(K, size(X, 2)); % 第一步：从样本中随机选一个作为第一个中心 centroids(1, :) = X(randi(n), :); % 计算每个样本到最近已有中心的距离平方 D = pdist2(X, centroids(1, :), 'euclidean').^2; for k = 2:K % 按距离平方的权重随机抽取 cumDist = cumsum(D) / sum(D); r = rand(); idx = find(cumDist >= r, 1, 'first'); centroids(k, :) = X(idx, :); % 更新最近中心距离 D = min(D, pdist2(X, centroids(k, :), 'euclidean').^2); end end

注意每次更新用的是min(D, 到新中心的距离)，代表“每个样本到当前已有中心的最近距离”。这个逻辑听起来简单，但写错成直接覆盖D就会导致后续选择全部偏向最新的中心，失去“铺满空间”的效果。我曾在一个项目里因为这一行写错，K-means++退化成近似随机初始化，实验结果惨不忍睹。

5.2 距离度量的选择：欧氏不是唯一答案

K-means教科书默认用欧氏距离，但我在处理文本TF-IDF向量时，效果一直不好。原因很好理解，文本向量的长度往往和文档篇幅强相关，欧氏距离会把“两篇内容相似但篇幅不同的文档”判得很远。换成余弦距离之后，聚出来的类别语义明显更合理。

但要注意一点：修改距离度量后，簇中心的更新公式也要变。余弦距离对应的合适中心不是算术平均向量，而是对簇内向量做归一化后求平均，再归一化。如果只知道改距离、不改更新规则，聚类很容易发散或者反复震荡。我在自己的代码里加了'cosine'分支，更新中心时现将簇内每个向量除以其模长，再取平均并归一化。这一步不写清楚，很多人会把欧氏版本硬套到余弦场景里，然后得出“余弦距离下K-means效果很差”的错误结论。

5.3 空簇问题只会在运行到一半的时候咬你

前面主函数里写了空簇处理，这里再展开说说这个坑为什么会发生。通常两种场景：第一种是初始化太差，比如某个离群样本单独成了初始中心，整个迭代过程没有任何样本靠它最近，这个簇永远是空的；第二种是数据集本身有极不平衡的分布，几个大簇把一个小簇“瓜分”了。处理空簇的方案有几种：

• 保留上一轮中心。实现最方便，但如果连续多轮都为空，这个中心就跟数据毫不相干，后续不会自己恢复。
• 从其他簇中随机拆分一些样本给它。实现稍复杂，但能帮助算法跳出局部最优。
• 用距离全局中心最远的样本重置该中心。这是我个人最推荐的做法，等于把空簇导向数据稀疏的区域，往往能找到之前遗漏的小簇。

% 空簇处理（替换上面主函数里的对应分支） for k = 1:K clusterPoints = X(labels == k, :); if isempty(clusterPoints) % 找到离全局中心最远的样本，用其初始化空簇 globalMean = mean(X, 1); distToGlobal = pdist2(X, globalMean); [~, farthestIdx] = max(distToGlobal); newCentroids(k, :) = X(farthestIdx, :); else newCentroids(k, :) = mean(clusterPoints, 1); end end

用“全局中心最远点”是因为空簇一旦分配到某个已有簇内部的样本，中心的更新会让这个样本重新归属到原来的簇，空簇问题还会复发；而最远点通常在大簇边缘，更可能站稳脚跟。

5.4 数据标准化与离群点影响：K-means的致命弱点

最后强调一个我在评审别人代码时反复指出的问题：很多人的数据根本没有标准化就直接聚类。K-means最小化的是“欧氏距离平方和”，如果第1个特征的量纲是0到1，第2个特征的量纲是0到1万，那么第2个特征对距离的贡献可以主导整个目标函数，第1个特征的信息完全被淹没。解决方法是zscore(X)，让每个特征均值0、方差1。这一步对聚类质量的影响往往比我调K值还大。

离群点方面，由于K-means是平方损失，离群点的距离平方则会被指数级放大，一个极端样本就能把一个簇的中心拉偏好几倍。预处理阶段应该先统一用箱线图或MAD（中位数绝对偏差）识别离群点，或者直接做截尾。如果用缩放后的数据，我一般会把超过3σ的样本标记出来，手动检查之后决定是否剔除。

6. 把代码整合成一套可复用的聚类分析工具箱

文章写到这里，核心代码和原理都讲完了。但我还留了一手：与其每次临时复制粘贴，不如把这些函数整合成一个结构清晰的分析工具箱。在实际项目里，我把myKmeans、initKmeansPP、四个指标函数、选K脚本归档为一个工具库，使用时只需要调用一个入口函数，从原始数据到指标报告一次性出结果。

6.1 主控脚本：一份数据跑出全部结论

function clusterReport(X, KRange) % clusterReport - 一站式聚类分析入口 % 输入: % X : n×d 原始特征矩阵（函数内部会做标准化） % KRange : 候选K值向量，例如 1:10 % 输出: % 控制台和图形窗口输出聚类报告 X = zscore(X); KRange = KRange(KRange >= 2); % K=1时内部指标无意义，跳过 n = length(KRange); sseAll = zeros(1, n); silAll = zeros(1, n); chAll = zeros(1, n); dbAll = zeros(1, n); dunAll = zeros(1, n); for i = 1:n K = KRange(i); rng(42); % 固定种子保证可复现 opts.InitMethod = 'kmeanspp'; [labels, centroids, sse] = myKmeans(X, K, opts); sseAll(i) = sse(end); silAll(i) = mean(silhouetteCoefficient(X, labels)); chAll(i) = calinskiHarabaszIndex(X, labels, centroids); dbAll(i) = daviesBouldinIndex(X, labels, centroids); if size(X, 1) <= 5000 % Dunn指数计算量大，限制样本量 dunAll(i) = dunnIndex(X, labels); else dunAll(i) = NaN; end end % 图形输出 figure; subplot(2,2,1); plot(KRange, sseAll, '-o'); title('SSE'); xlabel('K'); ylabel('SSE'); subplot(2,2,2); plot(KRange, silAll, '-o'); title('平均轮廓系数'); xlabel('K'); ylabel('Silhouette'); subplot(2,2,3); plot(KRange, chAll, '-o'); title('CH指数'); xlabel('K'); ylabel('CH'); subplot(2,2,4); plot(KRange, dbAll, '-o'); title('DB指数'); xlabel('K'); ylabel('DB'); % 控制台文本输出 fprintf('K\tSSE\t\tSilhouette\tCH\t\tDB\n'); for i = 1:n fprintf('%d\t%.2f\t%.4f\t\t%.2f\t\t%.4f\n', ... KRange(i), sseAll(i), silAll(i), chAll(i), dbAll(i)); end end

这个入口函数解决的问题是“报告里需要一张指标表”。以前每次都要复制五六段代码，换数据集时忘了改参数，现在一个函数搞定。里面的小细节也很实用：K=1时轮廓系数没有定义，所以直接过滤掉；Dunn指数计算量大，样本超过5000返回NaN，不硬算，也不让程序崩溃。

6.2 大数据量下的性能优化思路

当样本量超过十万量级时，K-means本身还能靠矩阵运算撑住，pdist2算全距离矩阵就不行了。优化方向有几个：

• 分块计算距离。把样本分成若干块，每块只与中心矩阵centroids计算距离，避免一次性生成n×K的大矩阵。K通常很小，所以n×K并不大，真正的大问题是pdist2(X, X)这种全距离矩阵（用于轮廓系数和Dunn指数）。样本一万时就是上亿规模，属于不可接受。
• 用抽样代替全量计算。先跑一遍K-means，然后每个簇抽取固定数量的样本计算轮廓系数近似值，代表整体。
• 换用更快的指标。CH和DB都只需要簇内距离或者中心距离，不需要全距离矩阵，在大样本场景下完全可以替代轮廓系数。

我在自己的项目里，数据量到5万行时，常规做法是用CH+DB选K，轮廓系数只在每个簇内抽500个样本估算，兼顾信息量和计算时间。

6.3 可视化：聚类结果不是一张散点图就结束了

聚类报告通常还需要可视化。PCA降到二维后画散点图是标配，但我会同时把每个簇中心也投影出来，用“×”标记，方便检查簇中心是否落在数据密集区域。还建议直接画轮廓系数的分布图，按簇分开展示每个样本的轮廓值，能立刻看出哪个簇里混入了边界样本。

% PCA降维后可视化 [coeff, score] = pca(X); proj = score(:, 1:2); % 前两个主成分 K = 3; [labels, centroids, ~] = myKmeans(X, K, optsWithKmeansPP); centroidProj = centroids * coeff(:, 1:2); figure; gscatter(proj(:,1), proj(:,2), labels); hold on; plot(centroidProj(:,1), centroidProj(:,2), 'kx', 'MarkerSize', 12, 'LineWidth', 2); hold off; title('PCA投影下的K-means聚类结果');

gscatter是MATLAB自带的分组散点图函数，一行就能按标签着色，配合簇中心标记，信息量比单色散点图高很多。这段代码也体现了为什么我在第一节说“MATLAB还是适合干这件事”：可视化、矩阵运算、内置数据集三方面配合得很顺手，你不需要为了画图再去学另一套工具链。

7. 如果你现在准备在自己的项目里用这套代码

建议的实践路径是：先跑一遍鸢尾花示例，把myKmeans和四个指标函数过一遍；然后换自己的数据集，先做标准化和离群点清洗，再用clusterReport画指标曲线；最后选一个稳定且可通过业务解释的K，用PCA可视化查看结果。我自己在这个框架上迭代了很长时间，它最核心的价值不是某一个函数的效率，而是让聚类过程里的每一步决策都可检查、可复现。

有一次带学生做实验，他用了这套代码，直呼“终于知道为什么报告里要放轮廓系数了”。其实算法本身不难，难的是把“公式为什么变成代码”“代码为什么变成报告”“报告为什么能支撑结论”这条链路打通。你如果也卡在这个链路上，不妨试试自己写一遍K-means，再把指标函数挂上去——很多之前想不明白的细节，会在第一个NaN出现在空簇中心的那一瞬间突然就懂了。