先说结论k-medoids这个算法凡是写论文或者做聚类对比实验的人迟早会想自己去实现一遍。我最近在做用户分群时就被k-means的均值中心坑了一回——一簇样本里头掺了两个离群点簇中心直接飘到异常值那边结果画出来的散点图简直没法看。后面换成k-medoids用medoid取代mean情况才好转。这篇文章我把我整理好的MATLAB实现放出来带中文注释从数据导入、PAM迭代到图形绘制每一步都讲清楚为什么这么写你拿去改改路径就能跑通自己的数据。先交代一下适用人群。如果你是在校学生需要一份能交作业、能复现论文的聚类代码如果你是工程师想快速验证k-medoids在自己的数据上是否比k-means更稳或者你就是想把PAM算法内部那套“试探替换”的机制彻底搞明白——这篇都能帮你省下不少时间。代码本身不复杂核心循环加起来不到三十行但我把容易出错的点、标准实现的取舍、以及数据导入和画图的坑都放在一起了比单纯甩一份源码有用得多。1. 为什么总有人要手写k-medoidsk-means的最大短板1.1 均值中心对离群点有多脆弱k-means聚类的中心是簇内样本的算术平均也就是centroid。均值这个东西有一个特点特别容易被极端值牵拉。假设一簇数据大部分集中在2到3附近但有一个离群点在20的位置按照最小化SSE的目标簇中心会被往离群点方向拉一大截。对这个簇里的正常样本来说它们距离中心的平均距离都会变大聚类结果就开始偏离直觉了。这个问题的本质在于k-means的目标函数。它试图最小化的其实是所有样本到它们所属簇中心的欧氏距离平方和平方项让离群点的误差贡献被进一步放大。有人会说那我做聚类之前先把离群点清洗掉不就行了问题是很多场景里“离群点”本身就是有业务含义的比如反欺诈场景里的异常交易、传感器数据里的瞬时毛刺、用户在某个时刻的极端行为。你把这些点直接删掉很可能会丢掉最想发现的那部分信息。我打个比方一个班同学的平均身高可能是175厘米但这175厘米可能对应班里任何一个真实同学也可能谁都不对应。k-medoids取的是“班里实际存在的某个同学”这个人的身高最能代表这一群人而不是一个虚构出来的数字。在聚类任务里用真实样本点当中心比用均值当中心要稳得多。1.2 medoid怎么解决这个问题medoid的定义是簇内“最接近所有其他样本”的那一个样本点。它不涉及平均而是直接比较簇内样本两两之间的距离。k-medoids的目标函数通常写成min Σ distance(x_i, medoid_of_its_cluster)跟k-means相比区别就一个字用真实样本点作为中心。出现极端样本时由于medoid必须是簇内实际存在的点它天然就是一种中位数性质的估计不会像均值那样被极端值大量牵拉。因此在包含噪声或离群点的场景中k-medoids的簇骨架往往更稳中心点的解释性也更强。比如落地到业务上你告诉运营“这批客户的核心画像就是ID 537这个人”比说“这批客户的平均特征是年龄35.2岁、消费频次7.8次”要直观得多。1.3 MATLAB自带kmedoids为什么还要手写MATLAB的Statistics Toolbox里确实有kmedoids函数一行[idx,C] kmedoids(X,k)就能跑。但我必须说自带函数适合“直接出结果”不适合以下几种情况你想用自定义距离度量比如DTW距离做时间序列聚类自定义距离矩阵直接喂给自带函数就很费劲。你想观察每次迭代的代价变化和交换路径自带函数给出的是最终结果中间的轨迹被封装死了。你想改初始化方式试试不同策略对聚类结果的影响自带函数开放不了这些自由度。你压根没装Statistics Toolbox。我自己平时做仿真需要输出中间可视化的收敛曲线自带函数给不了这些所以手写一份放在自己的工具箱里随时拿来扩展是很值得的。2. PAM算法核心medoid替换背后的代价计算逻辑2.1 初始化这步定基调PAMPartitioning Around Medoids是k-medoids最经典的实现。它的第一步要选出k个初始medoid。教材里的标准做法是贪心BUILD每一轮选一个新的medoid时计算当前所有候选点替换进去后总代价下降了多少选下降最多的那个点加进来。这个策略的效果好但写起来稍微绕一些。实际工程里大多数人直接用随机抽样randsample(N,k)逻辑简单、结果可复现。区别在于贪心BUILD得到的初始解通常总代价更低后面迭代次数更少随机抽样代码短但更容易掉进局部最优。我的建议是教学演示用随机因为逻辑直观正式做实验时多跑几次类似内置函数的Replicates20选代价最小的那一次结果来报告。2.2 分配步骤与总代价的定义选定k个medoid之后每个样本点归到离它最近的medoid这步没有任何争议。设样本到medoid的距离矩阵为D那么min(D, [], 2)取出的是每个样本到最近medoid的距离把这些距离加起来就是当前目标代价currentCost。注意一个写代码时的习惯问题这里用全量重算距离矩阵的方式每次交换尝试都重新算代价实现简单但会重复算很多距离。数据量小没关系数据量大了就必须优化这个后面专门讲。2.3 交换步骤替换谁、怎么判断换不换交换阶段是PAM的精华。标准做法是对每一个medoid点遍历所有非medoid候选点试探性地用候选点替换它重新计算全局总代价。如果新代价比当前低就接受这次替换否则放弃。这个“试探替换”的过程本质上是在解空间里做局部搜索每次只做一步最小的调整反复迭代直到找不到更优解。有个细节值得说明最严格的PAM会遍历所有非medoid点来搜索候选。我在下面代码里做了一点简化只在当前簇内部找候选点。为什么可以这么简化因为一个点如果跟某个簇的样本距离都很远它基本不可能成为那个簇的好medoid拿来换了大概率代价更高。当然这只是经验性的简化理论上会丢掉一些全局最优的候选。如果你要严谨版本把候选取成setdiff(1:N, medoidIdx)就行遍历范围从簇大小变成N-k代价会大不少。2.4 收敛判定要注意的细节PAM收敛的标准是在一整轮完整的交换扫描中没有任何一次交换能降低代价。注意代码里如果发现某次交换有用就立刻更新medoid并且继续往下走那么簇分配在扫描过程中会变化这跟“一轮扫描结束后再统一应用所有有效交换”的做法结果会有细微差别。但实践下来的差别通常很小用improved标志已经足够循环里直接break就行。收敛后输出的medoid就是当前目标函数下的局部最优解你可以对比一下每次迭代时的costHistory一般前面两三轮下降很快后面就稳定了。3. 带中文注释的MATLAB实现数据导入、PAM循环、画图一条龙3.1 数据导入的几种写法先给数据导入的代码。我平时收数据优先用readmatrix和readtable从R2019a开始这两个函数性能就很好别再用了csvread那种老接口。MATLAB内置的鸢尾花数据集很适合演示一共150条样本、4个特征、3个类别聚类结果容易看明白。%% 数据导入可选方案 % 方案1MATLAB内置数据集鸢尾花150x4 load fisheriris; rawX meas; % 150x4的double矩阵 trueLabel species; % 150x1的cell真实类别用于事后验证 % 方案2纯数值CSV第一行开始就是数据 % X readmatrix(myData.csv); % 方案3带表头的CSV或Excel % tab readtable(myData.csv); % X table2array(tab(:, 2:end)); % 去掉ID列/标签列按需选择这里有个特别常见的坑readtable读出来的东西是table类型不能直接丢给pdist2去算距离必须先用table2array转成double矩阵。如果表里混了字符串列table2array会生成cell数组这时候需要单独去掉文本列或者做标签编码。我见过太多人栽在这一步报错信息类似于“Undefined function pdist2 for input arguments of type table”。另外如果数据文件里存在缺失值最好在读入后先处理% 查看每一列的缺失数量 sum(ismissing(X)); % 删除包含缺失值的行或者用插值填充 X rmmissing(X);3.2 PAM主循环完整代码下面是核心部分PAM的分配、交换、收敛判断都在里面注释我写得很细。%% k-medoids聚类PAM算法主循环 - 带中文注释 clear; clc; close all; % 数据准备 load fisheriris; rawX meas; % 原始数据后面解释medoid时用 X zscore(rawX); % 标准化每列零均值、单位方差避免量纲主导距离 k 3; maxIter 100; % 最大迭代轮数防止死循环 N size(X, 1); % 样本总数 % 初始化 rng(42); % 固定随机种子结果可复现 medoidIdx randsample(N, k); % 随机选k个样本作为初始medoid costHistory zeros(maxIter, 1); % 记录每轮代价方便画收敛曲线 % PAM迭代 for iter 1:maxIter % ----- 分配步骤计算所有样本到当前medoid的距离 ----- D pdist2(X, X(medoidIdx, :), euclidean); [minD, assign] min(D, [], 2); % assign(i)是第i个样本的簇编号 currentCost sum(minD); % 目标函数值所有样本到最近medoid距离之和 costHistory(iter) currentCost; % ----- 交换步骤尝试在每个簇内部找一个更好的medoid ----- improved false; for j 1:k inCluster find(assign j); % 当前属于簇j的样本 for cand inCluster if cand medoidIdx(j) continue; end newMedoid medoidIdx; newMedoid(j) cand; % 用候选点替换第j个medoid Dnew pdist2(X, X(newMedoid, :), euclidean); newCost sum(min(Dnew, [], 2)); if newCost currentCost % 如果新代价更小就接受这次替换 medoidIdx newMedoid; currentCost newCost; improved true; end end end fprintf(第 %2d 轮完成总代价 %.4f\n, iter, currentCost); % ----- 收敛判定一轮扫描中没有任何交换被接受 ----- if ~improved break; end end % 最终聚类结果 D pdist2(X, X(medoidIdx, :), euclidean); [~, finalAssign] min(D, [], 2); fprintf(收敛于第 %d 轮medoid样本编号: %s\n, ... iter, mat2str(medoidIdx));这段代码有三个地方值得单独解释。第一pdist2(X, X(medoidIdx,:))返回的是一个N行k列的矩阵第i行第j列就是样本i到第j个medoid的距离。min(D, [], 2)这个写法初学者经常忘记第二个输出参数但它恰恰是关键——第二个输出assign记录了每个样本最近的是哪个medoid后面交换阶段要用它来定位每个簇的成员。如果不显式接收第二个输出MATLAB会自动丢掉它代码就会标错。第二交换阶段的候选范围。我上面用的是“只在当前簇内找候选点”的简化版本刚才第2.3节说过这个取舍。如果你做成严格PAM把inCluster换成setdiff(1:N, medoidIdx)就行但循环量会从簇大小变成N-k小数据无所谓大数据会明显拖慢速度。第三如果你没有Statistics Toolboxpdist2会报错。替代方案也很简单自己写一个算欧氏距离矩阵的函数就行function D myEuclideanDist(X, Y) % 自写欧氏距离矩阵避免依赖Statistics Toolbox % X: nxd, Y: mxd返回值是n x m n size(X, 1); m size(Y, 1); D zeros(n, m); for i 1:n D(i, :) sqrt(sum((X(i, :) - Y).^2, 2)); end end把代码里所有pdist2换成myEuclideanDist核心逻辑依然成立。3.3 聚类结果可视化代码接下来说图形绘制。我一般会画三类图分簇散点图、原始尺度下的medoid展示、以及收敛曲线。散点图注意把medoid突出标记出来这是中文注释版源码里最出彩的部分。%% 图形绘制1分簇散点图以第一、二列特征投影 figure(Color, w); colors [0.85 0.33 0.10; 0.00 0.45 0.74; 0.47 0.67 0.19; 0.49 0.18 0.56]; for j 1:k idx (finalAssign j); scatter(X(idx,1), X(idx,2), 40, colors(j,:), filled); hold on; end plot(X(medoidIdx,1), X(medoidIdx,2), kp, ... MarkerSize, 20, MarkerFaceColor, y); legend([arrayfun((j) sprintf(簇 %d, j), 1:k, ... UniformOutput, false), {Medoid}], Location, best); xlabel(特征1); ylabel(特征2); title(k-medoids聚类结果标准化空间); grid on; %% 图形绘制2原始尺度下的聚类结果 figure(Color, w); for j 1:k idx (finalAssign j); scatter(rawX(idx,1), rawX(idx,2), 40, colors(j,:), filled); hold on; end plot(rawX(medoidIdx,1), rawX(medoidIdx,2), kp, ... MarkerSize, 20, MarkerFaceColor, y); xlabel(花萼长度); ylabel(花萼宽度); title(原始尺度的聚类结果); grid on; %% 图形绘制3收敛曲线 costHistory(costHistory 0) []; % 去掉没被迭代到的位置 figure; plot(1:length(costHistory), costHistory, b-o, LineWidth, 1.5); xlabel(迭代轮次); ylabel(总代价); title(PAM收敛过程); grid on;绘制这里有三个容易被忽略的点。一是用循环scatter而不是gscatter是因为gscatter对颜色向量的处理和动态簇数没这么灵活。二是legend的句柄数量和标签数量必须一一对应代码里用arrayfun匹配簇数避免了把标签写死。三是散点图只能展示二维投影鸢尾花是四维数据平面图只能说明大致形状严谨判断还得靠下一节的PCA投影和轮廓系数。4. 实测跑一遍鸢尾花数据集从导入到出图4.1 输出长什么样怎么判断好不好上面代码直接用rng(42)时一般会在第4到第6轮收敛总代价从初始的某个值逐渐稳定下来我在普通笔记本上运行基本秒出结果。散点图上三类样本分得比较清楚两个medoid实际是三个但图里看有可能重叠会落在簇内比较靠中心的位置而不是被离群点拽走。如果你换成自己的数据只需要保证rawX是“每一行一个样本、每一列一个特征”的double矩阵其他代码都不用改动。判断聚类质量不能只看颜色分得开建议加一个Silhouette轮廓系数一行代码就出图figure; silhouette(X, finalAssign);轮廓系数的取值范围在-1到1之间越接近1说明簇内紧致、簇间分离越明显。鸢尾花数据集轮廓系数平均值在0.6左右就已经很理想了。如果这个值低于0.3说明数据可能并不适合做这个k值的聚类要么换个k要么换距离度量。4.2 高维数据先PCA降维再画图如果数据超过3维直接画第一二列特征很容易误导人。正确做法是先用PCA投影到二维再画散点图。PCA这一步不改变聚类结果只是帮你看清楚数据在高维空间里的整体结构。我经常是先算PCA再画图很多原本看起来重叠的簇在PCA平面上会明显分离出趋势。% 对标准化后的X做PCA用前两个主成分得分画图 [coeff, score, latent] pca(X); medoidScore X(medoidIdx, :) * coeff(:, 1:2); % medoid点映射到PC平面 figure(Color, w); for j 1:k idx (finalAssign j); scatter(score(idx,1), score(idx,2), 40, colors(j,:), filled); hold on; end plot(medoidScore(:,1), medoidScore(:,2), kp, ... MarkerSize, 22, MarkerFaceColor, y); xlabel(PC1); ylabel(PC2); title(PCA视图下的k-medoids聚类结果); grid on;4.3 和内置kmedoids的结果做个对比如果你装了Statistics Toolbox可以拿kmedoids对比一下[idx,c] kmedoids(X, 3, Replicates, 20);。大概率你会发现它返回的medoid编号和手写版不一样但轮廓系数和簇标签非常接近。原因很简单目标函数一致初始化随机、局部最优不同最终解落在不同的局部最优附近。这不代表哪一边实现错了聚类本身就是多解的。做论文时建议多初始跑几次选轮廓系数最高的那一次结果来报告这样比较公平。5. 几个中文MATLAB项目里容易翻车的细节5.1 中文注释的编码坑k-medoids的MATLAB源代码中文注释看着简单但实际踩坑的人不少。先说症状你在Windows中文版MATLAB里打开一个.m文件发现代码正常、中文全乱码。原因是编码不一致。R2020a之后MATLAB默认UTF-8R2020a之前的Windows中文版默认GBK。如果你的脚本是UTF-8拿到老版本MATLAB里打开所有中文注释就会变成乱码。反过来也一样。解决办法有几条按优先级排序写代码时统一在MATLAB编辑器里保存成UTF-8或者从VS Code等外部编辑器保存时显式选UTF-8。如果是Git协作在仓库根目录加一个.gitattributes文件写入*.m text eollf encodingutf-8可以避免团队成员之间编码互相污染。团队老项目用GBK那就全程GBK别混用。两种编码混在一个项目里才是最痛苦的。如果你发现注释乱码但代码能跑用edit打开看文件编码再统一转码。别用复制粘贴去救有时候会带进不可见字符。5.2 不标准化的后果你可能想不到k-medoids用距离驱动聚类如果特征量纲差异大量纲大的特征会完全主导距离计算。鸢尾花的花瓣长度和花萼宽度量级还比较接近但你在电商或金融数据里价格、次数、时间戳之间可能差好几个数量级。不做标准化聚类结果基本就是“按价格分堆”其他特征毫无贡献。标准化之后所有特征在距离计算里的权重才公平。还要提醒一点标准化之后你得到的是近似无量纲空间里的medoid。如果想解读业务含义直接用rawX(medoidIdx, :)取原始样本值不要反变换出一个不存在的“缩放后medoid”。因为medoid本来就是簇内真实样本原始尺度的特征值才是可解释的。5.3 距离度量怎么选才不坑pdist2默认欧氏距离但不同距离度量对聚类结果的影响可能非常大。我个人的经验是连续数值型、标准化之后分布比较规整的用欧氏距离没问题数据里偶尔有毛刺噪声的用曼哈顿距离cityblock会更稳一些文本向量或者方向比长度更重要的数据用余弦距离cosine。就一个实际例子来说我处理用户浏览行为序列时欧氏距离分出来的簇几乎都在比数值大小切到余弦距离之后用户的行为模式才被真正分开。在代码里最好用一个变量统一管理距离度量比如distMetric euclidean; % 可选 cityblock / cosine / squaredeuclidean D pdist2(X, X(medoidIdx, :), distMetric);这样调参的时候只改一行不用满文件找。6. 数据量变大以后PAM跑不动怎么办6.1 复杂度到底卡在哪里手写PAM最舒服的是小数据最怕的是大数据。它的主要瓶颈在交换阶段每一轮要尝试大约k*(N-k)次候选替换每次替换又要重新计算N*k的距离矩阵复杂度理论上可以到O(k * N^2 * d)。N1000的时候还行N10000的时候一轮就可能几十亿次距离计算MATLAB里跑完一轮要好几分钟。规范的PAM还有BUILD阶段仅构建初始解也要O(k * N^2 * d)所以数据量变大之后“逐点全量试探”的思路本身就要调整。6.2 实用的几条降复杂度路线第一条路线是换初始化。先用k-means跑一遍拿到簇心后在每个簇内找离簇心最近的样本点作为初始medoid。这个方法速度极快而且比随机初始化更能让后续交换快速收敛。我实测下来大多数数据上聚类结果和标准PAM非常接近工程上完全可接受。第二条路线是CLARA。它是PAM的抽样版本思路是随机抽几百个样本在样本上跑PAM然后把全部数据按最近的medoid分配多次抽样取最优结果。MATLAB里自己写CLARA大概二三十行代码就搞定适合N上万的场景。需要注意的是抽样数太少可能错过代表点一般至少要40 2*k个样本这是文献里常给的经验值。第三条路线是用并行或者分块。遍历候选点时可以用parfor但要注意MATLAB并行池的内存复制问题数据量很大的时候并行带来的收益会被通信开销抵消。我建议先把数据检查一遍如果确实是因为交换次数太多导致太慢优先考虑CLARA路线而不是硬上并行。我把不同数据规模下我的推荐方案整理成一个表方便你直接参考数据规模推荐方案说明N 1000手写PAM清晰可复现随便跑1000 ≤ N 5000手写PAM 固定种子跑多次还能接受注意时间5000 ≤ N 50000CLARA 或 k-means初始化medoid工程上推荐N ≥ 50000换k-means或Mini-batch聚类大数据场景再坚持PAM效率太低这张表是我自己的工程经验不是论文结论。如果你做科研需要严格的PAM保证那就降低N或者用经过优化的开源库实现MATLAB本身不太适合处理超大矩阵的迭代交换。对我来说手写k-medoids最大的收获其实不是代码能跑通而是彻底搞懂了“中心点的选择方式决定了聚类对离群点的态度”这件事。k-means在干净数据上又快又准但数据一旦脏起来k-medoids的鲁棒性就很值得用。后面如果时间充裕建议给上面的代码加一个distMetric参数再加一个自动k选择比如用轮廓系数最大化这套小工具在很多聚类实验里都能复用不用每次重新造轮子。