简介:面向机器学习初学者与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-means;m 取 2.5 以上,很多样本的最大隶属度会被压到 0.4 左右,看起来「每个簇都沾一点」,解释性反而变差。我的习惯是先固定 m=2,看最大隶属度的分布:如果 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(labels==j,:), 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(labels==others(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 左右,看不出边界
现象是 k=5 时每个样本对五个簇的隶属度都在 0.2 附近,最大隶属度没有区分度。原因通常是 m 太大、特征没标准化,或者簇数明显多于真实结构。先把 m 降到 2.0 甚至 1.7,检查zscore是否真的执行了,再看最大隶属度的直方图。如果仍然平均,减少 k 再跑一次。FCM 的目标函数在 m 较大时会把隶属度拉平,这不是算法坏了,而是参数在告诉你「这些簇之间没有明显边界」。
5.3 SOM 所有样本命中同一个神经元
现象是命中图只有一个格子里有数,其他全空。原因常见于学习率过大、训练轮数太少、权重初始化太集中,或者数据没有标准化导致某个维度主导距离。先把lr0从 0.5 降到 0.1,epochs翻倍,权重用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 扫,每次把参数和得分写进表格。复现清单我一般固定四样:原始数据哈希或文件名、标准化参数、随机种子、算法参数。少了任何一样,过两周回来都只能重新猜。
从那以后我每次跑聚类都先固定种子、保存归一化参数,再做参数扫描,不然连自己上周的图都复现不了。希望帮到你。
本文还有配套的精品资源,点击获取