K-means聚类收敛真相：局部最优、初始化与MATLAB避坑

简介围绕K-means聚类算法与MATLAB实现整理的一份入门型文档面向正在学习无监督学习、数据分析或需要快速复现聚类实验的高校学生与算法初学者重点解决聚类原理理解、内置函数调用和结果可视化等问题。压缩包共1个docx文件约15KB篇幅精炼适合作为课程笔记或实验参考随时查阅。内容从算法思想、距离度量与迭代流程讲起给出初始化数据、设定K值、调用kmeans函数、获取簇索引与质心、用scatter绘制聚类结果的完整示例并对代码逐段解析同时讨论K值预设、异常值敏感和初始化敏感等局限延伸介绍K-means、Bisecting K-means与Kernel K-means等改进思路。已有444人学习下载读者可据此快速掌握MATLAB聚类的核心调用方式理解参数含义与可视化表达并为后续调参与方法选择打下基础。1. 结论先行K-means 收敛的不是最优解而是离开局最近的不动点一段[idx, C] kmeans(data, 2)这样的 kmeans 聚类算法 matlab 代码敲下去屏幕上立刻给出干净的划分和两个质心看着就像标准答案。可只要把随机种子换一下同样七行数据、同样 K2簇的边界可能整体翻过去SSE 还更小。这个反直觉现象不是 bug它把 K-means 的底色暴露得很彻底算法从不搜索全局最优它只是在一个非凸目标上做交替下降停在哪取决于开局站在哪。搞懂这一点你才会明白为什么Start、Replicates、EmptyAction这些参数不是可选项而是算法本身的一部分。2. 目标函数与 Lloyd 迭代交替两步为什么单调下降还是会卡在局部极小2.1 把「聚类」翻译成一个最小化平方误差的优化问题K-means 的形式化非常朴素。给定 n 个 d 维样本和簇数 K它要找一个划分和一个质心集合让下面这个量最小J Σ(k1..K) Σ(x ∈ C_k) ||x - μ_k||²这就是常说的 SSESum of Squared Errors也叫惯性inertia。注意它把「簇内紧凑」定义成了平方距离和平方这个动作很关键它对大距离惩罚更重所以质心天然被拉向密集区也天然被远处的异常值拽得厉害。后面所有坑几乎都能从这个平方里找到源头。这个目标函数有个尴尬的性质如果允许任意划分它在欧氏空间里的一般情形是 NP-hard 的。也就是说除非 PNP否则没有多项式时间的精确解法。学界和工业界用的都是近似——精确求不行那就求一个「局部最优但够用」的解。理解 K-means本质上就是理解这个近似是怎么做的、误差从哪来。2.2 分配步与更新步两步各自的最优解是怎么推出来的Lloyd 算法的做法是把一个联合优化拆成两个能闭式求解的子问题轮流解。第一步固定质心 μ求最优划分。每个样本点之间独立要让 J 最小每个点只需要选离自己最近的那个质心idx_i argmin_k ||x_i - μ_k||²这一步没有近似是精确最优——因为各项独立、互不干扰。第二步固定划分求最优质心。对某个簇 C_k把 J 里与 μ_k 相关的项拿出来f(μ_k) Σ(x ∈ C_k) ||x - μ_k||²对 μ_k 求梯度令其为零得到Σ(x - μ_k) 0也就是μ_k (1/|C_k|) Σ x。这就是算术平均。二阶导是正定的所以均值是唯一极小点不是巧合是 L2 距离下解出来的唯一答案。两件事都是「给定对方就能算出精确最优」但合在一起就不行了分配是离散选择目标函数关于划分不连续、不凸。这正是交替优化最典型的处境——每一步都在下降但下降的方向被当前状态锁死。用一个手写循环把这两步摆出来比直接调内置函数看得清楚% 手写一版 Lloyd 迭代看清内置 kmeans 在循环什么 rng(7); % 锁种子保证这次讨论可复现 data [1 1; 1.5 2; 3 4; 5 7; 3.5 5; 4.5 5; 3.5 4.5]; K 2; C data(randperm(size(data,1), K), :); % 随机挑 K 个样本当初始中心 maxIter 100; tol 1e-6; prevSSE inf; for iter 1:maxIter % 分配步逐簇算平方距离取最小 D zeros(size(data,1), K); for k 1:K diff data - C(k, :); D(:, k) sum(diff.^2, 2); end [minD, idx] min(D, [], 2); % 更新步簇内取算术平均即平方误差下的最优中心 for k 1:K members data(idx k, :); if isempty(members) continue; % 空簇先跳过后面单独讲 end C(k, :) mean(members, 1); end SSE sum(minD); if abs(prevSSE - SSE) tol % 相对变化够小就停 break; end prevSSE SSE; endrng(7)决定初始中心落点这是整段代码里对结果影响最大的一个数tol是相对判据量级敏感SSE 本身就很大或很小时要相应调整maxIter是保险丝防止在浮点边界来回抖动。min(D, [], 2)里第二个参数2表示沿列方向取最小返回每个点对应的列号也就是簇编号。2.3 复杂度 O(n·K·d·t) 从哪来收敛判据又在测什么每一轮迭代的开销可以拆开算。分配步要对 n 个点、K 个中心、d 个维度两两算距离是 O(n·K·d)更新步只要扫一遍 n 个点求和是 O(n·d)通常被分配步压住。跑 t 轮就是 O(n·K·d·t)。这个式子解释了很多工程现象。n 到百万级、K 到几十、d 到几百时一轮就是几亿次浮点运算t 再取几十单机就吃紧了。所以 MATLAB 的实现在内部做了向量化和分块把||x - μ||²展开成||x||² - 2x·μ ||μ||²一次矩阵乘法搞定而不是真的写双重循环。第 4 章会把这步拆开看。收敛判据才是真正决定「停在哪」的东西。常见做法是看 SSE 的相对变化低于阈值就认为稳定了。要注意这个判据测的不是「离最优有多远」而是「这一轮比上一轮好了多少」。SSE 变化小可能意味着到了好解也可能意味着困在一个烂盆地里的平底上——这两种情况在数值上长得一模一样这就是 K-means 让人心里没底的地方。提示有限个划分状态加上 SSE 单调不增可以推出算法不会无限循环一定在有限步内停下。但「停下」和「停对」是两回事前者有数学保证后者没有。3. 初始化kmeans 的 D² 采样凭什么比随机撒点稳3.1 随机初始化的病根两个真簇被一个种子焊死把初始化想成往数据分布上撒 K 颗种子每一颗最终会长成一个簇。随机撒的问题在于它完全可能把两颗种子撒进同一个真簇里。设想数据其实是两团左团密集、右团稀疏。随机抽两个初始中心有一定概率两个都落在左团。接下来会发生什么左团被这两颗种子对半劈开形成两个簇右团没有任何种子靠近只能被就近归入其中一个左簇。迭代几轮之后两个质心会分别稳定在「左团左半边」和「左团右半边偏右」的位置右团那些点要么被硬塞进某个左簇要么把质心往右拽一点但拽不回来。SSE 明显偏大但它是局部最优算法照样收敛、照样返回结果一句警告都不给。这就是初始化敏感的本质不是算法不稳定是它压根没有跳出局部盆地的能力。Lloyd 迭代只会下降没有任何随机扰动或跳出机制。3.2 D² 采样把「选远一点」写成一个概率分布kmeans 的思路是把随机初始化改造成一种有偏采样让新中心倾向于落在离已有中心远的地方。具体规则是第一个中心在全体样本里均匀随机选此后每次选新中心时先算每个样本到「已选中心集合」的最近距离 D(x)然后按P(x) D(x)² / Σ D(x)²的概率抽取。分母是归一化常数分子里的平方是精髓——它把「远」这件事放大了。一个点到最近中心的距离是 1另一个是 3被选中的概率比不是 3 倍而是 9 倍。这种放大让采样强烈偏向那些还没被覆盖的区域。从效果上看它把「随机撒点」变成了「尽量铺开撒点」同时保留了随机性不是确定性地选最远点那样对异常值极其敏感。理论上能给出的保证是kmeans 的期望 SSE 不超过最优解的 8(ln K 2) 倍。这不算紧但它是第一个与 K 呈对数关系的乘性保证也是它被广泛内置的原因。不过要清醒常数 8 和 ln K 让它在 K 稍大时代价不小实际表现通常远好于这个上界但它仍然只是改善初始点不改变 Lloyd 迭代本身的局部性。3.3 MATLAB 的 Start 与 Replicates参数背后是两套治病的药MATLAB 把初始化的选择权交给了Start参数常见的取值有这么几类取值含义适用场景pluskmeans 采样较新版本的默认值通用数据分布不太规整时首选sample从数据里不放回随机抽 K 个点老版本默认行为快但敏感uniform在数据范围的超矩形内均匀撒点数据分布均匀、边界清晰时可用cluster先做分层抽样再按比例取点簇大小差异大时比纯随机稳数值矩阵直接指定 K×d 的初始中心做对比实验、复现历史结果Replicates是另一套药思路完全不同不改进单次初始化的质量而是跑多次每次用不同初始点最后返回sumd最小的那一组结果。它治的是「随机性」本身代价是时间线性增长。经验上Replicates取 10 到 20 已经能压掉大部分初始化翻车再多收益就很小了。% 两套药一起上kmeans 保证单次质量Replicates 抹平剩余随机性 opts statset(UseParallel, true); % 有并行工具箱时开启 [idx, C, sumd] kmeans(data, 3, ... Start, plus, ... % 或者传 cluster / 数值矩阵 Replicates, 20, ... % 跑 20 次取最优 MaxIter, 200, ... Options, opts);statset里的UseParallel在Replicates大于 1 时才有意义单次运行不会并行。如果环境里没装并行工具箱这个选项会被忽略而不是报错属于安静降级排查性能问题时容易看走眼。注意Replicates只保证在多次结果里挑最好的不保证这个「最好」就是全局最优。它提升的是概率不是确定性。4. 返回值解剖MATLAB 的 kmeans 算完到底交出了什么4.1 idx、C、sumd、D 四个输出各自的语义原文示例只取了idx和centroids其实完整签名能吐出四个东西每一个在排查问题时都用得上idxn×1 的簇编号向量取值 1..K。编号本身没有语义1 号簇不代表任何东西重跑一次编号可能整体换位。CK×d 的质心矩阵第 k 行是第 k 个簇的中心。sumdK×1 向量每个簇内所有点到该簇质心的距离之和。这是判断哪个簇「散」得最厉害的直接指标。Dn×K 矩阵每个点到每个质心的距离。想要做软分配、算置信度、画距离热力图都从这里取。sumd有个容易踩的语义陷阱它的度量单位跟着Distance参数走。默认Distance是sqeuclidean此时sumd是平方距离之和不是距离之和。有人拿sum(SSE)当「平均半径」用算出来的数对不上就是因为少了开方这一步。想看真正的欧氏距离和要么把Distance显式设成euclidean要么自己对sumd开方。% 默认 sqeuclidean 下sumd 是平方量纲 [idx, C, sumd, D] kmeans(data, 2); sse_default sum(sumd); % 平方距离总和 sse_euclid sum(sqrt(sumd)); % 换算成距离总和 fprintf(默认 SSE %.3f距离和 %.3f\n, sse_default, sse_euclid);如果后续要跨不同 K 值比较「聚类紧不紧」必须统一量纲否则 K 越大、簇越小、平方量纲下降得越快比较结果会系统性偏向大 K。4.2 距离度量默认是平方欧氏手写一版距离矩阵对照默认用平方欧氏而不是欧氏原因很实际平方省掉了每次比较前的开方而开方是单调的不改变 argmin 的结果。对聚类本身零影响对性能有正收益。但一旦你要复用D或者sumd量纲就得自己留神。MATLAB 内部算距离不走双重循环走的是矩阵展开% 用矩阵乘法一次算出 n×K 的平方距离矩阵 % 展开式||x - c||^2 ||x||^2 - 2·x·c^T ||c||^2 X2 sum(data.^2, 2); % n×1每个样本的模方 C2 sum(C.^2, 2); % 1×K每个质心的模方 D2 X2 C2 - 2 * (data * C); % n×K靠隐式扩展做加减 D2 max(D2, 0); % 浮点误差可能产生 -1e-16钳到 0 [~, idx] min(D2, [], 2); % 沿列取最小得到簇编号X2 C2那里用的是 MATLAB 的隐式扩展旧称广播n×1 加 1×K 自动铺成 n×K。max(D2, 0)这一行不是多余的当两个点数值上几乎重合时展开式会因为舍入误差算出极小的负数后面若对它开方就会得到复数。这类错误在小数据上基本碰不到上了十万级样本才开始零星出现属于典型的血泪经验。向量化的代价是内存。D2是 n×K 的稠密矩阵n100 万、K50 时就是 5000 万个双精度数约 400 MB还只是一份。真到这个量级就得改成分块计算每次只算几千行或者直接用kmedoids那种支持大数据的实现。4.3 空簇处理与 EmptyActionsingleton 是怎么把簇救回来的空簇是 Lloyd 迭代绕不开的边界情况某个簇在分配步之后一个点都没分到更新步里均值就成了 0/0质心变 NaN下一轮彻底失效。MATLAB 用EmptyAction参数处理常见取值有三种singleton默认把离自己所属质心最远的那个点强行挪过来作为这个空簇的新质心。本质是牺牲一个点换回簇的存活。error直接报错。适合做实验时确认「数据里是不是真有重复点」。drop删掉这个簇idx的取值会少于 K 个。后续代码若假设idx恰好有 K 个不同值就会崩。什么情况下会出空簇最常见的是数据里有大量完全相同的点。比如七个样本里有五个坐标一模一样K 设成 3 时这三个重合点很可能被分到同一簇另外两个簇分不到点。另一个来源是初始中心刚好落在数据范围外导致某个中心「谁也吸不到」。% 有重复点时默认的 singleton 会悄悄牺牲一个点 X [1 1; 1 1; 1 1; 1 1; 5 5; 9 9]; [idx, C] kmeans(X, 3, EmptyAction, singleton); tabulate(idx) % 看一眼三个簇各有多少点是否有一个只剩 1 个tabulate是个很轻的诊断手段比肉眼看散点图快得多。如果发现某个簇只有 1 个点基本可以断定是空簇被 singleton 救回来的这时候要回头查数据里是不是有重复行。提示去重unique(X, rows)之后再聚类可以一次性消掉大部分空簇同时也避免重复点人为放大某片区域的权重。5. 避坑与排查真实数据上最容易翻的五个车5.1 量纲不统一距离被大数量级特征垄断现象把身高厘米几十到两百和年收入元几万到几十万放进同一份数据做聚类结果永远按收入切身高完全不起作用。原因平方欧氏距离对每个维度一视同仁地求和。收入维度的差值动辄上万平方之后是 1e8 量级身高维度的差值也就几十平方后是 1e3 量级。两者差五个数量级身高那点贡献在总和里被彻底淹没等于没参与计算。解决聚类前先标准化。zscore按列减均值除标准差让每个维度的方差都是 1权重才对等。要是业务上认为某维度更重要就别用 zscore改用自定义权重矩阵或者手动缩放。X [height_cm, income_yuan]; range(X) % 先看清楚量纲差多少这一步别省 Z zscore(X); % 按列标准化方差归一到 1 [idx, C] kmeans(Z, 3, Replicates, 20); % 重要C 是标准化空间里的质心要还原回原量纲才能解释 C_orig C .* std(X, 0, 1) mean(X, 1);5.2 K 值拍脑袋肘部法失灵时用轮廓系数兜底现象画 SSE 随 K 变化的曲线想找那个「拐点」结果曲线平滑得像条直线根本没肘部。原因肘部法依赖曲线斜率突变数据分布本身没有明显层次时SSE 会随 K 近似平滑下降拐点纯靠肉眼猜不同人看出不同的 K。解决换成轮廓系数它对每个点的「簇内紧密度 vs 簇间分离度」打分取值 -1 到 1越大越好可以按数值直接比。Ks 2:6; s zeros(size(Ks)); for i 1:numel(Ks) [idx, ~] kmeans(Z, Ks(i), Replicates, 20); s(i) mean(silhouette(Z, idx)); end [best, pos] max(s); fprintf(K%d 时平均轮廓系数最高%.3f\n, Ks(pos), best);silhouette的复杂度是 O(n²)n 大的时候要抽样算。另外轮廓系数也偏好凸的、大小接近的簇数据形状离谱时它同样会骗人多指标交叉验证比单指标靠谱。5.3 异常值把质心拽走均值不是稳健统计量现象数据里混进几个极端值聚类完成后有一整个簇的质心落在荒郊野岭簇里只有那几个异常点加零星正常点。原因质心是均值均值对极端值没有抵抗力。一个坐标偏离一万倍的点和一百个正常点放在一格簇里均值会被它拉走接近一百个点的距离量级。解决先做离群点检测再聚类或者换用对异常值更稳的原型。kmedoids用中位数式的实际样本点当中心配合cityblock距离抗干扰能力明显强于均值。% 均值被拽走时换成基于绝对距离的中位数原型 [idx, C] kmedoids(Z, 3, Distance, cityblock, Replicates, 10);5.4 重复点与空簇NaN 质心从哪冒出来现象跑完聚类C里出现 NaN 行或者sumd里有 NaN。原因数据里有一批完全相同的样本初始化时这些点全被分进同一簇另一个簇分不到点更新步算出 0/0。解决第一去重第二显式指定EmptyAction别依赖默认行为第三跑之前查一下size(unique(X,rows),1)和原始行数的差距差距大就说明重复严重。5.5 图例硬编码与配色语义可视化层的低级错现象K 改成 3 之后图例还是只有Cluster 1和Cluster 2第三个簇在图上没有名字或者做多组对比实验时同一个真实簇在两张图里颜色不一样看图的人理解反了。原因原文示例里legend(Cluster 1, Cluster 2, Centroids)是写死的跟 K 没有任何联动。另外scatter(..., [], idx, filled)用数值索引当颜色MATLAB 会把它线性映射到当前 colormap 上K 一变映射区间就变同一个idx值对应的颜色跟着变。解决图例按 K 动态生成颜色显式指定别让 colormap 自己决定。colors lines(K); % 生成 K 个区分度稳定的颜色 figure; hold on; for k 1:K m idx k; scatter(data(m,1), data(m,2), 40, colors(k,:), filled); end scatter(C(:,1), C(:,2), 150, k, x, LineWidth, 2); labels arrayfun((k) sprintf(Cluster %d, k), 1:K, UniformOutput, false); legend([labels, {Centroids}], Location, best); hold off;arrayfun那行负责按 K 生成标签lines(K)保证颜色在 K 变化时仍然稳定可辨。这类问题不影响算法正确性但会直接毁掉结论的可信度做汇报材料时尤其致命。6. 设计取舍什么时候该把 K-means 整个换掉K-means 隐含了一组挺强的假设簇是凸的、近似等体积的球形、由均值代表合理、K 已知。这四个假设里破任何一个调参都救不回来只能换算法。我一般会先对着下面这张表过一遍再决定要不要在 Lloyd 迭代上继续耗时间。数据特征更合适的选择破的是哪条假设簇大小悬殊小簇被大簇吞掉高斯混合模型 / 谱聚类等体积假设月牙、环形、嵌套结构谱聚类 / 基于密度的方法凸性假设离群点占比不低k-medoids / 基于密度的方法均值可代表性K 完全未知层次聚类 轮廓系数K 已知假设高维稀疏文本、one-hot先降维或改用余弦距离欧氏距离在高维退化数据量到千万级Mini-batch 类实现全量迭代的复杂度余弦距离那一行值得多说一句。高维稀疏向量里两个几乎不相关的样本欧氏距离也会趋近于一个常数区分度被稀释得厉害余弦距离只看方向不看模长反而是这类数据的常规选择。这不是 K-means 本身的问题是「距离」这个前提在高维下失效了。验证上也别只盯着 SSE。同一个 K、同样的数据换三个随机种子跑如果三次的idx差异巨大说明数据本身没有强簇结构再精细的调参也是在拟合噪声。这时候更有价值的动作是回到特征工程而不是继续加Replicates。从那以后我每次跑聚类都强制走一遍「标准化 → 去重 → 多 K 轮廓系数 → 换种子复跑」这四步哪怕数据看着再干净也不例外。有一次偷懒跳过去重结果三个簇里有一个只含一个点差点把结论写进报告。希望帮到你。本文还有配套的精品资源点击获取