X-means聚类算法：用BIC自动选择K值，告别K-means猜K困境

简介X-means.zip是一套基于MATLAB实现的X-means聚类算法源码包面向需要自动确定聚类数量的算法研究者和数据分析人员。传统K-means需要预先指定聚类个数K而X-means借助贝叶斯信息准则BIC在迭代中尝试对每个聚类进行分裂并比较分裂前后的BIC值决定是否接受分裂从而自适应地确定最优聚类数量。压缩包共6个m文件大小仅4KB包含X-means主循环、标准K-means对照实现、BIC计算、两种初始中心选择策略和聚类中心扩展逻辑结构紧凑适合逐行阅读与调试。已有496人学习该压缩包。通过对比X-means与K-means的差异读者能深入理解BIC在聚类数选择中的权衡作用掌握不同初始中心选取方法对聚类结果的影响并可直接将代码迁移到自己的数据集完成更高效、更自动化的聚类任务。1. 为什么 X-means 能救回 K-meansK 值终于不用拍脑袋了做聚类分析的人迟早都会撞到同一个尴尬我用的不是 K-means而是“给 K-means 猜 K”。K 值拍脑袋定跑出来就得靠肘部法则、轮廓系数给结果圆场指标不好看就换一个 K 再跑——循环往复真正的业务问题反而被丢在一边。标题里的 X-means 是我在 MATLAB 里一直在用的改进 K-means 方案它把 K 值选择从人工试探变成算法内循环的一部分用二分分裂加 BIC 模型选择自动收敛到合理的簇数使用者只需要给它一个 maxK 上限。这篇文章用 MATLAB 手写一套可直接运行的 X-means讲清 BIC 怎么算、分裂怎么接、哪里容易翻车以及做完之后怎么验证结果适合正在做聚类分析、又想摆脱“试 K”循环的读者。2. K-means 的三个病灶与 X-means 的改进机制从 2-means 分裂到 BIC 裁决2.1 K-means 的三个老毛病K 未知、初始中心敏感、局部最优K-means 本身不是一个“完整”的聚类算法它只解决“给定 K 之后怎么把点分完”的问题。真实业务里最难的恰恰是 K 本身2000 条用户数据要分群老板问分几群没人知道答案。肘部法则看起来给了答案但换一个人画图“肘”的位置就可能变轮廓系数可以在 2 到 20 之间扫描一遍但这本质上就是把“猜 K”这件事用更耗时的办法又做了一遍。第二个毛病是初始中心敏感。同一个数据集、同一个 K换了随机种子可能收敛到完全不同的局部最优。MATLAB 的 kmeans 默认重复多次取最优缓解了一部分问题但“局部最优”这个结构性问题还在簇与簇之间如果存在重叠K-means 很容易在迭代中把两个真实簇合并成一个之后无论怎么迭代都回不去。第三个毛病是噪声样本会把中心拽离真实位置尤其当某个簇的样本量特别小时一两个离群点就足以让这个簇在下一轮迭代里被别的簇吞掉。这三个问题表面上是三件事本质上是一件事模型选择缺位。K 不知道中心初始化不对局部最优逃不出来都源于没有一个“判断当前聚类结果是否足够好”的准则。X-means 的做法是把“模型好不好”的评估直接放进循环里让算法自己决定要不要继续分。2.2 X-means 的核心机制每个簇尝试 2-means用 BIC 决定是否接受分裂X-meansPelleg Moore, 2000的思路非常工程化从 K2 开始跑一次标准 K-means 收敛然后对每一个簇单独做一次 2-means 分裂尝试用 BIC 比较“该簇保持一个高斯簇”与“该簇分裂成两个高斯子簇”哪个模型更优接受所有 BIC 显示更优的分裂合并所有中心后再跑一次全量 K-means重复直到没有分裂被接受或达到 maxK。关键在 BIC。BIC贝叶斯信息准则在聚类场景里可以写成BIC -2 * ln(L) p * ln(n)其中 L 是当前模型对数据取到的最大似然p 是自由参数个数n 是样本数。对一个簇来说如果假设簇内样本服从各向同性的高斯分布对数似然可以按簇展开ln(L) sum_i [ -n_i * d/2 * ln(2*pi) - n_i/2 * ln(sigma_i^2) - SSE_i / (2 * sigma_i^2) ]SSE_i 是第 i 个簇内所有点到中心的平方距离之和sigma_i^2 是该簇的方差估计。BIC 把似然增益和参数代价放在同一个天平上分裂成两个簇多出来的一个 d 维权值中心和一个方差参数共 d1 个自由参数必须用足够的似然增益来偿还否则 BIC 就会变大分裂被拒绝。这里有一个工程细节值得说X-means 的比较是“局部”的每个簇单独算自己的 BIC而不是每次分裂都重算全局 BIC。如果重算全局每次分裂都要把所有样本的似然重新加一遍复杂度是 O(ndK)循环几十次就非常难看。局部 BIC 只关心当前簇内部开销和这个簇的样本量线性相关代价是贪心决策可能错过全局最优但 X-means 本身就是贪心分裂策略局部比较和这个策略是自洽的。2.3 为什么选 BIC 而不是 AIC、肘部法则或轮廓系数AIC 的惩罚项是 2p不随样本量变化BIC 的惩罚项是 p*ln(n)在样本量稍大时明显更重。聚类数据动辄几千上万条用 AIC 做分裂判据模型几乎总是倾向于继续分最后 K 值一路冲到上限。BIC 的惩罚随样本量增长能在“数据量大”和“模型复杂”之间取得更好的平衡防止过分裂。肘部法则的问题是“肘”的定义依赖人的眼睛。不同人画同一张 SSE-K 曲线有人觉得 K3 是肘有人觉得 K4 才是而且在真实业务数据上曲线往往没有明显的肘只有一个缓坡。X-means 把“要不要分”转化成一个具体的数值比较省掉了这个主观环节。轮廓系数的问题在于计算量。它需要为每个样本计算它与同簇其他样本的平均距离、与最近邻簇所有样本的平均距离复杂度接近 O(n^2·d)。在 X-means 的每次分裂尝试里都跑一遍轮廓系数算法会慢几个数量级而且它优化的是“分离度”而不是模型结构最大化轮廓系数的 K 并不总是业务上想要的 K。轮廓系数适合作为跑完之后的验证指标不适合作为循环内的分裂判据。需要说明边界X-means 的高斯假设只是给 BIC 提供了似然标尺不要求数据真的服从高斯分布。任意形状的簇照样能跑只是 BIC 的统计意义会减弱数据分布极度偏斜时分裂决策可能失真。这一点在调参时心里要有数。3. 在 MATLAB 里从零手写 X-means主循环、BIC 子函数与验证代码3.1 先想清楚分工kmeans 管局部收敛computeBIC 管模型裁决拿到 X-means.zip 这类资源包时我一般不会直接 run 主函数而是先快速扫一眼它内部的“分裂接受准则”——网上流传的很多实现在这里偷懒有的用 SSE 差值代替 BIC有的干脆只要子簇不空就分裂。真正靠得住的实现必须满足两个条件BIC 算得对分裂保护做得到位。代码组织上我习惯拆成两个函数主函数 xmeans 负责整体循环和分裂调度子函数 computeBIC 负责单个簇的 BIC 计算。kmeans 本身由 MATLAB 的 Statistics and Machine Learning Toolbox 提供它做局部收敛比任何手写版本都可靠不值得自己造轮子。3.2 核心代码xmeans 主函数function [bestK, centers, labels, history] xmeans(X, maxK, opts) % X-means: 在 K-means 上自动选择簇数的 MATLAB 实现 % 输入: % X - n x d 样本矩阵 % maxK - 期望的最大簇数 % opts - 可选参数结构体, 支持字段: % .maxIter kmeans 最大迭代次数, 默认 100 % .splitBicGap 分裂需超过的 BIC 增益, 默认 0 % .minPoints 允许分裂的最小簇样本数, 默认 3 % .replicates kmeans 重复次数, 默认 5 % .verbose 打印迭代过程, 默认 false % 输出: % bestK - 最终簇数 % centers - bestK x d 的簇中心 % labels - n x 1 的标签向量 % history - 结构体数组, 记录每次迭代的 K / centers / labels if nargin 2, maxK 20; end if nargin 3, opts struct(); end maxK max(2, min(maxK, size(X,1))); % ---- 读取参数与默认值 ---- if isfield(opts, maxIter) maxIter opts.maxIter; else maxIter 100; end if isfield(opts, splitBicGap) splitBicGap opts.splitBicGap; else splitBicGap 0; end if isfield(opts, minPoints) minPoints opts.minPoints; else minPoints 3; end if isfield(opts, replicates) replicates opts.replicates; else replicates 5; end if isfield(opts, verbose) verbose opts.verbose; else verbose false; end % ---- 从 K2 启动 ---- [labels, centers] kmeans(X, 2, ... MaxIter, maxIter, Replicates, replicates, Start, plus); K 2; history(1) struct(K, K, centers, centers, labels, labels); % ---- 主循环 ---- for iter 1:maxIter improved false; newCenters []; newLabels zeros(size(X,1), 1); for i 1:K idx (labels i); Xi X(idx, :); ni size(Xi, 1); % 防退化: 样本过少的簇不参与分裂 if ni max(minPoints, 3) newCenters [newCenters; centers(i, :)]; newLabels(idx) size(newCenters, 1); continue; end % 簇内跑一次 2-means, 尝试分裂成两个子簇 [subLbl, subC] kmeans(Xi, 2, ... MaxIter, maxIter, Replicates, 1, Start, plus); parentBIC computeBIC(Xi, centers(i, :), ones(ni, 1)); childBIC computeBIC(Xi, subC, subLbl); % 分裂后 BIC 更低说明模型更好, 接受分裂 if parentBIC - childBIC splitBicGap improved true; globalIdx find(idx); for c 1:2 newCenters [newCenters; subC(c, :)]; newLabels(globalIdx(subLbl c)) size(newCenters, 1); end else newCenters [newCenters; centers(i, :)]; newLabels(idx) size(newCenters, 1); end end K size(newCenters, 1); % 用新中心集做一次全量 kmeans, 相当于精修 [labels, centers] kmeans(X, K, ... MaxIter, maxIter, Replicates, 1, Start, newCenters); history(end1) struct(K, K, centers, centers, labels, labels); if verbose fprintf([xmeans] iter %d - K %d\n, iter, K); end if ~improved break; end end bestK K; end这段代码里每一个簇的分裂尝试只做一次 2-means并且 Replicates 设成 1目的是让主循环的时间可控代价是子分裂结果对随机种子更敏感复现结果时需要固定 rng。分裂接受条件写的是 parentBIC - childBIC splitBicGap含义是“子模型的 BIC 比父模型小多少才值得分”默认 0 分不额外要求余量后续调参可以把这个值往正数调。newLabels 那两行的索引映射是典型的出错点subLbl 是 Xi 内部的标签长度 ni不能直接和全局 idx长度 n做逻辑与。先用 find(idx) 拿到这簇样本在原始数据里的全局下标再通过 subLbl c 取子集映射才是安全的。3.3 核心代码computeBIC 子函数function bic computeBIC(Xi, center, labels) % 计算单个簇(或一组样本)在高斯模型下的 BIC % Xi - m x d 样本矩阵 % center - K x d 的中心矩阵, K 是簇数 % labels - m x 1 标签, 取值为 1..K [m, d] size(Xi); K max(labels); if m K bic inf; % 样本不足以估计模型, 禁止分裂 return; end logL 0; for k 1:K Xk Xi(labels k, :); nk size(Xk, 1); if nk 1 bic inf; return; end mu center(k, :); % 当前簇中心 sqDist sum((Xk - mu).^2, 2); % 每个样本到中心的平方距离 varHat sum(sqDist) / nk; % 簇内方差的极大似然估计 varHat max(varHat, eps); % 保护: 防止零方差导致 log(0) logL logL - nk * d / 2 * log(2*pi) ... - nk / 2 * log(varHat) ... - sum(sqDist) / (2 * varHat); end % 自由参数个数: 每个簇一个 d 维权值中心 一个方差 p K * (d 1); bic -2 * logL p * log(m); endBIC 子函数里有三个细节值得注意。第一varHat 用的是极大似然估计分母是 nk 而不是 nk-1因为 BIC 里的模型是生成模型不是做无偏方差估计。第二p K * (d 1) 里的 d 是特征维数这个参数在后续调参会直接影响分裂倾向很多网上实现在这里写错要么漏乘 d要么写成样本数结果 BIC 惩罚项算出来完全不对。第三当某个簇只有一个样本时varHat 为 0log(0) 会算出 -Inf导致整个 BIC 变成 NaN所以我先用 max(varHat, eps) 打补丁再在调用端用 minPoints 做硬性保护——补丁是治标minPoints 才是治本。3.4 用一个三层高斯混合数据验证代码rng(2024); % 固定随机种子, 保证结果可复现 X [randn(300,2) * 0.3 [0, 0]; randn(300,2) * 0.3 [3, 0]; randn(300,2) * 0.3 [1.5, 2.5]]; opts struct(replicates, 3, verbose, true); [K, centers, labels] xmeans(X, 10, opts); figure; gscatter(X(:,1), X(:,2), labels); hold on; plot(centers(:,1), centers(:,2), kx, MarkerSize, 12, LineWidth, 2);这个实验数据刻意造得很“听话”三簇高斯方差一致中心分隔明显K3 是唯一合理的答案。如果 xmeans 输出 K3说明主循环和 BIC 子函数的工作正常如果输出 K2优先检查 BIC 里 p 的写法是否漏了 d如果输出 K4 或更大把 splitBicGap 调到 1 到 5 之间的正数再跑。gscatter 画出来之后可以直观看到中心点是否落在每个簇的几何中心附近——这能顺带验证精修那一步有没有把中心位置拉回正确的地方。4. 让 X-means 更稳的四个改进初始化、距离度量、过分裂保护与并行4.1 初始化改进k-means 与 Replicates 解决“每次结果不一样”MATLAB 的 kmeans 从 R2014a 开始支持 Startplus也就是 k-means 初始化。它比默认的随机采样给出的初始中心更分散能显著降低局部最优的概率。我代码里所有 kmeans 调用都加了 Start, plus包括簇内那次 2-means——子分裂的初始中心如果不好BIC 会基于一个糟糕的局部最优做裁决可能误判“不该分裂”。光有 k-means 还不够。主函数入口和簇内子分裂的初始中心仍然涉及随机性同一份数据跑两次可能得到不同的 K。最省事的做法是在调用端固定随机种子rng(42); [K, centers, labels] xmeans(X, 10, opts);固定种子保证“每次结果一样”但代价是只探索了一条路径。如果项目需要更可靠的结论把主循环里全量 kmeans 的 Replicates 从 1 调到 3 或 5每轮迭代都会用多个初始中心重新优化结果稳定性好很多只是耗时线性上升。我的经验是10 万条以内 Replicates3 是可接受的上限超了优先用固定种子而不是堆重复次数。4.2 距离度量改进余弦距离场景先做 L2 归一化X-means 的 BIC 子函数从欧氏距离出发方差估计和 SSE 都默认数据在欧氏空间里。做文本 TF-IDF、高维稀疏特征这类聚类时业务上更常用余弦距离但如果直接给 kmeans 换 Distance, cosineBIC 里 varHat 的含义就悬空了——余弦距离的“方差”和欧氏距离的“方差”不是同一个量强行混用会让 BIC 的分裂裁决失去统计意义。常见的做法是先对样本做 L2 行归一化再继续用 sqeuclidean。X X ./ vecnorm(X, 2, 2); % 每行除以 L2 范数, 样本落到单位球面归一化之后欧氏距离的平方与余弦距离只差一个常数倍数当 a 和 b 都是单位向量时||a - b||^2 2 - 2*cos(theta)。这意味着按 sqeuclidean 跑出来的簇结构和按余弦距离跑出来的一致而 BIC 里的高斯假设仍然说得通。代价是样本被投影到单位球面上原来的簇形状会被球面曲率稍微扭曲但只要簇不是特别狭长工程上影响很小。马氏距离我不建议在 X-means 循环里直接用它需要先估计全局协方差矩阵并做白化高维数据下协方差估计本身就不稳定白化后特征的含义也会变。真要处理各向异性很强的数据先白化一次再跑标准 X-means比给算法加一个马氏距离选项可靠得多。4.3 过分裂保护minPoints 和 splitBicGap 是防翻车的关键参数X-means 最容易被骂的一点是“分着分着就失控”K 值一路冲到 maxK分出来一堆一两个点的碎片簇。这不是算法本身的问题而是用户没有设置分裂保护。两个参数配合使用参数默认值设小了会怎样设大了会怎样minPoints3单点簇、碎片簇频繁出现小簇被直接忽略真实 K 被低估splitBicGap0轻微过拟合K 偏高K 停在很小的值甚至拒绝任何分裂maxIter100收敛前被截断空转浪费时间一般到不了这个上限replicates5局部最优概率上升时间成倍增加收益递减明显minPoints 我一般取 3 到 5。它保护的是“样本太少的分裂尝试”一个只有三五个点的簇分裂成两个子簇后每个子簇只有一两个点中心位置几乎完全被个别样本牵着走这种分裂即使 BIC 看起来更好也没有任何业务解释力。splitBicGap 是我调 K 值时第一个动的旋钮欠拟合就往负数调允许 BIC 小幅度恶化也接受分裂过拟合就往正数调要求分裂带来足够显著的似然增益。4.4 并行数据量上去之后才值得做的两件事对 kmeans 本身开并行是最低成本方案。MATLAB 的 kmeans 支持用 statset 配置并行多核机器上 Replicates 会被分散到多个 workeroptions statset(UseParallel, true); [labels, centers] kmeans(X, K, Options, options);这个做法只对“大样本 多次重复”的场景有明显收益10 万条以内基本感觉不到变化反而是并行池启动的 overhead 占了上风。另一个更高级的并行思路是在 parfor 里用不同的随机种子跑多个独立的 X-means收集 K 值分布用“多教投票”代替单次运行。parfor b 1:16 rng(b); % 每个 worker 用不同的种子 Klist(b) xmeans(X, 10, opts); end histogram(Klist);这个做法的代价是随机性变得不可精确复现——worker 端的随机数序列不完全由主线程的 rng 控制但用于“看 K 的分布范围是否稳定”已经足够。如果 K 在 3 和 5 之间跳说明数据本身没有清晰的簇数如果 K 每次都严格等于同一个值结构才算稳。5. X-means 避坑排查五个让聚类结果翻车的典型陷阱5.1 BIC 算出 NaN 或 -Inf程序静默停在 K2现象主循环跑完K 永远是 2控制台没有任何报错。原因某个簇样本数只有 1或者簇内所有点到中心的距离都是 0。前者让 varHat 为 0后者让 varHat 在取 max 之前就是 0log(0) 直接算出 -InfBIC 变成 NaN分裂条件永远不成立。解决computeBIC 开头保留 m K 的保护分支全量调用时强制 m size(Xi, 1) 而不是全局样本数varHat 用 max(varHat, eps) 兜底。真正根治还是在主函数里保留 minPoints 检查——靠补丁掩盖退化情况只是让代码“看起来没死”。5.2 maxK10K 一路冲到 10最后出现一两个点的碎片簇现象输出 K 等于 maxK其中一个簇只有 1 到 2 个点其他簇也很碎。原因BIC 惩罚被低估最常见的笔误是 p 只算了 K漏乘 (d 1)或者 splitBicGap 被设成了负数。另一种情况是数据本身噪声很大噪声点在欧氏空间里形成了一些“一对一”的假簇BIC 在样本足够多的局部会认为这些假簇值得保留。解决先把 splitBicGap 调到 2 到 10 之间的正数观察 K 是否回落。再用 minPoints5 兜底小于 5 个点的簇直接不参与分裂。要记住一个现实噪声数据里不存在真实的 K自动聚类算法在噪声上会造出看起来很合理、实际上没有任何可解释性的假象。5.3 同一份数据每次运行 K 值都不同没法向别人汇报现象第一次跑 K3第二次 K5第三次又回到 3。原因kmeans 是局部寻优算法X-means 的每个子分裂都依赖初始中心Replicates1 时随机性没有被平均掉几棵子树的贪心决策累积起来K 值就飘了。解决固定 rng 是第一层补救保证同一份数据在同一机器上可复现。第二层是接受事实并用统计方法量化不确定性Bootstrap 抽样跑 20 到 30 次看 K 的分布汇报时写“K 集中在 3 到 4众数为 3占比 70%”比硬写一个绝对 K 值诚实得多。除非业务明确需要一个整数否则不要为这一点随机性浪费时间。5.4 K 停在 2连已知有 3 类的数据都跑不出 3现象用 3.4 节的三簇高斯数据实验X-means 只给 2 个簇。原因如果数据肉眼可分问题几乎一定出在 BIC 惩罚项过大。高维数据下 p K * (d 1) 随维度线性增长样本量一大p*ln(n) 会吞掉分裂带来的似然增益让算法变得极其保守。另一个隐藏原因调用时把 maxK 传成了 1kmeans 报错后函数退回到默认值 2。解决先确认 maxK 至少大于预期 K。然后把 splitBicGap 改成 -1允许 BIC 轻微变差也接受分裂跑一次看分裂动作是否发生。再把数据降到 2 维做可视化——有些数据在二维投影里看起来可分在高维空间里实际上只有一个主要的分离方向K2 可能才是正确结论这时不要硬调参数。5.5 10 万乘 300 维的数据MATLAB 内存爆炸或慢到让人放弃现象代码在 kmeans 调用处长时间无响应或者直接报内存不足。原因主循环里反复切片 X(idx, :)每次切片都在内存里复制一份子矩阵如果还在 center 更新处写了三层 for 循环每层都对全量样本做距离计算内存和 CPU 都会被打爆。高维数据下 BIC 的 p*ln(n) 惩罚本来就偏大分裂决策也会失真。解决不要手写 kmeans 更新逻辑用 MATLAB 内置版本。大数据先 pca 降维到 30 到 50 维再跑 X-means把完整维度的 kmeans 留给最后精修。10 万条以上可以先用 2 万条子样本跑出 K再在全量数据上用固定 K 做一次 kmeans 分配——这是工程上的取舍不是为了追求“纯算法”。6. 验证 X-means 结果可信度的两个手段evalclusters 与 Bootstrap 稳定性6.1 用 evalclusters 先验证数据本身有多少簇X-means 跑出来的 K 值再漂亮也只是算法自己的判断作为工程师我会先用独立的指标工具对拍一次。MATLAB 的 evalclusters 可以直接对同一份数据跑 silhouette 和 Calinski-Harabasz 指标不需要自己写评价函数eva evalclusters(X, kmeans, silhouette, KList, 1:10); figure; plot(eva);把 evalclusters 给出的最优 K 和 X-means 给出的 K 放在一起看。两者一致时结论基本可靠不一致时优先怀疑 splitBicGap 的方向以及数据在高维空间里是否真的存在可分离结构。注意 evalclusters 的 KList 扫描本身就是一次“猜 K”它只用来交叉验证不替代 X-means 的自动选择。6.2 Bootstrap 稳定性抽样扰动下 K 是否守得住聚类结果的价值在于可重复抽掉 10% 的样本再跑一次结果就变那这个 K 就不值得写进报告里。我的习惯是对样本做 80% 有放回抽样跑 30 次 X-means把 K 值分布画成直方图Klist zeros(30, 1); rng(1); for b 1:30 idxs randsample(size(X,1), round(0.8 * size(X,1))); Klist(b) xmeans(X(idxs, :), 10, opts); end histogram(Klist, Normalization, probability);如果直方图在某个 K 值上的峰高超过 0.6说明数据结构稳定如果分布均匀地铺开说明数据本身没有强结构X-means 只是在噪声上做文章。这一步还能顺带看中心位移——每次运行的最佳簇中心之间的平均位移可以度量聚类的可重复性位移大说明初始化或局部最优问题没被完全压住。6.3 跑完先问三个问题再谈改进我自己在 X-means 上翻过不少车后来养成了一个习惯任何聚类算法跑完先问三个问题——K 值在 Bootstrap 抽样下稳不稳每个簇的簇内方差是不是同一个数量级最外层的语文描述能不能解释给一个不懂算法的人听。三个问题都能回答聚类结果才算数有一个回答不了就先回去调 splitBicGap 和 minPoints。X-means 不是一个黑匣子它把“猜 K”的负担从人移到了 BIC 上但模型选择是否合理最终还是要人来确认。希望这些拆解和踩坑记录帮到你。本文还有配套的精品资源点击获取