做聚类分析的时候k-means通常是新手接触的第一个算法但真正拿到带噪声的真实数据你会发现k-means经常被离群点带偏算出来的“中心”根本不是实际存在的样本。k-medoids是k-means的稳健替代方案它用簇内最靠近中心的实际样本点medoid来代表整个簇而不是用均值。最近我把这套算法用MATLAB完整实现了一遍代码带中文注释通过标识medoid直接确定每个样本的类别整个过程非常清晰。这篇文章就把实现细节、参数选择和踩过的坑一起分享出来适合正在做聚类分析、需要在MATLAB里快速落地k-medoids的读者参考。1. 为什么要用k-medoids从k-means的痛点说起1.1 k-means的均值陷阱k-means的核心思想是“聚类中心取均值”这个逻辑在数据分布比较均衡时表现不错一旦数据里出现离群点问题就来了。离群点距离正常样本很远均值会被它强行拉过去导致整个簇中心偏移最后聚类边界也跟着变形。我遇到过一组客户分群数据某个客户消费异常高硬生生把一整个簇的中心拽走分出来的群完全不符合业务直觉。根本原因在于均值本身不是一个“实际存在”的点它只是一个数学意义上的最小二乘解。当数据分布不是球形、噪声明显时均值对极端值太敏感。而k-medoids的做法是每个簇的中心必须是簇内某个真实样本称为medoid中心样本。这样中心一定落在数据区域内部不会凭空漂移到没有样本的地方。1.2 medoid的概念与稳定性优势medoid可以简单理解为“簇内最中心的那个样本”。判断标准是对于簇内所有样本点计算它到簇内其他样本的距离总和距离总和最小的那个样本就是medoid。这个定义决定了medoid天然是数据集中真实存在的点有业务含义也方便后续解释。比如做用户画像聚类k-means给出的“平均用户”可能是个既不出生也不消费的虚构用户而k-medoids选出的medoid会是某个真实存在的用户可以直接查看他/她的特征。这种可解释性在业务场景里特别重要。从稳定性角度说k-medoids对离群点和噪声的容忍度远高于k-means。因为medoid是在现有样本里挑出来的离群点很难“拖走”整个簇中心。代价是计算量变大每次更新medoid都要遍历簇内所有样本比较距离总和在大数据集上会明显变慢。2. 算法流程与MATLAB实现思路2.1 k-medoids核心步骤k-medoids的迭代流程和k-means很像但更新中心的方式不同。完整步骤如下确定聚类数k从数据集中随机选择k个样本作为初始medoid。计算每个样本到所有medoid的距离将其分配到距离最近的medoid所属的簇。在每个簇内重新寻找新的medoid。具体做法是遍历簇内每个样本点计算它到簇内其他样本的距离总和选取总距离最小的那个点作为新的medoid。重复步骤2和3直到medoid不再变化或者达到最大迭代次数。可以看到步骤2实际上就是“通过标识medoid确定类别”的过程。每个样本只需要比较它和k个medoid的距离哪个medoid离它最近就归属到哪个类别。最终输出的每个样本的类别标签完全由当前迭代中标识出来的medoid决定。2.2 代码框架设计距离矩阵、迭代更新、标识分类在MATLAB里实现k-medoids我建议按模块来写不要一股脑塞在一个脚本里。核心模块有这么几个距离矩阵计算函数输入数据集X和medoid集合输出每个样本到每个medoid的距离矩阵。这一步如果数据集不大可以直接用pdist2函数效率很高。初始化函数随机从样本中抽取k个行索引作为初始medoids。分配函数根据距离矩阵为每个样本找到最近的medoid得到类别标签。更新medoid函数对每个簇计算簇内所有样本两两之间的距离总和重新选出新的medoid。主程序控制迭代循环检查是否收敛输出最终类别标签和medoid位置。这种模块化设计的优点是方便调试比如你想换距离度量只需要改距离计算那一小块代码。我在实际写代码时习惯把距离矩阵先完整算一次再用索引方式查询避免在循环里反复调用pdist2能省不少事。3. 完整MATLAB源代码及中文注释3.1 主函数 kmedoids_demo.m 源码下面是一段可以独立运行的MATLAB源代码用中文注释详细说明了每一步的作用。整个代码包含了数据生成、k-medoids聚类、可视化三个部分你直接复制到MATLAB脚本文件里就能跑通。%% k-medoids聚类分析演示脚本 % 功能使用k-medoids方法对二维样本数据进行聚类 % 通过标识medoid确定每个样本的类别 % 作者分享者 % 版本MATLAB R2021a 及以上 clear; clc; close all; rng(42); % 固定随机种子保证结果可复现 %% 1. 生成测试数据3个高斯簇带少量噪声 % 这里生成了3个中心不同的高斯簇模拟实际聚类场景 n1 60; n2 50; n3 40; X1 randn(n1,2) * 0.6 [2, 2]; X2 randn(n2,2) * 0.5 [-2, 3]; X3 randn(n3,2) * 0.7 [1, -2]; % 添加几个离群点看看k-medoids的稳健性 outliers [8, 8; -6, -6; 5, -5]; X [X1; X2; X3; outliers]; true_label [ones(n1,1); 2*ones(n2,1); 3*ones(n3,1); 4*ones(3,1)]; % 绘制原始数据分布 figure; subplot(1,2,1); plot(X(:,1), X(:,2), k., MarkerSize, 10); title(原始数据分布); axis equal; grid on; %% 2. 调用k-medoids聚类函数 k 4; % 这里把离群点也算作一个簇可以直观看到medoid会选择到真实样本上 max_iters 100; [labels, medoid_idx, medoids] my_kmedoids(X, k, max_iters); %% 3. 可视化聚类结果并标注medoid subplot(1,2,2); colors lines(k); hold on; for i 1:k cluster_points X(labels i, :); plot(cluster_points(:,1), cluster_points(:,2), ., Color, colors(i,:), MarkerSize, 12); end % 把medoid用五角星标出来 plot(medoids(:,1), medoids(:,2), kp, MarkerSize, 20, MarkerFaceColor, y); legend_str cell(1, k); for i 1:k legend_str{i} [簇, num2str(i)]; end legend(legend_str, Location, best); title(k-medoids聚类结果五角星为medoid); axis equal; grid on; hold off; disp(聚类完成medoid对应的样本索引为); disp(medoid_idx);3.2 核心函数 my_kmedoids.m 的源码与讲解主脚本调用了子函数my_kmedoids这段代码是算法的核心。我在每一段都加了中文注释重点说明“通过标识medoid确定类别”的逻辑。function [labels, medoid_idx, medoids] my_kmedoids(X, k, max_iters) % MY_KMEDOIDS 使用k-medoids算法进行聚类 % 输入 % X - 数据矩阵每一行是一个样本每一列是一个特征 % k - 聚类簇数 % max_iters- 最大迭代次数 % 输出 % labels - N×1向量每个样本所属的类别编号1~k % medoid_idx- k×1向量medoid在X中的行索引 % medoids - k×2矩阵medoid的坐标 n size(X, 1); % 样本总数 dist_matrix pdist2(X, X, euclidean); % 预计算所有样本两两之间的距离 %% 1. 初始化随机选择k个样本作为初始medoid % 这里用随机选择更稳健的做法可以先用k-means的思想后面会提到 init_idx randperm(n, k); medoid_idx init_idx(:); %% 2. 迭代更新 for iter 1:max_iters % 计算每个样本到当前所有medoid的距离 % dist_to_medoids矩阵维度是 n×k第(i,j)个元素表示第i个样本到第j个medoid的距离 dist_to_medoids dist_matrix(:, medoid_idx); % 通过标识medoid确定类别 % 每个样本只取距离最小对应的那个medoid索引min函数返回两个输出 % 第二个输出就是最小值的列位置也就是类别编号 [~, labels] min(dist_to_medoids, [], 2); % 保存旧的medoid索引便于判断是否收敛 old_medoid_idx medoid_idx; %% 3. 在每个簇内重新寻找新的medoid % 簇内距离总和最小的样本点就是新的medoid for j 1:k % 找到当前属于第j个簇的所有样本索引 cluster_members find(labels j); if isempty(cluster_members) continue; % 如果某个簇为空保持原medoid不变实际很少见可加异常处理 end % 提取簇内样本两两之间的距离子矩阵 sub_dist dist_matrix(cluster_members, cluster_members); % 对子矩阵按行求和得到每个样本到簇内其他样本的距离总和 sum_dist sum(sub_dist, 2); % 距离总和最小的那个样本点就是新的medoid [~, min_pos] min(sum_dist); medoid_idx(j) cluster_members(min_pos); end % 收敛判断如果medoid索引不再变化直接退出循环 if isequal(old_medoid_idx, medoid_idx) % fprintf(已经收敛迭代次数%d\n, iter); break; end end % 返回medoid对应的坐标 medoids X(medoid_idx, :); end3.3 这段代码的几个关键设计点代码里最妙的部分是[~, labels] min(dist_to_medoids, [], 2);这一行。它直接利用了MATLAB的min函数在按行比较每个样本到所有medoid的距离后把最小值的列索引作为类别标签。Column index就是medoid的顺序编号比如第2列最小说明该样本离第2个medoid最近类别就是2。这样就把“通过标识medoid确定类别”落到了实处。另一个需要注意的是距离矩阵的预计算。我用pdist2(X, X, euclidean)一次性算出所有样本之间的两两距离之后每次迭代都只是做矩阵索引避免了重复计算。对于几千行数据这种写法效率足够高。如果数据量到了几万行全距离矩阵会占用大量内存到时候需要改成按需计算或者分块处理。4. 实操演示与结果解读4.1 跑通代码从生成数据到输出类别我把上一节的两段代码保存成kmedoids_demo.m和my_kmedoids.m放在同一目录下直接运行kmedoids_demo。数据生成了三个高斯簇和三个离群点聚类数k设为4这样离群点也被单独分出来便于观察medoid的标识效果。运行结束后命令窗口会输出medoid在原始数据中的索引。你可以打开变量查看会发现medoid一定是数据里的某个原始样本而不是像k-means那样得到一个坐标为小数的新点。比如运行一次可能得到medoid_idx [23; 76; 112; 148]代表第4个簇的medoid是第148个样本点。可视化图形里三个簇的颜色各不相同每个簇内的样本点落在一起而用黄色五角星标出的medoid恰好位于每个簇的密集区域中心。有意思的是那个离群点聚成的簇它的medoid就是离群点本身不是周围区域的平均值。这里能直观感受到k-medoids的稳健性。4.2 通过标识medoid确定类别的实际含义聚类完成后类别标签labels就是我们最终要的结果。它和第1节画图时输入的true_label不完全一致因为聚类本身不关心我们离散群点还是正常簇它只是根据距离把样本分组。这里labels的取值范围是1到k每个值对应一个medoid索引。“通过标识medoid确定类别”这句话可以这样理解聚类模型本质上就是一组medoid你拿到一个新的未知样本时只需要计算它到每个medoid的距离找最近的medoid编号就能判定它属于哪个类别。这比k-means需要在内存里保存均值中心更方便因为medoid可以直接指向某个已知样本做异常检测或者解释性分析时你可以直接围绕这个样本展开。我在实际项目中经常用这个特性来做数据清洗。例如聚类出用户分群后把每个簇的medoid打印出来人工一看就知道这个簇代表什么类型的用户群体。这样一来聚类结果不是黑盒业务同事也容易接受。4.3 不同距离度量的影响与选择我的代码里默认用的是欧氏距离你把pdist2(X, X, euclidean)换成cityblock、cosine或mahalanobis就能改变聚类行为。选距离度量时要结合数据特征。比如用户购买行为数据不同维度量纲差异大直接用欧氏距离会被数值大的维度主导这时候标准化的作用比换距离度量更明显。文本数据用余弦距离更合理因为只看方向、不看长度。k-medoids的优势在于只要距离定义好算法流程不需要改动它会自动根据距离矩阵寻找每个簇内最中心的样本。我建议在动手之前先画出数据分布或者做一下相关性分析判断距离度量的选择。如果数据维度高、无法直观观察可以多跑几个距离度量用轮廓系数辅助判断哪个结果更合理。5. 常见问题与排查技巧实录5.1 中文注释乱码的问题很多人在MATLAB 2023a或更早版本里写中文注释保存后重新打开会出现乱码。这是因为MATLAB默认编码和操作系统编码不一致造成的。我自己的解决办法是打开MATLAB的“预设项”在“常规”里找到“文本编码”直接改成UTF-8。改完以后再写中文注释就稳定了。如果已经有乱码文件可以先用纯文本编辑器把文件编码转成UTF-8再让MATLAB重新读取。需要注意转换前最好备份因为有些MATLAB版本对文件头有自己的处理方式转完可能注释还是乱这时可以让同事帮忙看下是不是编码被二次转换了。实在不行像我代码里一样用百分比注释配合简单英文说明至少能保证可读性。5.2 初始medoid的选择对结果的影响k-medoids对初始点很敏感随机初始化容易陷入局部最优。我在代码里用的randperm是基础做法但实际使用中如果你跑多次会发现最后聚类结果可能不同有时簇分配完全不一样。有两个改进方向。一个是多次随机初始化跑若干次后保留目标函数簇内总距离最小的那次结果。另一个是参考k-means的思想让初始medoid尽量分散。实现起来不复杂先随机选第一个medoid然后计算每个样本到当前medoid集合的最小距离用归一化后的距离作为概率权重按权重随机选下一个样本。这么做能显著提升稳定性。我在自己的代码里加了个简单的循环重复运行20次每次随机初始化后都记录最终的簇内距离总和最后挑最小的一次。在中小数据集上很实用虽然计算量大一点但值得。5.3 如何自动确定聚类数Kk值往往是最难定的参数。我的经验是先画出肘部图横轴是k从1到10纵轴是全部样本到各自簇medoid的距离总和。k增加时总距离会下降曲线会出现一个拐点这个拐点对应的k相对合理。另外可以用轮廓系数silhouette辅助判断。MATLAB自带silhouette函数输入样本和标签就能输出每个样本的轮廓值平均值越高说明聚类结构越清晰。我在演示中把k设成4正是因为原始数据里包含离群点但如果你只想聚类三个正常簇设k3跑出来的轮廓系数会更高。实际操作时我会结合业务需求来定k。比如业务方明确说想分3到5组那就分别跑k3、4、5比较轮廓系数和结果的可解释性而不是死磕一个“最优”数值。6. 个人实操体会与扩展建议我自己用了k-medoids差不多两年最大的感受是它比k-means“沉得住气”。数据里有少量异常值或者类别不平衡时medoid不会剧烈漂移业务解释也顺畅。但它的确比k-means慢不少尤其样本量上万每次更新medoid要计算簇内两两距离时间开销会明显上涨。如果数据量很大可以考虑先用k-means做粗聚类得到每个簇的中心点再在这些中心点附近寻找距离所有样本最近的真实样本作为候选medoid这样能大幅减少计算量。或者直接调用kmedoids函数——从R2018b开始MATLAB的Statistics and Machine Learning Toolbox里自带这个算法我的代码主要是帮助理解内部逻辑你在生产环境里完全可以优先使用官方函数。最后再分享一个小技巧如果你的数据维度很高做k-medoids前先用PCA降到二维或三维可视化看看真实分布这对选择k和距离度量特别有帮助。我就是在一次高维用户分群的项目里先降维到三维发现数据有明显长条分布于是把距离度量从欧氏距离换成了马氏距离聚类质量立刻提升。算法本身是死的你对数据的理解才是决定结果质量的关键。
