尧图网页设计 
简介面向Matlab用户的kMedoids聚类示例代码包聚焦于以真实样本作为聚类中心的鲁棒聚类方法。与K-means依赖质心不同kMedoids选择数据集中实际观测点作为medoids抗异常值干扰、结果更易解释对存在噪声或类别变量的数据尤为友好。包内仅含1个m脚本整体大小2KB可直接在Matlab中打开查看源码快速理解核心实现。脚本覆盖随机多维数据生成、调用内置kmedoids函数时对聚类数、距离度量、初始化方法、最大迭代次数与容差等关键参数的指定以及聚类结果查看等基本流程并可通过轮廓系数评估聚类质量。资源已有237人学习虽体量精简但对理解kMedoids原理、参数设置及扩展应用如市场细分、图像分割具有直接参考价值适合刚接触聚类分析或希望对比K-means差异的Matlab用户作为简洁起点。1. kMedoids 是什么聚类算法里那个不背均值包袱的鲁棒选项打开 Kmedios.rar解压出来通常是一堆 .m 文件和一份说明文档——这是 kMedoids 聚类算法在 MATLAB 里的实现也是很多模式识别课程作业里绕不开的一个坎。kMedoids 和 K-Means 长得很像但它在迭代中用真实样本点而不是均值当簇中心所以对离群点更鲁棒。这篇笔记直接讲三件事kMedoids 和 K-Means 到底差在哪、在 MATLAB 里怎么从零手写并跑通它、以及拿到 Kmedios.rar 这类代码后怎样验证它真的能用。适合正在做聚类相关课程设计或论文复现的人也适合初次接触这类资源包、想搞明白里面每行代码在干什么的 MATLAB 用户。2. 从 K-Means 到 kMedoids那个把均值换成样本点的关键差异2.1 K-Means 的均值中心为什么怕离群点一个具体例子先看一个最直白的例子。假设你拿到一维数据 [1, 2, 3, 8, 100]想分成两类。K-Means 的更新规则是每一轮把簇内所有点的均值当作新的簇中心。如果某个簇同时包含 8 和 100它的均值就是 54这显然不是任何真实样本而且和 8 的距离高达 46和 100 也有 46——这个虚拟中心根本没落在数据密集区。如果你把 100 单独视为离群点剩下的 [1, 2, 3, 8] 再分聚类结果就合理得多。问题出在均值的天性上它对极端值没有任何免疫力。一个量级特别大的观测点哪怕只有一个也会把簇中心拖向自己进而导致周围原本正常的点被错误归组。更麻烦的是K-Means 的目标函数是最小化点到中心的平方欧氏距离这是凸优化没错但离群点拉偏中心的能力同样得到强化——平方项给了大距离点更高的权重。而 kMedoids 把中心的候选限定在真实样本里。还是刚才这组数据分成两个簇时簇 {1,2,3} 的中心选 2簇 {8,100} 的中心选 8 或 100。无论选哪个中心都是数据集里实打实的观测值不会出现 54 这种人畜无害但毫无意义的数字。就这一个差别离群点对中心的影响就小得多——极端值最多让自己成为中心而不会把中心拖到一个不存在的坐标上。这里要说明一点kMedoids 并不保证聚类结果更好它保证的是中心可解释、结果对离群点更稳健。如果你的数据本身没有离群点且簇形状近似球形K-Means 的速度优势反而值得优先考虑。选型这件事继续往下说。2.2 kMedoids 的核心思想用真实样本做中心kMedoids 的算法骨架和 K-Means 几乎一样初始化 k 个中心点、把每个样本分到最近的中心、用簇内样本更新中心、重复直到收敛。唯一的区别在第 3 步——K-Means 直接计算簇内均值而 kMedoids 在簇内枚举每个候选样本挑到其他所有簇成员的总距离最小的那个点作为新中心。这个点就是 medoid也叫中心点或类中心。正式一点说kMedoids 要最小化的代价函数是每个样本都归到距离它最近的 medoid目标就是让全体样本到各自 medoid 的距离之和尽量小。算法通过迭代调整 medoid 集合来降低这个总代价。和 K-Means 一样这个代价是单调不增的所以算法保证收敛但不保证全局最优。业界最经典的实现是 Kaufman 和 Rousseeuw 提出的 PAMPartitioning Around Medoids算法它的流程分两步。BUILD 阶段先挑 k 个最分散的点作为初始 medoids避免随机初始化的偶然性SWAP 阶段每一轮遍历所有 (medoid, 非 medoid) 对计算交换后总代价的增量如果找到负增量就执行交换。这一轮全量扫描的计算量非常大SWAP 阶段的每轮迭代接近 O(k·(n−k)·n) 的量级所以 PAM 对大数据不友好。很多 MATLAB 代码习惯上把 PAM 的 SWAP 简化成只在簇内交换只尝试当前簇内的点而不试探跨簇的候选中心。这样做收敛更快但理论上可能错过更好的全局中心。我的建议是课程作业、小规模数据集、论文复现用简化版完全够生产环境对聚类质量要求高优先考虑 MATLAB 内置的完整 PAM 实现。另外BUILD 阶段虽然听起来很好手写实现中很多人直接跳过它用随机初始化。后果是初始点质量差时SWAP 阶段要多花很多轮次才能收敛。一个折中做法是执行 5 次随机初始化每次跑 20 轮保留代价最小的那次作为最终结果——用时间换稳定性对小数据集非常划算。2.3 选型决策什么时候用 kMedoids 而不是 K-Means到底什么场景值得牺牲速度去换 kMedoids 的鲁棒性我按经验排了一下。第一数据里明确有离群点或噪声。传感器数据、金融异常值、用户行为日志中那些抖一下的极端值会直接把 K-Means 的均值中心拉走。kMedoids 受的影响小因为离群点最多把自己标成一个 medoid不会成为幽灵中心。第二中心需要可解释。在推荐系统或客户分群场景里如果聚类中心对应着一位真实的用户画像、一条真实的交易记录业务方更容易理解这个簇代表谁。K-Means 的均值向量可能是现实中不存在的人。第三距离度量不是欧氏距离。K-Means 的均值在曼哈顿距离下没有闭合解kMedoids 只需要距离矩阵任何两个点之间能定义的距离都能用。kMedoids 因此常常和余弦距离、动态时间弯曲DTW配对用在文本和时序聚类上。反过来讲三种情况尽量别用 kMedoids数据量大比如超过 10 万条SWAP 阶段的复杂度会让你等到怀疑人生簇形状复杂kMedoids 和 K-Means 一样基于中心辐射假设对环形、嵌套簇都无能为力该换 DBSCAN 就换要聚类的维度特别高高维空间里最近邻的区分度会急剧下降这在避坑章单独说。最后补充一个经验值我在图像分割场景里发现kMedoids 作用于颜色直方图和纹理特征时经常比 K-Means 表现更好代价是每轮迭代多花几十秒。但这个时间成本换来的是分割区域更连续、边界更干净不少 MATLAB 图像处理大作业会用 kMedoids 替代 K-Means 来提升效果。3. MATLAB 手写 kMedoids从零实现一个可运行的聚类器3.1 最小可运行版本核心代码与逐段解释先上一个完整可跑的最小实现。我管它叫kmedoids_diy.m输入输出尽量向 MATLAB 内置kmedoids函数对齐方便你对比测试。function [idx, midx, total_cost] kmedoids_diy(X, k, max_iter) % kmedoids_diy 简易kMedoids聚类PAM简化版簇内交换策略 % 输入 % X : n×d 数据矩阵每行一个样本 % k : 聚类簇数必须大于0且小于n % max_iter : 最大迭代次数默认100 % 输出 % idx : n×1 簇标签值范围1~k % midx : k×1 每个簇的medoid在X中的行索引 % total_cost: 标量所有样本到其medoid的距离之和 if nargin 3 max_iter 100; end [n, ~] size(X); if k n error(k必须小于样本数n); end % 预计算n×n距离矩阵欧氏距离 D sqrt((sum(X.^2,2) sum(X.^2,2) - 2*(X*X))); D(1:n1:end) 0; % 对角元素清零避免自距干扰 % Step 1: 随机初始化k个不同的medoid索引 rng(42); % 固定种子保证结果可复现 midx randperm(n, k); % 迭代变量 old_cost inf; for iter 1:max_iter % Step 2分配: 每个样本归属最近的medoid [~, idx] min(D(:, midx), [], 2); % Step 3更新: 对每个簇尝试用簇内其他点替换当前medoid changed false; for j 1:k members find(idx j); if isempty(members) continue; % 空簇本迭代不处理留到避坑章讨论 end cur_m midx(j); % 当前簇的代价所有成员到当前medoid的距离和 cur_cost sum(D(members, cur_m)); % 遍历每个簇内样本找代价最小的候选中心 best_cost cur_cost; best_m cur_m; for t 1:length(members) cand members(t); if cand cur_m continue; % 自己换自己没有意义 end c_cost sum(D(members, cand)); % 候选点的总代价 if c_cost best_cost best_cost c_cost; best_m cand; end end if best_m ~ cur_m midx(j) best_m; % 执行替换 changed true; end end % 计算本轮总代价用于收敛判断 total_cost sum(min(D(:, midx), [], 2)); if ~changed || abs(old_cost - total_cost) 1e-6 break; end old_cost total_cost; end % 收敛后再做一次最终分配 [~, idx] min(D(:, midx), [], 2); total_cost sum(min(D(:, midx), [], 2)); end这段代码的核心逻辑和 K-Means 几乎同一个骨架差异只在更新中心这一步。先说距离矩阵的向量化写法sqrt((sum(X.^2,2) sum(X.^2,2) - 2*(X*X)))原理是欧氏距离展开式||a-b||² ||a||² ||b||² - 2ab一次矩阵乘法代替双重循环n 不超过一万时性能都还行。D(1:n1:end)0是把每个点和自己的距离清零避免数值误差带来的非零自距。min(D(:, midx), [], 2)这一步是整个算法的核心对第 i 个样本D(i, midx)给出它到所有 k 个 medoid 的距离取最小者所在列号就是它的簇标签。这里的向量化写法比for i1:n的循环快一个数量级建议不要改成循环。再说参数。max_iter100是一个偏保守的上限——PAM 类算法在中小数据集上通常 20 轮以内就能收敛100 轮基本不会触顶。真正的收缩点是随机种子rng(42)你把这个数字去掉每次跑出来的 medoid 初始位置都不同最终结果大概率也不同。这不是 bug而是 kMedoids 的非凸性质决定的避坑章会展开讲。3.2 用合成数据验证聚类效果从二维图上读结果光有函数不行得跑起来看效果。我一般用带离群点的二维高斯混合数据来验收算法——二维便于可视化离群点能检验算法鲁棒性。% 生成合成数据两个高斯簇 一个离群点 rng(10); X1 randn(40,2) * 0.6 repmat([2, 2], 40, 1); X2 randn(40,2) * 0.6 repmat([-2, -2], 40, 1); X [X1; X2; 10, 10]; % 最后一行是故意放的离群点 % 调用自定义kMedoids [idx, midx, cost] kmedoids_diy(X, 2, 100); % 可视化 figure; gscatter(X(:,1), X(:,2), idx, rb, xo); hold on; plot(X(midx,1), X(midx,2), kp, MarkerSize, 15, LineWidth, 2); title(kMedoids 聚类结果黑方块为medoid); legend(簇1, 簇2, Medoid, Location, best);跑完你会看到两个事实。第一离群点 (10,10) 会被独自分到另一个簇而两个高斯簇的分界线基本落在 (−2,−2) 到 (2,2) 之间的对角线上。第二medoid 点黑色方块一定是数据里真实存在的点而不是像 K-Means 那样输出一个坐标均值。这两个特征分别对应前面讲的鲁棒性和可解释性。值得试的对比实验把上面代码里的kmedoids_diy换成 MATLAB 内置的kmedoids函数两个结果在无离群点的数据上几乎一致。但在有离群点时自定义实现和内置函数通常都会正确隔离离群点而 K-Meanskmeans(X,2)则很可能把一个高斯簇劈成两半离群点拽着中心跑偏。这一步对比能帮你直观建立什么时候该信任 kMedoids的判断力。3.3 参数选择的几个要点标准化不是可选项。如果特征量纲差异大比如一列是 0 到 1另一列是 0 到 10000距离矩阵会被量纲大的特征主导。常规做法是X zscore(X)让所有特征方差为 1。这个操作对任何基于距离的聚类都适用不止 kMedoids。距离度量的选择要匹配数据语义。欧氏距离适合连续特征且各方向等权重的场景曼哈顿距离对坐标这类数据更自然文本 TF-IDF 向量适合余弦距离。MATLAB 内置kmedoids的Distance参数支持这些选项自定义实现里只需要替换距离矩阵的计算方式即可。空簇怎么处理。我的简化实现在更新步骤遇到空簇是直接跳过的这在 K-Means 里也有对应的经典问题。更稳妥的兜底策略是分配步骤结束后检查是否出现空簇如果出现就随机选一个非 medoid 样本补充进去。但这个策略要小心它可能打破算法已收敛的状态。4. 用好内置 kmedoids 与 Kmedios.rar下载代码之前先看这几点4.1 内置 kmedoids 函数与自定义实现的对比MATLAB 从 R2018a 开始在 Statistics and Machine Learning Toolbox 里提供了官方的kmedoids函数。接口非常干净% 基本用法X是数据矩阵k是聚类数 [idx, C] kmedoids(X, k); % 带更多输出sumd是每个簇的总距离D是样本到各medoid的距离 [idx, C, sumd, D, midx, info] kmedoids(X, k, Distance, sqeuclidean);其中C返回的是 k×d 的 medoid 坐标矩阵midx返回的是 medoid 在 X 中的行索引。info结构体里包含迭代次数、终止原因等诊断信息调试时多看一眼它能省不少时间。内置函数和自定义实现相比主要强在三点。第一内置的 SWAP 阶段是完整的 PAM 扫描搜索范围不只是当前簇内还会考虑所有非 medoid 候选理论上更容易逼近全局最优。第二Start参数允许你传入初始 medoid 矩阵k×d这对复用上一次结果继续迭代很有价值。第三内置函数针对大数据集提供了Algorithm, clara选项它是对行做抽样、在样本上跑 PAM能扛住几十万行的数据这是手写实现很难做到的。不过内置函数也有个短板它对输入矩阵的数据类型和维度有校验如果你传的是距离矩阵而不是原始数据需要显式换方案——内置kmedoids不支持直接传入成对距离矩阵这一点有时候让人挠头。4.2 Kmedios.rar 的典型结构与使用路径这类压缩包在资源分享站、GitHub 和网盘转存链接里非常常见。解压出来通常是三类文件一个主函数比如kmedoids.m或mykmedoids.m、一份测试脚本或 demo、以及一个说明文档。也有少数包里带示例数据或参考论文 PDF。拿到包的第一件事不是跑是看函数签名% 打开主函数文件后先看第一行function定义 % 常见签名有这几种形态 % [idx, C] kmedoids(X, k) % [idx, C, sumd] kmedoids(X, k, dist_type) % [label, center] KMEDOIDS(X, k, maxiter)签名决定了你调用它时传什么参数、返回什么数据结构。很多翻车现场都是因为没看签名拿着内置函数的调用方式去调一个自定义函数结果参数数量不匹配或输出错位。接下来是把它跑通的三个步骤。第一步确认依赖函数都在 MATLAB 路径上。压缩包里的 .m 文件可能散在不同子目录如果主函数调用了某个子函数而它不在路径上运行时会报 Undefined function。第二步用 demo 脚本里的数据跑一遍如果包里没带 demo你就用第 3 章那段合成数据生成代码自己构造数据。第三步把结果和内置kmedoids对比——同一个数据集、同一个 k 值两组输出的簇标签可能不同因为初始化和标签编号不同但total_cost应该非常接近。如果差得离谱说明下载的代码存在逻辑问题不值得继续投入。4.3 验证一个下载的 kMedoids 代码的 4 个步骤我在拿到任何一个网上下载的聚类代码时会做四件事来确认它没有暗伤。第一构建一个已知答案的数据集。生成 3 个远离的高斯簇标准差距大一点这样任何正常聚类算法都百分百能分开。如果下载的代码在这么简单的数据上都分不对那它显然有 bug。第二测试极端参数。k1时所有样本必须归同一簇kn时每个样本自成簇代价为 0样本点到自身的距离为 0。这两个边界测试能暴露索引越界和距离计算异常的问题。% 边界测试示例 [idx1, ~, cost1] kmedoids_diy(X, 1, 50); % 应全部归为1 [idx2, ~, cost2] kmedoids_diy(X, size(X,1), 50); % 每个样本一个簇 % 检查 cost2 是否约等于0数值误差范围内第三比较收敛代价。对同一个数据集跑kmeans和kmedoids观察 kMedoids 的代价是否低于 K-Means。这不是必然的因为目标函数不同但通常 kMedoids 的代价不会高于同规模聚类的 K-Means 太多。如果发现 kMedoids 代价高出几个量级多半是距离计算或分配逻辑出了问题。第四检查是否有编码问题。很多从国内资源站下载的 .m 文件是 GBK 编码在 UTF-8 环境用 MATLAB 直接打开中文注释会显示成乱码注释中包含引号时甚至报语法错误。解决办法是用文本编辑器把编码转成 UTF-8或者删掉注释只保留代码逻辑再保存。5. 避坑kMedoids 在 MATLAB 里最容易踩的 5 个坑5.1 距离矩阵算错导致聚类结果完全错乱现象聚类结果完全不符合直觉两个明显分开的簇被割裂开或者所有点都被归到同一个簇里。原因大多数手写实现用双重循环逐对计算距离稍不注意索引错位或者在对角元素清零时用了错误的线性索引就会让D(i,i)变成一个不该出现的干扰项。更隐蔽的错误是向量化公式里X*X的形状不匹配——X 是 n×d 时没问题但如果 X 里含有 NaN 或 Inf整个距离矩阵都会出问题。解决先把距离矩阵可视化。imagesc(D)能让你肉眼检查矩阵应该是对称的、对角线为 0、有离群点的数据集总会在某行某列出现明显的高亮。再用一个简单案例做验证[0,0; 1,0]两点之间的距离手算出应该是 1代码跑出来如果不是 1那就是距离计算本身的问题。5.2 随机种子影响结果同样代码两次跑出不同聚类现象同一份代码、同一个数据集第一次跑输出簇标签和代价与第二次完全不同甚至聚类结构都不一样。原因kMedoids 的初始中心是从样本中随机挑选的这是非凸优化问题——不同的初始点会落到不同的局部最优。特别是当你的数据有两个簇距离很近时初始化的微小差异可能把边界处的样本分到不同侧。解决固定随机种子rng(0)或rng(42)确保每次跑一致这是论文复现和课程作业提交的底线要求。如果你想评估这个数据集上 kMedoids 到底稳不稳可以跑 30 次不同的种子记录每次的代价看代价的方差——方差大说明聚类结构本身不够稳定下游分析要小心。内置函数也支持这个思路循环里每次调用前重新rng或删掉种子收集总代价。5.3 空簇问题现象某个 k 值下迭代结束后有一个簇的成员数等于 0medoid 成为一个无兵的将军。可视化时图例里出现空集。原因当初始 medoid 选得太靠近、而数据又有离群点分布时离群点可能单独占一个簇挤压其他簇的空间或者是某个 medoid 在更新过程中被换走之后没有新成员分给它。解决我一般用两招。第一招是初始化约束随机选初始 medoid 时要求两两之间的最小距离大于某个阈值比如所有样本距离的 10% 分位数从根源上降低初始中心扎堆的概率。第二招是运行后检查tabulate(idx)扫一眼成员数分布如果有空簇减少 k 值、更换种子或改用内置kmedoids的Start参数传入更好的初始中心。另外注意内置 kmedoids 遇到空簇时的行为可能和自定义实现不同它可能自动丢弃该簇中心并重建这会导致最终输出只有 k−1 个簇。如果你对簇数有刚性要求比如必须分 5 类空簇问题要在前置环节就规避掉。5.4 Distance 参数选错的连锁反应现象在某个数据集上用欧氏距离聚类效果不错换了个数据集后结果完全不可理喻但数据标准化似乎也没解决问题。原因kMedoids 的距离并不非得是欧氏距离。如果你的特征是词频、直方图或方向向量欧氏距离会把大量信息稀释掉。比如两篇文档各 1000 个词它们共享的词汇只有 20 个欧氏距离会很大但语义上它们可能非常相关。这时候用余弦距离才是正常的。解决根据数据语义选距离。MATLAB 内置kmedoids支持cityblock、cosine、correlation等选项。判断依据很简单算出来距离到底代表什么——对应到实际业务里这个距离大与小是否和人对事物的相似性判断一致。如果一致说明距离选对了如果不一致趁早换。5.5 高维数据的维度灾难现象数据维度从几十升到几千聚类结果开始变得随机代价曲线几乎分不出 k 值的拐点轮廓系数也接近 0。原因高维空间里所有点之间的距离都趋向于同量级欧氏距离的区分度急剧下降。这不是 kMedoids 的错是所有基于距离的聚类算法的通病。解决先降维再聚类。常见做法是 PCA 保留 90% 方差的前若干个主成分或者 t-SNE、UMAP 降到 2~3 维做可视化辅证。但注意——降维会丢失信息聚类结果要在原始维度上做代价评估不要在降维后的坐标上评估。另外选correlation或cosine这类对幅度不敏感的距离度量有时比降维更有效因为它把方向作为主信号天然忽略零均值的噪声维度。还有一个和维度灾难常常同时出现的坑是数据稀疏性。如果你做的是用户-物品偏好矩阵聚类每行只有少数非零元素直接在这个稀疏矩阵上算欧氏距离结果基本等于看谁和谁的零值坐标重合得多。这种场景正确的做法是先做 SVD 降维或使用correlation距离让聚类信号从共现模式里提取而不是从零值里提取。6. 选好 k 值与验证聚类质量收尾的两步实操6.1 肘部法画出代价曲线找拐点kMedoids 要求预先指定 k 值一个实用且直观的做法是肘部法对 k1 到 10 分别跑聚类记录每次的总代价然后画图。代价会随着 k 增加而下降但下降幅度在某个 k 值之后明显放缓那个拐点就是性价比最高的 k。% 肘部法示例 ks 1:10; costs zeros(size(ks)); for i 1:length(ks) [~, ~, costs(i)] kmedoids_diy(X, ks(i), 100); end plot(ks, costs, -o); xlabel(k); ylabel(total cost);注意一点肘部法的拐点在真实数据上不总是明显。如果曲线平滑下降没有拐点说明数据没有天然可分性这时别强行挑 k考虑先对数据做预处理或换聚类算法。6.2 轮廓系数一个本地就能算的质量指标轮廓系数不要求外部标签只计算每个样本和它所在簇的紧密度以及和其他簇的分离度。范围在 -1 到 1 之间接近 1 说明聚类结构清晰接近 0 说明样本在两个簇之间摇摆。轮廓系数的 MATLAB 实现可以直接用evalclusters这条路径但我更推荐手写一次因为能顺便加深对聚类的理解function s silhouette_diy(X, idx) % 简易轮廓系数欧氏距离 k max(idx); n size(X,1); s zeros(n,1); D pdist2(X, X); for i 1:n % 样本i到同簇其他样本的平均距离 a_i mean(D(i, idx idx(i) (1:n) ~ i)); if isempty(a_i), a_i 0; end % 样本i到最近异簇所有样本的平均距离 b_vals arrayfun((j) mean(D(i, idx j)), setdiff(1:k, idx(i))); b_i min(b_vals); s(i) (b_i - a_i) / max(a_i, b_i); end s mean(s); end轮廓系数大于 0.25 就说明聚类结构不是纯噪声了如果小于 0.1要么 k 没选对要么数据本身不适合做中心型聚类得回头检查特征工程。做轮廓系数评估时同样固定随机种子不同初始化的轮廓系数差异如果超过 0.05就多跑几次取平均。往深了说我养成的一个习惯是每次跑完 kMedoids除了看轮廓系数还会把紧密度最高的几个簇内样本单独拉出来看原始特征——聚类是工具不是答案真正交付给业务方的是这个簇代表什么含义。如果簇内样本在业务维度上完全说不过去哪怕轮廓系数是 0.8这个聚类结果也值不了几分钱。说到底算法就那几行循环真正值钱的是你对自己数据的理解。这次就聊到这儿希望帮到你。本文还有配套的精品资源点击获取 kMedoids聚类算法与MATLAB手写实现详解
简介面向Matlab用户的kMedoids聚类示例代码包聚焦于以真实样本作为聚类中心的鲁棒聚类方法。与K-means依赖质心不同kMedoids选择数据集中实际观测点作为medoids抗异常值干扰、结果更易解释对存在噪声或类别变量的数据尤为友好。包内仅含1个m脚本整体大小2KB可直接在Matlab中打开查看源码快速理解核心实现。脚本覆盖随机多维数据生成、调用内置kmedoids函数时对聚类数、距离度量、初始化方法、最大迭代次数与容差等关键参数的指定以及聚类结果查看等基本流程并可通过轮廓系数评估聚类质量。资源已有237人学习虽体量精简但对理解kMedoids原理、参数设置及扩展应用如市场细分、图像分割具有直接参考价值适合刚接触聚类分析或希望对比K-means差异的Matlab用户作为简洁起点。1. kMedoids 是什么聚类算法里那个不背均值包袱的鲁棒选项打开 Kmedios.rar解压出来通常是一堆 .m 文件和一份说明文档——这是 kMedoids 聚类算法在 MATLAB 里的实现也是很多模式识别课程作业里绕不开的一个坎。kMedoids 和 K-Means 长得很像但它在迭代中用真实样本点而不是均值当簇中心所以对离群点更鲁棒。这篇笔记直接讲三件事kMedoids 和 K-Means 到底差在哪、在 MATLAB 里怎么从零手写并跑通它、以及拿到 Kmedios.rar 这类代码后怎样验证它真的能用。适合正在做聚类相关课程设计或论文复现的人也适合初次接触这类资源包、想搞明白里面每行代码在干什么的 MATLAB 用户。2. 从 K-Means 到 kMedoids那个把均值换成样本点的关键差异2.1 K-Means 的均值中心为什么怕离群点一个具体例子先看一个最直白的例子。假设你拿到一维数据 [1, 2, 3, 8, 100]想分成两类。K-Means 的更新规则是每一轮把簇内所有点的均值当作新的簇中心。如果某个簇同时包含 8 和 100它的均值就是 54这显然不是任何真实样本而且和 8 的距离高达 46和 100 也有 46——这个虚拟中心根本没落在数据密集区。如果你把 100 单独视为离群点剩下的 [1, 2, 3, 8] 再分聚类结果就合理得多。问题出在均值的天性上它对极端值没有任何免疫力。一个量级特别大的观测点哪怕只有一个也会把簇中心拖向自己进而导致周围原本正常的点被错误归组。更麻烦的是K-Means 的目标函数是最小化点到中心的平方欧氏距离这是凸优化没错但离群点拉偏中心的能力同样得到强化——平方项给了大距离点更高的权重。而 kMedoids 把中心的候选限定在真实样本里。还是刚才这组数据分成两个簇时簇 {1,2,3} 的中心选 2簇 {8,100} 的中心选 8 或 100。无论选哪个中心都是数据集里实打实的观测值不会出现 54 这种人畜无害但毫无意义的数字。就这一个差别离群点对中心的影响就小得多——极端值最多让自己成为中心而不会把中心拖到一个不存在的坐标上。这里要说明一点kMedoids 并不保证聚类结果更好它保证的是中心可解释、结果对离群点更稳健。如果你的数据本身没有离群点且簇形状近似球形K-Means 的速度优势反而值得优先考虑。选型这件事继续往下说。2.2 kMedoids 的核心思想用真实样本做中心kMedoids 的算法骨架和 K-Means 几乎一样初始化 k 个中心点、把每个样本分到最近的中心、用簇内样本更新中心、重复直到收敛。唯一的区别在第 3 步——K-Means 直接计算簇内均值而 kMedoids 在簇内枚举每个候选样本挑到其他所有簇成员的总距离最小的那个点作为新中心。这个点就是 medoid也叫中心点或类中心。正式一点说kMedoids 要最小化的代价函数是每个样本都归到距离它最近的 medoid目标就是让全体样本到各自 medoid 的距离之和尽量小。算法通过迭代调整 medoid 集合来降低这个总代价。和 K-Means 一样这个代价是单调不增的所以算法保证收敛但不保证全局最优。业界最经典的实现是 Kaufman 和 Rousseeuw 提出的 PAMPartitioning Around Medoids算法它的流程分两步。BUILD 阶段先挑 k 个最分散的点作为初始 medoids避免随机初始化的偶然性SWAP 阶段每一轮遍历所有 (medoid, 非 medoid) 对计算交换后总代价的增量如果找到负增量就执行交换。这一轮全量扫描的计算量非常大SWAP 阶段的每轮迭代接近 O(k·(n−k)·n) 的量级所以 PAM 对大数据不友好。很多 MATLAB 代码习惯上把 PAM 的 SWAP 简化成只在簇内交换只尝试当前簇内的点而不试探跨簇的候选中心。这样做收敛更快但理论上可能错过更好的全局中心。我的建议是课程作业、小规模数据集、论文复现用简化版完全够生产环境对聚类质量要求高优先考虑 MATLAB 内置的完整 PAM 实现。另外BUILD 阶段虽然听起来很好手写实现中很多人直接跳过它用随机初始化。后果是初始点质量差时SWAP 阶段要多花很多轮次才能收敛。一个折中做法是执行 5 次随机初始化每次跑 20 轮保留代价最小的那次作为最终结果——用时间换稳定性对小数据集非常划算。2.3 选型决策什么时候用 kMedoids 而不是 K-Means到底什么场景值得牺牲速度去换 kMedoids 的鲁棒性我按经验排了一下。第一数据里明确有离群点或噪声。传感器数据、金融异常值、用户行为日志中那些抖一下的极端值会直接把 K-Means 的均值中心拉走。kMedoids 受的影响小因为离群点最多把自己标成一个 medoid不会成为幽灵中心。第二中心需要可解释。在推荐系统或客户分群场景里如果聚类中心对应着一位真实的用户画像、一条真实的交易记录业务方更容易理解这个簇代表谁。K-Means 的均值向量可能是现实中不存在的人。第三距离度量不是欧氏距离。K-Means 的均值在曼哈顿距离下没有闭合解kMedoids 只需要距离矩阵任何两个点之间能定义的距离都能用。kMedoids 因此常常和余弦距离、动态时间弯曲DTW配对用在文本和时序聚类上。反过来讲三种情况尽量别用 kMedoids数据量大比如超过 10 万条SWAP 阶段的复杂度会让你等到怀疑人生簇形状复杂kMedoids 和 K-Means 一样基于中心辐射假设对环形、嵌套簇都无能为力该换 DBSCAN 就换要聚类的维度特别高高维空间里最近邻的区分度会急剧下降这在避坑章单独说。最后补充一个经验值我在图像分割场景里发现kMedoids 作用于颜色直方图和纹理特征时经常比 K-Means 表现更好代价是每轮迭代多花几十秒。但这个时间成本换来的是分割区域更连续、边界更干净不少 MATLAB 图像处理大作业会用 kMedoids 替代 K-Means 来提升效果。3. MATLAB 手写 kMedoids从零实现一个可运行的聚类器3.1 最小可运行版本核心代码与逐段解释先上一个完整可跑的最小实现。我管它叫kmedoids_diy.m输入输出尽量向 MATLAB 内置kmedoids函数对齐方便你对比测试。function [idx, midx, total_cost] kmedoids_diy(X, k, max_iter) % kmedoids_diy 简易kMedoids聚类PAM简化版簇内交换策略 % 输入 % X : n×d 数据矩阵每行一个样本 % k : 聚类簇数必须大于0且小于n % max_iter : 最大迭代次数默认100 % 输出 % idx : n×1 簇标签值范围1~k % midx : k×1 每个簇的medoid在X中的行索引 % total_cost: 标量所有样本到其medoid的距离之和 if nargin 3 max_iter 100; end [n, ~] size(X); if k n error(k必须小于样本数n); end % 预计算n×n距离矩阵欧氏距离 D sqrt((sum(X.^2,2) sum(X.^2,2) - 2*(X*X))); D(1:n1:end) 0; % 对角元素清零避免自距干扰 % Step 1: 随机初始化k个不同的medoid索引 rng(42); % 固定种子保证结果可复现 midx randperm(n, k); % 迭代变量 old_cost inf; for iter 1:max_iter % Step 2分配: 每个样本归属最近的medoid [~, idx] min(D(:, midx), [], 2); % Step 3更新: 对每个簇尝试用簇内其他点替换当前medoid changed false; for j 1:k members find(idx j); if isempty(members) continue; % 空簇本迭代不处理留到避坑章讨论 end cur_m midx(j); % 当前簇的代价所有成员到当前medoid的距离和 cur_cost sum(D(members, cur_m)); % 遍历每个簇内样本找代价最小的候选中心 best_cost cur_cost; best_m cur_m; for t 1:length(members) cand members(t); if cand cur_m continue; % 自己换自己没有意义 end c_cost sum(D(members, cand)); % 候选点的总代价 if c_cost best_cost best_cost c_cost; best_m cand; end end if best_m ~ cur_m midx(j) best_m; % 执行替换 changed true; end end % 计算本轮总代价用于收敛判断 total_cost sum(min(D(:, midx), [], 2)); if ~changed || abs(old_cost - total_cost) 1e-6 break; end old_cost total_cost; end % 收敛后再做一次最终分配 [~, idx] min(D(:, midx), [], 2); total_cost sum(min(D(:, midx), [], 2)); end这段代码的核心逻辑和 K-Means 几乎同一个骨架差异只在更新中心这一步。先说距离矩阵的向量化写法sqrt((sum(X.^2,2) sum(X.^2,2) - 2*(X*X)))原理是欧氏距离展开式||a-b||² ||a||² ||b||² - 2ab一次矩阵乘法代替双重循环n 不超过一万时性能都还行。D(1:n1:end)0是把每个点和自己的距离清零避免数值误差带来的非零自距。min(D(:, midx), [], 2)这一步是整个算法的核心对第 i 个样本D(i, midx)给出它到所有 k 个 medoid 的距离取最小者所在列号就是它的簇标签。这里的向量化写法比for i1:n的循环快一个数量级建议不要改成循环。再说参数。max_iter100是一个偏保守的上限——PAM 类算法在中小数据集上通常 20 轮以内就能收敛100 轮基本不会触顶。真正的收缩点是随机种子rng(42)你把这个数字去掉每次跑出来的 medoid 初始位置都不同最终结果大概率也不同。这不是 bug而是 kMedoids 的非凸性质决定的避坑章会展开讲。3.2 用合成数据验证聚类效果从二维图上读结果光有函数不行得跑起来看效果。我一般用带离群点的二维高斯混合数据来验收算法——二维便于可视化离群点能检验算法鲁棒性。% 生成合成数据两个高斯簇 一个离群点 rng(10); X1 randn(40,2) * 0.6 repmat([2, 2], 40, 1); X2 randn(40,2) * 0.6 repmat([-2, -2], 40, 1); X [X1; X2; 10, 10]; % 最后一行是故意放的离群点 % 调用自定义kMedoids [idx, midx, cost] kmedoids_diy(X, 2, 100); % 可视化 figure; gscatter(X(:,1), X(:,2), idx, rb, xo); hold on; plot(X(midx,1), X(midx,2), kp, MarkerSize, 15, LineWidth, 2); title(kMedoids 聚类结果黑方块为medoid); legend(簇1, 簇2, Medoid, Location, best);跑完你会看到两个事实。第一离群点 (10,10) 会被独自分到另一个簇而两个高斯簇的分界线基本落在 (−2,−2) 到 (2,2) 之间的对角线上。第二medoid 点黑色方块一定是数据里真实存在的点而不是像 K-Means 那样输出一个坐标均值。这两个特征分别对应前面讲的鲁棒性和可解释性。值得试的对比实验把上面代码里的kmedoids_diy换成 MATLAB 内置的kmedoids函数两个结果在无离群点的数据上几乎一致。但在有离群点时自定义实现和内置函数通常都会正确隔离离群点而 K-Meanskmeans(X,2)则很可能把一个高斯簇劈成两半离群点拽着中心跑偏。这一步对比能帮你直观建立什么时候该信任 kMedoids的判断力。3.3 参数选择的几个要点标准化不是可选项。如果特征量纲差异大比如一列是 0 到 1另一列是 0 到 10000距离矩阵会被量纲大的特征主导。常规做法是X zscore(X)让所有特征方差为 1。这个操作对任何基于距离的聚类都适用不止 kMedoids。距离度量的选择要匹配数据语义。欧氏距离适合连续特征且各方向等权重的场景曼哈顿距离对坐标这类数据更自然文本 TF-IDF 向量适合余弦距离。MATLAB 内置kmedoids的Distance参数支持这些选项自定义实现里只需要替换距离矩阵的计算方式即可。空簇怎么处理。我的简化实现在更新步骤遇到空簇是直接跳过的这在 K-Means 里也有对应的经典问题。更稳妥的兜底策略是分配步骤结束后检查是否出现空簇如果出现就随机选一个非 medoid 样本补充进去。但这个策略要小心它可能打破算法已收敛的状态。4. 用好内置 kmedoids 与 Kmedios.rar下载代码之前先看这几点4.1 内置 kmedoids 函数与自定义实现的对比MATLAB 从 R2018a 开始在 Statistics and Machine Learning Toolbox 里提供了官方的kmedoids函数。接口非常干净% 基本用法X是数据矩阵k是聚类数 [idx, C] kmedoids(X, k); % 带更多输出sumd是每个簇的总距离D是样本到各medoid的距离 [idx, C, sumd, D, midx, info] kmedoids(X, k, Distance, sqeuclidean);其中C返回的是 k×d 的 medoid 坐标矩阵midx返回的是 medoid 在 X 中的行索引。info结构体里包含迭代次数、终止原因等诊断信息调试时多看一眼它能省不少时间。内置函数和自定义实现相比主要强在三点。第一内置的 SWAP 阶段是完整的 PAM 扫描搜索范围不只是当前簇内还会考虑所有非 medoid 候选理论上更容易逼近全局最优。第二Start参数允许你传入初始 medoid 矩阵k×d这对复用上一次结果继续迭代很有价值。第三内置函数针对大数据集提供了Algorithm, clara选项它是对行做抽样、在样本上跑 PAM能扛住几十万行的数据这是手写实现很难做到的。不过内置函数也有个短板它对输入矩阵的数据类型和维度有校验如果你传的是距离矩阵而不是原始数据需要显式换方案——内置kmedoids不支持直接传入成对距离矩阵这一点有时候让人挠头。4.2 Kmedios.rar 的典型结构与使用路径这类压缩包在资源分享站、GitHub 和网盘转存链接里非常常见。解压出来通常是三类文件一个主函数比如kmedoids.m或mykmedoids.m、一份测试脚本或 demo、以及一个说明文档。也有少数包里带示例数据或参考论文 PDF。拿到包的第一件事不是跑是看函数签名% 打开主函数文件后先看第一行function定义 % 常见签名有这几种形态 % [idx, C] kmedoids(X, k) % [idx, C, sumd] kmedoids(X, k, dist_type) % [label, center] KMEDOIDS(X, k, maxiter)签名决定了你调用它时传什么参数、返回什么数据结构。很多翻车现场都是因为没看签名拿着内置函数的调用方式去调一个自定义函数结果参数数量不匹配或输出错位。接下来是把它跑通的三个步骤。第一步确认依赖函数都在 MATLAB 路径上。压缩包里的 .m 文件可能散在不同子目录如果主函数调用了某个子函数而它不在路径上运行时会报 Undefined function。第二步用 demo 脚本里的数据跑一遍如果包里没带 demo你就用第 3 章那段合成数据生成代码自己构造数据。第三步把结果和内置kmedoids对比——同一个数据集、同一个 k 值两组输出的簇标签可能不同因为初始化和标签编号不同但total_cost应该非常接近。如果差得离谱说明下载的代码存在逻辑问题不值得继续投入。4.3 验证一个下载的 kMedoids 代码的 4 个步骤我在拿到任何一个网上下载的聚类代码时会做四件事来确认它没有暗伤。第一构建一个已知答案的数据集。生成 3 个远离的高斯簇标准差距大一点这样任何正常聚类算法都百分百能分开。如果下载的代码在这么简单的数据上都分不对那它显然有 bug。第二测试极端参数。k1时所有样本必须归同一簇kn时每个样本自成簇代价为 0样本点到自身的距离为 0。这两个边界测试能暴露索引越界和距离计算异常的问题。% 边界测试示例 [idx1, ~, cost1] kmedoids_diy(X, 1, 50); % 应全部归为1 [idx2, ~, cost2] kmedoids_diy(X, size(X,1), 50); % 每个样本一个簇 % 检查 cost2 是否约等于0数值误差范围内第三比较收敛代价。对同一个数据集跑kmeans和kmedoids观察 kMedoids 的代价是否低于 K-Means。这不是必然的因为目标函数不同但通常 kMedoids 的代价不会高于同规模聚类的 K-Means 太多。如果发现 kMedoids 代价高出几个量级多半是距离计算或分配逻辑出了问题。第四检查是否有编码问题。很多从国内资源站下载的 .m 文件是 GBK 编码在 UTF-8 环境用 MATLAB 直接打开中文注释会显示成乱码注释中包含引号时甚至报语法错误。解决办法是用文本编辑器把编码转成 UTF-8或者删掉注释只保留代码逻辑再保存。5. 避坑kMedoids 在 MATLAB 里最容易踩的 5 个坑5.1 距离矩阵算错导致聚类结果完全错乱现象聚类结果完全不符合直觉两个明显分开的簇被割裂开或者所有点都被归到同一个簇里。原因大多数手写实现用双重循环逐对计算距离稍不注意索引错位或者在对角元素清零时用了错误的线性索引就会让D(i,i)变成一个不该出现的干扰项。更隐蔽的错误是向量化公式里X*X的形状不匹配——X 是 n×d 时没问题但如果 X 里含有 NaN 或 Inf整个距离矩阵都会出问题。解决先把距离矩阵可视化。imagesc(D)能让你肉眼检查矩阵应该是对称的、对角线为 0、有离群点的数据集总会在某行某列出现明显的高亮。再用一个简单案例做验证[0,0; 1,0]两点之间的距离手算出应该是 1代码跑出来如果不是 1那就是距离计算本身的问题。5.2 随机种子影响结果同样代码两次跑出不同聚类现象同一份代码、同一个数据集第一次跑输出簇标签和代价与第二次完全不同甚至聚类结构都不一样。原因kMedoids 的初始中心是从样本中随机挑选的这是非凸优化问题——不同的初始点会落到不同的局部最优。特别是当你的数据有两个簇距离很近时初始化的微小差异可能把边界处的样本分到不同侧。解决固定随机种子rng(0)或rng(42)确保每次跑一致这是论文复现和课程作业提交的底线要求。如果你想评估这个数据集上 kMedoids 到底稳不稳可以跑 30 次不同的种子记录每次的代价看代价的方差——方差大说明聚类结构本身不够稳定下游分析要小心。内置函数也支持这个思路循环里每次调用前重新rng或删掉种子收集总代价。5.3 空簇问题现象某个 k 值下迭代结束后有一个簇的成员数等于 0medoid 成为一个无兵的将军。可视化时图例里出现空集。原因当初始 medoid 选得太靠近、而数据又有离群点分布时离群点可能单独占一个簇挤压其他簇的空间或者是某个 medoid 在更新过程中被换走之后没有新成员分给它。解决我一般用两招。第一招是初始化约束随机选初始 medoid 时要求两两之间的最小距离大于某个阈值比如所有样本距离的 10% 分位数从根源上降低初始中心扎堆的概率。第二招是运行后检查tabulate(idx)扫一眼成员数分布如果有空簇减少 k 值、更换种子或改用内置kmedoids的Start参数传入更好的初始中心。另外注意内置 kmedoids 遇到空簇时的行为可能和自定义实现不同它可能自动丢弃该簇中心并重建这会导致最终输出只有 k−1 个簇。如果你对簇数有刚性要求比如必须分 5 类空簇问题要在前置环节就规避掉。5.4 Distance 参数选错的连锁反应现象在某个数据集上用欧氏距离聚类效果不错换了个数据集后结果完全不可理喻但数据标准化似乎也没解决问题。原因kMedoids 的距离并不非得是欧氏距离。如果你的特征是词频、直方图或方向向量欧氏距离会把大量信息稀释掉。比如两篇文档各 1000 个词它们共享的词汇只有 20 个欧氏距离会很大但语义上它们可能非常相关。这时候用余弦距离才是正常的。解决根据数据语义选距离。MATLAB 内置kmedoids支持cityblock、cosine、correlation等选项。判断依据很简单算出来距离到底代表什么——对应到实际业务里这个距离大与小是否和人对事物的相似性判断一致。如果一致说明距离选对了如果不一致趁早换。5.5 高维数据的维度灾难现象数据维度从几十升到几千聚类结果开始变得随机代价曲线几乎分不出 k 值的拐点轮廓系数也接近 0。原因高维空间里所有点之间的距离都趋向于同量级欧氏距离的区分度急剧下降。这不是 kMedoids 的错是所有基于距离的聚类算法的通病。解决先降维再聚类。常见做法是 PCA 保留 90% 方差的前若干个主成分或者 t-SNE、UMAP 降到 2~3 维做可视化辅证。但注意——降维会丢失信息聚类结果要在原始维度上做代价评估不要在降维后的坐标上评估。另外选correlation或cosine这类对幅度不敏感的距离度量有时比降维更有效因为它把方向作为主信号天然忽略零均值的噪声维度。还有一个和维度灾难常常同时出现的坑是数据稀疏性。如果你做的是用户-物品偏好矩阵聚类每行只有少数非零元素直接在这个稀疏矩阵上算欧氏距离结果基本等于看谁和谁的零值坐标重合得多。这种场景正确的做法是先做 SVD 降维或使用correlation距离让聚类信号从共现模式里提取而不是从零值里提取。6. 选好 k 值与验证聚类质量收尾的两步实操6.1 肘部法画出代价曲线找拐点kMedoids 要求预先指定 k 值一个实用且直观的做法是肘部法对 k1 到 10 分别跑聚类记录每次的总代价然后画图。代价会随着 k 增加而下降但下降幅度在某个 k 值之后明显放缓那个拐点就是性价比最高的 k。% 肘部法示例 ks 1:10; costs zeros(size(ks)); for i 1:length(ks) [~, ~, costs(i)] kmedoids_diy(X, ks(i), 100); end plot(ks, costs, -o); xlabel(k); ylabel(total cost);注意一点肘部法的拐点在真实数据上不总是明显。如果曲线平滑下降没有拐点说明数据没有天然可分性这时别强行挑 k考虑先对数据做预处理或换聚类算法。6.2 轮廓系数一个本地就能算的质量指标轮廓系数不要求外部标签只计算每个样本和它所在簇的紧密度以及和其他簇的分离度。范围在 -1 到 1 之间接近 1 说明聚类结构清晰接近 0 说明样本在两个簇之间摇摆。轮廓系数的 MATLAB 实现可以直接用evalclusters这条路径但我更推荐手写一次因为能顺便加深对聚类的理解function s silhouette_diy(X, idx) % 简易轮廓系数欧氏距离 k max(idx); n size(X,1); s zeros(n,1); D pdist2(X, X); for i 1:n % 样本i到同簇其他样本的平均距离 a_i mean(D(i, idx idx(i) (1:n) ~ i)); if isempty(a_i), a_i 0; end % 样本i到最近异簇所有样本的平均距离 b_vals arrayfun((j) mean(D(i, idx j)), setdiff(1:k, idx(i))); b_i min(b_vals); s(i) (b_i - a_i) / max(a_i, b_i); end s mean(s); end轮廓系数大于 0.25 就说明聚类结构不是纯噪声了如果小于 0.1要么 k 没选对要么数据本身不适合做中心型聚类得回头检查特征工程。做轮廓系数评估时同样固定随机种子不同初始化的轮廓系数差异如果超过 0.05就多跑几次取平均。往深了说我养成的一个习惯是每次跑完 kMedoids除了看轮廓系数还会把紧密度最高的几个簇内样本单独拉出来看原始特征——聚类是工具不是答案真正交付给业务方的是这个簇代表什么含义。如果簇内样本在业务维度上完全说不过去哪怕轮廓系数是 0.8这个聚类结果也值不了几分钱。说到底算法就那几行循环真正值钱的是你对自己数据的理解。这次就聊到这儿希望帮到你。本文还有配套的精品资源点击获取