MATLAB手写K均值聚类：注释详细、可修改样本数据的完整实现

用MATLAB实现K均值聚类，几乎是每个接触机器学习的人都会做的一道基础练习。不过奇怪的是，网上能找到的版本要么注释少得可怜，复制下来只能干瞪眼；要么封装了一大堆高级函数，想改一个参数都得翻半天文档。这半年帮几个做课程设计的同学看过代码，发现大家的需求其实很统一：要一份原理不绕弯、注释写得明明白白、样本数据可以随手替换的“干净版K-means”。这篇文章就针对这个需求展开，给出完整代码、逐段讲透逻辑，再把换数据、调K值、做归一化这些实际操作中避不开的细节全部串一遍。适合正在学聚类算法、或者在课程设计/毕业设计里需要一份能读懂、能改、能跑的参考实现的人。

1. 这份代码解决什么问题，为什么要自己写

1.1 K均值聚类到底是干什么的

K均值聚类（K-means）是无监督学习里最经典、应用最广的算法之一。所谓“无监督”，就是拿到的数据没有标签，我们不知道哪些样本属于同一类，但希望算法能根据样本之间的距离关系，把相似的样本自动凑成一堆。它解决的问题可以概括成一句话：给定一堆样本点和一个整数K，把样本分成K个簇，使得同一个簇里的样本彼此尽量相近，不同簇的样本尽量远离。

举一个生活化的例子：假设你是一个快递配送站的站长，手上有几百个收件地址，需要安排K个配送员。你要做的就是把地址分成K片区域，让每个配送员负责一个区域，并且尽量让每个配送员跑的路程最短。K-means干的就是这件事，只是它落脚在数学上——那些“地址”变成了多维空间里的坐标点，“路程”变成了样本点到簇中心的欧氏距离。

这份代码解决的问题场景非常具体：给定一个N行D列的样本矩阵，每一行是一个样本点，每一列是一个特征维度，代码自动完成聚类，并把结果画出来。二维数据和三维数据可以直接可视化，高维数据则可以配合降维工具查看效果。对于初学聚类算法的人而言，它最大的价值在于完整、简单、可运行，并且所有注释都写在代码旁边，几乎没有“隐藏逻辑”。

1.2 为什么坚持手写而不是直接调内置函数

有个问题我几乎每次都会被问到：MATLAB的统计工具箱里不是有现成的kmeans函数吗？直接调用它不香吗？为什么要自己写？

答案是：没有错，内置kmeans函数功能更强、速度更快、还支持各种初始化策略，但它的内部实现对你来说是一个“黑盒”。对于正在学习算法原理的人来说，亲手把四步迭代写一遍，和直接调一次函数，收获完全不一样。自己写代码时，你会被迫搞清楚初始化怎么选、距离怎么算、中心点怎么更新、收敛条件怎么判断这几个核心环节；而直接调函数，你得到的只是一个结果。

另外还有两个现实原因。第一，很多同学的MATLAB是没有装统计工具箱或者机器学习工具箱的，而手写代码只需要基础MATLAB环境，完全没有附加工具箱依赖。第二，内置函数的封装程度太高，想在里面加一行“把每次迭代的误差记录下来”都得折腾半天，手写代码则可以随意插桩调试。这也正好对应标题里“注释详细、样本数据可修改”的定位——它是一份教学性质、可自由改造的基础实现。

2. 算法原理与代码设计思路

2.1 从原理层面拆解K-means的四步迭代

K-means的完整流程可以归纳为四步，循环往复直到收敛。

第一步是初始化。算法一开始并不知道K个簇的中心在哪里，所以需要随机选K个样本点作为初始簇中心。这里“随机”很关键，因为它直接影响最终结果的好坏，选得不巧可能陷入局部最优。第二步是分配。对于每个样本点，计算它到K个簇中心的距离，把它归到距离最近的那个簇。第三步是更新。每个簇的所有样本点取平均，得到新的簇中心。第四步是判断收敛。如果新旧簇中心的变化量小于某个阈值，说明聚类已经稳定，迭代停止；否则回到第二步继续。

整个过程的数学表达也不复杂。假设样本是x1到xN，簇中心是μ1到μK，每个样本到簇中心μk的距离用欧氏距离来度量：d(xi, μk) = ||xi - μk||²。样本xi会被分配到使这个距离取最小值的那个簇。更新簇中心时，μk就取簇k内所有样本的均值，即μk = (1/|Ck|) * Σxi∈Ck xi。所谓“收敛”，就是SSE，也就是所有样本到各自簇中心距离平方之和，不再明显下降。

这段逻辑对应的代码就是主程序里的for循环。把四步拆清楚之后再去看代码，你会发现每条语句都有明确归属，完全没有“不明觉厉”的片段。

2.2 代码模块的划分与变量规划

写这份代码之前，我先把结构分成六大块，这也是我建议初学者拿到一份代码后最先看的部分：

1. 清空环境并准备数据：定义data矩阵，N和D自动读取
2. 参数区：K值、最大迭代次数、收敛阈值
3. 初始化簇中心：用randperm随机抽取K个样本
4. 主迭代循环：分配样本、更新中心、检查收敛
5. 可视化：用散点图画出簇归属和最终中心点
6. 控制台输出：打印簇中心和标签

变量设计上，我用data存放样本矩阵，用N表示样本个数，D表示特征维数；center是K行D列的簇中心矩阵，label是N行1列的标签向量，distMatrix是N行K列的距离矩阵。变量名尽可能做到“看名知义”，这样不管是谁拿到代码，都不需要来回翻注释才能理解某个变量在干嘛。

这种模块化设计还有一个额外好处：改动数据、调整K值时，不需要动主迭代逻辑，边界清晰。这也正是“样本数据可以进行修改”的工程基础。

3. 核心代码逐段解析

3.1 完整代码参考

先给整份代码，后面我会逐段拆开讲关键逻辑。

% ===================================================================== % 基本的K均值(K-means)聚类算法 % 数据格式：data为N行D列矩阵，N为样本个数，D为特征维数 % 输出：簇标签label(N行1列)，簇中心center(K行D列) % ===================================================================== clc; clear; close all; %% 1. 准备样本数据 % 每行是一个样本点，这里给的是二维坐标，方便画图观察效果 % 可以看到样例数据大致分成三堆：左下角、右上角、右下角 data = [ 1.0, 1.0; 1.2, 0.8; 1.5, 1.3; 0.8, 1.2; 5.0, 4.0; 5.3, 4.2; 4.8, 3.8; 5.5, 4.5; 9.0, 2.0; 9.3, 2.4; 8.8, 1.8; 9.5, 2.2; ]; [N, D] = size(data); % N是样本个数，D是特征维数 %% 2. 聚类参数设置 K = 3; % 聚类簇数，可修改 maxIter = 100; % 最大迭代次数，防止代码不收敛时陷入死循环 tol = 1e-4; % 收敛阈值：中心点变化量小于该值时停止迭代 %% 3. 初始化簇中心 % 从所有样本中随机抽取K个，作为初始簇中心 rng(0); % 固定随机种子，保证每次运行结果一致；想体验随机性可注释掉 initIdx = randperm(N, K); center = data(initIdx, :); %% 4. K-means主迭代 label = zeros(N, 1); % 预分配标签向量 for iter = 1:maxIter % ---- 步骤1：计算每个样本到每个簇中心的欧氏距离 ---- distMatrix = zeros(N, K); for k = 1:K diff = data - center(k, :); % 每个样本与第k个中心的坐标差 distMatrix(:, k) = sqrt(sum(diff.^2, 2)); % 按行平方求和再开根号 end % ---- 步骤2：把每个样本分配到距离最近的簇 ---- [~, label] = min(distMatrix, [], 2); % label中存的是最近簇的编号 % ---- 步骤3：将每个簇的样本均值作为新的簇中心 ---- for k = 1:K members = (label == k); % 逻辑索引，找出属于簇k的样本 if sum(members) > 0 % 新中心 = 该簇所有样本点在各维度上的平均值 center(k, :) = mean(data(members, :), 1); else % 空簇处理：保持原中心不变，避免代码直接崩溃 end end % ---- 步骤4：判断收敛 ---- if iter > 1 change = norm(center - oldCenter, 'fro'); % 中心矩阵变化量 if change < tol fprintf('第%d次迭代后收敛\n', iter); break; end end oldCenter = center; % 保存本次中心，供下次比较 end %% 5. 可视化聚类结果 figure; hold on; grid on; % 预设12种颜色，超过12个簇时自动循环取色 colors = [1 0 0; 0 0 1; 0 1 0; 1 0 1; 0 1 1; 1 1 0; ... 0.5 0 0; 0 0 0.5; 0 0.5 0; 0.5 0.5 0; 0 0.5 0.5; 0.5 0 0.5]; for k = 1:K members = data(label == k, :); scatter(members(:, 1), members(:, 2), 60, colors(mod(k-1,12)+1, :), 'filled'); end % 画出最终簇中心，用黑色五角星突出显示 plot(center(:, 1), center(:, 2), 'kp', 'MarkerSize', 14, 'MarkerFaceColor', 'k'); xlabel('特征1'); ylabel('特征2'); title(['K-means聚类结果，K = ' num2str(K) ', 迭代 ' num2str(iter) ' 次']); hold off; %% 6. 输出结果 disp('最终簇中心：'); disp(center); disp('每个样本的簇标签：'); disp(label);

3.2 关键代码段的逐行逻辑

先看初始化部分。randperm(N, K)是从1到N的整数序列里随机抽取K个不重复的值，返回值是一个包含K个元素的向量。用这K个下标去索引data矩阵，就得到了K个初始簇中心。这里有一个容易忽略的细节：k-means的初始中心一般是从样本点里抽，而不是在整个特征空间里随机生成坐标。原因在于，从样本中抽取能让中心点落在数据分布范围内，减少初始中心远离数据群的风险。

再看距离计算。我用了两层循环，外层遍历K个中心，内层直接从data矩阵减掉center(k,:)。这里利用了MATLAB的矩阵广播机制：data - center(k,:)会得到一个N行D列的矩阵，每一行都是“当前样本减去第k个中心”的差向量。sum(diff.^2, 2)表示按行求和，得到每个样本点到该中心的距离平方；sqrt开根号以后就是欧氏距离。为什么不直接用pdist2这类现成距离函数？因为手写循环能让你明确看到“每个距离是多少”，也更容易在调试时输出中间结果。

min函数的用法需要单独说一句。[~, label] = min(distMatrix, [], 2)中，第一个输出是每行的最小值，这里用~丢弃；第二个输出是每行最小值的列索引，也就是最近的簇编号。中间的[]是固定写法，表示在第二维方向求最小值。这个语法在MATLAB里新手很容易漏掉，少了[]或2都会报错或者得到完全不同的结果。

更新中心点的代码里，members = (label == k)是一个逻辑向量，成员为1说明样本属于簇k。data(members, :)通过逻辑索引取出这些样本，mean(..., 1)对每一列求均值，即取所有坐标的平均值。这一行是整个算法数学上最直观的落点——簇中心就是簇内样本的几何中心。

收敛判断部分用了一个辅助变量oldCenter。每次迭代开始时把中心存下来，迭代结束后计算新旧中心的Frobenius范数。这个范数可以理解为“中心矩阵所有元素变化的平方和再开根号”，用它来衡量整体变化幅度非常合适。当变化量小于tol时，说明聚类已经稳定，跳出循环。

空簇处理是我特意加进去的。K-means在迭代过程中可能因为初始中心选得不好或者某个簇只有一两个样本，导致某个簇在分配阶段没有任何成员。如果没有空簇保护，mean(data(members, :), 1)运算后中心会变成NaN，后面的距离计算全部崩掉。这里选择“保持原中心不动”，是最简单稳妥的处理方式。

4. 样本数据的修改方法与结果解读

4.1 三种方式把自己的数据放进去

标题里专门提到“聚类的样本数据可以进行修改”，如果你把代码下载下来，最自然的改动方式就是直接编辑脚本里的data矩阵。原来的样例数据是12行2列，你完全可以改成自己的几十行甚至几百行数据。需要注意一个约定：每一行必须是一个样本，每一列是一个特征。千万别搞成“每一列是一个样本”的矩阵排列方式，否则聚类结果会完全乱套。

如果数据量比较大，或者数据存在Excel表格里，就直接改读取方式。下面两种替换方式最常见：

% 从Excel读取，数据从第一个sheet的A1开始，不要包含文本表头 data = xlsread('mydata.xlsx'); % 从CSV读取，同样要求第一行就是纯数值 data = readmatrix('mydata.csv');

注意，如果用readmatrix读取带表头的CSV，需要额外指定跳过表头行数。比如表头占了第一行，就写data = readmatrix('mydata.csv', 'NumHeaderLines', 1);。还有一个高频坑：Excel读取时如果表格里有空行、空列或者中文字段，读进来以后data里会出现NaN，聚类代码跑起来会直接报错。稳妥的做法是先对data做一次清洗：data(isnan(data)) = [];再重塑成矩阵，或者干脆在Excel里把数据整理成一整块矩形区域。

修改完数据以后，K值也需要跟着改。样例数据肉眼就能看出三堆，所以K=3是合理的；如果换成自己的数据却忘记把K一起改掉，聚类结果画出来就会非常怪异。我的习惯是先用网格扫描一遍K，画一个SSE折线图再定K值，这个流程放在后面讲。

4.2 聚类参数的三个关键选择技巧

K值怎么选，是K-means里最经典的问题，没有标准答案，但有一个非常实用的经验方法叫肘部法则。它的思路很简单：从小到大依次尝试K=1、2、3……每次记录总误差平方和SSE，然后画出SSE随K变化的曲线。随着K增大，SSE一定会下降，因为簇变多了，样本离中心的整体距离会越来越小。但曲线在某个K附近会出现一个明显的拐点，像手肘一样，拐点之后的下降速度骤减，这个拐点对应的K通常就是比较合适的取值。

要计算SSE，你可以在代码里加一句统计：

% 在迭代结束后计算SSE（误差平方和） SSE = 0; for k = 1:K members = data(label == k, :); diff = members - center(k, :); SSE = SSE + sum(sum(diff.^2)); end disp(['SSE = ', num2str(SSE)]);

试着把这个SSE分别算出来，画一条折线，肉眼选点。这个方法虽然朴素，但在实践里特别抗打。

其次是迭代次数与收敛阈值的关系。maxIter设100对大多数数据都够，但如果你的样本有上万个点、维度又高，距离计算量大，迭代一次可能比较慢，maxIter可以放宽到500。tol决定聚类的精度，太小会导致迭代次数增加，但一般不需要低于1e-6；太大则可能在中心点还没稳定时就提前退出，聚类最终看起来比较粗糙。

最后是随机种子。rng(0)这一行把随机数生成器固定下来，效果是每次运行得到一模一样的初始中心，因此结果完全可复现。如果你是做实验、写报告、或者想把聚类结果给导师复核，固定种子会省掉很多“为什么我跑出来和你不一样”的麻烦。但如果你是在探索数据、想观察不同初值对结果的影响，把rng(0)注释掉即可。

4.3 从聚类结果里能读出什么

运行完代码，你会在图上看到几坨颜色不同的点，每个点属于一个簇，中心是一个黑色五角星。此时不要只感叹“效果不错”，要养成读图的习惯。第一，看同一簇内的点是否彼此靠近、不同簇之间是否明显分开。如果簇与簇之间大面积重叠，说明特征选择不好或者K值偏大，聚类没有实际意义。第二，看簇中心的实际位置代表什么。在客户分群场景中，中心点各维度的坐标值就是该类客户在各指标上的平均水平，比如“簇1是低消费高频次客户，簇2是高消费低频次客户”这种结论就是从中心坐标推出来的。第三，看每个簇的样本数是否均衡。如果某个簇只有一两个样本，可能是离群点干扰，也可能是K取大了。

另外有一点我要特别提醒：K-means对特征的量纲非常敏感。举个例子，如果数据有两个特征，年龄范围是0到100，收入范围是20000到500000，直接聚类的话，距离贡献几乎全部由收入决定，年龄等于被无视了。这种情况下聚类效果很差。正确的做法是先归一化，让每个特征的尺度一致。MATLAB里最省事的归一化是用zscore：

data_std = zscore(data); % 每一列减均值、除标准差

归一化之后再跑K-means，通常聚类曲线会顺眼很多。不过也有反例：如果某个特征的数值差异本身就是业务上最关心的信息，强行归一化反而会丢失区分度。要不要归一化，取决于你想让哪个维度的差异主导聚类结果。

5. 常见问题与排查技巧实录

5.1 运行报错与异常结果速查表

这些坑几乎每个人都会踩，我整理成一张表，遇到对应问题可以直接翻：

5.2 局部最优问题与“多次运行取最优”策略

你可能会遇到一种情况：同样的数据、同样的K，只是把rng(0)注释掉再跑几次，聚类结果画出来完全不一样，有时候明显合理，有时候离奇。这不是代码有bug，而是K-means从随机初始中心出发，本质是在最小化一个有多个“坑”的目标函数，算法很容易掉进局部最优解，而不是全局最优解。

MATLAB自带的kmeans函数默认就做了这件事：'Replicates'参数控制重复实验次数。手写版本想达到类似效果，代码也不复杂，核心思路是：对多个随机初始中心分别跑完整聚类，记录每次的SSE，最后取SSE最小的那次结果作为最终输出。可以这样改造：

bestSSE = inf; for r = 1:10 % 重复10次 % 这里重新随机初始化center，然后完整跑一遍主循环 % 用4.2节的方法计算本次的SSE if SSE < bestSSE bestSSE = SSE; bestLabel = label; bestCenter = center; end end % 最终用bestLabel和bestCenter输出结果

这个方案在实际操作中非常管用。虽然多跑了九次迭代，但对中小规模样本完全无压力，得到的聚类稳定性会明显提升。

5.3 数据预处理与距离度量的进阶提醒

如果发现数据本身就是数量级相差很大，先从归一化入手，而不是换算法。归一化之后依然聚类不好，再考虑数据是否真的存在簇结构。K-means隐含假设每个簇是各向同性的“球形团”，如果数据是细长的环形、螺旋形，K-means怎么跑都分不对，这时候应该换方法。

另外，这篇代码用的是欧氏距离，这是K-means的默认配置。如果你的特征包含类别型取值，比如“红绿蓝”“高中低”，欧氏距离意义不大，得先做数值化编码，或者换用其他聚类算法。在MATLAB里如果确实想换距离度量，可以用pdist2(data, center, 'cityblock')这种写法替代手写的欧氏距离，但要意识到，换了距离之后，中心点更新公式也需要跟着调整，不能继续无脑求均值，否则收敛性和结果含义都会出问题。

我自己在使用这份代码时，还经常把每次迭代的中心点变化量打印出来看。往往刚开始几次迭代变化幅度很大，后面逐渐缩小，最后小于阈值跳出。这是一个直观的“观察算法收敛”的过程，比单纯看最终结果更有助于理解K-means的迭代本质。

最后再分享一个小技巧。如果你要处理的是带分类标签的公开经典数据集，比如Iris鸢尾花数据，先把真实标签去掉跑K-means，再把聚类结果和真实标签做对比，计算一下正确匹配的样本占比。这样一来，手写代码的价值就不仅是跑通一个Demo，而是能帮你直观理解“无监督聚类”和“有监督分类”之间的差距。这个过程放到课程设计里，是很加分的扩展实验方向。