平时用MATLAB做聚类分析绕不开k-means但一旦数据里混了几个离群点k-means的均值中心就会被拽得七荤八素。这时候该换k-medoids了。我在实际项目里经常碰到这种场景传感器数据偶尔跳一个异常值用户行为数据带点噪声标签用medoid簇内真正存在的样本代替均值聚类中心就不会被离群点带偏。这篇博文就把我整理的k-medoids聚类MATLAB源代码、数据导入和图形绘制整套流程放出来全程带中文注释拿过去就能改能跑适合正在做数据分析、模式识别实验或者毕业论文需要聚类对比的读者。我会把每一段代码的来龙去脉、参数选择、画图细节都讲透尽量减少你踩坑的时间。1. k-medoids聚类原理与选型分析1.1 为什么用medoid替代centroidk-means和k-medoids的差异本质上是“平均数”和“中位数”的差异。k-means每个簇的中心是簇内所有样本的算术平均可能计算出一个人造坐标k-medoids的中心则是簇内某个真实存在的样本点该点是簇内所有样本到其他样本总距离最小的那个。这个区别带来的直接影响就是鲁棒性k-means对离群点极度敏感一个偏离很远的点就能把centroid拉过去导致几个簇被错误合并k-medoids不会因为中心必须选在真实数据上离群点通常孤立成簇或对中心选择影响有限。从距离度量上说k-medoids对距离函数的要求更宽松。k-means在欧氏距离下有闭式解求均值但如果你需要做曼哈顿距离、余弦距离甚至自定义相似度k-means很难推广——均值在非欧空间没有定义。k-medoids只要求“能算样本两两距离”然后枚举簇内样本选最优中心所以处理分类变量、混合类型数据、缺失值较多的数据时反而更方便。当然代价是计算量PAMPartitioning Around Medoids的经典实现复杂度大约是O(k(n-k)²)比k-means重不少。下面会讲怎么在MATLAB里平衡效果和效率。1.2 算法核心流程与收敛性理解k-medoids的核心步骤只有四步但每一步都有细节讲究初始化从n个样本里随机选k个作为初始medoid或者用启发式方法选分散的样本。分配把每个样本分配到距离最近的medoid所在簇。更新对每个簇遍历簇内所有样本选一个能使簇内总距离样本到中心距离之和最小的样本作为新medoid。收敛判断如果medoid集合不再变化或者代价函数的变化量小于阈值就停止迭代。值得注意的收敛性问题是k-medoids的目标函数是“所有样本到其medoid距离的总和”这个函数在枚举更新下是单调递减的所以算法保证收敛到局部最优但可能是局部最优而非全局最优。初始点选得不好很容易收敛到差的结果。这也是为什么测试代码时要多跑几次随机初始化或者用上一节说的k-means式启发式初始化。1.3 适用场景什么样的数据值得用k-medoidsk-medoids不是要完全取代k-means。经验上这几类场景优先选k-medoids数据存在明显离群点且离群点不是你需要单独聚出来的噪声类而是混在正常样本附近。距离度量是曼哈顿距离或自定义相似度无法直接求均值。聚类中心后续需要解读而你的“中心”必须是真实对象。比如客户聚类希望中心对应一个真实客户画像商品聚类希望中心是某件实际商品。样本量不太大基本在几千到一两万量级k也不太大。如果你有百万级样本建议先抽样或者用Mini-batch策略否则迭代速度会很难看。2. MATLAB数据导入与预处理实战2.1 从Excel、CSV、TXT导入数据MATLAB读取外部数据的主要函数有readmatrix、readtable、xlsread老版本、load。其中readmatrix在R2019a之后是首选它自动识别数值和文本返回纯数值矩阵省去table2array转换。我平时用得最多的代码是% 读取CSV文件跳过第一行表头只取数值列 data readmatrix(iris_data.csv, NumHeaderLines, 1); % 或者读取Excel指定工作表 data readmatrix(dataset.xlsx, Sheet, Sheet1, Range, A2:C150);如果数据里第一列是样本ID或标签文本需要分离data_all readtable(dataset.xlsx); labels data_all{:, 1}; % 文本标签 features data_all{:, 2:end}; % 数值特征 features varfun(str2double, features); % 如果需要转换从txt导入时经常遇到分隔符问题readmatrix默认按逗号但如果txt是用空格或Tab分隔要指定Delimiter, \t或Delimiter, 。踩过多次坑建议先打开txt预览一下第一行确认分隔符和表头位置再写导入代码。2.2 数据标准化聚类前必须做的事k-medoids的距离计算完全依赖量纲。如果第一个特征范围是0到1第二个特征范围是0到10000第二个特征会主导整个距离第一个特征等于没参与。所以在聚类前我一般做z-score标准化也就是每个特征减去均值除标准差mu mean(data); sigma std(data); data_std (data - mu) ./ sigma;还有一种min-max归一化把数据缩放到[0,1]区间data_min min(data); data_max max(data); data_norm (data - data_min) ./ (data_max - data_min);两种方法选哪个如果你的数据分布近似高斯z-score更稳如果有硬边界min-max更直观。需要注意标准化必须在划分训练集和测试集之前或之后聚类是无监督任务不存在训练测试划分所以直接对全部数据标准化即可。但要保存mu和sigma以便新样本加入时用同样的参数变换。2.3 缺失值与异常值的处理策略缺失值在MATLAB里通常表现为NaN。直接用含NaN的数据算距离结果会全部变成NaN。处理方案分几种如果缺失比例小于5%可以直接删除对应行data_clean data(~any(isnan(data), 2), :);用列均值或中位数填充for j 1:size(data, 2) col data(:, j); col(isnan(col)) median(col(~isnan(col))); data(:, j) col; end如果你的特征是分类变量可以用众数填充不过k-medoids本身不太吃这个亏因为它能处理非欧距离。异常值的处理反而有讲究k-medoids本身抗离群点所以不建议轻易删除异常值删了反而会丢失信息。我一般先画出箱线图或做z-score检测看一下异常点占比如果异常点是某种测量错误就剔掉如果是真实分布的一部分就保留。保留的情况正好发挥k-medoids的鲁棒性优势。3. 核心源代码实现与中文注释详解3.1 主函数结构设计与参数定义我写的k-medoids实现主函数输入是特征矩阵X和簇数k输出是簇分配标签idx、medoid索引、每轮代价函数历史。函数头长这样function [idx, medoid_idx, cost_history] myKmedoids(X, k, maxIter, distType) % myKmedoids 基于PAM思路的k-medoids聚类实现 % 输入 % X - 样本特征矩阵每一行是一个样本每一列是一个特征 % k - 聚类簇数 % maxIter - 最大迭代次数默认100 % distType - 距离度量类型euclidean或manhattan默认euclidean % 输出 % idx - 每个样本所属簇的标签n行1列 % medoid_idx - 每个簇的medoid在原始数据中的行号 % cost_history - 每轮迭代的代价函数值用于绘制收敛曲线这个函数接口是通用的后续你换数据集、换簇数都不需要改动主循环。参数默认值用nargin判断在MATLAB里比较常见if nargin 3 maxIter 100; end if nargin 4 distType euclidean; end3.2 初始化策略随机选择与启发式初始化PAM算法的原始初始化是随机抽k个样本做medoid但这样在多峰数据上容易掉进局部最优。我在代码里加了两个选项随机初始化和基于最大距离的初始化类似k-means的思路。第二种做法是第一个medoid随机选后续每个medoid选择离已有medoid最远的样本。这个策略能让初始中心尽量分散收敛稳定不少。n size(X, 1); if strcmpi(initMethod, random) medoid_idx randperm(n, k); elseif strcmpi(initMethod, maxdist) % 第一个中心随机 medoid_idx randi([1, n], 1, 1); % 后续中心选择距离已有中心最远的样本 for j 2:k D pdist2(X, X(medoid_idx, :), distType); minDist min(D, [], 2); [~, idx_new] max(minDist); medoid_idx [medoid_idx; idx_new]; end medoid_idx medoid_idx; end需要说明maxdist初始化会比随机初始化多一次n×k的距离计算但换来的是迭代次数明显减少整体耗时反而更低。对于n在几万以下的数据这个开销完全可以接受。3.3 计算距离矩阵与簇分配分配阶段需要计算每个样本到k个medoid的距离。用MATLAB内置的pdist2一次性算距离矩阵又快又简洁D pdist2(X, X(medoid_idx, :), distType); [~, idx] min(D, [], 2);pdist2的distType支持euclidean、squaredeuclidean、cityblock曼哈顿、cosine、correlation等。如果你需要自定义距离就自己写一个函数句柄传入pdist2或者干脆用三重循环算距离矩阵。但三重循环在MATLAB里很慢能向量化尽量向量化。这里有个细节pdist2对大矩阵比较吃内存n×k的矩阵还好但当n到了几十万就要考虑分块计算。我还没在这个代码里做分块如果你的数据量特别大建议把距离计算改成for循环按块处理避免内存爆掉。3.4 更新medoid簇内枚举与总距离最小化更新阶段是k-medoids和其他算法最不同的地方。对每一个簇取出簇内所有样本计算簇内样本两两之间的距离再选出“到其他样本总距离最小”的那个样本作为新medoid。这里不能用pdist2对整个簇矩阵做两步操作直接算平方距离矩阵再求和for j 1:k cluster_idx find(idx j); if isempty(cluster_idx) % 如果某个簇为空重新随机指定一个样本 cluster_idx randi([1, n], 1, 1); idx(cluster_idx) j; end % 簇内样本的特征矩阵 cluster_X X(cluster_idx, :); % 簇内两两距离矩阵欧氏距离的平方也等价 Dc pdist2(cluster_X, cluster_X, distType); totalDist sum(Dc, 2); % 每个样本到簇内其他样本的总距离 [~, bestLocalIdx] min(totalDist); medoid_idx(j) cluster_idx(bestLocalIdx); end这里有一个容易踩的坑sum(Dc, 2)会把样本到自己距离0也算进去但0不影响最小值比较所以没问题。另一个坑是如果簇内有重复样本完全一样的行pdist2会产生0距离可能导致totalDist偏小但不影响选出正确medoid。3.5 收敛判断与代价函数计算收敛判断我用了双条件medoid不再变化或者达到最大迭代次数。每次迭代计算一次代价函数cost sum(min(D, [], 2)); % 所有样本到最近medoid的距离和如果新老medoid集合完全一致说明已经稳定直接break。如果代价函数的变化小于某个阈值比如1e-6也可以提前终止。我在代码里用一个isequal判断if isequal(medoid_idx, prev_medoid_idx) break; end prev_medoid_idx medoid_idx; cost_history [cost_history; cost];有人会问为什么不用代价函数变化阈值因为medoid是离散选择代价变化不一定是平滑下降的可能出现两个不同medoid集合代价几乎相同。用medoid索引变化的判断更直觉更稳。3.6 完整源代码中文注释版下面给出完整的函数实现可以直接保存为myKmedoids.m使用function [idx, medoid_idx, cost_history] myKmedoids(X, k, maxIter, distType, initMethod) % myKmedoids 基于PAM思路的k-medoids聚类实现 % 输入 % X - 样本特征矩阵每一行是一个样本 % k - 聚类簇数 % maxIter - 最大迭代次数默认100 % distType - 距离度量euclidean或manhattan默认euclidean % initMethod - random或maxdist默认random % 输出 % idx - 每个样本所属簇的标签 % medoid_idx - 每个簇的medoid在原始数据中的行号 % cost_history - 每轮迭代的代价函数值 % 参数默认值处理 if nargin 3 maxIter 100; end if nargin 4 distType euclidean; end if nargin 5 initMethod random; end n size(X, 1); % 初始化medoid if strcmpi(initMethod, random) medoid_idx randperm(n, k); elseif strcmpi(initMethod, maxdist) medoid_idx randi([1, n], 1, 1); for j 2:k D pdist2(X, X(medoid_idx, :), distType); minDist min(D, [], 2); [~, idx_new] max(minDist); medoid_idx [medoid_idx; idx_new]; end medoid_idx medoid_idx; end cost_history []; prev_medoid_idx medoid_idx; for iter 1:maxIter % 分配每个样本归属最近的medoid D pdist2(X, X(medoid_idx, :), distType); [~, idx] min(D, [], 2); % 更新对每个簇选择总距离最小的样本作为新medoid for j 1:k cluster_idx find(idx j); if isempty(cluster_idx) cluster_idx randi([1, n], 1, 1); idx(cluster_idx) j; end cluster_X X(cluster_idx, :); Dc pdist2(cluster_X, cluster_X, distType); totalDist sum(Dc, 2); [~, bestLocalIdx] min(totalDist); medoid_idx(j) cluster_idx(bestLocalIdx); end % 计算当前代价 D pdist2(X, X(medoid_idx, :), distType); cost sum(min(D, [], 2)); cost_history [cost_history; cost]; % 收敛判断medoid集合是否变化 if isequal(medoid_idx, prev_medoid_idx) break; end prev_medoid_idx medoid_idx; end end这段代码加起来不到60行核心逻辑是清晰的。如果你要对比实验只需要把初始化的随机种子固定rng(42)之类就能复现结果。3.7 复杂度分析与调优技巧k-medoids的每次迭代复杂度是O(knF)其中F是特征维度。更新阶段每个簇内算pdist2矩阵的开销是O(n^2/F)整体在k个簇上更接近O(n^2)因此样本量是主要瓶颈。实测数据n10000、k5、特征维度10欧氏距离迭代30轮在我的笔记本上耗时大约8秒。这个量级日常实验完全没问题。如果想提速有一个替代方案更新阶段不需要计算簇内所有样本两两距离只要在簇内每个样本到其他样本的距离时用距离公式直接算不要缓存整矩阵。但MATLAB向量化之后一次性算pdist2往往比循环快这个取舍要在效率上自己测一下。另一个方案是每轮分配只用部分样本估算也就是mini-batch代价是收敛会抖一点但大样本时值得。4. 图形绘制与分类结果可视化4.1 二维散点图绘制与medoid标记聚类做完不画图等于白做。我的标准可视化代码是这样的figure; gscatter(X(:,1), X(:,2), idx); hold on; plot(X(medoid_idx, 1), X(medoid_idx, 2), kx, MarkerSize, 12, LineWidth, 2); hold off; legend(簇1, 簇2, 簇3, Medoid); title(k-medoids聚类结果);gscatter会自动给不同类别分配颜色省去手动配色。plot用黑色叉号标记medoid和散点区分明显。这里有个小细节如果X的特征维度不止2画图前先做主成分分析降维[coeff, score] pca(X); X2D score(:, 1:2); % 取前两个主成分 gscatter(X2D(:,1), X2D(:,2), idx);降维后的可视化只用于展示聚类过程仍用原始高维特征。4.2 轮廓图评估聚类质量聚类质量不能光靠肉眼轮廓系数silhouette是常用指标。MATLAB内置silhouette函数直接能用figure; silhouette(X, idx, distType); title(K-medoids聚类轮廓图);轮廓系数范围是[-1,1]越接近1说明样本离自己簇内近、离别的簇远。绘制出的图如果大部分样本的轮廓值在0.5以上说明聚类结构清晰如果很多值接近0甚至负数说明有样本被分错了位置。可以在不同k值下跑多次观察轮廓均值的变化趋势来选k。这个步骤建议固定随机种子再比较否则每次结果可能不同。4.3 代价函数收敛曲线每次运行保存的cost_history可以直接画收敛曲线判断算法是否稳定figure; plot(1:length(cost_history), cost_history, -o, LineWidth, 1.5); xlabel(迭代轮次); ylabel(代价函数值); title(k-medoids收敛曲线); grid on;正常情况下代价函数单调下降并趋于平稳。如果曲线出现反复震荡多半是初始化太差或者k值不合适。如果迭代20轮代价还在明显下降说明maxIter设置小了建议调大。4.4 中文注释乱码与绘图中文显示问题这里专门说一下MATLAB中文显示的两个坑。第一个是源代码里的中文注释保存后乱码常见于直接用记事本保存的.m文件MATLAB默认按系统编码读取一旦文件编码是UTF-8而系统区域是GBK中文注释就全乱。解决办法在MATLAB编辑器里设置文件编码为UTF-8或直接在preferences中调整另外新建脚本用MATLAB编辑器保存而不是外部编辑器。如果代码已经乱码了用MATLAB的edit菜单里重新以正确编码打开一次就能恢复。第二个是绘图的标题和图例中文显示成方框这是字体问题。绘图之前加一行set(0, DefaultAxesFontName, SimHei);或者用set(gca, FontName, SimHei)。Windows一般有SimHei黑体和Microsoft YaHei微软雅黑macOS用PingFang SC。加上这个设置标题里的中文才能正常显示。5. 常见问题与排查经验速查5.1 聚类结果不稳定怎么办k-medoids对初始值敏感每次运行结果可能都不一样。解决思路有三个固定随机种子做复现运行前执行rng(42)保证每次结果一致。多次运行取最优循环运行比如20次保留代价最小的结果。改用maxdist初始化稳定性明显好于随机初始化。如果多次运行后代价差异仍然很大说明数据本身聚类结构不清晰或者k值选得不合适。这时候先画一下数据分布或者尝试不同k值。5.2 数据导入报错与乱码问题汇总数据导入最常见的几个报错我在表里汇总一下问题原因解决办法readmatrix读CSV后全是NaN文件里有文本列或表头没跳过用readtable读再用table2array或者readmatrix指定NumHeaderLines中文表头导入后乱码CSV编码和系统编码不一致用UTF-8编码保存CSV或用fopen指定编码读取Excel里数字带千分位逗号变成文本Excel单元格格式是文本在Excel里先转成数值格式再导出xlsread提示找不到文件路径含中文改工作目录到数据目录用相对路径读取矩阵维度不一致数据里有空行或空列导入前删除空行空列检查最后一行是否完整以我自己的经验CSV导入踩坑率最高。有些软件导出CSV时会带BOM头MATLAB某些版本读BOM会出现第一个变量名多出字符。解决办法是用fopenfgetl读第一行手动清洗或者用readtable的PreserveVariableNames选项。实在不行就用Excel另存为xlsx格式readmatrix对付xlsx更省心。5.3 medoid更新后簇变空的情况在更新步骤中如果某个簇的成员在重新分配后为空需要处理否则下一步会报错。我的代码里已经加了随机填充逻辑if isempty(cluster_idx) cluster_idx randi([1, n], 1, 1); idx(cluster_idx) j; end更严谨的做法是找到当前距离矩阵中离所有medoid最远的样本把它强制归入空簇。但随机填充在多数情况下也能用因为下一步的分配会重新调整归属。还有一种情况簇内只有一个样本它的medoid就是本身距离为0更新后不变正常。5.4 高维数据下的可视化与评估技巧高维数据没法直接画散点图我的做法是先画距离矩阵热力图。把样本按聚类标签排序后用imagesc绘制样本间距离矩阵能直观看到块状结构。代码片段Dfull pdist2(X, X, distType); [~, order] sort(idx); imagesc(Dfull(order, order)); colormap(jet); colorbar;对角线附近的方块颜色越深说明簇内紧凑方块之间的边界颜色越亮说明簇间分离明显。这个方法比散点图更能反映真实聚类质量。另外还可以用簇间距离/簇内距离的比值Davies-Bouldin指数做定量评估。6. 完整项目文件结构与调用示例6.1 文件清单与各模块职责一个完整的k-medoids实验项目我建议这么组织文件project/ ├── myKmedoids.m % 主聚类函数 ├── loadData.m % 数据导入封装 ├── visualizeResult.m % 可视化封装 ├── evaluateCluster.m % 轮廓系数等评估 ├── runExperiment.m % 主脚本调用所有模块 └── data/ └── iris.csv % 实验数据把导入、聚类、可视化、评估拆成独立函数的好处是换数据集时只需要改loadData.m换算法时只改runExperiment.m里的调用可视化可以复用。6.2 完整调用示例鸢尾花数据集下面是一个完整的演示脚本直接用runExperiment.m跑通整个流程%% 数据导入 data readmatrix(data/iris.csv, NumHeaderLines, 1); features data(:, 1:4); % 标准化z-score mu mean(features); sigma std(features); features_std (features - mu) ./ sigma; %% 聚类运行 rng(42); k 3; [idx, medoid_idx, cost_history] myKmedoids(features_std, k, 100, euclidean, maxdist); %% 结果展示 fprintf(聚类标签分布\n); tabulate(idx); fprintf(Medoid样本行号); disp(medoid_idx); %% 可视化 visualizeResult(features_std, idx, medoid_idx, cost_history);其中visualizeResult函数里画三张图散点图、代价收敛曲线、轮廓图。这套流程在新数据集上只需要改数据和k值其他代码不用动。6.3 后续扩展k值选择、并行与对比实验如果你要把这个代码用在正式实验里可以扩展三个方向自适应k值写一个循环从k2到k10跑聚类每次计算轮廓系数均值画折线图选拐点。k-medoids的轮廓系数计算可以直接复用silhouette函数代价不太高。并行加速不同k值的聚类任务相互独立用parfor替代for循环跑多次实验。注意parfor里不能共用随机流需要每个worker独立设置随机种子。与k-means对比在相同数据和相同k值下跑MATLAB内置kmeans函数比较代价、迭代次数、轮廓系数。这个对比表格放在论文里是很有说服力的实验数据。我个人在实际操作中最深的一个体会是k-medoids的代码逻辑不难难在调试时对“初始化敏感”有心理预期。第一次跑出不太好的结果不要急着怀疑代码先加rng固定种子、换maxdist初始化、多跑几次取最优。另外写代码时务必把中文注释的编码问题提前处理掉否则图里和注释里一堆乱码会浪费你半小时排查时间。最后再分享一个小技巧画聚类图时如果两类样本重叠严重可以在gscatter里加上透明度参数Alpha, 0.5点重叠的地方颜色变深比默认散点更容易看出边界走势。这套代码和思路我一直在用希望也能帮你省事。