MATLAB聚类源码：Kmeans、FCM、SOM与层次聚类对比

简介面向机器学习初学者与MATLAB实践者的聚类算法代码合集聚焦无监督分类分析系统整理FCM模糊C均值、K-means、Kmedoid等经典方法的实现与示例适合课程实验、算法对比和科研数据预处理。压缩包共74个文件约2.05MB以40个m脚本为核心辅以gif、png、jpg演示图txt说明、htm与css页面以及mat数据、pdf文档兼顾算法运行、结果展示和资料查阅。已有4888人学习参考基础较好。读者可获得可直接运行的聚类流程包括数据标准化与反标准化、隶属度与聚类中心迭代、中心类算法对比、扩展聚类与有效性评价、PCA投影可视化并通过Demos示例理解参数调整、收敛判断及不同算法在模糊边界、非球形分布数据上的适用差异。资源目录按算法实现、演示脚本、评价指标与可视化辅助分区便于按主题检索和二次修改。1. 从无标签散点开始四种聚类源码先解决哪几个重复动作手上有一批没有标签的二维或高维样本想快速比较 Kmeans、FCM、SOM 神经网络聚类和层次聚类的效果四种聚类算法源代码及示例代码就是为这个场景准备的。它给出 MATLAB 可运行脚本和核心函数覆盖硬划分、模糊隶属度、竞争学习和树状切分适合做模式识别前的簇结构探查、图像分割粗分、故障状态划分和客户分群的从业者。真正省时间的地方不是算法本身而是同一套数据接口、同一套评估口径和同一套结果落盘方式让你不用每换一种算法就重写一遍预处理和可视化。2. K-means 与 FCM硬划分和隶属度矩阵怎么写进同一套 MATLAB 脚本2.1 硬聚类与模糊聚类的选型边界K-means 的目标很直接把 n 个样本分到 k 个簇最小化簇内平方误差。每个样本只属于一个簇输出是标签向量和簇中心。它的优点是对球形簇、尺度接近的数据收敛快缺点是初始化敏感、对噪声和离群点不客气。FCM 则在目标函数里引入隶属度矩阵 U样本以 u_ij 的隶属度同时属于多个簇目标函数写成 sum sum u_ij^m ||x_i - c_j||^2约束是每个样本对全部簇的隶属度之和为 1。m 是模糊指数m 越接近 1隶属度越接近 0/1 的硬划分m 越大隶属度越平均簇的边界越模糊。选型上如果样本边界清晰、簇数已知、数据量偏大K-means 更稳如果样本存在过渡状态比如设备从正常到异常的中间退化段或者图像里相邻组织没有硬边界FCM 的隶属度矩阵更有解释空间。常见做法是先用 K-means 跑一个基线再用 FCM 看隶属度分布如果大量样本的最大隶属度都在 0.5 附近说明簇结构本身不清晰不要急着调 m先检查特征和簇数。MATLAB 里kmeans和fcm都有现成函数但源码包通常会自带实现原因是工具箱依赖和参数默认值不一致。下面按自写核心循环的方式拆方便你改成自己的距离度量。2.2 K-means 核心代码与参数先构造一份带轻微重叠的三簇二维数据标准化后跑自写 K-means。标准化这一步不要省否则量纲大的特征会直接主导欧氏距离。rng(7); X [randn(80,2)*0.7 [0 0]; randn(80,2)*0.8 [3 1]; randn(80,2)*0.6 [1 3]]; X zscore(X); % 按列标准化保存 mu 和 sigma 以便后续新样本对齐 k 3; maxIter 200; tol 1e-6; [n, d] size(X); idx randperm(n, k); % 随机选 k 个样本做初始中心 C X(idx, :); labels zeros(n, 1); for iter 1:maxIter D pdist2(X, C); % n x k 距离矩阵 [~, labels] min(D, [], 2); Cnew zeros(k, d); for j 1:k if any(labels j) Cnew(j, :) mean(X(labels j, :), 1); else Cnew(j, :) X(randi(n), :); % 空簇重新播种避免中心留在原地 end end if norm(Cnew - C, fro) tol C Cnew; break; end C Cnew; end fprintf(迭代次数: %d\n, iter); fprintf(簇中心:\n); disp(C);这段代码的逻辑是「分配—更新—收敛判断」三步循环。pdist2(X,C)计算每个样本到每个中心的欧氏距离min(D,[],2)取最近中心作为硬标签更新阶段对每个簇取均值空簇用随机样本重新播种这是实际数据里必须处理的分支。maxIter设 200 足够覆盖大多数二维和中等维数据tol用 1e-6 控制中心矩阵的 Frobenius 范数变化数据量很大时可以放宽到 1e-4 减少空转。如果直接用 MATLAB 的kmeans建议显式写出重复次数和初始化方式[labels, C, sumd] kmeans(X, k, ... Distance, sqeuclidean, ... Start, plus, ... Replicates, 10, ... MaxIter, 200);Start,plus是 k-means 初始化Replicates,10表示用不同初始中心跑 10 次取最优能显著降低「一次跑出两个簇、另一次跑出三个簇」的玄学。sumd是每个簇内样本到中心的距离平方和后面做评估时可以直接用。2.3 FCM 核心代码与参数FCM 的核心是隶属度矩阵更新。下面这段代码用向量化方式写注意距离平方为 0 时要用 eps 兜底否则会出现 Inf 或 NaN。m 2.0; % 模糊指数常用 1.5 到 2.5 maxIter 150; tol 1e-5; U rand(n, k); U U ./ sum(U, 2); % 每行隶属度和为 1 for iter 1:maxIter Um U .^ m; C (Um * X) ./ sum(Um, 1); % 按隶属度加权更新中心 D2 pdist2(X, C).^2; D2(D2 eps) eps; % 避免除零 Unew (D2 .^ (-1/(m-1))) ./ ... sum(D2 .^ (-1/(m-1)), 2); obj sum(sum((Unew .^ m) .* D2)); if norm(Unew - U, fro) tol U Unew; break; end U Unew; end [~, labels] max(U, [], 2); fprintf(FCM 迭代次数: %d, 目标函数: %.6f\n, iter, obj);m是最需要试的参数。m 取 1.5 时隶属度已经比较硬接近 K-meansm 取 2.5 以上很多样本的最大隶属度会被压到 0.4 左右看起来「每个簇都沾一点」解释性反而变差。我的习惯是先固定 m2看最大隶属度的分布如果 80% 以上样本的最大隶属度大于 0.7说明簇结构清楚如果大量样本集中在 0.3 到 0.5先别调 m优先检查特征是否重复、簇数是否偏多。D2 .^ (-1/(m-1))是 FCM 隶属度更新公式的向量化写法等价于对每个样本计算1 / sum((d_ij/d_il)^(2/(m-1)))。目标函数obj用来观察单调下降趋势正常情况下应该逐步减小或基本稳定如果来回震荡通常是 m 太小或数据里有重复点。MATLAB 的 Fuzzy Logic Toolbox 提供fcm调用形式大致是[centers,U] fcm(X,k,options)options 里依次放 m、最大迭代次数和收敛阈值。源码包如果自带实现建议对照一次中心更新公式确认它用的是距离平方还是距离两者在 m 不等于 2 时结果会有差异。2.4 示例数据与运行结果检查跑完 K-means 和 FCM至少检查三件事簇中心是否稳定、目标函数是否收敛、标签和隶属度是否一致。轮廓系数可以用 MATLAB 的silhouette也可以自写后面第 4 章会给自写版本。先看一个快速检查脚本sK silhouette(X, labels_kmeans); sF silhouette(X, labels_fcm); fprintf(K-means 平均轮廓系数: %.4f\n, mean(sK)); fprintf(FCM 平均轮廓系数: %.4f\n, mean(sF)); figure; subplot(1,2,1); gscatter(X(:,1), X(:,2), labels_kmeans); title(K-means 硬划分); subplot(1,2,2); gscatter(X(:,1), X(:,2), labels_fcm); title(FCM 最大隶属度标签);轮廓系数接近 1 表示簇内紧、簇间远接近 0 表示簇边界重叠负值表示样本可能分错。二维示例上 K-means 和 FCM 的平均轮廓系数通常不会差太多FCM 的优势体现在隶属度矩阵你可以把最大隶属度低于 0.6 的样本单独标出来这些就是边界样本。实际项目里我会把这类样本交给人工复核或下游规则而不是硬塞进某个簇。提示K-means 和 FCM 的簇编号没有对应关系比较标签时不要直接算准确率先做簇匹配或用调整兰德指数。3. SOM 神经网络聚类竞争学习、邻域更新与命中图3.1 SOM 为什么在无标签数据里还能用SOM 是一种竞争学习神经网络结构上是一张低维网格每个网格节点有一个与输入同维的权重向量。训练时对每个样本找到距离最近的节点称为最佳匹配单元 BMU然后更新 BMU 及其邻域内的权重让网格在输入空间中展开。它和 K-means 的关键差别是K-means 只更新获胜中心SOM 还更新获胜节点的邻居因此保留了输入空间的拓扑关系。这个性质让 SOM 在无标签数据里既能做聚类也能做降维可视化特别适合先看「数据大概分几片」而不是一上来就定 k。SOM 的代价是参数更多网格尺寸、学习率、邻域半径、训练轮数都会影响结果。网格太小所有样本挤到少数节点网格太大节点稀疏命中图全是空格。常见做法是先用 5×5 或 6×6 网格跑一遍看命中分布再决定要不要放大。3.2 网络初始化与训练参数下面是一个简化版 SOM 训练循环输入还是前面的二维标准化数据。权重张量 W 的尺寸是 gridX×gridY×d。rng(11); gridX 5; gridY 5; W rand(gridX, gridY, d); % 随机初始化权重 lr0 0.5; % 初始学习率 sigma0 max(gridX, gridY) / 2; epochs 200; for ep 1:epochs lr lr0 * exp(-ep / epochs); sigma sigma0 * exp(-ep / epochs); order randperm(n); % 每个 epoch 打乱样本顺序 for ii 1:n x X(order(ii), :); best inf; br 1; bc 1; for a 1:gridX for b 1:gridY w squeeze(W(a,b,:)); dist sum((w - x).^2); if dist best best dist; br a; bc b; end end end for a 1:gridX for b 1:gridY h exp(-((a-br)^2 (b-bc)^2) / (2*sigma^2)); w squeeze(W(a,b,:)); W(a,b,:) w lr * h * (x - w); end end end endlr0控制每一步权重向样本靠拢的幅度0.5 对二维小数据够用特征维度高或样本噪声大时降到 0.1 到 0.2 更稳。sigma0是初始邻域半径取网格最大边长的一半随着 epoch 指数衰减后期基本只更新 BMU 附近。epochs设 200 是保守值样本多时可以到 500但要注意学习率衰减曲线别让后期学习率接近 0 还在空转。训练完后把每个样本映射到 BMU统计每个节点的命中数hit zeros(gridX, gridY); bmuIdx zeros(n, 1); for ii 1:n x X(ii, :); best inf; br 1; bc 1; for a 1:gridX for b 1:gridY w squeeze(W(a,b,:)); dist sum((w - x).^2); if dist best best dist; br a; bc b; end end end hit(br, bc) hit(br, bc) 1; bmuIdx(ii) sub2ind([gridX gridY], br, bc); end figure; imagesc(hit); colorbar; title(SOM 命中图);命中图上一个节点对应一个 BMU 编号可以直接把bmuIdx当簇标签用也可以对权重矩阵再做一次 K-means 或层次聚类把相邻节点合并成更大的簇。前者适合看局部密度后者适合得到少量可解释的簇。3.3 示例代码逐段拆解初始化阶段用rand给权重小随机值不要全部设成 0否则所有节点对同一个样本的距离相同BMU 选择会退化成固定顺序。训练阶段每个 epoch 打乱样本顺序是为了避免样本输入顺序影响拓扑展开这一点和随机梯度下降的常识一致。BMU 搜索用双重循环写目的是让逻辑清楚数据量大时可以改成张量运算把 W 拉平成(gridX*gridY)×d的矩阵用pdist2一次算完。邻域更新里的h是高斯邻域函数(a-br)^2(b-bc)^2是网格上的平方距离sigma越小邻域越窄。学习率和邻域半径都按exp(-ep/epochs)衰减这是最常用的指数衰减简单且不容易在后期震荡。参数上网格尺寸建议从 5×5 起步样本量几百到几千都可以先跑如果命中图里超过一半节点是空的把网格缩到 4×4 或 3×3如果大量样本集中在两三个节点网格放大到 8×8 再看。学习率 0.5 适合小数据0.1 适合高维或噪声数据。邻域半径的初始值不要小于 1否则训练早期就退化成只更新 BMU拓扑保持效果会变差。3.4 聚类结果可视化SOM 最实用的可视化是 U-matrix 和命中图。U-matrix 的每个节点值是该节点权重与邻居权重的平均距离距离大的地方通常是簇边界。Umat zeros(gridX, gridY); for a 1:gridX for b 1:gridY neigh []; for da -1:1 for db -1:1 if da 0 db 0 continue; end aa a da; bb b db; if aa 1 aa gridX bb 1 bb gridY neigh [neigh; squeeze(W(aa,bb,:))]; end end end w0 squeeze(W(a,b,:)); Umat(a,b) mean(sqrt(sum((neigh - w0).^2, 2))); end end figure; subplot(1,2,1); imagesc(Umat); colorbar; title(U-matrix); subplot(1,2,2); imagesc(hit); colorbar; title(命中图);U-matrix 里亮色带通常对应簇与簇之间的过渡区命中图里密集区对应样本集中的节点。把两张图叠在一起看能快速判断「数据是分成两片还是三片」。这里没有直接给出最终簇标签因为 SOM 的簇数不是预设的需要你根据命中图和 U-matrix 做一次切分。常见做法是对 BMU 坐标或权重矩阵再做 K-means簇数从 2 试到 6结合轮廓系数选。提示SOM 训练结果对随机初始化敏感正式跑之前固定rng种子并把网格尺寸、学习率、epoch 数一起记下来否则第二次跑出来的命中图可能完全不一样。4. 层次聚类与统一评估第四种算法怎么接进同一套输出接口4.1 层次聚类接入理由与距离选择层次聚类不要求预先给定簇数它先构造一棵树再按切割高度或最大簇数切分。对样本量不大、想保留层次结构的场景很合适比如先看样本之间是怎么合并的再决定切几刀。MATLAB 的linkage支持ward、average、complete、single等方法距离度量支持euclidean、cosine、correlation等。ward适合欧氏距离且希望簇内方差小average对噪声更稳complete倾向于紧凑簇。我的默认组合是ward euclidean先标准化再跑。Z linkage(X, ward, euclidean); figure; dendrogram(Z, 0); % 0 表示显示全部叶节点 T cluster(Z, maxclust, 3); % 按最大簇数切分Z是 (n-1)×3 的连接矩阵每行记录合并的两个簇编号和合并距离。dendrogram(Z,0)在样本多时会很密可以改成dendrogram(Z,30)只显示最后 30 次合并。cluster(Z,maxclust,3)按最大簇数切也可以换成cutoff, 阈值按距离切。注意层次聚类的时间和内存复杂度是 O(n^2)样本上万时不要直接跑全量先抽样或先用 K-means 粗分。4.2 统一函数签名与输出接口四种算法输出不一致K-means 给标签和中心FCM 给隶属度和中心SOM 给 BMU 命中层次聚类给树和切分标签。比较之前先统一成labels、centers、score三个输出。下面是一个批量脚本框架核心是把每个算法的结果塞进同一个结构体。methods {kmeans, fcm, som, hier}; results struct(method, {}, labels, {}, centers, {}, sil, {}, seconds, {}); for mi 1:numel(methods) t0 tic; switch methods{mi} case kmeans [labels, C] kmeans(X, 3, Start, plus, Replicates, 10); case fcm [labels, C] run_fcm(X, 3, 2.0, 150, 1e-5); case som [labels, C] run_som_cluster(X, 5, 5, 200); case hier Z linkage(X, ward, euclidean); labels cluster(Z, maxclust, 3); C zeros(3, size(X,2)); for j 1:3 C(j,:) mean(X(labelsj,:), 1); end end results(mi).method methods{mi}; results(mi).labels labels; results(mi).centers C; results(mi).sil mean_silhouette(X, labels); results(mi).seconds toc(t0); endrun_fcm和run_som_cluster分别封装前面的 FCM 循环和 SOM 命中统计返回最大隶属度标签或 BMU 标签以及对应的中心。层次聚类的中心用簇内均值补出来方便统一画图。mean_silhouette是自写函数下面给实现。这个框架的好处是后面换数据、换 k、换参数只需要改methods和参数表不用重复写评估代码。4.3 评估指标轮廓系数与 Davies-Bouldin轮廓系数对每个样本计算 a到同簇其他样本的平均距离和 b到最近其他簇的平均距离再取(b-a)/max(a,b)最后对所有样本取平均。自写版本如下注意排除样本自身。function s mean_silhouette(X, labels) n size(X, 1); s zeros(n, 1); for i 1:n same labels labels(i); same(i) false; if sum(same) 0 s(i) 0; continue; end a mean(pdist2(X(i,:), X(same,:))); others unique(labels(labels ~ labels(i))); b inf; for j 1:numel(others) b min(b, mean(pdist2(X(i,:), X(labelsothers(j),:)))); end s(i) (b - a) / max(a, b); end s mean(s); end这个实现是 O(n^2) 的二维几千样本还能跑上万样本要改用抽样或近似算法。Davies-Bouldin 指数用簇内散度和簇间距离的比值值越小越好适合和轮廓系数一起看。实际比较时不要只看一个指标轮廓系数偏好凸簇DBI 对球形簇更敏感SOM 的 BMU 标签在轮廓系数上可能吃亏但命中图的可解释性更好。算法需要预设簇数输出适合场景主要代价K-means是标签、中心球形簇、样本量大初始化敏感FCM是隶属度、中心边界模糊、需要过渡信息m 参数敏感SOM否但网格尺寸要定BMU 命中、权重拓扑探索、可视化训练参数多层次聚类否切割时定树、切分标签小样本、层次结构O(n^2) 内存4.4 批量跑示例脚本与结果落盘批量跑完把结果保存成.mat和 CSV后面复现或写报告直接取。writetable对字符数组和数值混合的表格支持很好注意列名不要用中文特殊符号。methodNames {results.method}; silVals [results.sil]; secVals [results.seconds]; T table(methodNames, silVals, secVals, ... VariableNames, {method, silhouette, seconds}); writetable(T, cluster_compare.csv); save(cluster_results.mat, results, X);落盘时把X和results一起保存避免后面想复现图却找不到当时的数据。seconds字段不是用来比谁快而是提醒你层次聚类在样本量上去后时间会陡增。批量脚本里每次循环前固定rng否则 K-means 和 SOM 的对比没有意义。结果表里如果某个算法的轮廓系数明显低先看它的簇数是不是和别的方法一致簇数不同时轮廓系数不能直接横比。提示统一评估前先确认所有算法用的是同一份标准化数据。层次聚类对距离尺度敏感K-means 对特征量纲敏感FCM 对 m 敏感三者不在同一预处理下比较结论很容易翻车。5. 避坑与排查MATLAB 聚类示例跑不通时的五个翻车现场5.1 每次 K-means 结果不一样簇编号还跳变现象是同一份数据跑两次簇的形状差不多但标签 1 和标签 2 对调画图颜色全变。原因有两个初始中心随机选以及不同簇的编号本身没有语义。解决分两步第一步固定随机种子并在kmeans里加Replicates第二步比较结果时不要直接比标签用簇中心匹配或调整兰德指数。代码上至少写成rng(42)再调用批量实验里每个重复也用独立种子不要把种子设在循环外面导致所有重复其实只跑了一次初始化。5.2 FCM 隶属度全是 0.2 左右看不出边界现象是 k5 时每个样本对五个簇的隶属度都在 0.2 附近最大隶属度没有区分度。原因通常是 m 太大、特征没标准化或者簇数明显多于真实结构。先把 m 降到 2.0 甚至 1.7检查zscore是否真的执行了再看最大隶属度的直方图。如果仍然平均减少 k 再跑一次。FCM 的目标函数在 m 较大时会把隶属度拉平这不是算法坏了而是参数在告诉你「这些簇之间没有明显边界」。5.3 SOM 所有样本命中同一个神经元现象是命中图只有一个格子里有数其他全空。原因常见于学习率过大、训练轮数太少、权重初始化太集中或者数据没有标准化导致某个维度主导距离。先把lr0从 0.5 降到 0.1epochs翻倍权重用rand而不是全零再检查X每列的均值和方差。如果网格是 5×5 而样本只有几十个也可能出现大量空节点这时候缩小网格比调学习率更直接。5.4 层次聚类内存爆掉或树状图糊成一团现象是样本量到几千以上时linkage报内存不足或者dendrogram画出来全是黑线。原因是层次聚类要存 O(n^2) 的距离矩阵树状图叶节点太多也无法阅读。解决方式是先抽样到几百到一千条或者先用 K-means 粗分到几十个微簇再对微簇中心做层次聚类。树状图用dendrogram(Z, 30)限制显示层数完整树仍然保留在Z里需要切分时用cluster即可。5.5 轮廓系数很高但业务上没法用现象是某个 k 的轮廓系数最高但簇里混进了完全不同的样本。原因是轮廓系数只看几何紧密度不看特征语义和业务约束。比如客户分群时收入维度和消费频次维度标准化后权重相同但业务上收入差异可能更重要。解决方式是把轮廓系数当筛选条件而不是决策条件配合簇中心解释、稳定性检验和人工抽样复核。我的血泪经验是先用轮廓系数排除明显坏的 k再用业务规则检查每个簇的中心是否说得通两步都过才采纳。6. 换成自己的数据参数扫描、种子固定与结果复现把自己的 CSV 或 MAT 数据接进来时先做三件事按列zscore并保存mu和sigma检查缺失值和重复点确认距离度量与标准化方式匹配。欧氏距离配标准化是默认组合余弦距离则不要先做中心化。参数扫描不要手写十几个脚本用一个循环把 k、m、网格尺寸串起来每次固定种子记录轮廓系数和运行时间。rng(2024); kList 2:8; best struct(sil, -inf, k, NaN, labels, []); for k kList for rep 1:10 rng(rep); % 每个重复独立种子保证可复现 [labels, C] kmeans(X, k, Start, plus, Replicates, 1); s mean_silhouette(X, labels); if s best.sil best.sil s; best.k k; best.labels labels; best.centers C; end end end fprintf(最佳 k%d, 轮廓系数%.4f\n, best.k, best.sil);这段代码的关键是「每个重复独立种子」而不是在循环外设一次rng。如果只在最外层设种子十次重复其实共享同一串随机数K-means 的初始化多样性会被吃掉。FCM 的 m 可以按 1.5、1.8、2.0、2.2、2.5 扫SOM 的网格按 4×4、5×5、6×6 扫每次把参数和得分写进表格。复现清单我一般固定四样原始数据哈希或文件名、标准化参数、随机种子、算法参数。少了任何一样过两周回来都只能重新猜。从那以后我每次跑聚类都先固定种子、保存归一化参数再做参数扫描不然连自己上周的图都复现不了。希望帮到你。本文还有配套的精品资源点击获取