做聚类分析的时候我最早用的也是k-means毕竟它简单、跑得快MATLAB里一行kmeans就能出结果。但后来处理一批含离群点的客户分群数据时均值中心被几个极端样本拉得偏得离谱同一个簇里的样本被切得七零八落。那时候我才意识到距离均值这种“虚拟中心”在某些场景下确实不够稳。换用k-medoids之后聚类中心必须是真实存在的样本点离群点的影响一下子小了很多。这篇文章就围绕我常用的k-medoids MATLAB源码展开把数据导入、中文注释写法、聚类结果绘图这三个环节完整拆开讲代码都是可以直接拿回去改用的。整套代码是模块化写的主函数负责聚类迭代距离计算单独拆出来绘图部分独立成小节。这样做的原因是实际项目中数据格式、聚类个数、可视化需求都在变如果所有逻辑揉在一个脚本里每次改动都要来回翻找特别浪费时间。下面我按照从整体设计到细节实现再到踩坑记录的顺序来写不管你之前用没用过k-medoids都能照着落地。1. 为什么用k-medoids而不是k-means核心差异必须搞清1.1 “用样本代表簇”带来的鲁棒性提升k-means和k-medoids最大的区别就在“中心”的定义上。k-means计算每一簇的均值作为中心这个均值是特征空间里的一个虚拟点很可能在原始数据中完全不存在而k-medoids必须在每一簇内挑一个真实样本出来当代表这个代表就是medoid。挑的标准不是随机选一个而是在簇内计算所有样本两两之间的距离选出到其他样本总距离最小的那个样本。为什么这个差异在很多场景下很关键我举一个生产环境里的例子。当时做的是设备传感器数据聚类传感器偶发漂移会产生一些离群点。用k-means时均值中心直接被漂移点拽过去正常簇的边界跟着变形。换成k-medoids后离群点距离其他样本远它在簇内作为候选中心的总距离很大基本永远不会被选中聚类中心就能稳定落在样本密集区域。所以在数据含噪声、聚类结果要求中心可解释为“典型样本”的场合k-medoids是比k-means更合理的选择。从数学上看k-medoids优化的目标函数是“所有样本到其所在簇medoid的绝对距离之和”通常用曼哈顿距离或欧氏距离。由于中心取自样本点目标函数对离群点的敏感性天然低于k-means。代价是计算复杂度比k-means高不少这是下面要说的另一个话题。1.2 PAM算法思路与我的代码架构选择k-medoids最经典的实现是PAMPartitioning Around Medoids它的核心流程分两步也被称为构建阶段和交换阶段。构建阶段先选k个初始medoid可以是随机挑也可以用某种启发式交换阶段则逐个尝试用簇内非medoid样本替换当前medoid看目标函数是否下降能下降就保留替换。完整的PAM每次迭代要遍历所有可能的交换计算开销很大复杂度接近O(k(n-k)^2)。我实际写源码时没有把交换阶段做得这么彻底而是采用了一种“批量更新”策略每次迭代先按当前medoid分配样本然后在每个簇内部找出总距离最小的样本作为新medoid。这种方式本质上是在一个簇内做局部寻优虽然不保证全局最优但收敛快、编码简单在绝大多数中小规模数据集上效果和标准PAM非常接近。如果追求绝对精确可以在初始化时多跑几次随机种子选目标函数最小的那次结果。代码架构上我分成四个模块数据导入与预处理模块负责读取xlsx/csv/txt做标准化距离计算与分配模块计算样本与medoid距离分配簇标签medoid更新模块在每个簇内重新选举代表样本可视化模块绘制聚类散点图和轮廓图这样的分层让我在换数据集、换距离度量时不用动主循环逻辑只改对应模块即可。下面每一部分我都会贴出关键源码并解释为什么这么写。2. 数据导入与预处理MATLAB读取外部数据最容易踩坑2.1 三种常见格式的导入代码与适用场景MATLAB导入数据的方式很多早期版本有xlsread、csvread、load新版本统一推荐readmatrix和readtable。我在代码里按数据来源分成三种方式用注释标明读者按需取消注释即可% 方式1Excel表格最常用适合带特征名 % 注意readmatrix默认把第一行作为数值如果第一行是文本列名会丢弃 data readmatrix(iris.xlsx); % 方式2CSV逗号分隔文件 % data readmatrix(data.csv); % 方式3纯文本数据文件空格或制表符分隔 % data load(data.txt);如果你手里的数据第一行是特征名称比如“身高,体重,年龄”直接readmatrix会报错或者跳过非数值。这种情况我更推荐readtable它能保留列名后续做图形标注很方便T readtable(data.csv); % 自动识别列名 X table2array(T(:, 2:end)); % 去掉第一列样本ID转成double矩阵还有一个容易被忽略的问题文件路径。路径里如果带中文、空格或特殊字符MATLAB的读取函数经常莫名其妙报“文件不存在”或“无法打开”。解决办法是不依赖当前工作目录直接用绝对路径data readmatrix(D:\myproject\data\iris.xlsx);如果已经用cd切换了工作目录用相对路径也行但一定要确认当前目录下有对应文件。我的习惯是脚本开头统一判断文件是否存在if ~exist(iris.xlsx, file) error(文件不存在请检查路径%s, fullfile(pwd, iris.xlsx)); end2.2 标准化不是可选项是聚类前的必选项在代码中我专门加了一步zscore标准化。很多初学者直接拿原始数据丢进聚类算法结果量纲大的特征比如收入、购买金额完全主导了距离计算量纲小的特征比如年龄段、评分几乎不起作用。k-medoids用的是距离度量这个问题尤其突出。% 标准化每个特征减去均值除以标准差 X zscore(X);为什么要用zscore而不是最大最小归一化因为zscore对离群点更稳健它不会把极值强行压缩到固定区间而是根据数据的整体分布来放缩。如果你的数据已经有明确业务含义且所有特征量纲一致可以跳过标准化。但我在多数场景下都建议做哪怕只是简单试验也能避免很多“聚类结果看起来很奇怪”的问题。标准化之后需要保留原始数据的均值mu和标准差sigma因为后面绘制图形或者做业务解释时要把聚类中心还原回原始尺度mu mean(X_original); sigma std(X_original);如果你用的是readmatrix导入的矩阵建议在标准化前先复制一份原始数据否则你可能发现绘图时坐标轴全变成均值为0、标准差为1的范围不利于业务人员理解。3. 核心源代码逐段解析带中文注释的k-medoids实现3.1 主函数结构与参数设计我把整个聚类过程封装成一个函数这样可以直接在脚本里调用也可以放在循环里做多次随机初始化筛选。函数签名如下function [idx, medoids, obj_history] kmedoids_custom(X, k, max_iter, dist_metric) % 自定义k-medoids聚类函数 % 输入: % X: n行d列数值矩阵n为样本数d为特征数 % k: 聚类数 % max_iter: 最大迭代次数默认100 % dist_metric: 距离度量方式默认euclidean % 输出: % idx: n行1列整数向量样本所属簇编号(1到k) % medoids: k行d列矩阵每一行是一个真实样本点 % obj_history: 每次迭代的总体目标函数值用于观察收敛 % 说明: % 整个过程使用medoid簇内真实样本代替均值中心 % 对离群点有更好的鲁棒性。max_iter设100次其实在大多数情况下用不到因为批量更新策略收敛很快通常十几次迭代目标函数就不再下降了。但我仍保留这个参数防止个别数据分布下循环不终止。距离度量我会单独支持欧氏距离和曼哈顿距离两种。欧氏距离适合连续型特征曼哈顿距离对量纲和离群点更鲁棒。具体实现上不必手写距离公式直接调用MATLAB自带pdist2函数它支持多种距离度量还支持GPU加速大数据场景可用。function D compute_distance(X, center, dist_metric) % 计算样本矩阵X到单个中心center的距离向量 % 返回n行1列距离值 switch dist_metric case euclidean D pdist2(X, center, euclidean); case manhattan D pdist2(X, center, manhattan); otherwise error(不支持的距离度量); end end3.2 迭代主循环分配样本与更新medoid这部分是整个代码的核心。先随机从样本中选择k个作为初始medoid然后循环执行“分配-更新”两个步骤。初始选择要注意随机种子我习惯用rng(default)固定种子这样每次运行结果可复现。如果要做对比实验可以在外层再套一层循环跑20次取目标函数最小的结果。n size(X, 1); rng(default); init_idx randperm(n, k); % 随机选k个不同样本下标 medoids X(init_idx, :); % 初始化medoid obj_history zeros(max_iter, 1); for iter 1:max_iter % 步骤1计算每个样本到k个medoid的距离矩阵 D zeros(n, k); for j 1:k D(:, j) compute_distance(X, medoids(j, :), dist_metric); end % 分配每个样本选择距离最近的medoid的簇 [~, idx] min(D, [], 2); % 步骤2在每个簇内重新选择medoid new_medoids zeros(k, size(X, 2)); for j 1:k cluster_samples X(idx j, :); % 取当前簇所有样本 if isempty(cluster_samples) % 如果这个簇没有分配到任何样本空簇 % 随机选一个样本重新作为中心防止算法崩溃 new_medoids(j, :) X(randi(n), :); else % 计算簇内两两距离矩阵 D_in pdist2(cluster_samples, cluster_samples); % 行和表示当前样本到其他所有簇内样本的总距离 total_dist sum(D_in, 2); % 找总距离最小的样本作为新的medoid [~, best_idx] min(total_dist); new_medoids(j, :) cluster_samples(best_idx, :); end end % 判断是否收敛新旧medoid完全相同则停止 if isequal(medoids, new_medoids) medoids new_medoids; break; end medoids new_medoids; % 计算本次迭代的目标函数值所有样本到其簇medoid距离之和 D_final zeros(n, 1); for j 1:k members idx j; if any(members) D_final(members) D(members, j); end end obj_history(iter) sum(D_final); end % 截断为空迭代之后的部分 obj_history obj_history(1:iter);这一步里的“更新medoid”我用的是簇内距离矩阵行和最小法。假设簇内有m个样本我计算一个m乘m的距离矩阵然后对每一行求和行和最小的那个样本就是距离其他所有样本最近的样本即当前簇的medoid。这个做法比标准PAM的交换法简单很多但效果已经足够好尤其是当簇内样本分布呈凸形时。3.3 图形绘制让聚类结果一眼看懂的中文标注图形绘制是很多人容易忽略的一环其实它对快速验证聚类效果非常重要。二维数据可以直接画散点图每个簇用不同颜色medoid用叉号突出显示。如果特征超过两个可以选前两个主成分作为坐标轴来画或者用pairs画矩阵图。我在这里给出最常用的二维散点图代码并且特别处理了中文注释乱码问题。% 绘制聚类散点图 figure; gscatter(X(:,1), X(:,2), idx); hold on; plot(medoids(:,1), medoids(:,2), kx, MarkerSize, 15, LineWidth, 2); legend(簇1, 簇2, 簇3, medoid中心, Location, best); xlabel(特征1标准化后); ylabel(特征2标准化后); title(k-medoids聚类结果); set(gca, FontName, SimHei); % 设置黑体防止中文乱码 grid on; hold off;这里有几个细节值得说。gscatter是MATLAB自带的按分组上色散点图函数它需要第三个参数是分组向量我们的idx正好符合要求。plot中的kx表示黑色叉号MarkerSize需要设大一点否则medoid中心不够明显。legend如果不设置MATLAB会自动用“数据1”“数据2”这种名称完全看不出哪个对应哪个所以最好手动指定。还有一个更专业的可视化方式——轮廓图silhouette。轮廓图可以直观反应每个样本在簇内的紧密程度和簇间分离程度取值范围从-1到1值越大说明聚类效果越好。绘制起来也超简单figure; silhouette(X, idx); title(k-medoids聚类轮廓图); set(gca, FontName, SimHei);轮廓图对判断k值选得是否合理非常有帮助。如果大量样本的轮廓值接近0甚至为负说明这些样本处于簇边界或者可能放错了簇这时候就要考虑增大k或检查数据预处理。3.4 结果输出与保存聚类完成后最好把结果保存成文件方便后续分析和别人复现。我用writetable把样本ID、簇标签和原始特征都存成一张表T table(original_data, idx, VariableNames, {SampleID, Cluster}); writetable(T, clustering_result.xlsx);MATLAB的table支持中文列名但写Excel时如果有中文列名会遇到编码问题更稳妥的做法是用英文字段名存文件写报告时再人工映射。我在保存前总是会检查idx中的每个簇是否有足够样本如果某个簇只有一两个样本可能是k选大了或者数据分布本身就不适合当前k值。4. 实操过程与常见问题实测效果和排错经验4.1 用带离群点的二维数据对比k-means与k-medoids为了验证代码效果我构造了一组二维数据三个簇分别生成100个样本每个簇内部服从高斯分布再额外追加10个离群点分布在三个簇之外。分别用MATLAB自带kmeans和我写的kmedoids_custom跑聚类初始中心都固定随机种子。结果非常直观k-means的第一个簇中心明显被离群点拉走导致边界样本被错误划分而k-medoids的中心仍然停留在原始簇密集区准确率明显更高。我把这个过程也写进了一个简单脚本每次更换数据集时只需要改文件路径和k值就能复用。判断聚类效果时除了看准确率有真实标签的情况还要看目标函数的变化曲线。我在主函数中返回了每次迭代的目标函数值obj_history画成折线图能看到收敛过程。通常前三次迭代目标函数就下降明显后面逐渐趋于平稳。4.2 MATLAB中文注释乱码的完整解决办法很多读者下载别人源码后打开发现中文注释全部变成“涓冨瓧”之类的乱码。这基本不是代码问题而是文件编码不一致导致的。老版MATLAB2016之前默认用系统本地编码新版默认UTF-8如果你的源码编辑环境用了某种编码换台电脑就会乱。我自己的解决经验有三条首先统一使用新版MATLABR2017b以上在“预设项常规字符编码”中把默认字符集设为UTF-8。其次写注释时避免使用生僻字和特殊符号尽量用简体中文常用字。最后如果已经乱码可以在命令窗口执行feature(DefaultCharacterSet,UTF-8)然后重新打开文件大概率能恢复。如果代码文件本身已经损坏最笨但有效的方式是用纯文本编辑器Notepad或VS Code重新打开并另存为UTF-8格式再拷回MATLAB。这个方法我帮别人解决过好多次基本能救回来。4.3 数据导入失败的常见坑和排查逻辑数据导入这块遇到的报错多是文件路径、数据类型和表头问题。下面列一个我实际工作中整理的速查表报错场景原因解决方式“无法读取文件”或“文件不存在”路径含中文/空格或文件不在工作目录用绝对路径或先cd到目标目录readmatrix报错包含非数值列第一行是文本列名改用readtable或用data readmatrix(file, NumHeaderLines, 1)数据显示NaN或错误值单元格有空值或文本夹带逗号用readtable后填缺失值T rmmissing(T)导入后矩阵维度与期望不符文件末尾有汇总行或表头有多行增加参数readmatrix(file, NumHeaderLines, 2)数据量太大导入很慢文件超100MB且格式为xlsx转成csv或用parquet格式或用readmatrix的分块导入还有一个非常容易被坑的点中文列名在Excel里看起来很完美导入到MATLAB后列名变成Var1、Var2因为readtable默认会把无法识别为合法变量名的中文字段替换掉。所以我在数据预处理阶段会把表头改名成英文T.Properties.VariableNames {ID, Feature1, Feature2};4.4 k值选择和初始medoid随机性问题k-medoids和k-means一样都需要事先指定k。我在实际项目中一般用三种方法交叉验证一是肘部法画目标函数随k变化的折线寻找拐点二是轮廓系数法对不同k值计算平均轮廓系数取最大值三是业务法直接按业务可解释性定簇数。肘部法的MATLAB实现可以复用我写的目标函数。比如对k从2到10循环调用kmedoids_custom记录每次的obj_history最后一项然后plot出来ks 2:10; total_dists zeros(size(ks)); for i 1:length(ks) [~, ~, obj] kmedoids_custom(X, ks(i)); total_dists(i) obj(end); end plot(ks, total_dists, o-); xlabel(k); ylabel(总体距离); title(肘部图);需要注意的是由于初始medoid是随机选择的每次运行得到的最终目标函数可能略有差异。特别是当数据本身簇结构不明显时局部最优问题会被放大。我建议对每个k运行5到10次取目标函数最小值作为该k的代表值这样选出来的k更稳。很多期刊论文里也会注明“重复30次取最优”也是为了避免初始值带来的偏差。4.5 运行效率优化从PAM到CLARA的一个思路如果你的数据样本量在几千到几万之间我上面的批量更新方案跑起来还行。但如果到十万以上pdist2距离矩阵会直接耗尽内存。一种可行的替代是只对每簇采样一部分样本来计算候选medoid而不是全部样本。这其实就是CLARA算法的思想先对大样本集进行多次抽样对每个抽样子集运行PAM最后在所有子集结果里选目标函数最小的medoid集合。我在自己的代码里留了一个可选参数sample_size当簇内样本数超过该阈值时只随机抽取sample_size个样本参与medoid选举其余样本只做分配。这个改动在百万级样本上能显著降低时间和内存开销同时聚类质量和全量计算差距很小。如果以后遇到更大规模的数据建议转向Apache Spark或Python里的scikit-learn-extra实现MATLAB更适合教学验证和中小型项目。4.6 整套代码的注释风格与维护建议写中文注释其实不只是给人看还是帮助自己两个月后快速回忆的关键。我注意到很多网上的源代码注释写得太简略比如“更新中心”四个字就算完事但没说明为什么更新、更新原则是什么。我的注释习惯是“做什么 为什么 输入输出是什么”每段核心代码前都会写一两句背景。比如在更新medoid那几行我不仅写了“找簇内样本点为新的medoid”还补了一句“不能用均值替代否则退化为k-means”。这种注释对新手非常友好也避免自己以后误改。还有一点文件名最好能体现算法和用途比如kmedoids_custom.m而不是新建文档.m。函数内变量命名尽量语义化比如cluster_samples、total_dist、best_idx比用a、b、tmp好得多。这对代码维护的长期价值远大于写注释本身。最后分享一个小技巧在脚本里如果调试时想看看每次迭代的中间结果可以在循环里加上pause(1)配合drawnow这样你能在图上看到聚类中心逐步移动的过程。我第一次把这段代码跑起来时看着叉号在图上慢慢挪向簇中心那种理解算法执行过程的直观感远远强过只读文字说明。把这份代码吃透之后再回头去看k-means你会发现两者之间的差异远不止“中心怎么定义”这么简单它牵扯到鲁棒性、计算代价、场景适用性一系列选择。我做聚类项目时现在默认先跑一次k-medoids如果结果和k-means差异不大才考虑用k-means换效率如果差异明显那就说明数据里很可能存在离群点或者簇形状不规整这时k-medoids反而更可靠。这些经验都是踩过坑才总结出来的希望能帮你少走弯路。