
用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 代码模块的划分与变量规划写这份代码之前我先把结构分成六大块这也是我建议初学者拿到一份代码后最先看的部分清空环境并准备数据定义data矩阵N和D自动读取参数区K值、最大迭代次数、收敛阈值初始化簇中心用randperm随机抽取K个样本主迭代循环分配样本、更新中心、检查收敛可视化用散点图画出簇归属和最终中心点控制台输出打印簇中心和标签变量设计上我用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值也需要跟着改。样例数据肉眼就能看出三堆所以K3是合理的如果换成自己的数据却忘记把K一起改掉聚类结果画出来就会非常怪异。我的习惯是先用网格扫描一遍K画一个SSE折线图再定K值这个流程放在后面讲。4.2 聚类参数的三个关键选择技巧K值怎么选是K-means里最经典的问题没有标准答案但有一个非常实用的经验方法叫肘部法则。它的思路很简单从小到大依次尝试K1、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里最省事的归一化是用zscoredata_std zscore(data); % 每一列减均值、除标准差归一化之后再跑K-means通常聚类曲线会顺眼很多。不过也有反例如果某个特征的数值差异本身就是业务上最关心的信息强行归一化反而会丢失区分度。要不要归一化取决于你想让哪个维度的差异主导聚类结果。5. 常见问题与排查技巧实录5.1 运行报错与异常结果速查表这些坑几乎每个人都会踩我整理成一张表遇到对应问题可以直接翻问题现象可能原因解决方案报错“索引超出矩阵维度”K值大于样本数N或者data是行向量检查K取值确认data每行是一个样本某个簇始终没有成员初始中心选得太集中或K值偏大换随机种子重跑减小K检查数据分布每次运行聚类结果不同初始中心随机陷入不同局部最优固定rng种子多次运行取SSE最小的结果聚类效果很差、簇边界奇怪数据未归一化量纲差异大先zscore归一化再聚类迭代到maxIter仍未收敛收敛阈值太小或maxIter不足适当增大maxIter检查距离计算是否正确读取Excel后程序报NaN错误数据中有空行、空列、文本在Excel里整理成纯数值矩形区域或读取后清理NaN可视化时只有两个簇显示颜色矩阵行数小于K索引越界扩展颜色矩阵或参考代码中用mod循环取色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而是能帮你直观理解“无监督聚类”和“有监督分类”之间的差距。这个过程放到课程设计里是很加分的扩展实验方向。