简介这份MATLAB编写的GMM代码包聚焦高斯混合模型在图像背景减除中的应用适合计算机视觉学习者、图像处理研究者及需要快速实现前景检测的开发者。压缩包共7个文件以5个m脚本、1个eps图和1个mat数据文件为主体整体约40KBm脚本覆盖模型初始化、EM迭代、高斯概率计算、结果绘制与演示入口。目前已有373人学习下载。代码结构清晰注释完整便于二次开发。通过阅读这些代码可深入理解GMM利用多个高斯分量对像素分布建模、通过EM算法迭代更新参数并设置阈值区分背景与前景的完整流程配套数据与结果图可直接运行验证直观展示分割效果方便对照调试。在监控视频运动目标检测、智能交通等场景中具备实用参考价值也便于在此基础上进行在线学习和自适应阈值优化。1. GMM 的 MATLAB 实现手写 EM 与 fitgmdist 该怎么选GMM高斯混合模型在 MATLAB 里最常见的实现方式有两种直接调 fitgmdist或者从零手写 EM 迭代。前者一行代码出结果后者能让你在参数翻车时知道该查哪里。这套代码包把两条路线都打包好了包含手写 EM 训练函数、fitgmdist 标准调用模板、BIC 选 K 脚本、协方差椭圆与概率密度曲面可视化以及异常检测和图像分割的扩展 Demo。适合正在做课程设计、论文实验的初学者也适合在正式项目里需要做密度估计或软聚类的工程师。下面从模型公式走到每一步代码再落到参数怎么设、报错怎么查。2. 高斯混合模型与 EM 迭代从责任度公式到矩阵化更新2.1 三个核心参数GMM 到底在拟合什么一个 K 分量的混合高斯模型概率密度是 K 个高斯分布按权重叠加p(x) Σ w_k · N(x | μ_k, Σ_k)其中 k 从 1 到 K。w_k 是第 k 个分量的权重满足 w_k ≥ 0 且所有权重之和等于 1μ_k 是 d 维均值向量Σ_k 是 d×d 协方差矩阵。在 MATLAB 的数据结构里一组完整参数对应三个变量权重 w 是 K×1 向量均值 mu 是 K×d 矩阵协方差 Sigma 是 d×d×K 的三维数组。用 fitgmdist 训练完这三个量分别挂在 gm.ComponentProportion、gm.mu 和 gm.Sigma 上。需要留意的是参数数量。以 full 协方差为例一个 K 分量、d 维数据的 GMM 总共有 (K−1) K·d K·d·(d1)/2 个自由参数。这个数字在选 K 和选协方差类型时是核心依据K 每加 1就要多付 d·(d1)/2 个协方差参数样本量不够时模型必然过拟合这也是后面 BIC 选 K 的逻辑起点。为什么不能用 K-means 代替因为 K-means 只给硬标签每个样本要么属于这簇要么属于那簇而 GMM 输出的是软标签也就是每个样本属于每个分量的后验概率。这个后验概率在做异常检测、不确定性估计时比硬标签有用得多。代价是 GMM 没有闭式解必须用 EM 迭代逼近这也是这份代码包最核心的部分。2.2 E 步与 M 步先算责任度再更新参数EM 的核心工具是责任度记作 γ(i,k)表示样本 i 由第 k 个高斯生成的概率。E 步就是把它算出来γ(i,k) w_k · N(x_i | μ_k, Σ_k) / Σ_j w_j · N(x_i | μ_j, Σ_j)。M 步用这些责任度加权更新参数N_k 等于所有 γ(i,k) 求和w_k N_k / nμ_k 等于 γ 加权平均Σ_k 等于 γ 加权的叉积和再除以 N_k。这套公式写成 MATLAB 最顺的方式是矩阵化。下面这份代码是这个资源包里的核心函数完整实现了从初始化到收敛判断的整个流程function [mu, Sigma, w, llh] gmm_em(X, K, maxIter, tol) % 手写 GMM-EM 训练函数 % 输入: % X - n x d 数据矩阵, n 个样本, d 维特征 % K - 高斯分量个数 % maxIter - 最大迭代轮数 % tol - log-likelihood 变化量阈值, 低于它提前停止 % 输出: % mu - K x d 均值矩阵 % Sigma - d x d x K 协方差数组 % w - K x 1 权重向量 % llh - 每轮迭代的 log-likelihood 序列 [n, d] size(X); % ---------- 初始化: 用 K-means 给一个合理的起点 ---------- [idx, C] kmeans(X, K, Replicates, 5, MaxIter, 200); mu C; % C 是 K x d, 直接作为均值初值 Sigma zeros(d, d, K); for k 1:K Sigma(:, :, k) cov(X(idx k, :)) 1e-6 * eye(d); end w accumarray(idx, 1, [K 1]) / n; % 每个簇的样本比例做权重初值 llh zeros(maxIter, 1); for iter 1:maxIter % E 步 gamma zeros(n, K); for k 1:K gamma(:, k) w(k) * mvnpdf(X, mu(k, :), Sigma(:, :, k)); end gamma gamma ./ sum(gamma, 2); % 行归一化, 得到责任度 % M 步 Nk sum(gamma, 1); % 每个分量的有效样本数 w Nk / n; % 权重更新 for k 1:K mu(k, :) (gamma(:, k) * X) / Nk(k); % 均值更新 Xc X - mu(k, :); % 中心化 Sigma(:, :, k) (Xc .* gamma(:, k)) * Xc / Nk(k); Sigma(:, :, k) Sigma(:, :, k) 1e-6 * eye(d); % jitter 保正定 end % 监控: 计算本轮 log-likelihood L zeros(n, K); for k 1:K L(:, k) w(k) * mvnpdf(X, mu(k, :), Sigma(:, :, k)); end llh(iter) sum(log(sum(L, 2))); if iter 1 abs(llh(iter) - llh(iter - 1)) tol llh llh(1:iter); break; end end代码里三个容易忽略的点值得单独说。第一kmeans 的 Replicates 参数设为 5意思是 kmeans 自己从 5 个不同起点跑 5 遍取最优这一步能显著降低后续 EM 陷进局部最优的概率代价只是几十毫秒。第二M 步里(Xc .* gamma(:, k)) * Xc / Nk(k)这行是整段代码最精华的矩阵化写法先把中心化矩阵按列乘上责任度再左乘转置一步得到加权的叉积和。它等价于对每个样本乘上责任度后做外积累加避免了在 d 维空间里写双重循环。注意这种写法依赖 R2016b 起的隐式扩展老版本会报维度不匹配改成Xc .* repmat(gamma(:, k), 1, d)即可。第三每次更新后加 1e-6 单位阵是对奇异协方差的兜底后面第 4 章会讲它和 fitgmdist 的 RegularizationValue 的对应关系。关于 log-likelihood 的计算sum(log(sum(L, 2)))是对每个样本取 log 后求和这是 EM 的收敛判据。理论上 EM 保证它单调不降如果你看到它下降基本可以断定数值上出了问题——E 步分母出现 0 或者 mvnpdf 下溢第 5 章会专门处理这种情况。2.3 手写 EM 和 fitgmdist两条路线的边界既然 MATLAB 自带 fitgmdist为什么这份资源还要保留手写版本我的判断是两者用途不同不是替代关系。fitgmdist 是个黑匣子给一个 X 和一个 K 就出结果内部自动做了初始化、正则、复跑但在以下场景它帮不了你你要固定某个分量的均值或协方差比如先验知道某个簇的中心位置你要在每一轮迭代后插入自定义监控或约束你要验证自己的新算法在线 GMM、变分 GMM 之类时需要绝对可靠的 baseline。反过来手写版本在数据量大、维度高的时候性能不如 fitgmdist。fitgmdist 在多轮迭代里做了不少工程优化比如对初始化起点做筛选、对协方差分解做缓存这些手写代码照抄成本很高。所以我的建议是学习、调试、验算法用手写版正式实验和数据生产用 fitgmdist但要把手写版的 log-likelihood 终值拿来做对拍确认两个结果一致我习惯管这叫双跑验证。对拍时注意一个坑fitgmdist 对象返回的 NlogL 是负对数似然数值越小越好手写函数返回的 llh 是正的 log-likelihood数值越大越好。比较之前记得把其中一个取负不然差出几个量级都不知道是代码写错了还是在比两个相反的数。3. 把 GMM 跑通合成数据、训练对拍与可视化的完整流程3.1 准备数据用 mvnrnd 造三类形状不同的高斯样本验证聚类算法最忌讳用完美的球形数据那样 K-means 都能做对GMM 的优势体现不出来。我一般会造三个形状差异明显的簇一簇是旋转的椭圆协方差一簇是标准圆一簇是稍微带倾斜的窄椭圆而且每簇样本数不一样用来考验 EM 对权重不对称的恢复能力。rng(42); % 固定随机种子, 保证结果可复现 % 簇1: 旋转椭圆, 协方差有非对角项 X1 mvnrnd([0 0], [1.2 0.5; 0.5 0.8], 300); % 簇2: 标准圆, 对角且各项相同 X2 mvnrnd([4 4], [0.6 0; 0 0.6], 400); % 簇3: 窄椭圆, 带轻微旋转 X3 mvnrnd([-2 3], [1.5 -0.4; -0.4 0.5], 300); X [X1; X2; X3];mvnrnd 的第一个参数是均值向量第二个是协方差矩阵第三个是样本数。簇 1 的协方差矩阵 [1.2 0.5; 0.5 0.8] 非对角意味着两个特征线性相关密度等高线是倾斜的椭圆这种结构只有 full 协方差的 GMM 才能还原。簇 2 的协方差是 0.6 的纯对角对应圆形。数据拼成 X 之后n1000d2。这里为什么要固定 rng(42)两个原因一是你调参时要能确认结果变好是参数改对了而不是随机数变了二是后面要和 fitgmdist 的对拍保持同一份输入。血的教训不固定种子你永远说不清是模型改进了还是运气改进了。3.2 训练fitgmdist 一行接住所有参数手写版负责对拍% fitgmdist 训练: full 协方差 1e-5 正则 5 次复跑 gm fitgmdist(X, 3, ... CovarianceType, full, ... RegularizationValue, 1e-5, ... Replicates, 5, ... Options, statset(MaxIter, 500, TolFun, 1e-6)); % 手写 EM 训练, 300 轮上限, 1e-6 收敛阈值 [mu_my, Sigma_my, w_my, llh_my] gmm_em(X, 3, 300, 1e-6); % 对拍: fitgmdist 给的是负对数似然, 取负号和手写版比较 disp([fitgmdist 负对数似然: , num2str(gm.NlogL)]); disp([手写 EM 负对数似然: , num2str(-llh_my(end))]);fitgmdist 的训练结果都以属性形式挂在 gm 对象上gm.mu 是 K×d 均值gm.Sigma 是 d×d×K 协方差数组gm.ComponentProportion 是 K×1 权重gm.NlogL 是负对数似然。注意 Options 里的 TolFun 控制的是对数似然变化量的收敛阈值MaxIter 是迭代上限这两个参数在生产环境里建议直接放宽到 500 和 1e-6默认值在数据不干净时经常提前停下。对拍结果一般两类如果手写版的负对数似然比 fitgmdist 大 1e-3 以上先别怀疑手写代码查查 fitgmdist 是不是在初始化阶段选了更差的起点。把 Replicates 加到 10 再比一次。如果两者相等或差在 1e-6 量级说明实现没问题可以放心用手写版做定制。由于手写版加的 jitter 是 1e-6 而 fitgmdist 是 1e-5终值差一点点属于正常量级一致就行。3.3 可视化聚类标签、协方差椭圆与密度曲面训练完先看硬聚类结果用 cluster 函数给每个样本分配它后验概率最大的分量label cluster(gm, X); % 返回 n x 1 的标签向量 figure; gscatter(X(:,1), X(:,2), label); title(GMM 聚类结果 (K3, full 协方差)); hold on; for k 1:gm.NumComponents plot_ellipse(gm.mu(k, :), gm.Sigma(:, :, k), 2); end hold off;椭圆用 plot_ellipse 辅助函数画逻辑是对协方差做特征分解eig 得到特征向量方向 V 和特征值 Dsqrt(D) 就是椭圆长短轴半径scale 取 2 对应二维高斯约 95% 的置信区域。这份代码包里已经把这个函数单独拆出来了function plot_ellipse(mu, Sigma, scale) % 画高斯分布的协方差椭圆 % mu: 1 x d 均值; Sigma: d x d 协方差; scale: 半径倍数 [V, D] eig(Sigma); theta linspace(0, 2*pi, 100); xy [scale * sqrt(D(1,1)) * cos(theta); scale * sqrt(D(2,2)) * sin(theta)]; pts V * xy; % 把单位椭圆旋到协方差的方向 plot(mu(1) pts(1,:), mu(2) pts(2,:), k-, LineWidth, 1.2); end这套可视化能一眼看出三个东西椭圆是否覆盖对应簇的数据、椭圆间有没有过度重叠、以及有没有空分量。如果某个椭圆里几乎没有点基本就是 K 设大了或者那个分量在 EM 中被淘汰这就是第 5 章要展开的坑。密度曲面是另一种视角适合展示 GMM 作为密度估计器的能力% 在数据范围内铺 80 x 80 的网格 gx linspace(min(X(:,1))-1, max(X(:,1))1, 80); gy linspace(min(X(:,2))-1, max(X(:,2))1, 80); [GX, GY] meshgrid(gx, gy); Z reshape(pdf(gm, [GX(:) GY(:)]), size(GX)); figure; surf(GX, GY, Z, EdgeColor, none); xlabel(x1); ylabel(x2); zlabel(密度); view(2); % 俯视角度, 密度越高颜色越暖pdf(gm, points) 接收 n×d 的查询点矩阵返回每个点上的混合密度。网格采样时注意范围要稍微外扩否则边界处的密度被截断曲面看起来像被切了一刀。4. 参数调优BIC 选 K、协方差结构与正则项的三板斧4.1 用 BIC/AIC 选 K别让模型参数超过信息量K 是 GMM 里最重要的超参数。选小了欠拟合两个真实簇被并成一个选大了出现过拟合加空分量。BIC 和 AIC 是两条现成的信息准则fitgmdist 把两个值都算好了直接读取 gm.BIC 和 gm.AIC 就行。% 扫 K1:8, 每个 K 用 3 次复跑取最优, 避免初始化偏差 Ks 1:8; bic zeros(size(Ks)); aic zeros(size(Ks)); for i 1:numel(Ks) gm_tmp fitgmdist(X, Ks(i), ... RegularizationValue, 1e-5, ... Replicates, 3); bic(i) gm_tmp.BIC; aic(i) gm_tmp.AIC; end figure; plot(Ks, bic, o-, Ks, aic, s--); xlabel(K); ylabel(信息准则值); legend(BIC, AIC);选 K 的标准是看曲线拐点K 从 1 开始信息准则值快速下降到某个 K 之后下降变缓甚至回升那个位置就是合适的候选。BIC 对参数数量的惩罚比 AIC 重因此 BIC 选出来的 K 通常比 AIC 小两类结果差 1 到 2 属于正常。单独看数值大小没有意义BIC 是相对量只在同一份数据的 K 之间比较才有意义。这里有个细节容易踩每个 K 都要开 Replicates否则某个初始化运气差K4 的 BIC 可能虚高曲线出现假拐点。我一般扫 K 时 Replicates 开到 3 到 5BIC 曲线稳定后再对候选 K 做一次大复跑确认。4.2 协方差类型full、diagonal 与共享协方差怎么选fitgmdist 的 CovarianceType 控制每个分量的协方差是 full 还是 diagonalSharedCovariance 控制所有分量是否共用同一个协方差矩阵。三种组合的自由参数规模和适用场景差别很大参数设置协方差自由参数规模适用场景典型风险CovarianceTypefullK·d(d1)/2各簇形状、取向完全不同小样本易奇异、易过拟合CovarianceTypediagonalK·d特征近似独立或维度高样本多丢失簇内特征相关性SharedCovariancetruefulld(d1)/2全体共用各簇大小形态接近只有位置不同假设太强容易欠拟合% 三种配置一次都跑出来, 对比它们的 BIC 和聚类效果 gm_diag fitgmdist(X, 3, CovarianceType, diagonal, Replicates, 5); gm_tied fitgmdist(X, 3, SharedCovariance, true, Replicates, 5); gm_full fitgmdist(X, 3, CovarianceType, full, Replicates, 5); % 用各自的 BIC 判断哪个结构更适合这份数据 [gm_diag.BIC, gm_tied.BIC, gm_full.BIC]在第 3 章的合成数据上full 的 BIC 应该明显低于 diagonal因为簇 1 和簇 3 都有非对角协方差diagonal 把相关性全部抹掉了。反过来如果你处理的是高维稀疏特征比如 200 维的词频向量d(d1)/2 高达两万哪怕样本量十万也拟合不稳这时 diagonal 就是唯一现实的选择。共享协方差是这三者里最容易被忽视的当你有先验知识说这几个簇只是位置不同、形状相同时它能用最少参数拿到接近 full 的效果而且永远不会奇异。判断用哪种协方差还有一个实用信号full 训练报奇异错误而 diagonal 不报说明问题出在协方差自由度太高先切 diagonal 稳一手再用第 4.3 节的正则处理。4.3 RegularizationValue给协方差上保险的正确姿势% 同样的 K5, K 偏大, full 协方差容易退化 gm_no fitgmdist(X, 5, Replicates, 3); % 可能直接报错 gm_ok fitgmdist(X, 5, RegularizationValue, 1e-4, Replicates, 3);RegularizationValue 的作用是在每次 M 步更新协方差时给对角线注入一个常数让矩阵强制正定。对应到手写代码里就是 M 步里那句Sigma 1e-6 * eye(d)。它解决的是一类很典型的问题某个分量有效样本太少或特征之间存在近似线性关系导致协方差矩阵奇异或半正定接下来的 mvnpdf 就全崩了。这个参数值的量级是经验活多少带点玄学成分。太小1e-8挡不住奇异太大1e-1会把所有分量的协方差都拉向单位阵聚类边界被压成规则形状彻底失真。我一般从 1e-6 开始按 10 倍步长试到模型不再报错的最小值为止常见落脚点在 1e-5 到 1e-3。注意正则项会真实抬高 BIC因为它改变了似然所以做 BIC 选 K 的时候对每个 K 用同一个正则值横向比较才公平。5. GMM 调试避坑五个让 log-likelihood 翻车的典型问题5.1 报错 SIGMA must be symmetric positive definite现象手写版在迭代过程中断报 Sigma 非正定fitgmdist 干脆拒绝训练或者跑完后查协方差矩阵有负特征值log-likelihood 变成 NaN。原因排查下来九成是两类情况。一是 K 设得比真实结构多某个分量在 M 步只拿到两三个样本协方差由两三个点估计出来必然奇异二是原始特征里有近似线性相关的列比如一列是另一列的 2 倍协方差矩阵秩亏。解决先加 jitter手写版在 M 步加 1e-6 单位阵fitgmdist 用 RegularizationValue如果加了正则还报错就降 K 或降维。我的经验是加了正则后一半问题直接消失另一半是你把 K 设过头了先跑一遍 BIC 再回来调。5.2 log-likelihood 不升反降甚至出现 NaN现象打印 llh 序列不是理论上的单调递增中间某轮突然掉下去或者直接出现 NaN。原因这是典型的数值下溢。E 步里每个样本都要算 K 个高斯密度当样本离某个高斯中心有几十个标准差远时mvnpdf 返回 0分母归一化时出现 0/0log 里吃到 0 就出 -Inf。数据没做标准化时最爱出这种事。解决标准做法是在 E 步换成对数域计算用 log-sum-exp 技巧。先写一个数值稳定的对数密度函数再替换原来的 E 步function logp logmvnpdf(X, mu, Sigma) % 数值稳定的多元高斯对数密度 % 用 Cholesky 分解避免显式求逆, 离均值很远也不会下溢为 0 [d, ~] size(Sigma); [L, p] chol(Sigma); if p 0 error(Sigma 非正定, 请检查正则项); end Xc X - mu; % n x d 中心化矩阵 U Xc / L; % 三角回代, 等价于 Xc * inv(L) logp -0.5 * d * log(2*pi) - sum(log(diag(L))) - 0.5 * sum(U.^2, 2); end对应的稳定 E 步是% 数值稳定版 E 步: 先算对数再减去行最大值, 防止指数爆炸 loggamma zeros(n, K); for k 1:K loggamma(:, k) log(w(k)) logmvnpdf(X, mu(k, :), Sigma(:, :, k)); end m max(loggamma, [], 2); % 每行最大值 gamma exp(loggamma - m); % 减去最大值后再 exp gamma gamma ./ sum(gamma, 2);logmvnpdf 用 Cholesky 分解有三个好处不用显式求逆更快也更稳sum(log(diag(L)))一步拿到 log|Sigma| 的一半检测到非正定矩阵时立刻报错让问题暴露在训练早期而不是等到 NaN 扩散到全是再回头查。5.3 空分量EM 把某个高斯饿死了现象训练完打印 gm.ComponentProportion发现某个权重是 1e-10 量级对应分量的均值漂移到某个数据孤岛或者和其他分量的均值几乎重合。原因EM 迭代在最大化似然的过程中会把责任度太小的分量权重不断压低直到归零。这个行为在数学上没毛病但如果你业务上明确要求 K 个分量比如要做 5 类图像分割空分量就是问题。常见诱因是 K 设得偏大或者初始化阶段起点太差。解决第一选择是回到 BIC 确认这个 K 是否合理第二是增加 Replicates 次数让 fitgmdist 从多个起点里挑似然最高的结果第三手写版可以加一个保护逻辑每轮检查 Nk如果某个分量有效样本数小于 n 的千分之一就把它重新初始化到当前 log-likelihood 最低的几个样本附近恢复它的权重。这相当于人为给空分量续命。5.4 聚类结果每次运行都不一样现象同一份数据不固定 rng跑两次得到的 mu 和 label 完全不同甚至分量顺序互相交换。原因两个问题叠在一起。第一EM 是局部优化算法高斯混合的似然面存在多个局部最优不同初始化会收敛到不同结果。第二即便似然相同分量编号本身可以任意置换同一个模型把簇 1 和簇 2 的编号换一下是数学上完全等价的另一个解这叫标签置换。解决可复现性靠固定种子把 rng(seed) 存成脚本最顶部的常量每个实验一个唯一 seed模型质量靠 Replicatesfitgmdist 会自动保留对数似然最高的那次结果跨模型比较分量时不能直接用第 k 个分量要先按均值大小排序或者做一次标签匹配再对齐。5.5 高维数据跑不动还频繁奇异现象d 超过 30 以后手写 EM 一轮迭代肉眼可见地慢fitgmdist 也经常在训练中途报错。原因GMM 的参数规模和维度是平方关系full 协方差在 d50 时每个分量要估 1275 个参数几千样本根本喂不饱mvnpdf 每轮要对每个分量做一次 d×d 分解计算量随维度平方增长。解决先对数据做 PCA 或标准化降到 5 到 20 维再上 GMM维度实在要保留时先用 diagonal 协方差跑一轮确定合理 K再在低维空间换 full如果只是做聚类不需要密度估计直接降级到 K-means别用 GMM 硬扛。这个坑在高维文本向量和基因表达数据里出现频率最高属于走完弯路才想明白的问题。6. 进阶GMM 做异常检测与像素级图像分割的两个落地技巧6.1 用密度阈值做异常检测GMM 训练完成后pdf 函数直接给出每个样本的混合密度。正常数据的密度集中在某个区间异常点通常落在密度极低的区域。关键是阈值怎么定不要看密度绝对值它对维度和数据尺度很敏感要看训练集密度的分位数。% 只用正常样本训练 gm fitgmdist(X_normal, 4, RegularizationValue, 1e-5, Replicates, 5); % 训练集和测试集各自的密度 p_train pdf(gm, X_normal); p_test pdf(gm, X_test); % 阈值取训练集密度的 1% 分位数, 低于它就是异常 thresh prctile(p_train, 1); isAnomaly p_test thresh;这个阈值的意思很直观训练集里最不像正常的 1% 样本被划为异常边界。和单纯算马氏距离比GMM 密度考虑了每个簇的形状和权重能处理多模态的正常数据。一份正常的工业传感器数据常常本身就有三四个工作状态K-means 或单高斯都搞不定GMM 可以。6.2 像素级图像分割RGB 三通道直接喂给 GMM图像分割是 GMM 聚类最直观的应用。把每个像素当成一个三维特征向量聚出来的每个簇就是一类颜色区域。img imread(demo.png); [h, w, c] size(img); pixels double(reshape(img, h*w, c)) / 255; % 归一化到 [0,1] gm fitgmdist(pixels, 5, ... CovarianceType, diagonal, ... % RGB 通道近似独立, 速度更快 RegularizationValue, 1e-4, ... Replicates, 2); label cluster(gm, pixels); seg reshape(label, h, w); figure; imagesc(seg); axis image;这里用 diagonal 协方差就够了RGB 三个通道之间的相关性在自然图像里并不强而且像素量动辄几十万full 协方差的迭代成本高到不划算。归一化到 [0,1] 是必须的不然 RGB 的数值量级差异会让协方差被某一个通道主导。如果图像很大可以先对像素做随机抽样训练再把全部像素喂给 cluster这样训练时间能砍掉大半。从那以后我每次在项目里用 GMM 之前都会强制走一遍固定随机种子标准化数据扫 BIC 定 K最后设好正则再开跑。这套流程已经把我在 EM 上踩过的坑都挡在了门外。直接下载这份代码包把 gmm_em.m、logmvnpdf.m、plot_ellipse.m 和示例脚本从头跑一遍比只看公式更能建立手感。希望帮到你。本文还有配套的精品资源点击获取