简介:模糊C均值聚类(FCM)算法因其能刻画数据点对多个类别的模糊隶属关系,在模式识别与数据挖掘中常被用于处理边界不清或含噪声的数据。这份MATLAB源码包面向需要快速上手FCM的科研人员与学生,包含完整的算法主函数与运行脚本,并附有Iris等测试数据,可直接演示聚类效果与优缺点。包内共5个文件,以.m源码和.txt数据说明文件为主,整体仅206KB,轻量易部署。已有3500余人学习下载,内容覆盖模糊因子m、隶属度迭代更新、类中心重算等关键环节,并支持结合Davies-Bouldin指数等内部指标评估结果,适合用于课程实验、论文复现及算法对比。通过运行示例可直观感受FCM对非球形簇的适应能力,同时理解其对初始中心和参数设置的敏感性。
1. 模糊C均值聚类:先承认数据边界是模糊的,再谈聚类
做客户分群时,最常遇到的问题是“这个用户到底属于哪一群”,但真实答案往往是他在两个群之间都有归属。图像分割里的边缘像素同样棘手,它既像前景又像背景,用 K-means 做非此即彼的硬判断,效果会很生硬。模糊C均值聚类(FCM)允许一个样本同时以不同隶属度属于多个类别,因此对噪声、重叠分布和异常值更稳健。这个算法由 Bezdek 在 1973 年系统化,把经典 K-means 的目标函数改写成带模糊指数 m 的隶属度加权形式。源码包里的 FCMCluster.m 是核心聚类函数,FCMmain.m 是主脚本,iris.txt 则提供了可直接验证效果的 Iris 数据,适合想快速复现算法、理解参数影响的人。
2. FCM 的目标函数与隶属度更新公式:模糊C均值聚类的数学骨架
2.1 从硬划分到模糊隶属度
K-means 的核心是给每个样本指派唯一的簇标签,目标函数是组内平方和。FCM 把这种硬指派改成软指派,引入一个 n×c 的隶属度矩阵 U,其中 n 是样本数,c 是聚类数。矩阵里的元素 u_{ij} 表示第 j 个样本属于第 i 个簇的程度,取值范围在 0 到 1 之间,并且对任意一个样本 j,所有类别上的隶属度之和为 1。这个约束是 FCM 与 K-means 最根本的差异:一个样本可以同时属于多个簇,只是程度不同。
在 MATLAB 实现里,U 的布局需要先确认。我习惯把 U 排成 n 行 c 列,每一行是某个样本对所有簇的隶属度,每一列对应一个簇中心。有的老代码会把 U 存成 c 行 n 列,读源码时如果发现维度对不上,先转置再往下看。FCM 的目标函数把 K-means 的 0/1 指派权重替换成 u_{ij}^m,m 被称为模糊指数,它直接控制隶属度曲线的软硬程度,也是模糊C聚类算法里最需要调的核心参数。
2.2 拉格朗日乘子法与两个迭代公式
在约束 Σ_{i=1}^c u_{ij}=1 下最小化目标函数,标准的做法是用拉格朗日乘子法。对每个样本引入一个乘子 λ_j,把约束写进目标函数,再对 u_{ij} 和 v_i 分别求偏导并令导数为零,就能得到 FCM 的两个核心迭代公式。
隶属度更新公式是:
u_{ij} = 1 / Σ_{k=1}^{c} ( ||x_j - v_i|| / ||x_j - v_k|| )^{2/(m-1)}
簇中心更新公式是:
v_i = Σ_{j=1}^{n} u_{ij}^m x_j / Σ_{j=1}^{n} u_{ij}^m
第二个公式看起来复杂,但本质上就是对所有样本做加权平均,权重是隶属度的 m 次方。m 越大,权重越平滑;m 越小,距离近的样本越能主导中心位置。两个公式交替执行,直到簇中心的变化量小于预设阈值或达到最大迭代次数。分母里的 2/(m-1) 值得注意,当 m 接近 1 时指数非常大,距离差异会被急剧放大,所以 m 越接近 1,FCM 的行为就越像 K-means。
2.3 FCMCluster.m 核心循环:向量化实现与防除零
下面这段代码对应 FCMCluster.m 的核心迭代部分,我按最常见的 MATLAB 实现风格给出:
function [U, V, obj_fcm] = FCMCluster(X, c, m, max_iter, epsilon) % X: n x d 数据矩阵 % c: 聚类数 % m: 模糊指数 % max_iter: 最大迭代次数 % epsilon: 中心变化阈值 n = size(X, 1); V = X(randperm(n, c), :); % 随机选 c 个样本做初始中心 obj_fcm = zeros(max_iter, 1); for iter = 1:max_iter D = pdist2(X, V); % n x c 距离矩阵 D(D < eps) = eps; % 避免除零 invD = 1 ./ D; U = invD .^ (2 / (m - 1)); % 按距离倒数计算隶属度 U = U ./ sum(U, 2); % 行归一化,满足隶属度和为 1 V_new = (U .^ m)' * X ./ sum(U .^ m, 1)'; V_new(isnan(V_new)) = V(isnan(V_new)); % 空簇保护 obj_fcm(iter) = sum(sum((U .^ m) .* (pdist2(X, V_new) .^ 2))); if norm(V_new - V, 'fro') < epsilon V = V_new; break; end V = V_new; end obj_fcm = obj_fcm(1:iter); end代码里 pdist2 用于计算任意样本到任意簇中心的欧氏距离,输出是一个 n×c 矩阵。D(D<eps)=eps 是必须的,因为如果一个样本恰好和某个初始中心重合,距离为 0,invD 会变成 Inf,最终导致 U 出现 NaN。V_new 的更新式里,(U .^ m)' * X先对所有样本的加权特征求和,再除以sum(U .^ m, 1)'即该簇的隶属度总和。空簇保护把分母为 0 时产生的 NaN 替换成旧中心,让后续迭代有机会把样本重新拉回来。
正常情况下 obj_fcm 应该单调下降。如果迭代到一半目标函数上升,问题多半出在初始化:随机选中的中心离其他样本太远,某个簇在迭代初期就几乎失去所有隶属度。遇到这种情况,可以多跑几次,或者直接在初始化时用 K-means 的结果作为 FCM 起点。
2.4 模糊指数 m 的边界行为与建议
m 的选取直接影响聚类质量,这是模糊C均值聚类算法最被人诟病的缺点之一。下面的表格总结了不同 m 值下的表现:
| m 取值 | 聚类行为 | 工程建议 |
|---|---|---|
| m → 1 | 逼近 K-means 硬聚类 | 很少使用,结果与 K-means 几乎一致 |
| m = 1.5 | 隶属度较锐利 | 簇重叠明显时可用 |
| m = 2 | 最常用默认值 | 大部分论文和 MATLAB 示例的起点 |
| m = 2.5 | 边界更模糊,对噪声更钝 | 数据噪声大时尝试,但簇心会互相牵制 |
| m > 3 | 所有隶属度趋近 1/c | 通常不建议,中心退化成全局均值 |
从工程角度看,m=2 是最安全的起点,但不应该不做验证就直接用。我一般会在固定 c 和初始化的前提下,把 m 从 1.5 到 2.5 按 0.1 步长扫一遍,同时观察目标函数下降曲线和最终簇中心距离。如果某个 m 值附近出现目标函数突然上升或簇中心大幅跳变,说明该区域的数值稳定性已经变差,应该避开。
提示:不确定 m 时,先在 m=2 跑通流程,再对比 1.8 和 2.2 的结果。如果三个配置给出的硬标签差异小于 5%,说明数据本身对 m 不敏感,可以放心使用默认可调低迭代次数。
3. 用 FCMmain.m 跑通 iris.txt:MATLAB 下的模糊聚类参数设置与坑
3.1 读取数据与标准化
iris.txt 是经典的 Iris 数据集,包含 150 个样本、4 个特征和 3 个品种。源码包把 iris.txt 和 FCMmain.m 放在同一目录,因此直接用 load 就能读进来。主脚本里常见的读取方式是:
data = load('iris.txt'); X = data(:, 1:end-1); % 前四列是花萼长度、花萼宽度、花瓣长度、花瓣宽度 true_label = data(:, end); % 如果文件最后一列是类别编号则取出 X = zscore(X); % 标准化,避免花瓣长度的量纲压过花萼宽度zscore 会把每个特征变成均值为 0、方差为 1,这对基于欧氏距离的 FCM 很重要。如果不做标准化,花瓣长度在数值上远大于花萼宽度,距离计算会被它主导,聚类结果几乎等价于只用这一个特征。如果 iris.txt 里没有真实标签,就把 true_label 那行注释掉,只保留 X 继续跑。
3.2 FCMmain.m 的调用方式与输出含义
FCMmain.m 是主脚本,它的职责是设置聚类数 c、模糊指数 m、迭代次数和终止阈值,然后调用 FCMCluster.m。一个可以直接运行的结构如下:
clc; clear; close all; data = load('iris.txt'); X = data(:, 1:end-1); X = zscore(X); c = 3; % Iris 有三个品种 m = 2; % 模糊指数 max_iter = 200; epsilon = 1e-5; [U, V, obj_fcm] = FCMCluster(X, c, m, max_iter, epsilon); [~, label] = max(U, [], 2); % 每个样本取隶属度最大的类别作为硬标签 disp('前 10 个样本的隶属度向量:'); disp(U(1:10, :)); figure; plot(obj_fcm, 'o-'); xlabel('迭代次数'); ylabel('目标函数值');隶属度矩阵 U 的每一行是某个样本对 3 个簇的归属程度,行和恒为 1。max(U,[],2) 会把模糊结果转成硬标签,这个硬标签只用于后续评估,如果业务上需要软分类,直接保留 U 即可。obj_fcm 记录了每次迭代的目标函数值,正常情况下会单调下降并在 10 到 20 轮内收敛。如果看到折线呈锯齿状,优先检查距离矩阵里是否出现了 NaN,或者 m 是否设置得过小导致数值震荡。
3.3 四个关键参数:c、m、max_iter、epsilon
FCMmain.m 的核心参数只有四个,理解了它们的作用,整个算法就可以迁移到其他数据上:
| 参数 | 默认值 | 作用 | 建议范围 |
|---|---|---|---|
| c | 3 | 聚类数,需要先验或有效性指标确定 | 2~sqrt(n) |
| m | 2 | 模糊指数,控制隶属度软硬程度 | 1.5~2.5 |
| max_iter | 100 | 最大迭代次数,防止不收敛死循环 | 100~500 |
| epsilon | 1e-5 | 中心变化的 Frobenius 范数阈值 | 1e-6~1e-3 |
c 的确定不依赖 matlab 优化工具箱,FCM 的迭代本身只需要矩阵运算。c 太小会把不同类别强行压到一起,太大则会出现簇中心扎堆,后文会提到用 Davies-Bouldin 指数来选。m 影响的是边界样本,对全局结构影响没那么大。max_iter 在数据量大时可以调到 300,因为每次迭代都要算 n×c 的距离矩阵,迭代次数设置得过大会浪费计算资源。epsilon 越小结果越精确,但对初始化越敏感,工程上 1e-5 已经足够。
3.4 从模糊结果到硬标签:可视化与准确率评估
Iris 数据自带标签,可以用来评估模糊C均值聚类效果。先画出前两个特征和簇中心:
figure; gscatter(X(:,1), X(:,2), true_label); hold on; plot(V(:,1), V(:,2), 'kx', 'MarkerSize', 14, 'LineWidth', 2); legend('真实类别', '簇中心');聚类得到的簇标签编号和真实类别编号不一定对应,直接算准确率会偏低,所以需要做一个标签匹配:
conf = confusionmat(true_label, label); [~, perm] = max(conf, [], 2); % 每一行真实类别被映射到哪个聚类标签 mapped_label = zeros(size(label)); for i = 1:c mapped_label(label == i) = perm(i); end acc = sum(mapped_label == true_label) / length(true_label); fprintf('匹配后准确率: %.2f%%\n', acc * 100);confusionmat 得到 3×3 计数矩阵,max 的第二个输出把每个真实类别映射到出现次数最多的聚类簇。这个贪心匹配在类别数少时足够用,类别多了可以用匈牙利算法做全局最优匹配。需要提醒的是,聚类准确率不能完全代表 FCM 本身的性能,真实标签只是评价基准,在实际应用中我们往往更关心 U 矩阵能否揭示重叠区域的分布。
3.5 常见坑:CFM.txt、CMF.txt 和空簇
压缩包里的 CFM.txt 和 CMF.txt 通常是聚类过程中导出的中心矩阵或隶属度矩阵。运行type CFM.txt就能看到内容,如果是 c×4 的数值,那是最终簇中心;如果是 n×c 的概率数值,那就是最终隶属度矩阵。很多 MATLAB 老代码会把中间结果写进文本文件,方便在 Python 或 Origin 里二次绘图。
最常见的坑是空簇。随机初始中心如果选到孤立点,它可能在整个迭代过程中没有一个样本真正靠近,导致该簇隶属度趋近于 0。FCMCluster.m 里的空簇保护会把 NaN 替换成旧中心,让后续迭代有机会“救回来”。另一个坑是 m 设置过大,比如 m=5,此时所有隶属度都接近 1/c,簇中心会被拉向全局均值,聚类结果几乎失去区分度。遇到这种结果,先看目标函数曲线是否平缓,再检查 m 是否超出合理范围。对比 kmeans聚类算法matlab 自带的 kmeans,FCM 的时间开销明显更高,但边界样本的软归属信息是 K-means 给不了的。
4. 从 FCMmain.m 向实用化:初始化修正、有效性指标与加速
4.1 多次启动选最优目标函数
FCM 对初始簇中心敏感,这是它最主要的短板。随机初始化一次可能陷入局部极小值,不同运行得到的簇中心差异很大。常见做法是把 FCMmain.m 里单次调用改成多次启动:
best_obj = inf; for r = 1:20 [U_tmp, V_tmp, obj_tmp] = FCMCluster(X, c, m, max_iter, epsilon); if obj_tmp(end) < best_obj best_obj = obj_tmp(end); U = U_tmp; V = V_tmp; end end这不会增加多少代码量,但能把局部最优带来的结果方差按数量级压低。如果还想更稳,可以把初始中心替换成 K-means++ 的输出,让初始中心彼此拉开距离,减少空簇出现的概率。
4.2 用 Davies-Bouldin 指数确定聚类数 c
确定 c 是模糊C聚类算法里比调 m 更麻烦的事情。手写 Davies-Bouldin 指数并不难,它衡量簇内离散度与簇间中心距离的比值,越小越好:
[~, hard_label] = max(U, [], 2); S = zeros(c, 1); for i = 1:c pts = X(hard_label == i, :); if ~isempty(pts) S(i) = mean(sqrt(sum((pts - V(i,:)).^2, 2))); end end M = pdist2(V, V); db = 0; for i = 1:c ratios = zeros(1, c-1); k = 1; for j = 1:c if i ~= j ratios(k) = (S(i) + S(j)) / M(i, j); k = k + 1; end end db = db + max(ratios); end db = db / c;对 c=2 到 8 各跑一轮多次启动,然后取 DB 最小的 c。MATLAB 自带的 evalclusters 能算 Calinski-Harabasz 等指标,但手写这段不需要额外工具箱,并且可以直接复用 FCMCluster.m 输出的 U 和 V。
4.3 向量化与降采样:百万级数据的生存之道
FCM 每一步都要对每个样本计算到所有簇中心的距离,复杂度是 O(n·c·d·iter)。当 n 到百万、d 上百时,pdist2 一次性生成 n×c 矩阵会直接占满内存。我一般会先把特征用 PCA 降到 10 维以内,再做一个无放回抽样,比如只抽 5 万条训练 FCM,然后把簇中心拿到全量数据上做一次最近邻分配。具体做法是在 FCMmain.m 读取数据之后插入一行X = X(randperm(n, min(n, 50000)), :);,其余逻辑完全不用改。这样既保留模糊聚类发现重叠结构的能力,又避免每次迭代都在全量数据上算距离。当样本量达到百万级时,任何矩阵形式的距离计算都会成为瓶颈,提前抽样的收益远比调整 m 和 epsilon 来得明显。
本文还有配套的精品资源,点击获取