k-medoids聚类MATLAB实现：抗离群点聚类源代码与可视化全流程

平时用MATLAB做聚类分析，绕不开k-means，但一旦数据里混了几个离群点，k-means的均值中心就会被拽得七荤八素。这时候该换k-medoids了。我在实际项目里经常碰到这种场景：传感器数据偶尔跳一个异常值，用户行为数据带点噪声标签，用medoid（簇内真正存在的样本）代替均值，聚类中心就不会被离群点带偏。这篇博文就把我整理的k-medoids聚类MATLAB源代码、数据导入和图形绘制整套流程放出来，全程带中文注释，拿过去就能改能跑，适合正在做数据分析、模式识别实验，或者毕业论文需要聚类对比的读者。我会把每一段代码的来龙去脉、参数选择、画图细节都讲透，尽量减少你踩坑的时间。

1. k-medoids聚类原理与选型分析

1.1 为什么用medoid替代centroid

k-means和k-medoids的差异，本质上是“平均数”和“中位数”的差异。k-means每个簇的中心是簇内所有样本的算术平均，可能计算出一个人造坐标；k-medoids的中心则是簇内某个真实存在的样本点，该点是簇内所有样本到其他样本总距离最小的那个。这个区别带来的直接影响就是鲁棒性：k-means对离群点极度敏感，一个偏离很远的点就能把centroid拉过去，导致几个簇被错误合并；k-medoids不会，因为中心必须选在真实数据上，离群点通常孤立成簇或对中心选择影响有限。

从距离度量上说，k-medoids对距离函数的要求更宽松。k-means在欧氏距离下有闭式解（求均值），但如果你需要做曼哈顿距离、余弦距离甚至自定义相似度，k-means很难推广——均值在非欧空间没有定义。k-medoids只要求“能算样本两两距离”，然后枚举簇内样本选最优中心，所以处理分类变量、混合类型数据、缺失值较多的数据时反而更方便。当然代价是计算量：PAM（Partitioning Around Medoids）的经典实现复杂度大约是O(k(n-k)²)，比k-means重不少。下面会讲怎么在MATLAB里平衡效果和效率。

1.2 算法核心流程与收敛性理解

k-medoids的核心步骤只有四步，但每一步都有细节讲究：

1. 初始化：从n个样本里随机选k个作为初始medoid，或者用启发式方法选分散的样本。
2. 分配：把每个样本分配到距离最近的medoid所在簇。
3. 更新：对每个簇，遍历簇内所有样本，选一个能使簇内总距离（样本到中心距离之和）最小的样本作为新medoid。
4. 收敛判断：如果medoid集合不再变化，或者代价函数的变化量小于阈值，就停止迭代。

值得注意的收敛性问题是：k-medoids的目标函数是“所有样本到其medoid距离的总和”，这个函数在枚举更新下是单调递减的，所以算法保证收敛到局部最优，但可能是局部最优而非全局最优。初始点选得不好，很容易收敛到差的结果。这也是为什么测试代码时要多跑几次随机初始化，或者用上一节说的k-means++式启发式初始化。

1.3 适用场景：什么样的数据值得用k-medoids

k-medoids不是要完全取代k-means。经验上这几类场景优先选k-medoids：

• 数据存在明显离群点，且离群点不是你需要单独聚出来的噪声类，而是混在正常样本附近。
• 距离度量是曼哈顿距离或自定义相似度，无法直接求均值。
• 聚类中心后续需要解读，而你的“中心”必须是真实对象。比如客户聚类希望中心对应一个真实客户画像，商品聚类希望中心是某件实际商品。
• 样本量不太大，基本在几千到一两万量级，k也不太大。如果你有百万级样本，建议先抽样或者用Mini-batch策略，否则迭代速度会很难看。

2. MATLAB数据导入与预处理实战

2.1 从Excel、CSV、TXT导入数据

MATLAB读取外部数据的主要函数有readmatrix、readtable、xlsread（老版本）、load。其中readmatrix在R2019a之后是首选，它自动识别数值和文本，返回纯数值矩阵，省去table2array转换。我平时用得最多的代码是：

% 读取CSV文件，跳过第一行表头，只取数值列 data = readmatrix('iris_data.csv', 'NumHeaderLines', 1); % 或者读取Excel指定工作表 data = readmatrix('dataset.xlsx', 'Sheet', 'Sheet1', 'Range', 'A2:C150');

如果数据里第一列是样本ID或标签（文本），需要分离：

data_all = readtable('dataset.xlsx'); labels = data_all{:, 1}; % 文本标签 features = data_all{:, 2:end}; % 数值特征 features = varfun(@str2double, features); % 如果需要转换

从txt导入时经常遇到分隔符问题，readmatrix默认按逗号，但如果txt是用空格或Tab分隔，要指定'Delimiter', '\t'或'Delimiter', ' '。踩过多次坑，建议先打开txt预览一下第一行，确认分隔符和表头位置再写导入代码。

2.2 数据标准化：聚类前必须做的事

k-medoids的距离计算完全依赖量纲。如果第一个特征范围是0到1，第二个特征范围是0到10000，第二个特征会主导整个距离，第一个特征等于没参与。所以在聚类前我一般做z-score标准化，也就是每个特征减去均值除标准差：

mu = mean(data); sigma = std(data); data_std = (data - mu) ./ sigma;

还有一种min-max归一化，把数据缩放到[0,1]区间：

data_min = min(data); data_max = max(data); data_norm = (data - data_min) ./ (data_max - data_min);

两种方法选哪个？如果你的数据分布近似高斯，z-score更稳；如果有硬边界，min-max更直观。需要注意：标准化必须在划分训练集和测试集之前或之后？聚类是无监督任务，不存在训练测试划分，所以直接对全部数据标准化即可。但要保存mu和sigma，以便新样本加入时用同样的参数变换。

2.3 缺失值与异常值的处理策略

缺失值在MATLAB里通常表现为NaN。直接用含NaN的数据算距离，结果会全部变成NaN。处理方案分几种：

• 如果缺失比例小于5%，可以直接删除对应行：

  data_clean = data(~any(isnan(data), 2), :);

• 用列均值或中位数填充：

  for j = 1:size(data, 2) col = data(:, j); col(isnan(col)) = median(col(~isnan(col))); data(:, j) = col; end

• 如果你的特征是分类变量，可以用众数填充，不过k-medoids本身不太吃这个亏，因为它能处理非欧距离。

异常值的处理反而有讲究：k-medoids本身抗离群点，所以不建议轻易删除异常值，删了反而会丢失信息。我一般先画出箱线图或做z-score检测看一下异常点占比，如果异常点是某种测量错误就剔掉，如果是真实分布的一部分就保留。保留的情况正好发挥k-medoids的鲁棒性优势。

3. 核心源代码实现与中文注释详解

3.1 主函数结构设计与参数定义

我写的k-medoids实现，主函数输入是特征矩阵X和簇数k，输出是簇分配标签idx、medoid索引、每轮代价函数历史。函数头长这样：

function [idx, medoid_idx, cost_history] = myKmedoids(X, k, maxIter, distType) % myKmedoids 基于PAM思路的k-medoids聚类实现 % 输入： % X - 样本特征矩阵，每一行是一个样本，每一列是一个特征 % k - 聚类簇数 % maxIter - 最大迭代次数，默认100 % distType - 距离度量类型，'euclidean'或'manhattan'，默认'euclidean' % 输出： % idx - 每个样本所属簇的标签，n行1列 % medoid_idx - 每个簇的medoid在原始数据中的行号 % cost_history - 每轮迭代的代价函数值，用于绘制收敛曲线

这个函数接口是通用的，后续你换数据集、换簇数都不需要改动主循环。参数默认值用nargin判断，在MATLAB里比较常见：

if nargin < 3 maxIter = 100; end if nargin < 4 distType = 'euclidean'; end

3.2 初始化策略：随机选择与启发式初始化

PAM算法的原始初始化是随机抽k个样本做medoid，但这样在多峰数据上容易掉进局部最优。我在代码里加了两个选项：随机初始化和基于最大距离的初始化（类似k-means++的思路）。第二种做法是：第一个medoid随机选，后续每个medoid选择离已有medoid最远的样本。这个策略能让初始中心尽量分散，收敛稳定不少。

n = size(X, 1); if strcmpi(initMethod, 'random') medoid_idx = randperm(n, k)'; elseif strcmpi(initMethod, 'maxdist') % 第一个中心随机 medoid_idx = randi([1, n], 1, 1); % 后续中心选择距离已有中心最远的样本 for j = 2:k D = pdist2(X, X(medoid_idx, :), distType); minDist = min(D, [], 2); [~, idx_new] = max(minDist); medoid_idx = [medoid_idx; idx_new]; end medoid_idx = medoid_idx'; end

需要说明，maxdist初始化会比随机初始化多一次n×k的距离计算，但换来的是迭代次数明显减少，整体耗时反而更低。对于n在几万以下的数据，这个开销完全可以接受。

3.3 计算距离矩阵与簇分配

分配阶段需要计算每个样本到k个medoid的距离。用MATLAB内置的pdist2一次性算距离矩阵，又快又简洁：

D = pdist2(X, X(medoid_idx, :), distType); [~, idx] = min(D, [], 2);

pdist2的distType支持'euclidean'、'squaredeuclidean'、'cityblock'（曼哈顿）、'cosine'、'correlation'等。如果你需要自定义距离，就自己写一个函数句柄传入pdist2，或者干脆用三重循环算距离矩阵。但三重循环在MATLAB里很慢，能向量化尽量向量化。

这里有个细节：pdist2对大矩阵比较吃内存，n×k的矩阵还好，但当n到了几十万就要考虑分块计算。我还没在这个代码里做分块，如果你的数据量特别大，建议把距离计算改成for循环按块处理，避免内存爆掉。

3.4 更新medoid：簇内枚举与总距离最小化

更新阶段是k-medoids和其他算法最不同的地方。对每一个簇，取出簇内所有样本，计算簇内样本两两之间的距离，再选出“到其他样本总距离最小”的那个样本作为新medoid。这里不能用pdist2对整个簇矩阵做两步操作，直接算平方距离矩阵再求和：

for j = 1:k cluster_idx = find(idx == j); if isempty(cluster_idx) % 如果某个簇为空，重新随机指定一个样本 cluster_idx = randi([1, n], 1, 1); idx(cluster_idx) = j; end % 簇内样本的特征矩阵 cluster_X = X(cluster_idx, :); % 簇内两两距离矩阵（欧氏距离的平方也等价） Dc = pdist2(cluster_X, cluster_X, distType); totalDist = sum(Dc, 2); % 每个样本到簇内其他样本的总距离 [~, bestLocalIdx] = min(totalDist); medoid_idx(j) = cluster_idx(bestLocalIdx); end

这里有一个容易踩的坑：sum(Dc, 2)会把样本到自己距离0也算进去，但0不影响最小值比较，所以没问题。另一个坑是如果簇内有重复样本（完全一样的行），pdist2会产生0距离，可能导致totalDist偏小，但不影响选出正确medoid。

3.5 收敛判断与代价函数计算

收敛判断我用了双条件：medoid不再变化，或者达到最大迭代次数。每次迭代计算一次代价函数：

cost = sum(min(D, [], 2)); % 所有样本到最近medoid的距离和

如果新老medoid集合完全一致，说明已经稳定，直接break。如果代价函数的变化小于某个阈值，比如1e-6，也可以提前终止。我在代码里用一个isequal判断：

if isequal(medoid_idx, prev_medoid_idx) break; end prev_medoid_idx = medoid_idx; cost_history = [cost_history; cost];

有人会问：为什么不用代价函数变化阈值？因为medoid是离散选择，代价变化不一定是平滑下降的，可能出现两个不同medoid集合代价几乎相同。用medoid索引变化的判断更直觉更稳。

3.6 完整源代码（中文注释版）

下面给出完整的函数实现，可以直接保存为myKmedoids.m使用：

function [idx, medoid_idx, cost_history] = myKmedoids(X, k, maxIter, distType, initMethod) % myKmedoids 基于PAM思路的k-medoids聚类实现 % 输入： % X - 样本特征矩阵，每一行是一个样本 % k - 聚类簇数 % maxIter - 最大迭代次数，默认100 % distType - 距离度量，'euclidean'或'manhattan'，默认'euclidean' % initMethod - 'random'或'maxdist'，默认'random' % 输出： % idx - 每个样本所属簇的标签 % medoid_idx - 每个簇的medoid在原始数据中的行号 % cost_history - 每轮迭代的代价函数值 % 参数默认值处理 if nargin < 3 maxIter = 100; end if nargin < 4 distType = 'euclidean'; end if nargin < 5 initMethod = 'random'; end n = size(X, 1); % 初始化medoid if strcmpi(initMethod, 'random') medoid_idx = randperm(n, k)'; elseif strcmpi(initMethod, 'maxdist') medoid_idx = randi([1, n], 1, 1); for j = 2:k D = pdist2(X, X(medoid_idx, :), distType); minDist = min(D, [], 2); [~, idx_new] = max(minDist); medoid_idx = [medoid_idx; idx_new]; end medoid_idx = medoid_idx'; end cost_history = []; prev_medoid_idx = medoid_idx; for iter = 1:maxIter % 分配：每个样本归属最近的medoid D = pdist2(X, X(medoid_idx, :), distType); [~, idx] = min(D, [], 2); % 更新：对每个簇选择总距离最小的样本作为新medoid for j = 1:k cluster_idx = find(idx == j); if isempty(cluster_idx) cluster_idx = randi([1, n], 1, 1); idx(cluster_idx) = j; end cluster_X = X(cluster_idx, :); Dc = pdist2(cluster_X, cluster_X, distType); totalDist = sum(Dc, 2); [~, bestLocalIdx] = min(totalDist); medoid_idx(j) = cluster_idx(bestLocalIdx); end % 计算当前代价 D = pdist2(X, X(medoid_idx, :), distType); cost = sum(min(D, [], 2)); cost_history = [cost_history; cost]; % 收敛判断：medoid集合是否变化 if isequal(medoid_idx, prev_medoid_idx) break; end prev_medoid_idx = medoid_idx; end end

这段代码加起来不到60行，核心逻辑是清晰的。如果你要对比实验，只需要把初始化的随机种子固定（rng(42)之类），就能复现结果。

3.7 复杂度分析与调优技巧

k-medoids的每次迭代复杂度是O(knF)，其中F是特征维度。更新阶段每个簇内算pdist2矩阵的开销是O(n^2/F)，整体在k个簇上更接近O(n^2)，因此样本量是主要瓶颈。实测数据：n=10000、k=5、特征维度10，欧氏距离，迭代30轮在我的笔记本上耗时大约8秒。这个量级日常实验完全没问题。

如果想提速，有一个替代方案：更新阶段不需要计算簇内所有样本两两距离，只要在簇内每个样本到其他样本的距离时，用距离公式直接算，不要缓存整矩阵。但MATLAB向量化之后，一次性算pdist2往往比循环快，这个取舍要在效率上自己测一下。另一个方案是：每轮分配只用部分样本估算，也就是mini-batch，代价是收敛会抖一点，但大样本时值得。

4. 图形绘制与分类结果可视化

4.1 二维散点图绘制与medoid标记

聚类做完不画图等于白做。我的标准可视化代码是这样的：

figure; gscatter(X(:,1), X(:,2), idx); hold on; plot(X(medoid_idx, 1), X(medoid_idx, 2), 'kx', 'MarkerSize', 12, 'LineWidth', 2); hold off; legend('簇1', '簇2', '簇3', 'Medoid'); title('k-medoids聚类结果');

gscatter会自动给不同类别分配颜色，省去手动配色。plot用黑色叉号标记medoid，和散点区分明显。这里有个小细节：如果X的特征维度不止2，画图前先做主成分分析降维：

[coeff, score] = pca(X); X2D = score(:, 1:2); % 取前两个主成分 gscatter(X2D(:,1), X2D(:,2), idx);

降维后的可视化只用于展示，聚类过程仍用原始高维特征。

4.2 轮廓图评估聚类质量

聚类质量不能光靠肉眼，轮廓系数（silhouette）是常用指标。MATLAB内置silhouette函数直接能用：

figure; silhouette(X, idx, distType); title('K-medoids聚类轮廓图');

轮廓系数范围是[-1,1]，越接近1说明样本离自己簇内近、离别的簇远。绘制出的图如果大部分样本的轮廓值在0.5以上，说明聚类结构清晰；如果很多值接近0甚至负数，说明有样本被分错了位置。可以在不同k值下跑多次，观察轮廓均值的变化趋势来选k。这个步骤建议固定随机种子再比较，否则每次结果可能不同。

4.3 代价函数收敛曲线

每次运行保存的cost_history可以直接画收敛曲线，判断算法是否稳定：

figure; plot(1:length(cost_history), cost_history, '-o', 'LineWidth', 1.5); xlabel('迭代轮次'); ylabel('代价函数值'); title('k-medoids收敛曲线'); grid on;

正常情况下代价函数单调下降并趋于平稳。如果曲线出现反复震荡，多半是初始化太差或者k值不合适。如果迭代20轮代价还在明显下降，说明maxIter设置小了，建议调大。

4.4 中文注释乱码与绘图中文显示问题

这里专门说一下MATLAB中文显示的两个坑。第一个是源代码里的中文注释保存后乱码，常见于直接用记事本保存的.m文件，MATLAB默认按系统编码读取，一旦文件编码是UTF-8而系统区域是GBK，中文注释就全乱。解决办法：在MATLAB编辑器里设置文件编码为UTF-8，或直接在preferences中调整，另外新建脚本用MATLAB编辑器保存而不是外部编辑器。如果代码已经乱码了，用MATLAB的edit菜单里重新以正确编码打开一次就能恢复。

第二个是绘图的标题和图例中文显示成方框，这是字体问题。绘图之前加一行：

set(0, 'DefaultAxesFontName', 'SimHei');

或者用set(gca, 'FontName', 'SimHei')。Windows一般有SimHei（黑体）和Microsoft YaHei（微软雅黑），macOS用'PingFang SC'。加上这个设置，标题里的中文才能正常显示。

5. 常见问题与排查经验速查

5.1 聚类结果不稳定怎么办

k-medoids对初始值敏感，每次运行结果可能都不一样。解决思路有三个：

1. 固定随机种子做复现：运行前执行rng(42)，保证每次结果一致。
2. 多次运行取最优：循环运行比如20次，保留代价最小的结果。
3. 改用maxdist初始化，稳定性明显好于随机初始化。

如果多次运行后代价差异仍然很大，说明数据本身聚类结构不清晰，或者k值选得不合适。这时候先画一下数据分布，或者尝试不同k值。

5.2 数据导入报错与乱码问题汇总

数据导入最常见的几个报错我在表里汇总一下：

以我自己的经验，CSV导入踩坑率最高。有些软件导出CSV时会带BOM头，MATLAB某些版本读BOM会出现第一个变量名多出字符。解决办法是用fopen+fgetl读第一行手动清洗，或者用readtable的'PreserveVariableNames'选项。实在不行，就用Excel另存为xlsx格式，readmatrix对付xlsx更省心。

5.3 medoid更新后簇变空的情况

在更新步骤中，如果某个簇的成员在重新分配后为空，需要处理，否则下一步会报错。我的代码里已经加了随机填充逻辑：

if isempty(cluster_idx) cluster_idx = randi([1, n], 1, 1); idx(cluster_idx) = j; end

更严谨的做法是：找到当前距离矩阵中离所有medoid最远的样本，把它强制归入空簇。但随机填充在多数情况下也能用，因为下一步的分配会重新调整归属。

还有一种情况：簇内只有一个样本，它的medoid就是本身，距离为0，更新后不变，正常。

5.4 高维数据下的可视化与评估技巧

高维数据没法直接画散点图，我的做法是先画距离矩阵热力图。把样本按聚类标签排序后，用imagesc绘制样本间距离矩阵，能直观看到块状结构。代码片段：

Dfull = pdist2(X, X, distType); [~, order] = sort(idx); imagesc(Dfull(order, order)); colormap(jet); colorbar;

对角线附近的方块颜色越深说明簇内紧凑，方块之间的边界颜色越亮说明簇间分离明显。这个方法比散点图更能反映真实聚类质量。另外还可以用簇间距离/簇内距离的比值（Davies-Bouldin指数）做定量评估。

6. 完整项目文件结构与调用示例

6.1 文件清单与各模块职责

一个完整的k-medoids实验项目，我建议这么组织文件：

project/ ├── myKmedoids.m % 主聚类函数 ├── loadData.m % 数据导入封装 ├── visualizeResult.m % 可视化封装 ├── evaluateCluster.m % 轮廓系数等评估 ├── runExperiment.m % 主脚本：调用所有模块 └── data/ └── iris.csv % 实验数据

把导入、聚类、可视化、评估拆成独立函数的好处是：换数据集时只需要改loadData.m，换算法时只改runExperiment.m里的调用，可视化可以复用。

6.2 完整调用示例：鸢尾花数据集

下面是一个完整的演示脚本，直接用runExperiment.m跑通整个流程：

%% 数据导入 data = readmatrix('data/iris.csv', 'NumHeaderLines', 1); features = data(:, 1:4); % 标准化z-score mu = mean(features); sigma = std(features); features_std = (features - mu) ./ sigma; %% 聚类运行 rng(42); k = 3; [idx, medoid_idx, cost_history] = myKmedoids(features_std, k, 100, 'euclidean', 'maxdist'); %% 结果展示 fprintf('聚类标签分布：\n'); tabulate(idx); fprintf('Medoid样本行号：'); disp(medoid_idx); %% 可视化 visualizeResult(features_std, idx, medoid_idx, cost_history);

其中visualizeResult函数里画三张图：散点图、代价收敛曲线、轮廓图。这套流程在新数据集上只需要改数据和k值，其他代码不用动。

6.3 后续扩展：k值选择、并行与对比实验

如果你要把这个代码用在正式实验里，可以扩展三个方向：

1. 自适应k值：写一个循环从k=2到k=10跑聚类，每次计算轮廓系数均值，画折线图选拐点。k-medoids的轮廓系数计算可以直接复用silhouette函数，代价不太高。
2. 并行加速：不同k值的聚类任务相互独立，用parfor替代for循环跑多次实验。注意parfor里不能共用随机流，需要每个worker独立设置随机种子。
3. 与k-means对比：在相同数据和相同k值下跑MATLAB内置kmeans函数，比较代价、迭代次数、轮廓系数。这个对比表格放在论文里是很有说服力的实验数据。

我个人在实际操作中最深的一个体会是：k-medoids的代码逻辑不难，难在调试时对“初始化敏感”有心理预期。第一次跑出不太好的结果不要急着怀疑代码，先加rng固定种子、换maxdist初始化、多跑几次取最优。另外，写代码时务必把中文注释的编码问题提前处理掉，否则图里和注释里一堆乱码会浪费你半小时排查时间。最后再分享一个小技巧：画聚类图时，如果两类样本重叠严重，可以在gscatter里加上透明度参数（'Alpha', 0.5），点重叠的地方颜色变深，比默认散点更容易看出边界走势。这套代码和思路我一直在用，希望也能帮你省事。