别再只用K-Means了！用MATLAB手把手教你搞定更抗噪的K-Medoids聚类（附完整代码）
超越K-Means：用MATLAB实战K-Medoids聚类算法解决噪声数据难题
当你的数据集里混入了异常值，K-Means的表现往往会让你失望——那些偏离群体的数据点像磁铁一样把聚类中心拽离合理位置。这时候，K-Medoids算法就该登场了。与K-Means不同，K-Medoids选择实际存在的样本点作为聚类中心（medoids），这种特性让它对噪声和异常值有着天然的抵抗力。
1. 为什么K-Means在噪声数据上会失败？
K-Means算法通过计算簇内数据点的均值来确定新的聚类中心，这个看似合理的策略在面对异常值时却成了致命弱点。想象一下，在一个客户消费数据集中，99%的客户月消费在1000-5000元之间，但有少数几个极端客户消费高达50万元。这些异常值会显著拉高均值，导致聚类中心严重偏离大多数数据点的真实分布。
K-Means的核心问题：
- 对异常值极度敏感
- 聚类中心可能落在没有实际数据点的位置
- 使用欧氏距离，对数据分布形状有特定假设
相比之下，K-Medoids选择实际存在的数据点作为中心，避免了均值带来的偏差。在同一个客户消费数据集中，即使存在极端高消费客户，K-Medoids也会选择一个最能代表群体的实际客户作为中心点，而不是计算可能没有实际意义的平均值。
2. K-Medoids算法原理深度解析
K-Medoids算法的核心思想可以概括为：选择能够最小化簇内不相似度的实际数据点作为代表。这种策略让它特别适合处理以下场景：
- 数据包含噪声和异常值
- 距离度量不是欧氏距离（如曼哈顿距离）
- 需要聚类中心有实际意义（如零售选址问题）
2.1 算法步骤详解
K-Medoids的标准实现（PAM算法）包含以下关键步骤：
- 初始化：随机选择k个数据点作为初始medoids
- 分配阶段：将每个非medoid点分配到最近的medoid所在的簇
- 交换阶段：尝试用非medoid点替换当前medoid，计算总代价变化
- 迭代：重复步骤2-3，直到medoids不再变化或达到最大迭代次数
% PAM算法核心交换阶段伪代码 for each medoid m for each non-medoid point o temp_medoids = medoids; temp_medoids(m) = o; total_cost = calculate_total_cost(temp_medoids); if total_cost < current_cost medoids = temp_medoids; current_cost = total_cost; end end end2.2 代价函数设计
K-Medoids算法的优化目标是最小化总绝对误差（Total Absolute Error）：
TAE = Σ Σ d(p, m_i) i p∈C_i其中d是距离度量（常用曼哈顿距离），m_i是第i个簇的medoid，C_i是第i个簇的所有点。
距离度量选择对比：
| 度量类型 | 公式 | 对异常值敏感度 | 适用场景 |
|---|---|---|---|
| 欧氏距离 | √(Σ(x_i-y_i)²) | 高 | 球形分布数据 |
| 曼哈顿距离 | Σ | x_i-y_i | |
| 余弦相似度 | (x·y)/( | x |
3. MATLAB实战：从数据准备到结果可视化
让我们通过一个完整的MATLAB示例，演示如何处理含噪声数据的聚类问题。
3.1 生成模拟数据
首先创建包含三个簇和随机噪声的测试数据：
% 生成三个清晰簇 cluster1 = randn(100,2) * 0.5 + repmat([2,2],100,1); cluster2 = randn(80,2) * 0.6 + repmat([6,3],80,1); cluster3 = randn(120,2) * 0.4 + repmat([4,7],120,1); % 添加噪声点 noise = [rand(10,2)*10; [8,0.5]; [1,8]; [0.5,0.5]; [9,9]]; % 合并数据 X = [cluster1; cluster2; cluster3; noise]; % 可视化原始数据 figure; scatter(X(:,1), X(:,2), 'filled'); title('含噪声的原始数据分布');3.2 实现K-Medoids算法
下面是完整的K-Medoids算法MATLAB实现：
function [medoids, clusters] = kmedoids(X, k, max_iter) % 输入参数： % X - 数据矩阵(n×d) % k - 簇数量 % max_iter - 最大迭代次数 [n, d] = size(X); % 1. 随机初始化medoids rng(42); % 设置随机种子保证可重复性 medoid_indices = randperm(n, k); medoids = X(medoid_indices, :); % 2. 初始化变量 clusters = zeros(n, 1); prev_medoids = zeros(size(medoids)); iter = 0; % 3. 主循环 while ~isequal(medoids, prev_medoids) && iter < max_iter prev_medoids = medoids; % 3.1 分配点到最近的medoid distances = pdist2(X, medoids, 'cityblock'); % 使用曼哈顿距离 [~, clusters] = min(distances, [], 2); % 3.2 更新每个簇的medoid for i = 1:k cluster_points = X(clusters == i, :); if ~isempty(cluster_points) % 计算所有点作为medoid的代价 intra_distances = pdist2(cluster_points, cluster_points, 'cityblock'); total_cost = sum(intra_distances, 2); % 选择代价最小的点作为新medoid [~, min_idx] = min(total_cost); medoids(i, :) = cluster_points(min_idx, :); end end iter = iter + 1; fprintf('迭代 %d: medoids变化量 = %.4f\n', iter, norm(medoids - prev_medoids)); end end3.3 结果可视化与分析
运行算法并可视化结果：
% 运行K-Medoids k = 3; [medoids, clusters] = kmedoids(X, k, 100); % 可视化聚类结果 figure; colors = ['r', 'g', 'b', 'c', 'm', 'y', 'k']; for i = 1:k scatter(X(clusters==i,1), X(clusters==i,2), 'filled', 'MarkerFaceColor', colors(i)); hold on; end % 标记medoids scatter(medoids(:,1), medoids(:,2), 200, 'k', 'p', 'filled'); title('K-Medoids聚类结果'); legend('簇1', '簇2', '簇3', 'Medoids'); hold off;提示：在实际应用中，建议多次运行算法（不同随机初始化）并选择代价最小的结果，以减少局部最优的影响。
4. 高级技巧与实战建议
4.1 如何确定最佳K值
确定合适的簇数量k是聚类分析中的关键挑战。以下是几种实用方法：
- 肘部法则：绘制不同k值对应的总代价曲线，选择拐点
- 轮廓系数：计算每个点的轮廓系数并取平均值
- Gap统计量：比较实际数据与参考分布的聚类质量差异
% 肘部法则实现示例 k_range = 1:8; costs = zeros(size(k_range)); for i = 1:length(k_range) [~, ~, distances] = kmedoids(X, k_range(i), 100); costs(i) = sum(min(distances, [], 2)); end figure; plot(k_range, costs, '-o'); xlabel('簇数量 k'); ylabel('总代价'); title('肘部法则确定最佳k值');4.2 处理大规模数据的优化策略
标准PAM算法的时间复杂度为O(k(n-k)²)，对于大规模数据效率较低。可以考虑以下优化：
- CLARA算法：在数据子集上应用PAM，选择最佳结果
- 并行计算：利用MATLAB的并行计算工具箱加速距离计算
- 近似算法：如FastPAM等改进版本
算法变体对比：
| 算法 | 时间复杂度 | 适用数据规模 | 精度 |
|---|---|---|---|
| PAM | O(k(n-k)²) | 小规模(<1000) | 高 |
| CLARA | O(ks² + k(n-k)) | 中大规模 | 中 |
| FastPAM | O(n²) | 中等规模 | 高 |
| BanditPAM | O(n log n) | 大规模 | 较高 |
4.3 常见问题排查
当K-Medoids表现不佳时，可以检查以下方面：
- 距离度量选择不当：尝试不同的距离度量（如欧氏距离、曼哈顿距离、余弦距离）
- 异常值过多：考虑先进行异常值检测和过滤
- 数据尺度不一致：对特征进行标准化处理
- 初始medoids选择不佳：尝试k-means++风格的初始化
% 改进的medoids初始化（类似k-means++） function medoids = kmedoids_plusplus(X, k) [n, ~] = size(X); medoids = zeros(k, size(X,2)); % 1. 随机选择第一个medoid medoids(1,:) = X(randi(n),:); for i = 2:k % 2. 计算每个点到最近medoid的距离 distances = pdist2(X, medoids(1:i-1,:), 'cityblock'); min_distances = min(distances, [], 2); % 3. 按距离平方的概率选择下一个medoid probabilities = min_distances.^2 / sum(min_distances.^2); medoids(i,:) = X(find(rand <= cumsum(probabilities), 1), :); end end5. K-Medoids在实际项目中的应用案例
5.1 零售店铺选址分析
在连锁零售店铺选址中，K-Medoids可以帮助识别具有相似特征的区域中心。与K-Means相比，它的优势在于：
- 选择的中心点是实际存在的位置（如现有商圈）
- 对异常区域（如机场、火车站等特殊商圈）不敏感
- 结果更容易解释和落地
% 店铺选址数据示例 load('store_locations.mat'); % 加载经纬度数据 % 使用K-Medoids寻找区域中心 k = 5; % 计划开设5家新店 [medoids, ~] = kmedoids(store_locations, k, 50); % 在地图上可视化 figure; geoscatter(store_locations(:,1), store_locations(:,2), 20, 'b', 'filled'); hold on; geoscatter(medoids(:,1), medoids(:,2), 200, 'r', 'p', 'filled'); title('零售店铺选址分析 - K-Medoids聚类结果'); legend('现有店铺', '推荐新店位置');5.2 医疗图像分析
在医疗图像分割中，K-Medoids可用于识别组织样本中的典型细胞形态：
- 选择最具代表性的细胞作为medoids
- 排除异常细胞（如病变细胞）对聚类中心的干扰
- 结果更符合医生的诊断习惯（基于实际样本）
实际项目经验：在处理病理切片图像时，我们发现K-Medoids选择的代表性细胞比K-Means的虚拟中心更能帮助医生理解算法决策过程。特别是在排除异常细胞影响方面，K-Medoids的表现明显优于K-Means。