尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

MATLAB k-means聚类实战:从k-means++初始化到质量评估

MATLAB k-means聚类实战:从k-means++初始化到质量评估 简介一套MATLAB实现k-means聚类算法的示例代码包面向需要快速上手无监督学习、数据分组或k-means初始化技巧的初学者和科研人员。压缩包共2个文件均为.m脚本整体仅2KB包含一个演示主程序和一个辅助函数前者展示数据加载、调用kmeans函数执行聚类以及结果输出的完整流程后者用于配合计算距离或归类操作便于读者在MATLAB中直接运行并观察每一步的中间结果。包内代码简洁注释清晰适合作为课堂作业、实验报告或项目起步的参考模板。资源已有1071人学习下载运行示例后既能掌握kmeans与init,kmeanspp参数的组合用法也能深入理解聚类中心初始化对最终划分质量的影响。1. 从一次失败的聚类说起初始化决定了 k-means 的上限先抛一个反直觉的结论同样的数据、同样的 k 值k-means 跑两次可能得到完全不同的结果而且第二次的结果未必比第一次好。问题不在迭代次数而在初始聚类中心的选择。传统 k-means 随机选中心一旦初始点扎堆迭代很容易收敛到局部最优簇与簇之间的边界明显偏离真实分布。k-means 的出现就是为了解决这个问题它不是改迭代策略而是改初始化策略——让初始中心尽可能分散从源头压低局部最优的概率。这份 k-means.zip 压缩包恰好把这两条路线都覆盖了k-means demo.m走的是 MATLAB 内置kmeans函数的完整演示流程zixie.m是手写实现适合想搞懂每个细节的人。后续所有代码和参数说明都基于 R2023b 版本验证旧版本差异不大个别参数名需要对照文档确认。如果你正在做数据聚类、图像分割或特征分组这篇文章能让你既会用函数也能在结果不理想时定位到具体环节。2. 算法核心与 MATLAB 内置 kmeans 函数的参数边界2.1 从 Lloyd 算法到 k-means 的完整迭代逻辑k-means 的迭代过程可以拆成四步初始化 k 个中心、计算每个样本到所有中心的距离、把样本分配到最近的中心、重新计算每个簇的均值作为新中心。重复后两步直到中心不再变化或达到最大迭代次数。这个流程本质上是 Lloyd 算法MATLAB 内置函数底层用的就是它。关键点在于距离度量。默认是欧氏距离适合凸形簇但离群点对均值影响极大。实际项目中我一般先做标准化否则量纲大的特征会主导距离计算导致聚类结果偏向某些维度。% 手写核心迭代逻辑对应 zixie.m 的简化版 function [idx, C] my_kmeans(X, k, max_iters) [n, d] size(X); % 随机选 k 个样本作为初始中心 rng(42); init_idx randperm(n, k); C X(init_idx, :); for iter 1:max_iters % 计算每个样本到所有中心的欧氏距离 dist pdist2(X, C, euclidean); % 分配样本到最近中心 [~, idx] min(dist, [], 2); % 更新中心 for j 1:k C(j, :) mean(X(idx j, :), 1); end end endpdist2是 MATLAB 的距离计算函数euclidean指定欧氏距离也可以换成cityblock或cosine。min(dist, [], 2)表示按行取最小值返回的idx就是每个样本的簇编号。更新中心时用mean对每个簇内样本取均值这就是 Lloyd 算法中的「质心更新」。这个实现没有处理空簇的情况如果某个簇在迭代中没有任何样本mean会得到NaN后续计算直接崩掉。2.2 kmeans 函数的核心参数与选型建议MATLAB 内置kmeans远比我上面手写的版本健壮它默认处理了空簇、收敛判断和并行计算。最常用的调用方式是[idx, C] kmeans(X, k)但仅用默认参数很难发挥这个函数的全部能力。参数可选值默认值适用场景initplus即 k-means、sample、kmeans或矩阵plus默认plus就是 k-means无需手动指定MaxIter正整数100数据量大或簇形状复杂时调大比如 500Replicates正整数1多次运行取最优避免局部最优一般设 5-10Displayoff、final、iteroff调试时看迭代过程用finalDistancesqeuclidean、cityblock、cosinesqeuclidean高维文本数据用cosine有离群点用cityblock% 实际项目中的推荐配置 [idx, C, sumd] kmeans(X, 5, ... init, plus, ... % k-means 初始化 MaxIter, 500, ... % 最大迭代次数 Replicates, 10, ... % 重复 10 次取最优 Distance, sqeuclidean); % 平方欧氏距离Replicates参数不是简单跑 10 次取平均而是每次用不同的随机初始中心最终返回目标函数组内平方和最小的那次结果。sumd返回每个簇内样本到中心的距离之和这个值可以用来做后续的质量评估。提示init参数如果传一个 k×d 的矩阵可以手动指定初始中心这在某些业务场景下很有用——比如按已知的分群人工指定初始点然后让算法微调。2.3 内置函数与手写实现的差异在哪里zixie.m的价值不在性能而在可读性。它把 k-means 的每个环节暴露出来方便你在中间加调试代码、输出每次迭代的可视化结果。但它缺少几个关键机制空簇处理、收敛阈值判断、多次运行取最优。内置kmeans默认的收敛阈值是1e-4即中心变化小于该值就停止迭代而手写版通常只靠max_iters控制。另一个差异在距离计算方式。内置函数在sqeuclidean下直接算平方距离省去开方运算性能更好。手写版用pdist2后还要再平方虽然结果等价但多余的开方在百万级数据点上会明显拖慢速度。如果你要处理大规模数据建议直接用内置函数如果只是学习原理zixie.m更合适。3. k-means 初始化概率分布如何化解局部最优困局3.1 为什么随机初始化会失败随机初始化选中心时每个点被选中的概率是相等的。如果数据分布不均匀比如三个簇大小差异明显随机选到的中心很可能全部落在最大的簇里。迭代开始后这些挤在一起的初始中心会把大簇切成几块而真正的小簇完全没有竞争到中心最终收敛到一个看起来合理但违背真实结构的局部最优。k-means 的核心思想是第一个中心随机选之后每个新中心优先选在离已有中心远的位置。这样做的好处是初始中心覆盖整个数据空间迭代起点更接近全局最优附近。MATLAB 从 2015a 开始plus作为默认初始化方法底层实现就是 k-means 的变体。3.2 概率加权采样的完整 MATLAB 实现k-means 的第二步是关键——每个样本点被选为中心的概率与该点到最近已有中心的距离平方成正比。这种加权抽样保证了远的点有更大概率被选中但又保留一定随机性避免每次都选到离群点。% k-means 初始化的手写实现 function C kmeans_pp_init(X, k) [n, ~] size(X); C zeros(k, size(X, 2)); % 第一步随机选第一个中心 rng(42); C(1, :) X(randi(n), :); for t 2:k % 计算每个点到最近中心的距离平方 dist pdist2(X, C(1:t-1, :), euclidean); min_dist min(dist, [], 2); prob min_dist.^2 / sum(min_dist.^2); % 按概率加权随机选择下一个中心 cum_prob cumsum(prob); r rand(); idx find(cum_prob r, 1); C(t, :) X(idx, :); end end这里的核心逻辑在prob min_dist.^2 / sum(min_dist.^2)先取每个样本到最近中心的距离平方后再归一化成概率分布。cumsum生成累积概率find(cum_prob r, 1)找到随机数r落在哪个区间对应选哪个点。这种写法比直接调用datasample更透明也方便调试。注意距离平方可能导致概率分布极度不均匀某个离群点会被反复选中影响稳定性。常见的做法是在平方前加上很小的平滑项或者直接用距离而非距离平方。MATLAB 内置实现用的是距离平方但不加平滑——这意味着如果数据里有极端离群点k-means 可能选到它作为中心但这个中心在后续迭代中会快速被拉回正常区域。3.3 用可视化验证 k-means 的效果只看代码不够直观我一般会在二维数据上把初始中心和最终中心都画出来对比。% 生成三个高斯簇 rng(1); X [randn(100,2)*0.5 [2,2]; randn(100,2)*0.5 [-2,2]; randn(100,2)*0.5 [0,-2]]; k 3; % 分别用随机初始化和 k-means 初始化 init_random X(randperm(size(X,1), k), :); init_plus kmeans_pp_init(X, k); % 画初始中心分布 figure; scatter(X(:,1), X(:,2), 10, filled); hold on; scatter(init_random(:,1), init_random(:,2), 100, r, x); scatter(init_plus(:,1), init_plus(:,2), 100, g, o); legend(数据点, 随机初始化, k-means初始化);通常你会看到随机初始化的三个中心可能都落在同一片区域而 k-means 的中心会分散开分别靠近三个簇的质心附近。这个可视化能直观解释为什么后续迭代中 k-means 的收敛速度更快——起点更接近最终解需要的迭代次数自然更少。4. 聚类质量评估肘部法则、轮廓系数与 k 值选择4.1 困惑点k 值不是算法算出来的k-means 需要用户预先指定 k这个参数不属于算法自动优化范畴。常见的误区是直接用最大轮廓系数选 k但这并不适合所有场景——轮廓系数偏向紧凑的球形簇对拉长的簇分布会给出偏低分数。工程上更稳妥的做法是把业务需求作为第一约束肘部法则和轮廓系数作为辅助判断。肘部法则的原理是随着 k 增大组内平方和within-cluster sum of squares, WSS单调下降但下降速度在某个 k 值后明显放缓这个拐点就是合理的 k。轮廓系数则衡量每个样本与自己簇内样本的相似度高于与其他簇样本的程度取值范围是 [-1, 1]越接近 1 越好。4.2 MATLAB 实现同步计算 WSS 和轮廓系数% 用内置 kmeans 函数测试多个 k 值 X load(your_data.mat); % 替换成你的数据 X X.data; k_list 2:8; wss zeros(length(k_list), 1); sil_scores zeros(length(k_list), 1); for i 1:length(k_list) k k_list(i); % 每个 k 值跑 5 次取最优 [idx, ~, sumd] kmeans(X, k, Replicates, 5, MaxIter, 500); wss(i) sum(sumd); % 总 WSS % 随机抽样 1000 个样本计算轮廓系数加速计算 s silhouette(X, idx); sil_scores(i) mean(s); end % 绘制双轴图对比 figure; yyaxis left; plot(k_list, wss, -o); ylabel(WSS); yyaxis right; plot(k_list, sil_scores, -s); ylabel(轮廓系数); xlabel(k);silhouette函数直接接收数据矩阵和聚类索引返回每个样本的轮廓系数。数据量大时全量计算很慢我的做法是用randsample抽取一部分样本算轮廓系数得到的均值作为整体近似。WSS 和轮廓系数同时画在一张图上能更清晰看到拐点和峰值是否一致。提示肘部法则的「拐点」有时并不明显尤其是数据本身没有天然簇结构时WSS 曲线可能是一条平滑下降的曲线。这时别硬找拐点结合业务判断或者改用更严格的稳定性检验。4.3 Gap Statistic比肘部法则更可靠的补充方法肘部法则的缺陷在于它没有比较基准——WSS 下降多少算「明显放缓」没有客观标准。Gap Statistic 通过比较真实数据的 WSS 与均匀分布数据的 WSS 来定义 gap 值取 gap 最大值对应的 k。MATLAB 没有内置实现需要手写但逻辑不复杂。% Gap Statistic 的核心逻辑简化版 function [best_k, gap] gap_statistic(X, max_k, B) [n, ~] size(X); wss_obs zeros(1, max_k); wss_ref zeros(max_k, B); for k 1:max_k % 计算真实数据的 WSS [~, ~, sumd] kmeans(X, k, Replicates, 3); wss_obs(k) sum(sumd); % 生成 B 次均匀分布参考数据 for b 1:B X_ref unifrnd(min(X), max(X), n, size(X, 2)); [~, ~, sumd_ref] kmeans(X_ref, k, Replicates, 1); wss_ref(k, b) sum(sumd_ref); end end % gap 参考数据的期望 WSS - 真实 WSS l mean(log(wss_ref), 2); gap l - log(wss_obs); [~, best_k] max(gap); end注意这里对 WSS 取了对数原因是参考数据的 WSS 波动较大对数变换可以稳定方差。参考数据用unifrnd在真实数据的取值范围内生成均匀分布——这种做法对多维数据效果一般因为均匀分布和真实分布可能差异很大。更严谨的替代方案是 PCA 白化后在主成分空间内均匀采样但实现复杂度会高很多。Gap Statistic 的计算量是 B 倍于普通 k-meansB 一般取 20 左右。如果数据量大这个方法的耗时可能让人难以接受我的建议是中小数据集用 Gap Statistic 做交叉验证大数据集直接看肘部曲线加业务判断。5. 两个隐藏技巧空簇兜底与高维数据的降维预处理5.1 空簇问题的正确处理方式k-means 迭代过程中可能出现某个簇没有分配到任何样本的情况。内置kmeans默认会终止并把该簇标记为空但这不是我们要的结果。常见做法是重新初始化该簇中心或者把离其他中心最远的样本点设为该簇中心。后者更推荐因为它不会破坏已有结构。% 空簇重初始化策略 function C handle_empty_clusters(X, idx, C, k) for j 1:k if sum(idx j) 0 % 找到离所有现有中心最远的点 dist pdist2(X, C, euclidean); min_dist min(dist, [], 2); [~, far_idx] max(min_dist); C(j, :) X(far_idx, :); end end end这段逻辑把空簇的中心重置到离当前所有中心最远的样本点上这样下一轮迭代大概率会有样本被分配过来。为什么不直接随机选因为随机选可能又选到已有中心附近导致空簇反复出现。最远点策略在增量式聚类中也常用思路一致。5.2 高维数据先降维还是先聚类高维数据对 k-means 最致命的冲击是维度灾难——高维空间中所有点对之间的距离趋同距离度量区分度急剧下降。我的经验是先用 PCA 降到 10-50 维观察可解释方差如果前两个主成分能解释超过 60% 的方差直接可视化后用二维或三维的 k-means 就足够了。% PCA 降维后聚类 [coeff, score, ~, ~, explained] pca(X); cum_var cumsum(explained); % 选择保留 90% 方差的维度 d find(cum_var 90, 1); X_pca score(:, 1:d); [idx, C_pca] kmeans(X_pca, k, Replicates, 5);pca返回的explained是每个主成分的方差百分比cumsum累加后找到第一个超过 90% 的位置作为保留维度。需要注意的是PCA 是线性降维如果数据分布呈流形结构应该考虑 t-SNE 或 UMAP——但 t-SNE 的结果受困惑度等超参数影响很大且每次运行结果不稳定用于聚类前的预处理反而可能引入噪声。一个容易踩的坑是降维后再聚类聚类的簇边界已经不在原始特征空间里了无法直接把簇解释到原始特征上。如果业务上需要知道每个簇的「特征画像」我一般会在聚类完成后对每个簇单独计算原始特征的中位数而不是均值——中位数对离群点更鲁棒得到的画像更接近典型样本。5.3 轮廓系数异常时的排查路径如果某个 k 值下轮廓系数极低甚至有大量负值不要急着换 k先检查数据预处理是否正确。最常见的问题是标准化漏了某个分类变量或稀疏特征导致距离被少数维度主导。其次是存在明显的离群点轮廓系数对离群点极其敏感一个离群点就可能拉低整个簇的均值建议在聚类前用 z-score 检测并处理。当数据中存在明显非凸的簇形状时k-means 天然无法正确分簇此时无论怎么调 k 和初始化都无济于事。与其硬调参数不如换个算法——DBSCAN 基于密度而非距离均值对任意形状的簇都有效但它的参数eps和minPts对结果影响极大需要事先估计密度。谱聚类能处理任意形状但复杂度高适合中小规模数据。k-means 的场景是凸形簇、中等规模、快速迭代理解这个边界比记住更多参数更有价值。本文还有配套的精品资源点击获取
返回列表