
简介面向机器学习初学者、课程实验与科研复现人员这份MATLAB聚类算法代码合集围绕无监督分类任务整理FCM模糊C均值、K-means、Kmedia等经典方法的实现与调用示例并附带FuzzyClusteringToolbox相关函数便于理解隶属度更新、聚类中心迭代与收敛判断。压缩包共74个文件、约2.05MB以40个m脚本为核心另有17个gif、3个png展示聚类过程与界面素材6个txt、2个htm及2个css补充说明1个mat、1个pdf和1个wk1可用于数据实验与文档查阅。已有4888人学习下载。读者可获得可直接运行的聚类主函数、数据标准化与反标准化脚本、聚类有效性评估代码及PCA、投影等示例配合示例数据完成算法对比、结果可视化和参数调试进而辨析不同算法在模糊边界、凸形分布和非球形数据上的适用差异也方便快速搭建课程作业或项目实验流程。1. 四种聚类算法摆在一起为什么“多出来的三个算法”不是备胎而是校准你直觉的标尺跑过聚类的人大多遇到过这样的场景同一份数据Kmeans 跑出来是三个簇换成 FCM 就裂成五个再用神经网络里的自组织映射跑一遍边界又是另一副样子。多数人的第一反应是“算法不稳定”但真正反直觉的结论是——聚类结果的差异很少来自算法本身的随机性而是每种算法对“簇应该长什么形状”内置了完全不同的几何假设。在 MATLAB 里把 Kmeans、FCM、层次聚类和 SOM 神经网络这四种聚类算法放进同一套示例代码最大的价值不是多给三个备选方案而是逼着你在动手前先回答一句我的数据是球状的、边界重叠的还是需要保留拓扑关系的把这个问题想清楚后面四份源码才读得下去。2. Kmeans 的硬划分为什么 Lloyd 迭代一定收敛却常常停在错误的局部最优Kmeans 的整个过程可以压缩成一句话在“固定中心算归属”和“固定归属算中心”之间反复横跳。这个朴素到近乎粗暴的交替恰恰是它能在几十行内跑完、又能在大部分球状数据上给出稳定结果的原因。但要真正读懂这份源码得先明白它在优化一个什么样的目标以及这个目标天然带着什么缺陷。2.1 目标函数与两步交替每次迭代到底在下降什么Kmeans 优化的目标函数是簇内平方误差和记作 JJ Σ_{k1..K} Σ_{x∈C_k} ||x - μ_k||²其中 C_k 是第 k 个簇的样本集合μ_k 是该簇的中心。算法的迭代过程只有两步分配步固定所有 μ_k把每个样本 x 归到距离最近的中心等价于让该项 ||x - μ_k||² 最小。更新步固定所有归属把每个 μ_k 重算成簇内样本的均值等价于让该簇的平方误差和对 μ_k 求导为零。两步各自都让 J 单调不增而 J 有下界 0所以迭代必然收敛。这里的关键细节是收敛的是“不再下降”不是“降到全局最小”。分配步本质上是一个离散优化J 关于归属是非凸的所以 Lloyd 迭代只保证落在一个局部极小值上。换句话说Kmeans 给出的是一个“自洽”的解不是“最优”的解。读源码时要特别留意初始化部分——初始化决定了你最终落在哪个局部极小值里这也是后面 k-means 存在的原因。2.2 距离度量就是几何假设欧氏距离下Kmeans 只能切出凸簇很多人把 Kmeans 当成“万能聚类”直到在环形数据或月牙形数据上遇到问题。根本原因不在实现而在欧氏距离本身在欧氏距离下每个中心的“势力范围”是一个 Voronoi 胞腔而 Voronoi 胞腔一定是凸的。这意味着 Kmeans 的决策边界是若干条直线高维下是超平面的组合它天生只能切出凸的、接近球状的簇。一旦数据的真实簇是环形、螺旋形、或者两个簇被一条细长的“桥”连在一起Kmeans 的切法就会把同一个环切成两半或者把桥两端的点强行归到同一侧。这不是参数没调好是几何假设和数据结构根本不匹配。示例代码里如果提供了“月牙形测试数据”的生成脚本读的时候要顺手把 Kmeans 和后面 FCM 的结果叠在一张图上对比视觉上的差异比看数字直观得多。提示判断数据适不适合 Kmeans最省事的办法是先做一次 PCA 降到二维画散点图。如果簇看起来是一团一团的Kmeans 大概率没问题如果是弯月、环形、嵌套结构直接跳到层次聚类或 SOM 那两节。2.3 MATLAB 源码拆解初始化、空簇处理与收敛判据的工程细节示例代码里 Kmeans 的核心实现通常是这个结构下面这段是带 k-means 初始化和空簇处理的版本function [labels, centers, J_hist] kmeans_core(X, K, maxIter, tol) % X: N x D 样本矩阵N 个样本D 维特征 % K: 聚类数maxIter: 最大迭代次数tol: 中心点移动的收敛阈值 [N, D] size(X); % 1. k-means 初始化第一个中心随机后续按距离平方概率采样 centers zeros(K, D); centers(1, :) X(randi(N), :); for k 2:K d2 min(pdist2(X, centers(1:k-1, :)).^2, [], 2); prob d2 / sum(d2); centers(k, :) X(find(rand cumsum(prob), 1), :); end J_hist zeros(maxIter, 1); for iter 1:maxIter % 2. 分配步每个样本归到最近的中心 dist pdist2(X, centers); [~, labels] min(dist, [], 2); % 3. 更新步按当前归属重算均值 newCenters centers; for k 1:K idx labels k; if any(idx) newCenters(k, :) mean(X(idx, :), 1); else % 空簇处理把离所有中心最远的样本拎出来当新中心 [~, farIdx] max(min(pdist2(X, centers), [], 2)); newCenters(k, :) X(farIdx, :); end end % 4. 记录目标函数并判断收敛 J_hist(iter) sum(min(pdist2(X, newCenters).^2, [], 2)); if norm(newCenters - centers, fro) tol centers newCenters; break; end centers newCenters; end J_hist J_hist(1:iter); end几个参数和逻辑值得单独拎出来说。tol控制的是中心点矩阵的 Frobenius 范数变化量而不是目标函数的变化量这样做的好处是阈值含义稳定——不管数据量纲怎么变中心移动小于tol就停。空簇处理是很多教学代码会省掉的部分但真实数据里一旦某个中心“抢不到”样本不处理就会产生 NaN整条链路直接断掉。把最远点拉过来做新中心是一个常见但不算最优的补救办法示例代码里如果用的是随机重选也要能看懂两种做法的差别。J_hist留了每次迭代的目标函数值后面做收敛曲线、判断是否震荡时都用得上建议不要删。3. FCM 模糊 C 均值隶属度矩阵为什么能缓解“边界样本”问题Kmeans 最让人不放心的部分是它对边界样本的态度一个点离两个中心差不多远它也会被硬塞进其中一个而且这个归属在下一次迭代里会被当成既成事实继续放大。FCM模糊 C 均值换了一种思路——不给非黑即白的归属而是给每个样本对每个簇一个 0 到 1 之间的隶属度让“不确定”本身成为算法的一部分。关键词里把 FCM 和 Kmeans 并列正是因为它们在结构上同源、在边界处理上互为对照一起读最能看出设计动机。3.1 隶属度矩阵把“非此即彼”变成“灰度可调”FCM 的目标函数在 Kmeans 的基础上改了一处平方误差不再按归属累加而是按隶属度的 m 次方加权J Σ_{i1..N} Σ_{j1..C} u_ij^m · ||x_i - c_j||² 约束对每个样本 iΣ_j u_ij 1且 u_ij ∈ [0, 1]这里 u_ij 是第 i 个样本对第 j 个簇的隶属度c_j 是簇中心m 是模糊指数。约束“每行隶属度之和为 1”是关键它保证了一个样本的“归属总量”守恒不会出现所有隶属度一起趋近 0 的退化情况。Kmeans 可以看作 FCM 在 m→1 时的极限隶属度退化成 0/1 的硬分配。所以读 FCM 源码时可以先把它当成“Kmeans 加了权重”的版本去理解再把权重和约束单独拿出来看。3.2 m 参数到底在控制什么为什么经验值取 2m 是模糊指数它决定了隶属度分布的“陡峭程度”。m 越接近 1隶属度越接近 0/1FCM 退化为 Kmeansm 越大隶属度分布越平缓一个样本对多个簇都给出中等偏高的隶属度簇的边界越模糊。那为什么绝大多数示例代码默认 m2原因不是理论推导出来的最优值而是实践经验m2 时隶属度的区分度和抗噪性之间平衡得比较好。m 取 1.1 会让隶属度几乎退化成硬划分FCM 相对 Kmeans 的优势消失m 取 4 以上又会让所有样本对所有簇的隶属度都挤在 1/C 附近去模糊化之后簇的边界几乎随机。所以这个参数与其说是“调出来的”不如说是“不能乱动”的经验值——源码里如果把它写成可配参数注释里最好标明常用区间是 [1.5, 2.5]。注意m 参数和样本量纲是耦合的。如果数据没有归一化距离平方的量级很大m 的作用会被放大表现为隶属度整体偏平。归一化之后再调 m才是有效的调参顺序。3.3 MATLAB 实现去模糊化与收敛判据FCM 的迭代公式可以直接从拉格朗日乘子法推出来落到代码里是两个交替步骤function [U, centers, objHist] fcm_core(X, C, m, maxIter, tol) % X: N x DC: 聚类数m: 模糊指数常用 2 [N, D] size(X); % 1. 隶属度矩阵初始化每行归一化避免某列全零 U rand(N, C); U U ./ sum(U, 2); objHist zeros(maxIter, 1); for iter 1:maxIter U_old U; % 2. 更新聚类中心隶属度的 m 次方作为权重 Um U .^ m; centers (Um * X) ./ sum(Um, 1); % 3. 更新隶属度距离越小隶属度越高 dist pdist2(X, centers); dist max(dist, eps); % 防止除零 ratio dist .^ (-2 / (m - 1)); U ratio ./ sum(ratio, 2); % 4. 目标函数与收敛判断 objHist(iter) sum(sum((U .^ m) .* (dist .^ 2))); if norm(U - U_old, fro) tol break; end end objHist objHist(1:iter); % 5. 去模糊化取每行最大隶属度作为硬标签 [~, labels] max(U, [], 2); end这里有两个容易看漏的细节。第一收敛判据用的是隶属度矩阵的 Frobenius 范数变化量而不是目标函数值——因为 FCM 的目标函数在迭代后期变化非常小用目标函数做判据容易过早停止而隶属度矩阵的变化能更敏感地反映状态是否稳定。第二dist max(dist, eps)这一行不是可有可无的防御性代码当某个样本恰好和某个中心重合时距离为 00^(-2/(m-1))在 MATLAB 里会得到 Inf整行隶属度都会坏掉。示例代码里如果做了这个保护说明作者至少在这个点上实际跑过。最后的max(U, [], 2)就是去模糊化把软隶属度转成硬标签方便和 Kmeans 的结果做对比评估。4. 层次聚类与 SOM 神经网络一个用树状结构看全局一个用拓扑保持做降维前两种算法都在“先定簇数、再迭代”的框架里打转层次聚类和 SOM 跳出了这个框架。层次聚类不预设 K而是先建一棵完整的合并树再在树上切一刀SOM 则借了神经网络的结构用竞争学习和邻域更新让高维数据在低维网格上保留拓扑关系。这两种算法和前面两种放在一起覆盖的正是“簇数未知”和“需要保留邻近关系”这两类真实场景。4.1 层次聚类链接准则的选择决定了簇的形状凝聚层次聚类的流程很直接把每个样本当成一个簇反复合并距离最近的两个簇直到只剩一个簇整个过程记录成一棵 dendrogram。真正影响结果的是“簇与簇之间的距离怎么定义”也就是链接准则。单链接取两个簇中最近点对的距离全链接取最远点对平均链接取所有点对的平均。单链接最大的问题是链式效应只要两个簇之间有一串距离都很近的点它就会把它们串成一个细长的簇哪怕这个簇在几何上根本不应该存在。平均链接和 Ward 方法在实际数据上更稳示例代码里如果提供了可切换的linkageType参数建议至少把单链接和平均链接各跑一遍亲眼看一次链式效应长什么样。function labels hierarchical_core(X, K, linkageType) % linkageType: single | complete | average | ward Y pdist(X); % 成对距离 Z linkage(Y, linkageType); % 建合并树 labels cluster(Z, maxclust, K); % 在树上切出 K 个簇 end这段代码短到几乎不需要解释但cluster(Z, maxclust, K)的切树逻辑值得注意它切的是“让簇数不超过 K 的最大不一致系数”位置而不是简单地在第 N-K 次合并处切两者在噪声数据上会得到不同结果。读源码时如果发现示例里换成了cophenet做树质量评估那是更严谨的做法。4.2 SOM 神经网络竞争学习与邻域衰减如何保留拓扑SOM自组织映射是这四种算法里唯一带“神经网络”标签的。它的结构是一张低维网格通常二维网格上每个节点挂一个和输入同维的权重向量。训练时对每个输入样本做两件事先找权重向量离它最近的节点获胜节点然后让获胜节点和它邻域内的节点一起向这个样本靠拢。关键在于邻域半径和学习率都随训练进度衰减——早期大范围调整让网格整体铺开后期小范围微调让局部收敛。这种“先粗后细”的调度正是它能保留拓扑关系的原因输入空间里相近的样本会激活网格上相近的节点。function [W, bmu] som_core(X, gridSize, maxEpoch, eta0, sigma0) % gridSize: 网格边长如 5 表示 5x5eta0: 初始学习率sigma0: 初始邻域半径 [N, D] size(X); W rand(gridSize, gridSize, D); [gx, gy] meshgrid(1:gridSize, 1:gridSize); totalIter N * maxEpoch; t 0; for epoch 1:maxEpoch for i randperm(N) t t 1; % 1. 找获胜节点 diff W - reshape(X(i, :), 1, 1, D); d2 sum(diff .^ 2, 3); [~, winIdx] min(d2(:)); [wx, wy] ind2sub([gridSize, gridSize], winIdx); % 2. 邻域函数与学习率随时间指数衰减 sigma sigma0 * exp(-t / totalIter); eta eta0 * exp(-t / totalIter); h exp(-((gx - wx).^2 (gy - wy).^2) / (2 * sigma^2)); % 3. 获胜节点及其邻域向输入靠拢 W W eta * h .* (reshape(X(i, :), 1, 1, D) - W); end end % 4. 用最终权重给样本打标签 bmu zeros(N, 1); for i 1:N diff W - reshape(X(i, :), 1, 1, D); d2 sum(diff .^ 2, 3); [~, bmu(i)] min(d2(:)); end endsigma0一般取网格边长的一半左右eta0取 0.5 上下两者都按exp(-t/totalIter)衰减保证前期粗调、后期细调。如果示例代码里把衰减函数写成线性而不是指数收敛轨迹会更“硬”但对小数据集够用。最后一步的bmu是每个样本的获胜节点编号可以拿它当簇标签用也可以进一步对网格节点做层次聚类得到更少的簇。4.3 四个算法的统一接口与对比示例把四种算法放在同一个脚本里调用时接口统一比算法本身更影响可读性。示例代码一般会约定输入都是N×D的X和簇数K层次聚类和 SOM 可以传也可以不传输出都是长度为 N 的labels外加一个可选的簇中心或权重矩阵。这种约定让“换算法跑同一份数据、用同一个指标评估”变成一行代码的事也是后面验证章节能顺利展开的前提。5. 排查与避坑MATLAB 里跑聚类最容易翻车的五个细节聚类代码跑通不难跑出可信结果不容易。下面这五条是实际调试时最容易撞上的每一条都按“现象 → 原因 → 解决”拆开写方便对照自己的运行结果。5.1 Kmeans 两次运行结果完全不同现象同一份数据、同一个 K跑两次kmeans_core得到的标签几乎对不上簇中心也差很远。原因初始中心是随机选的不同的初始化会落到不同的局部极小值而 Kmeans 没有任何机制保证跳出局部最优。解决固定随机种子rng(42)只解决可复现性真正降低方差的是 k-means 初始化或者干脆跑多次取目标函数最小的那次。示例代码里如果提供了nInit参数就是干这个的别把它当摆设。5.2 FCM 的隶属度矩阵几乎全是 1/C现象去模糊化之后簇的边界和 Kmeans 几乎一样甚至更乱隶属度矩阵每行都在 1/C 附近。原因m 取太大了或者数据没有归一化导致距离平方量级过大两者都会让dist^(-2/(m-1))的区分度被压平。解决先把 m 固定在 2检查数据是否做了标准化如果仍然偏平再往 1.5 调。调参顺序一定是先归一化、再调 m反过来做等于在错的基础上调错。5.3 SOM 训练后期震荡权重不收敛现象损失曲线前期下降后期在某个值附近来回跳最终权重矩阵看起来像随机噪声。原因学习率或邻域半径没有衰减或者衰减得太慢导致训练后期胜者邻域仍然在整体移动把已经收敛的局部结构又打散。解决确认eta和sigma都按exp(-t/totalIter)或线性衰减到接近 0如果数据集很小totalIter要相应调小否则衰减在训练结束前根本走不完。5.4 层次聚类出现一条细长的“链”现象dendrogram 上出现一个非常长、非常细的簇把本应分开的两团数据串在一起。原因用了单链接single链式效应把距离上“首尾相连”的点串成一条。解决换成average或ward如果数据噪声大先做一次粗粒度的 Kmeans 或 SOM 降采样再对簇中心做层次聚类链式效应会明显减弱。5.5 四个算法结果差异巨大却找不到谁对现象同一份数据Kmeans 给 3 个簇FCM 给 5 个层次聚类切出 4 个SOM 打出来 6 个不知道信谁。原因四个算法优化的目标函数根本不同簇数也不是同一个含义直接比标签数量没有意义。解决先用轮廓系数或 DB 指数在同一套标签上做横向比较下一章展开再结合数据的几何形状判断——球状数据上 Kmeans 的指标通常最高边界重叠的数据上 FCM 更稳拓扑结构明显的数据上 SOM 的邻域关系更有解释力。6. 设计取舍与验证用轮廓系数和 DB 指数把“感觉对”变成“证据对”四种算法放在一起最终要回答的问题是这份数据到底该用哪个。凭肉眼观察散点图只能排除明显不合适的选项真正做取舍时需要一套能在同一尺度上比较的指标。常用的两个是轮廓系数Silhouette和 Davies-Bouldin 指数前者越大越好、上限为 1后者越小越好、下限为 0两者都不依赖真实标签适合无监督场景下的横向对比。算法时间复杂度几何假设关键参数更适合的场景KmeansO(N·K·D·iter)凸的、球状簇K、初始化次数大样本、簇形状规则FCMO(N·K·D·iter)允许边界重叠C、m边界模糊、需要软归属层次聚类O(N²logN) 到 O(N³)依赖链接准则K、linkageType小样本、需要看合并过程SOMO(N·D·grid²·epoch)保留拓扑关系gridSize、eta0、sigma0需要降维可视化、邻近关系重要轮廓系数的计算逻辑是对每个样本算它到同簇其他点的平均距离 a再算它到最近的其他簇的平均距离 b单个样本的轮廓系数是(b-a)/max(a,b)最后对所有样本取均值。下面这段是示例代码里常见的实现function s silhouette_score(X, labels) % 返回所有样本轮廓系数的均值范围 [-1, 1] N size(X, 1); s zeros(N, 1); for i 1:N sameIdx labels labels(i); sameIdx(i) false; if ~any(sameIdx) % 单点簇轮廓系数记为 0 s(i) 0; continue; end a mean(pdist2(X(i, :), X(sameIdx, :))); otherLabels unique(labels(~sameIdx)); b inf; for k 1:numel(otherLabels) idx labels otherLabels(k); b min(b, mean(pdist2(X(i, :), X(idx, :)))); end s(i) (b - a) / max(a, b); end s mean(s); end第 5 行那个sameIdx(i) false是必须的否则样本会把自己算进同簇平均距离里a 被系统性拉低轮廓系数整体偏高。单点簇的处理也要留意b有定义但a无定义标成 0 是最常见的约定但如果你在对比四种算法时发现某个算法的单点簇特别多这个 0 会稀释整体均值需要单独看簇大小分布。DB 指数的思路更“簇级”对每个簇算它和别的簇的相似度簇内散度之和除以簇间中心距离取每个簇最相似的那个再对所有簇取平均。它不像轮廓系数那样逐样本算所以在样本量很大时更快但对簇形状的假设更粗。实际做对比时我一般两个都算轮廓系数看整体分离度DB 指数看簇间最坏情况两个指标结论一致时才敢下判断。验证这件事最容易犯的错是拿同一份数据既调 K 又看指标指标自然一路向好。比较稳妥的做法是留一份固定的评估集或者至少固定随机种子后只改算法、不改 K先看四种算法在同一个 K 下的相对排序再单独讨论 K 怎么选。K 的选择本身也是一个取舍轮廓系数在 K 增大时常常单调上升不能只看最大值要结合 dendrogram 或者业务上对簇数的约束一起判断。写到这里我想起早期做项目时的一次教训当时为了省事四种算法共用一份没归一化的数据结果 FCM 和 SOM 的结果完全不可信还以为是算法实现有问题翻了两天源码才发现是量纲在作怪。从那以后我每次跑聚类之前都强制走一遍“归一化 → 画二维散点 → 四算法同 K 对比指标”的流程再决定往哪个方向深挖。希望这份资源里的源码和示例能帮你把这套流程跑顺少走那两天的弯路。本文还有配套的精品资源点击获取