MATLAB实现k-medoids聚类：PAM算法与medoid标识详解

如果你点进这篇文章估计多半是正在被聚类分析某个环节卡住了。要么是数据里离群点太多跑完k-means一看中心点飘得离谱要么是老板丢给你一堆客户特征让你“分个群看看”但你心里清楚普通聚类出来的“虚拟均值中心”根本没法跟业务解释。这个项目标题里写的k-medoids方法、MATLAB源代码、中文注释、通过标识medoid确定类别恰好就是针对这些痛点的方案。文章会把原理、完整代码、注释逻辑和踩坑经验全部拆开讲适合正在用MATLAB做数据分析的工程师、本科生研究生以及所有需要在业务里做聚类但受不了k-means中心点“不是真实样本”的人。k-medoids的核心就是每个类别用一个真实存在的样本点作为代表这个代表点叫做medoid。对MATLAB用户来说这套算法代码量不大但非常考验对索引、距离矩阵、空簇处理这些细节的把控。这文章里我给的源代码带完整中文注释直接粘到编辑器里就能跑还会重点解释怎么用medoid标识去确定样本类别。1. k-medoids是什么为什么要用它替代k-means1.1 k-medoids与k-means的差异k-means大概是全世界用得最多的聚类算法思路很直接随机初始化K个中心把每个样本分给最近的中心然后重新计算每个簇的均值作为新中心重复到收敛。问题就出在“均值”这两个字上均值是一个数学构造出来的点它不一定存在于你的原始数据里。一旦数据里有几个离群点均值会被拉偏聚类边界就跟着变形。k-medoids的做法完全不同。它的代表点必须是从原数据里挑出来的一个实际样本这个实际样本就叫medoid。更新阶段不是算均值而是在每个簇内部遍历所有点找出“到本簇其他点距离之和最小”的那个点来当新medoid。这样做有两个明显好处第一对离群点和噪声更鲁棒因为离群点很难被选成medoid也不会因为参与均值计算而污染簇心第二结果是可解释的比如你聚类客户得到的每个medoid就是一条真实的客户记录你可以直接看这条记录长什么样然后说“这个群的代表客户是消费频次高、客单价中等的那类人”。k-meansk-medoids代表点是均值不属于原始数据代表点是真实样本叫medoid对离群点敏感对离群点鲁棒计算量相对小每次更新要遍历簇内点计算量更大结果解释性弱能直接拿出代表样本1.2 PAM算法流程与medoid的概念经典实现k-medoids的算法叫PAMPartitioning Around Medoids。它的流程比k-means多一个内层循环但不难理解。第一步随机选K个样本作为初始medoid。第二步分配阶段计算每个样本到K个medoid的距离就近归入对应的簇。第三步更新阶段对每个簇遍历簇内每一个样本计算它到同簇所有其他样本的距离总和谁的总和最小谁就是这个簇的新medoid。第四步重新分配计算目标函数如果所有medoid都不再变化或达到最大迭代次数停止。这里面的目标函数一般写成所有样本到它所属簇medoid的距离之和也就是总绝对误差。因为用的是绝对距离而不是平方距离所以离群点对大数值的惩罚不像k-means那么夸张这就是它更稳的原因。一句话概括medoid的作用它既是一个“簇中心”又是一个可追溯的真实样本ID。你在MATLAB里拿到medoid在原始矩阵中的行号就等于拿到了这个类别的“身份证”。2. 项目设计思路如何用MATLAB实现并“标识medoid确定类别”2.1 整体函数设计标题里特别强调“通过标识medoid确定类别”这句话翻译成代码逻辑就是每个样本的类别由它离哪个medoid最近来决定。所以我设计实现时没有把类别当成一个孤立变量而是让类别标签和medoid索引强绑定。我实现的时候把它封装成一个函数kmedoidsClustering(X, K, ...)输入数据矩阵和聚类数返回结构体。结构体里有几个关键字段idx是每个样本的类别标签取值范围1到KmedoidIdx是每个类别medoid在原始数据矩阵中的行号。使用者拿到medoidIdx再用X(medoidIdx(k), :)就能提取出第k类的代表样本。为什么要用结构体而不是返回多个零散变量因为这个项目的使用场景通常是一套完整分析流程后面要接可视化、业务解读、报告输出把结果打包成一个结构体后续代码会清爽很多。另外我在函数里支持了Replicates参数也就是多个随机起点跑多次取最优这个对k-medoids尤其重要因为PAM本质是一个局部搜索算法初始点选得不好很容易陷进局部最优。2.2 距离矩阵计算细节项目里我同时支持了欧氏距离和曼哈顿距离默认欧氏。距离计算是整个算法最频繁的部分如果每次迭代现算现用数据量大一点就慢得没法看。所以我的做法很简单在进入迭代循环之前用一次矩阵运算把样本两两之间的距离矩阵D一次性算出来之后所有分配和更新操作都直接查这张表。欧氏距离矩阵用了一个常用技巧D(i,j)^2 ||x_i||^2 ||x_j||^2 - 2 * x_i * x_j用矩阵一次算出来比双重循环快非常多。计算过程中要注意浮点误差可能导致括号里出现一个极小的负数开根号之前先max(x,0)保护一下。曼哈顿距离我直接用了pdist2(X, X, cityblock)MATLAB自带实现性能稳定不用自己写循环。距离矩阵还有一个细节对角线元素要强制置0。虽然理论上样本到自己的距离就是0但浮点运算下可能出现极小的非零值影响后续min操作的判断。我在代码里用D(1:n1:end) 0一次性把对角线清零。2.3 分配与更新步骤的核心逻辑分配阶段代码里我用min(D(:, medoidIdx), [], 2)得到每个样本最近的medoid编号这里D(:, medoidIdx)是n行K列的矩阵第i行第k列就是样本i到第k个medoid的距离min函数返回的第一参数就是列下标也就是类别标签。更新阶段是PAM最核心的部分。对第k个簇我取出簇内所有样本的下标然后用D(clusterPoints, clusterPoints)截出簇内样本两两之间的距离子矩阵对每一行求和和最小的那个位置就对应新medoid。这里求和的意义是“如果选这个点当medoid整个簇的总距离是多少”所以取最小值就是找最优代表点。我看到网上很多实现会在更新时只比较替换前后的目标函数变化量这当然也行但对新手不友好。我采用的是直接遍历比总距离理解起来更直观代码也更短。代价是每个簇内部要做一次全对全距离求和复杂度大约是O(n²/K)样本量过万后会明显变慢这个问题我在后面“常见问题”章节专门讲优化思路。空簇处理也在这部分做了。如果某个簇在分配阶段一个点都没分到find(idx k)会返回空数组。这种情况不做防护的话下一步取clusterPoints(1)就直接越界报错。我的处理是从未被选作medoid的样本里随机挑一个补位并打印警告。虽然不优雅但至少不会让程序崩掉大多数情况下后续迭代能自动恢复。3. 完整MATLAB源代码与运行说明3.1 完整源代码下面是核心函数的完整代码。注释我写得比较细尤其是索引操作的部分几乎每行都有说明。function result kmedoidsClustering(X, K, varargin) % KMEDOIDSCLUSTERING 基于PAM思想的k-medoids聚类算法 % 输入: % X - n行d列数据矩阵每行是一个样本点 % K - 聚类簇数必须为正整数且小于等于样本数 % 可选参数键值对: % Distance - 距离度量可选 euclidean(默认) 或 manhattan % MaxIter - 最大迭代次数默认100 % Replicates - 随机初始化次数每次独立运行取最优结果默认1 % Verbose - 是否打印每次迭代信息默认false % 输出: % result - 结构体包含: % .idx - n行1列类别标签取值1~K % .medoidIdx - K行1列第k个簇的medoid在X中的行号 % .cost - 最终目标函数值总绝对距离 % .iterUsed - 实际迭代次数 % .history - 每次迭代的目标函数值用于观察收敛过程 % ---------- 检查输入 ---------- [n, d] size(X); % n是样本数d是特征维度 if K 1 || K n error(K必须满足 1 K 样本数); end if d 2 error(数据至少需要两列特征才能进行聚类); end % ---------- 解析可选参数 ---------- p inputParser; addRequired(p, X); addRequired(p, K); addParameter(p, Distance, euclidean, ... (x) any(strcmpi(x, {euclidean, manhattan}))); addParameter(p, MaxIter, 100, (x) isnumeric(x) isscalar(x) x 0); addParameter(p, Replicates, 1, (x) isnumeric(x) isscalar(x) x 1); addParameter(p, Verbose, false, islogical); parse(p, X, K, varargin{:}); opts p.Results; % ---------- 多次随机初始化保留最优结果 ---------- bestCost inf; for rep 1:opts.Replicates [idx, medoidIdx, cost, iterUsed, history] ... runSinglePam(X, K, opts.Distance, opts.MaxIter, opts.Verbose); if cost bestCost bestCost cost; result.idx idx; result.medoidIdx medoidIdx; result.cost cost; result.iterUsed iterUsed; result.history history(1:iterUsed); % 只保留实际迭代次数的记录 end end end % function [idx, medoidIdx, cost, iterUsed, history] ... runSinglePam(X, K, distType, maxIter, verbose) % 单次PAM运行不对外公开 [n, ~] size(X); % 随机初始化K个medoidrandperm无放回抽取保证不重复 medoidIdx randperm(n, K); medoidIdx medoidIdx(:); % 统一转成列向量避免版本差异 % 一次性算好距离矩阵后续所有查找都走这张表 D computeDistMatrix(X, distType); % 预分配内存 history zeros(maxIter, 1); for iter 1:maxIter % 1. 分配步骤每个样本归到距离最近的medoid所在簇 % D(:, medoidIdx)是n行K列第i行第k列样本i到第k个medoid距离 % min返回的col是每行最小距离对应的列下标正好就是簇编号 [~, assign] min(D(:, medoidIdx), [], 2); idx assign; % 2. 更新步骤在每个簇内部选总距离最小的点作为新medoid newMedoidIdx zeros(K, 1); for k 1:K clusterPoints find(idx k); if isempty(clusterPoints) % 空簇保护从未选中的点里随机补一个 used newMedoidIdx(newMedoidIdx ~ 0); unused setdiff(1:n, used); if isempty(unused) % 极端情况所有点都已被选为medoid保留原值 newMedoidIdx(k) medoidIdx(k); else newMedoidIdx(k) unused(randi(length(unused))); end continue; end % 取簇内所有点两两之间的距离子矩阵按行求和得到每个点的总距离 % 总距离最小的那个点就是当前簇的最优medoid subD D(clusterPoints, clusterPoints); totalDist sum(subD, 2); [~, bestPos] min(totalDist); newMedoidIdx(k) clusterPoints(bestPos); end % 3. 计算本次迭代后的目标函数值总绝对距离 cost 0; for k 1:K clusterPoints find(idx k); if ~isempty(clusterPoints) cost cost sum(D(clusterPoints, newMedoidIdx(k))); end end history(iter) cost; if verbose fprintf(迭代 %d目标函数 %.4f\n, iter, cost); end % 4. 收敛判断比较medoid集合是否一致 % 用排序后的集合比较避免簇编号互换导致的误判 if isequal(sort(newMedoidIdx), sort(medoidIdx)) medoidIdx newMedoidIdx; break; end medoidIdx newMedoidIdx; end iterUsed iter; % 5. 收敛后重新按最终medoid计算一次类别标签确保对应关系准确 [~, idx] min(D(:, medoidIdx), [], 2); end % function D computeDistMatrix(X, distType) % 计算样本间两两距离矩阵D(i,j)是第i个样本到第j个样本的距离 % 核心思路一次矩阵运算代替双重循环 [n, ~] size(X); if strcmpi(distType, euclidean) % 欧氏距离的矩阵化展开 % ||a-b||^2 ||a||^2 ||b||^2 - 2*a*b % 其中||a||^2是逐行元素平方再求和 XX sum(X .* X, 2); % n×1列向量每个样本的模长平方 D sqrt(max(bsxfun(plus, XX, XX) - 2 * (X * X), 0)); % max(...,0)防止浮点误差导致负数开根号 elseif strcmpi(distType, manhattan) % 曼哈顿距离直接用MATLAB内置函数性能稳定 D pdist2(X, X, cityblock); end % 对角线强制置0避免浮点误差影响min操作 D(1:n1:end) 0; end3.2 调用方式与参数说明这个函数用起来很简单。主脚本里先生成示例数据再调用最后画图。注意一点如果不同特征的量纲差异很大聚类结果会被数值大的特征主导。比如做客户分群时消费金额可能几千消费频次只有个位数这时候不标准化频次特征基本等于被忽略。正确做法是先做z-score标准化再聚类。% 示例主脚本 % 生成三簇高斯分布数据并混入少量离群点 rng(123); % 固定随机种子保证结果可复现 X [randn(80, 2) [3, 3]; % 第一簇中心在(3,3) randn(80, 2) [-3, -2]; % 第二簇中心在(-3,-2) randn(80, 2) [8, -6]; % 第三簇中心在(8,-6) randn(5, 2) * 5]; % 离群点故意放大方差 % 特征标准化消除量纲影响 X zscore(X); % 跑k-medoidsK3欧氏距离最多跑10个随机初始化保留最优 result kmedoidsClustering(X, 3, ... Distance, euclidean, MaxIter, 100, Replicates, 10, ... Verbose, true); % 查看结果 disp(result.idx); % 每个样本的类别标签 disp(result.medoidIdx); % 每个簇的medoid在X中的行号 fprintf(总距离: %.4f迭代次数: %d\n, result.cost, result.iterUsed); % 可视化用不同颜色区分簇黑色三角标出medoid figure; gscatter(X(:, 1), X(:, 2), result.idx, rgbcmyk, ., 14); hold on; plot(X(result.medoidIdx, 1), X(result.medoidIdx, 2), ... k^, MarkerSize, 16, LineWidth, 2.5); legend off; axis equal; title(k-medoids聚类结果黑色三角为medoid);如果不想传入那么多参数可以直接result kmedoidsClustering(X, 3)其余全部走默认值。实测下来默认100次最大迭代对大多数中小规模数据都足够收敛Replicates建议至少设5可以有效避开局部最优。3.3 结果可视化与解读运行上面的代码会得到一张散点图不同颜色代表不同簇黑色三角代表每个簇的medoid。你会注意到一个很有意思的特点medoid落在真实数据点上而且通常处在簇的“视觉中心”附近但又不会被离群点拉着跑。对比一下同样是K3的k-means结果两者在分布规则的数据上差别不大但一旦数据里混入离群点k-means的一个簇心会被离群点拽走而k-medoids的medoid依然稳如泰山。再解释一下“通过标识medoid确定类别”的含义。当你跑完聚类得到result之后如果想判断一个新样本属于哪个簇不需要重新跑算法只需要计算这个新样本到K个medoid的距离取最近的medoid编号作为类别即可。换句话说medoid就是类别的“标识物”类别标签不是独立存在的而是由medoid的行号所定义的。在代码里result.medoidIdx就承担了这个标识功能。4. 常见问题与排查技巧实录4.1 中文注释乱码问题这个项目是全中文注释那大概率会遇到一个很烦的问题在别人电脑上或者换了MATLAB版本之后打开代码一看中文全变成了乱码。这个问题不是代码错误而是历史版本的编码差异。旧版MATLAB比如R2020b及更早默认用系统本地编码保存.m文件在中文Windows上通常是GBK而新版MATLAB从某个版本开始默认用UTF-8。两个编码打开同一份文件就容易乱。解决办法有两个路径。如果代码需要在新版本MATLAB里长期维护用编辑器打开乱码文件后选择“另存为”在编码下拉框选UTF-8保存一遍以后再打开就正常了。如果你手头的脚本要发给别人跑建议在文件头部加一行注释说明编码或者干脆把中文注释集中在代码顶部避免在关键算法行里夹中文变量名之类的写法。4.2 聚类结果不稳定k-medoids非常依赖初始medoid的选取。不同随机起点可能得到不同的局部最优解这是PAM这类算法的本性不是代码bug。我在函数里用Replicates参数来解决这个问题本质上就是“多跑几次取最优”和MATLAB自带的kmeans函数里Replicates参数思路一致。如果你发现多次运行结果差异还是很大可以从两个方向排查。第一确认K值是否合理K设置不合理时本来就很难稳定。第二看数据里是否存在大量离群点离群点太多会导致某些簇只靠一两个离群样本撑着这种情况下建议先做离群点剔除或改用带噪声过滤的聚类方法。另外固定随机种子rng(123)可以让结果可复现这在写报告和调试时是必须养成的习惯。4.3 K值怎么选K值选择是个老生常谈又绕不开的问题。k-medoids没有返回“最优K”的功能需要自己选。我在实际项目里最常用的两个方法手肘法和轮廓系数。手肘法就是跑不同K值记录每个K对应的目标函数值也就是result.cost然后画一条折线。随着K增大cost一定会下降但下降速度会在某个点突然变缓这个拐点也就是折线上的“手肘”就是比较合理的K。轮廓系数法用MATLAB内置的silhouette(X, idx)可以算但注意它接受的idx是类别标签也就是result.idx取值范围要在1到K之间这个我们正好满足。% 快速跑K2到K8观察手肘 Ks 2:8; costs zeros(size(Ks)); for i 1:length(Ks) tmp kmedoidsClustering(X, Ks(i), Replicates, 5); costs(i) tmp.cost; end plot(Ks, costs, o-); xlabel(K); ylabel(总距离);这里的思路是K每增加1总距离一定会降一点但如果新增一个簇只是把原有簇硬拆开距离下降幅度会很小表现在图上就是曲线变平。选变平之前那个拐点位置既保证簇内紧凑又不至于过度细分。4.4 大数据集太慢怎么办PAM算法本身计算量不小尤其是更新阶段每个簇都要算一次簇内全对全距离。我实测过样本量在几千左右完全没压力但上了几万就会明显感觉到每次迭代变慢十万级别基本不建议跑经典PAM了。如果数据量确实大优先考虑三条优化路径。第一条换用距离矩阵的近似比如只保留每个样本最近的几个邻居距离其他用一个大数填充。第二条改用CLARA算法思路先对数据做随机采样在样本子集上跑k-medoids再用样本子的medoid去给全体数据分配类别。第三条如果项目允许直接换用k-means或者DBSCAN毕竟k-medoids的鲁棒性优势主要体现在中小规模数据和高维离群点上数据量上去了计算成本会快速抵消它的优点。另外一个容易忽略的性能瓶颈是computeDistMatrix里bsxfun那行。虽然做了矩阵化但在超大矩阵上它仍然会一次性生成n行n列的距离矩阵内存占用是O(n²)。一万个样本算下来就是上亿个浮点数几十个GB内存直接卡死。所以上面说的采样思路不仅节约时间更是为了节约内存。5. 从代码到业务的落地经验5.1 实际业务中的典型使用场景项目跑通了代码没问题接下来更重要的问题是这套东西到底能用在哪。我接触过的真实案例里最有代表性的是客户分群和供应链异常检测。客户分群场景下特征一般有消费频次、客单价、最近一次消费时间、品类偏好等等。聚类完成后k-means给出的中心点是一个“平均人”而k-medoids给出的medoid是一个真实的客户样本。你可以直接调出这条客户的档案跟业务同事说“你看这一群的代表客户就是这位消费频次高但客单价低偏好日用百货。”这种沟通方式业务同事能听懂也愿意配合你后续制定运营策略。这就是medoid标识类别的可解释性红利。另一个场景是设备或传感器异常分析。工业数据里经常有大量噪声点用k-means做聚类时异常点会污染簇心可能把一个本来正常的簇拉到偏向异常的位置。k-medoids对离群点天然鲁棒因为medoid是从真实点位中选出来的异常点很难被选为代表。这样聚类出来之后你直接看每个簇的medoid对应的设备参数就能理解这个簇的“正常运行模式”长什么样离medoid特别远的点就是潜在异常候选。5.2 我的一些使用心得与后续扩展方向代码写完只是开始真正用的时候有几个细节值得注意。一是特征标准化一定要做而且是每次聚类前都要做哪怕上一次已经做过。因为只要数据源稍微更新几条整体均值和方差就会变标准化结果也随之变化。二是多跑几次Replicates不丢人反而能帮你确认结果是否稳定如果10次跑出来medoid完全一致说明这个聚类结构比较可靠如果每次都不一样那你要先检查K值是否合理。三是在写报告时把result.medoidIdx和result.cost都记录下来一个代表结果结构一个代表结果质量下次对比不同方案时非常有用。后续如果要扩展比较推荐的方向是轮廓系数评估和可视化升级。轮廓系数可以定量评估每个样本分得好不好而不只是看总距离。可视化方面可以结合gscatter把原始数据、medoid、以及每个样本到medoid的距离用线条画出来一张图就能说明聚类结构的紧密度。甚至还能把medoid提取出来做下一轮分析的输入比如对medoid再做一次层次聚类得到“类别的类别”这在客户画像分层里很有用。最后再提一个我踩过好几次的坑千万记得在代码开头rng固定随机种子否则你分析到凌晨终于得到一版满意的结果第二天一打开MATLAB重新运行聚出来的结果完全变了那种绝望我不想再经历第二次。