k-medoids聚类算法MATLAB实现与代码详解：从原理到可视化

看到这份标题，第一反应是很典型的“课程作业/入门项目”配置——k-medoids 聚类算法的 MATLAB 源代码，数据导入和画图部分带中文注释。别看它结构简单，这套东西背后把聚类算法的几个关键环节全串起来了：数据读入、距离计算、中心点更新、结果可视化。我已经帮多位同学把这个模型从“能跑通”调到“能讲清”，今天就把这份代码的完整拆解、踩坑记录和可复现步骤一起整理出来。

先说清楚适用范围：这份代码适合入门机器学习、模式识别课程设计，或者需要在 MATLAB 里快速做无监督聚类分析的开发者。如果你的数据量在几千条以内、维度不超过 5 维，这份代码的训练速度和内存占用都能接受。要处理百万级样本，那就得换 Mini-Batch 方案或上 Python 版本了。

1. 项目背景与 k-medoids 算法原理拆解

1.1 为什么选 k-medoids 而不是 k-means

k-means 是大多数人第一个接触的聚类算法，它的核心计算是“簇内均值”。均值这个东西有两个天然缺陷：第一，它会被离群点带偏。比如某个簇里混入一个极端值，簇中心会被硬生生拉过去，导致边界划分错位。第二，均值点是一个虚拟点，不一定对应真实存在的样本，后续分析时想“看看这个簇的代表长什么样”，得额外找最近的真实点，很麻烦。

k-medoids 的做法完全不同：它强行要求每个簇的中心必须是簇内真实存在的样本点，这个点就叫 medoid（中心样本）。每次迭代只在“现有样本”里选代表，天然免疫离群值，而且聚类结果可以直接解释成“这个点就是本簇最典型的用户/商品/案例”。

如果你处理的是金融交易记录、传感器异常数据、图像特征向量这类离群值占比高的数据，k-medoids 的稳定性优势非常明显。当然代价也有——每次更新中心点要做多次相异度比较，计算量比 k-means 大一截。

1.2 PAM 算法的两步迭代逻辑

MATLAB 源码里最常见的 k-medoids 实现是 PAM（Partitioning Around Medoids）算法，它的思路拆开就两步：

构建阶段的目的是选初始中心。做法是遍历所有样本，选出一个对象作为第一个中心，使得它到其余所有样本的相异度总和最小，然后反复挑选能使“整体相异度增量”最小的新样本加入中心集，直到攒够 k 个中心。

交换阶段则是循环优化：对每一对“当前中心点 m_i”和“非中心点 h”，计算把 m_i 换成 h 之后代价的变化量。如果代价减小，就执行替换，重复这个操作直到任何替换都无法再降低代价。

这里的“代价”通常定义为：每个普通样本到其所属最近中心点的距离之和。代码里核心循环就是围绕这个代价函数转的。理解了这两个阶段，你就知道了源码里 200 行左右的逻辑主线是什么——它不像神经网络那样有复杂的反向传播，所有核心都集中在“替换中心点是否划算”这一个判断上。

提示：PAM 算法对初始中心比较敏感，不同初始点导致最终聚类结果差异明显。代码里如果没做多次随机初始化取最优的操作，建议自己加上，后面我会给具体做法。

2. 代码框架与数据导入部分的实现解读

2.1 完整代码模块划分

整份源码按功能可以切成四个模块，从上到下依次是：数据准备模块、算法初始化模块、PAM 迭代主循环模块、结果可视化模块。数据准备和可视化是注释最密集的区段，这是有意的安排——这两块高度依赖外部环境（文件路径、坐标系、数据范围），不写清楚使用者根本没法改。

数据准备模块的核心任务就一句话：把存储在外部文件里的矩阵读进 MATLAB 工作区，并做初步清洗。这里涉及三个问题：文件格式、路径写法、空值处理。源码里默认支持 .xlsx、.txt 和 .mat 三种格式，用一行 switch 分支控制。

2.2 数据加载区注释翻译与设计要点

先看导数据这一段，典型的源码注释和实现逻辑大致是这样：

% 导入数据部分 % 请根据你的实际数据文件格式修改文件名 file_name 和路径 file_path % 数据格式要求：每一行是一个样本，前 d 列是特征，最后一列可以是标签（可省略） file_path = './data/'; file_name = 'iris_data.xlsx'; data = readmatrix([file_path, file_name]); % 如果文件最后一列是标签，先剥离标签列，仅保留特征矩阵 if size(data, 2) > size(data, 2) - 1 % 占位写法，实际以代码为准 X = data(:, 1:end-1); true_labels = data(:, end); else X = data; end

设计要点有三个。第一是readmatrix这个函数，它在较新版本的 MATLAB（R2019b 之后）里能自动识别数值和文本混合的 Excel 表格，老代码里常用的xlsread已经标了 deprecated 状态，新写代码建议直接上readmatrix。第二是“剥离标签列”这个操作，很多公开数据集最后带真实标签，聚类是无监督算法不需要标签参与训练，但后面对比聚类准确率时要用，所以注释里专门提示用户保留。第三是./data/这种相对路径写法，它假设你当前工作目录和控制脚本在同一级，避免了绝对路径在换电脑后失效的坑。

这里面最容易被忽略的是数据标准化。如果不同特征列的数值范围差异很大（比如一列是 0~1，另一列是 0~10000），距离计算会被大数值列主导，聚类效果几乎必然不理想。源码注释里如果没有提醒这一步，我强烈建议你导入数据后加一行：

% 标准化处理：每个特征列减去均值并除以标准差，消除量纲影响 X = (X - mean(X)) ./ std(X);

把标准化注释也写进中文注释里，之后调参时不容易忘。

3. 核心聚类过程的 MATLAB 实现

3.1 初始化策略

PAM 的初始化有两种常见写法。第一种是随机抽 k 个样本做 medoid，实现最简单，但结果方差大。第二种是“最大最小法”：第一个中心选离所有样本最远的点，之后每个新中心都选“离已有中心最近距离最大”的样本，这样中心散布较均匀，比纯随机稳定得多。

源码里通常提供的是随机初始化，注释里会提示运行多次取最优。我的建议是直接改成最大最小法，代码量增加不超过十行：

% 初始化 medoid：采用最大最小距离策略，使初始中心尽可能分散 rng(42); % 固定随机种子，保证结果可复现 medoid_idx = zeros(k, 1); % 先选距离矩阵总和最大的样本作为第一个中心（简化为随机选第一个） medoid_idx(1) = randi(n); for i = 2:k % 计算每个非中心样本到已有最近中心的距离，取最小值中的最大值 dist_to_center = min(dist_mat(:, medoid_idx(1:i-1)), [], 2); [~, new_idx] = max(dist_to_center); medoid_idx(i) = new_idx; end

固定随机种子（比如rng(42)）是个好习惯。它能保证你每次运行结果一致，方便调试和比较实验——尤其在做参数调优时，不确定的初始点会让前后结果没法公平对比。

3.2 分配与更新循环

主循环的代码逻辑是容易绕晕的地方，我把伪代码结构抽出来看：

主循环分为两个阶段：

• 分配阶段：根据当前 medoid 集，把每个样本归到距离最近的 medoid 所属簇。这一步用距离矩阵的行比较实现。源码里最常见的做法是先预计算一个n x k的矩阵，第 i 行存的是样本 i 到每个 medoid 的距离，然后取最小值所在列索引。

• 更新阶段：对每个簇，遍历簇内所有样本，计算把它设为新 medoid 后簇内总距离的变化量。若总距离更小，则替换。这个阶段最耗时，因为涉及“簇内每个点都试一遍”的穷举式搜索。

更新阶段的计算量能占到整个程序的 70% 以上。如果样本量大，这一段建议改成向量化实现——用矩阵运算替掉内层循环，而不是一个一个算距离：

% 更新 medoid 的向量化近似写法 for c = 1:k members = (assignments == c); cluster_data = X(members, :); % 计算簇内两两距离矩阵，选行和最小的样本作为新 medoid intra_dist = pdist2(cluster_data, cluster_data, 'squaredeuclidean'); [~, new_medoid_in_cluster] = min(sum(intra_dist, 2)); new_idx = find(members); medoid_idx(c) = new_idx(new_medoid_in_cluster); end

这样写的好处有二：一是pdist2内部用 BLAS 优化过的矩阵乘法算距离，比自己写双重 for 循环快不少；二是行和最小的样本理论上就是到簇内所有样本距离和最小的点，正符合 medoid 的定义。代价是额外开了一整块n^2的内存，数据量上万时留意内存占用。

3.3 距离矩阵计算的向量化技巧

距离度量直接决定聚类形态。源码里常见写法是支持三种距离：欧氏距离、曼哈顿距离、余弦距离。欧氏适合连续实数特征；曼哈顿（L1）在高维稀疏数据上更稳，它不像欧氏那样对大数值差敏感；余弦适合文本向量或用户行为向量这类重视方向、不重视绝对大小的数据。

三种距离的向量化 MATLAB 写法供参考：

% 欧氏距离：||a-b||^2 = ||a||^2 + ||b||^2 - 2*a*b' dist_mat = sqrt(sum(X.^2, 2) + sum(X.^2, 2)' - 2 * (X * X')); % 曼哈顿距离：按列分别计算绝对差再求和（用循环实现） % 余弦距离：1 - 向量夹角的余弦值 dist_cos = 1 - (X * X') ./ (sqrt(sum(X.^2, 2)) * sqrt(sum(X.^2, 2))');

欧氏距离的展开式是一个很常用的技巧，原理是二次展开避免双循环，网上搜“距离矩阵向量化”几乎都会看到这一行。注意X * X'得到的矩阵要小心对角线数值误差——因为浮点计算可能导致对角线出现 -1e-15 这种极小负数，开方前先max(0, dist_sq)钳一下。

4. 可视化模块：画图部分的注释与配色设计

4.1 二维散点图绘制

聚类结果的可视化通常针对二维数据，源码里的画图部分是把每个簇画成不同颜色散点，中心点用大号星标标出。核心用scatter函数，分散参数设置如下：

• S（点大小）：设为 30，太小会被整体淹没，太大在样本密集处会粘连。
• filled：填充散点颜色，比空心圆清晰很多。
• MarkerFaceAlpha：设置透明度为 0.6，样本重叠多时能看出密度层次。
• LineWidth和MarkerEdgeColor：中心点要突出，边框加粗、颜色用黑边。

画图部分的主要代码结构类似这样：

% 画图部分 % 画出聚类结果，不同簇用不同颜色区分，medoid 用星号突出标记 figure; hold on; for c = 1:k scatter(X(assignments == c, 1), X(assignments == c, 2), 30, colors(c, :), 'filled', 'MarkerFaceAlpha', 0.6); end scatter(X(medoid_idx, 1), X(medoid_idx, 2), 200, 'k', 'p', 'LineWidth', 2); title(['k-medoids 聚类结果，k = ', num2str(k)]); xlabel('特征 1'); ylabel('特征 2'); grid on; hold off;

这里有个注释要特别提醒：colors这个变量是代码内部定义的一个k x 3矩阵，每行是一个 RGB 三元组。如果你的 k 大于 7，直接用 MATLAB 默认的lines或parula色带抽取颜色即可，否则手工写 10 种以上的 RGB 值很容易撞色。

4.2 中心点突出显示与图例

画中心点有个细节：星标'p'只有 200 point size 时才足够醒目，但也会遮挡周围样本。实测最好的做法是分两层画——先画一个带透明度的浅色大圆作为光环，再在上面画实心星标。这样既突出中心又不挡信息：

% 画中心点光环：半透明的大空心圆，让星标更突出且不遮挡附近样本 scatter(X(medoid_idx, 1), X(medoid_idx, 2), 400, colors(c, :), 'o', 'LineWidth', 1.5, 'MarkerFaceAlpha', 0.15); scatter(X(medoid_idx, 1), X(medoid_idx, 2), 200, 'k', 'p', 'LineWidth', 2);

图例方面，如果是人工标注簇名称的场景可以直接用legend加字符串元胞数组。但聚类结果通常没有先验名称，建议在标题里写清楚是第几次运行、k 值和距离度量，图例反而可以省略。可视化最重要的目标是让读者一眼判断“簇分得开吗”“边界在哪”，而不是精确读出每个点的编号。

高维数据没法直接画二维散点，这时候可以在画图模块前加一步 PCA 降维，保留前两主成分再画，注释里记得说明这属于“降维可视化，不等同于聚类输入”。

5. 常见问题排查与调试笔记

5.1 问题速查表

我在实际使用中把这些年跑 k-medoids 遇到过的问题做成速查表，基本覆盖了新手会踩的坑：

5.2 防止局部最优的三个操作习惯

PAM 是局部搜索算法，理论上不保证全局最优，实操层面有三个有效习惯：

第一个是多次随机重启。同一个 k 值多跑 20 次，每次用不同随机种子，记录目标函数（总距离和）最小的那一次作为最终结果。代价是耗时变成 20 倍，但聚类质量在大部分数据集上明显稳定。

第二个是“先用 k-means 找中心位置，再替换成最近的真实点”作为初始化。具体做法是把 k-means 跑一遍，得到 k 个均值中心，然后每个均值找离它最近的原始样本作为 k-medoids 的初始中心。这个方案在 sklearn 的KMedoids里叫k-medoids++，实测收敛一致性大幅提高。

第三个是兜底策略：朴素地检查每次迭代是否有中心点重复。如果第 i 个和第 j 个簇的中心相同，那就是初始化失败或者 k 选大了，果断换初始化重新跑，别让代码带着明显错误继续跑完所有迭代。

5.3 中文注释和显示乱码的处理

这份源码最大的特色就是中文注释，所以乱码问题必须单独讲。MATLAB 的编辑器默认按本地编码（通常是 GBK）打开文件，但源码文件保存成 UTF-8 时，中文字符在旧版本里就可能显示成乱码。

好用的处理步骤是：先用记事本打开源码文件，另存为“带有 BOM 的 UTF-8”编码，再用 MATLAB 打开重试。更省事的办法是在 MATLAB 命令行执行feature('DefaultCharacterSet', 'UTF-8')切换默认编码，但这只对当前会话有效，下次启动要重新执行。

另外建议在源码这样写注释：

% 功能：对输入样本矩阵 X 执行 k-medoids 聚类 % 输入：X —— n x d 样本矩阵，n 为样本数，d 为特征维度 % k —— 聚类簇数，须为小于 n 的正整数 % 输出：medoid_idx —— k x 1 向量，中心样本的行索引 % assignments —— n x 1 向量，每个样本所属的簇编号 % total_cost —— 数值，最终簇内总相异度（越小说明聚类越紧密）

把人话写进注释，几个月后回看既有“现在在干什么”也有“为什么这么干”。这也是我改别人代码时最看重的一点——不是注释越多越好，而是关键决策处必须写清动机。总结一句：这份代码值得留作基础模板，把距离度量参数扩展到三种，初始化换成最大最小距离，迭代过程改成向量化，再做中文注释，就是一个能放进项目文档和课程报告的完整方案。