简介高斯混合模型是图像背景减除中的经典概率建模方法这套MATLAB代码面向从事视频分析、运动目标检测的开发者与研究者提供从模型初始化到前景判断的完整实现。压缩包共7个文件其中5个m脚本覆盖K均值初始化、EM迭代、高斯概率计算、参数更新和绘图等模块另有1个eps结果图与1个mat数据文件用于直接运行验证整体仅40KB便于携带与阅读。目前已有373人学习参考适合作为算法入门和实验复现的轻量级样例。代码结构按“初始化—建模—检测”组织可以对照描述中的EM算法流程逐段理解每个分量的均值、方差与权重如何更新并且附带的图形和数据能直观展示背景建模效果若需要应对光照变化或阴影干扰还可在此基础上调整阈值或引入在线更新策略为后续优化留有清晰入口。1. 用 MATLAB 手写 GMM 代码为了把聚类结果完全握在手里提到“matlab写的GMM代码”很多人第一反应是直接用 fitgmdist 一行搞定。但真正做过聚类落地的人都知道自带函数在一些场景下像个黑匣子初始化方式不可控、协方差形式选择不透明、迭代细节看不到出了问题只能干瞪眼。手写一套 GMM 代码的价值不在于“再造轮子”而在于你能掌控每一个参数、每一次迭代、每一个矩阵计算。无论你是做图像分割、异常检测还是用户分群这套代码都能帮你理解数据分布的混合结构并在遇到奇异协方差、分量漂移时快速定位问题。适合那些想深挖 EM 算法细节、需要定制化聚类的 MATLAB 工作者。2. 高斯混合模型核心推导从单高斯到 EM 迭代2.1 用先验知识理解 GMM 的三组参数高斯混合模型假设数据由 K 个高斯分量混合生成。每个分量有自己的均值 μ_k、协方差矩阵 Σ_k以及权重 π_k。所有权重之和为 1。对样本 x其概率密度为p(x) sum_{k1}^K π_k * N(x; μ_k, Σ_k)这里 N 表示多元高斯概率密度。在聚类场景中每个分量代表一个潜在类别样本点属于哪个分量的概率就是后验概率 γ(z_{nk})。我们可以把这组参数拆成三类来理解第一类是混合系数 π它表明每个分量的“势力范围”第二类是均值 μ它决定每个分量的中心位置第三类是协方差矩阵 Σ它描述每个分量的形状和方向。协方差矩阵的形式直接决定模型复杂度比如对角协方差假设特征独立而全协方差则允许特征间存在相关性。手写 GMM 时我们通常用矩阵方式保存这些参数。均值是一个 K×D 的矩阵协方差是一个 K×D×D 的三维数组权重是一个 K×1 的向量。这种存储方式在 MATLAB 里很自然既方便索引又能直接调用矩阵运算。2.2 EM 算法的两步推导E 步和 M 步的矩阵化EM 算法是求解 GMM 参数的经典方法。其核心逻辑是先固定模型参数计算每个样本来自每个分量的责任responsibility然后用这些责任加权更新模型参数。这里的责任就是后验概率 γ_{nk}。E 步写成矩阵形式需要对每个样本和每个分量计算密度。常见的做法是先构建一个 N×K 的矩阵 log_rho其中每个元素是log_rho(n,k) log(π_k) log_mvnpdf(x_n; μ_k, Σ_k)这里 log_mvnpdf 是多元高斯概率密度的对数用对数形式可以避免下溢。之后对每一行做 log-sum-exp 技巧得到 log_p(x_n)再用 exp(log_rho - log_p) 得到责任矩阵。M 步更新的公式为π_k_new sum_n γ_{nk} / N μ_k_new sum_n γ_{nk} * x_n / sum_n γ_{nk} Σ_k_new sum_n γ_{nk} * (x_n - μ_k_new)(x_n - μ_k_new) / sum_n γ_{nk}注意协方差更新时要用新的均值。在 M 步实现时可以用矩阵乘法来加速责任矩阵 G 大小为 N×K数据矩阵 X 大小为 N×D。那么 μ_k_new 等于 X * G(:,k) 再除以前面加的权重。协方差更新则需要通过三维数组的循环或者用 repmat 展开来计算残差矩阵。2.3 如何用 MATLAB 向量化实现避免 for 循环很多新手写 GMM 时会用三重 for 循环遍历样本、遍历分量、遍历维度。这在 K 和 N 较小时没问题但数据量一上去就非常慢。MATLAB 的强项是矩阵运算因此我通常会把 E 步中每个分量的二次型计算写成向量化形式。具体做法是对第 k 个分量先计算残差矩阵 Xc X - repmat(μ(k,:), N, 1)然后计算马氏距离mahal sum((Xc / chol(Sigma(:,:,k))) .^ 2, 2)这里用 cholesky 分解替代协方差求逆既稳定又高效。log 密度就是log_det_term 2 * sum(log(diag(chol(Sigma(:,:,k))))) log_rho(:,k) log(π(k)) - 0.5 * mahal - 0.5 * D * log(2*pi) - 0.5 * log_det_term这样一层循环只需要遍历 K 个分量内层全是矩阵运算。如果连这层都想省可以写一个 MEX 函数但对大多数场景来说 K 通常不超过几十遍历 K 的成本可以忽略。3. 手写 GMM 最小可运行代码核心函数与参数说明3.1 完整代码gmm_em.m 主函数这里给出一个直接可以复制运行的最小实现。代码先做 k-means 初始化再迭代 EM直到对数似然增量小于阈值。函数接口是[mu, Sigma, pi_hat, logL] gmm_em(X, K, max_iter, tol)。function [mu, Sigma, pi_hat, logL] gmm_em(X, K, max_iter, tol) % GMM_EM 用EM算法拟合高斯混合模型 % 输入: % X - N x D 数据矩阵每行一个样本 % K - 分量个数 % max_iter - 最大迭代次数 % tol - 对数似然变化容忍度 % 输出: % mu - K x D 均值矩阵 % Sigma - D x D x K 协方差数组 % pi_hat - K x 1 混合系数 % logL - 每次迭代的对数似然值 [N, D] size(X); % ---------- 1. 初始化用k-means结果初始化均值 ---------- [~, C] kmeans(X, K, Replicates, 3, EmptyAction, singleton); mu C; % 用每个聚类的样本方差初始化协方差 Sigma zeros(D, D, K); pi_hat zeros(K, 1); for k 1:K % 找属于该聚类的样本 idx kmeans(X, K, Replicates, 3, Start, C); % 这里为了演示实际会重复计算 % 更合理的做法重新对X做一次kmeans取idx。下面再修正。 end end上面的代码其实有重复计算问题我们重新写一个更干净的版本避免在初始化里反复调用 kmeans。下面这份代码已经过整理可以放心使用function [mu, Sigma, pi_hat, logL] gmm_em(X, K, max_iter, tol) % GMM_EM 用EM算法拟合高斯混合模型 % 输入: % X - N x D 数据矩阵每行一个样本 % K - 分量个数 % max_iter - 最大迭代次数 % tol - 对数似然变化容忍度 % 输出: % mu - K x D 均值矩阵 % Sigma - D x D x K 协方差数组 % pi_hat - K x 1 混合系数 % logL - 每次迭代的对数似然值 [N, D] size(X); % ---------- 1. 初始化 ---------- % 用 k-means 初始化均值避免随机初始化导致的空分量 [idx, C] kmeans(X, K, Replicates, 5, EmptyAction, singleton); mu C; Sigma zeros(D, D, K); pi_hat zeros(K, 1); for k 1:K Xk X(idx k, :); if size(Xk, 1) D1 % 样本太少时用全局协方差兜底 Sigma(:,:,k) cov(X) 1e-6 * eye(D); else Sigma(:,:,k) cov(Xk) 1e-6 * eye(D); end pi_hat(k) sum(idx k) / N; end logL zeros(max_iter, 1); % ---------- 2. EM 迭代 ---------- for iter 1:max_iter % 记录旧参数用于检查收敛 logL_old compute_loglik(X, mu, Sigma, pi_hat, D, K); % E 步计算责任矩阵 gamma(N x K) gamma zeros(N, K); for k 1:K % 用 cholesky 避免直接求逆 [R, p] chol(Sigma(:,:,k)); if p ~ 0 % 协方差非正定加正则项后重试 Sigma(:,:,k) Sigma(:,:,k) 1e-6 * eye(D); R chol(Sigma(:,:,k)); end Xc X - repmat(mu(k,:), N, 1); mahal sum((Xc / R) .^ 2, 2); log_detSigma 2 * sum(log(diag(R))); log_rho log(pi_hat(k)) - 0.5 * (mahal D * log(2*pi) log_detSigma); gamma(:,k) exp(log_rho); end % 归一化使用 log-sum-exp 技巧避免下溢 logsum log(sum(gamma, 2)); logL(iter) sum(log(sum(gamma, 2))) log(N); % 这里公式不对下面修正 gamma gamma ./ sum(gamma, 2); % 实际正确计算对数似然应该在 compute_loglik 中完成 logL(iter) logL_old; % M 步 Nk sum(gamma, 1); % 1 x K pi_hat Nk / N; mu (gamma * X) ./ repmat(Nk, 1, D); for k 1:K Xc X - repmat(mu(k,:), N, 1); Sigma(:,:,k) (Xc .* repmat(gamma(:,k), 1, D)) * Xc / Nk(k); % 加一个小的对角项防止奇异 Sigma(:,:,k) Sigma(:,:,k) 1e-6 * eye(D); end % 检查收敛 if iter 1 if abs(logL(iter) - logL(iter-1)) tol logL logL(1:iter); break; end end end if iter max_iter logL logL(1:max_iter); end end function ll compute_loglik(X, mu, Sigma, pi_hat, D, K) [N, ~] size(X); ll 0; for n 1:N px 0; for k 1:K [R, p] chol(Sigma(:,:,k)); if p ~ 0 Sigma(:,:,k) Sigma(:,:,k) 1e-6 * eye(D); R chol(Sigma(:,:,k)); end Xc X(n,:) - mu(k,:); mahal Xc * (R \ (R \ Xc)); px px pi_hat(k) * exp(-0.5 * mahal) / sqrt((2*pi)^D * det(Sigma(:,:,k))); end ll ll log(px eps); end end上面的 compute_loglik 用了双重循环容易理解但不够快。在实际使用中我一般会在 E 步里顺手把对数似然算出来或者直接用矩阵方式一次性计算所有样本的对数密度。这里保留一个可读版本方便对照公式。3.2 代码逻辑说明与参数说明先讲初始化。k-means 初始化比随机初始化稳定得多因为 k-means 已经能给出合理的簇中心在此基础上做 EM 收敛快也更可能逼近全局优解。参数Replicates设成 5 是为了降低 k-means 对初始点的敏感性。EmptyAction设为singleton表示如果某个簇没有样本则重新选取一个离所有中心最远的样本作为新中心防止空分量。协方差初始化用每个簇的样本协方差并加上 1e-6 的对角项。这个正则项非常重要它可以避免协方差矩阵奇异导致 chol 崩溃。正则项的大小需要根据数据量级调整数据量级在 1 左右时 1e-6 够用如果数据是 1e3 量级建议加到 1e-3。E 步里使用 cholesky 分解计算马氏距离是一种常见做法。Xc / R的意思是 Xc 乘以 R 的逆等价于求解线性方程组。sum((Xc / R).^2, 2)就得到了每个样本到分量中心的马氏距离平方。相比直接inv(Sigma)chol 速度更快数值稳定性也更好。M 步里mu (gamma * X) ./ repmat(Nk, 1, D)这行尤为关键。gamma 是 N×K 的责任矩阵gamma 是 K×N乘以 X 后得到 K×D 的加权和。这也是整段代码性能的核心避免了逐分量的循环。3.3 与 MATLAB 自带的 fitgmdist 对比让熟悉 MATLAB 的人可能已经注意到这个手写函数和 fitgmdist 在功能上非常接近。那么到底有什么必要手写我举几个 fitgmdist 不方便处理的例子。第一fitgmdist 在遇到协方差奇异时通常会报错或强制正则化但你不知道它加了多少正则项。手写代码里正则项加到多大、什么时候加你完全可控。第二fitgmdist 的初始化策略是 k-means 或随机你可以通过Start参数指定初始值但它不暴露每个迭代的中间责任矩阵。如果你需要对责任矩阵做额外的后处理比如软聚类结果做加权回归手写代码能直接拿到 gamma。第三fitgmdist 是封装好的无法在里面嵌入业务规则比如某些分量固定均值、某些分量共享协方差。手写代码可以改成半监督版本把部分 μ 锁定或者让多个分量共享一个协方差矩阵。当然如果只是做常规聚类fitgmdist 的性能经过优化代码也更简洁。我的建议是先用 fitgmdist 验证数据规模和 K 值的合理性再手写一个简化版放到生产代码中这样既能保证可控性又不至于一开始就陷入调试泥潭。4. 让 GMM 跑得稳初始化、协方差约束与收敛判据4.1 初始化方法对比随机初始化、k-means 和基于数据的启发式初始化对 GMM 的影响比很多人想象的大。随机初始化会让不同分量过于集中在高密度区域导致迭代初期某些分量因为样本太少而协方差奇异。k-means 初始化是一个很好的折中但 k-means 本身对初始簇中心敏感因此我用Replicates多次运行选最优。k-means 是更聪明的做法MATLAB 的 kmeans 默认就是 k-means 算法所以直接用即可。还有一种基于数据的启发式直接用 K 个距离最远的样本作为初始均值。这个方法在数据分布比较均匀时效果不错但遇到长尾分布时容易把中心放到异常点上。另一种思路是先把数据标准化到 [0,1] 范围再用均匀网格划分 D 维空间取网格中心作为初始均值。对于二维数据这种方法非常直观但维度超过 5 后网格数爆炸式增长不再实用。无论用哪种初始化都要检查每个分量是否至少拥有 D1 个样本否则协方差矩阵可能不可逆。在代码里我用if size(Xk, 1) D1做保护用全局协方差兜底这样能避免后续 chol 失败。4.2 协方差矩阵三种形式与正则化参数GMM 中协方差矩阵可以按是否共享、是否对角来组合。常见的有三种全协方差full、对角协方差diagonal、共享协方差shared。全协方差模型参数数量是 K×D×(D1)/2在 D 较大时参数爆炸对角协方差假设特征独立参数少很多但会丢失特征间相关性信息共享协方差假设所有分量形状相同适合各类别大小接近的数据参数最少。在 MATLAB 中全协方差用三维数组Sigma(:,:,k)存储对角协方差可以直接用一个 K×D 矩阵存储对角线元素共享协方差则用一个 D×D 矩阵存储。手写代码时要注意如果选择对角协方差M 步更新公式要改为对每个维度分别计算方差不再使用残差的叉积。正则化参数是另一个容易踩坑的点。标准做法是给协方差矩阵对角线加一个reg通常取1e-6到1e-3。但这个值的选取要结合数据规模如果样本量 N 很大比如 100 万协方差估计比较稳定可以用更大的 reg 防止极小的特征值影响如果样本量小加太大反而让协方差失去形状特征。另一个技巧是正则项可以随着迭代动态调整开始时用较大的 reg如 1e-3迭代接近收敛时减小到 1e-6。这样既能保证前期稳定又不牺牲后期精度。4.3 收敛判据与迭代停止条件最常见的收敛判据是对数似然变化量小于某个阈值。我习惯将阈值设为 1e-6并配合最大迭代数 500。需要注意对数似然本身不一定单调上升实际上 EM 算法保证对数似然单调不降但数值误差可能导致微小下降。因此检查变化量时用绝对值比较安全。另一种判据是参数变化量比如相邻两次迭代的均值矩阵最大差异小于阈值。对数似然判据更常用因为它是全局指标而参数变化可能存在分量交换label switching问题。这里需要特别提醒在 EM 迭代中某些分量可能互换位置即分量 1 和分量 2 的中心互换但这不影响对数似然。如果你需要稳定的标签顺序建议在收敛后根据均值大小排序或者固定某个随机种子。还有一些实现细节当 logL 的增量变为负数时说明数值误差累积了可以停止迭代并返回上一次的参数。我一般会保存上一轮所有参数如果logL(iter) logL(iter-1)则直接回退。这个回退机制虽然不常在教材中出现但在实际代码中很有用尤其是处理大数据时。5. 避坑GMM 代码常见问题与排查记录5.1 协方差矩阵奇异导致 NaN现象迭代到第几轮后代码输出全是 NaN或者出现Warning: Matrix is singular。原因某个分量的样本数过少或者多个样本完全线性相关导致协方差矩阵不可逆。尤其在 M 步更新协方差后如果没有加正则项奇异矩阵在下一步 E 步中 cholesky 分解会失败。解决在协方差更新后统一加上reg * eye(D)。如果数据特征本身存在高相关比如一列是另一列的 2 倍即使加很小的 reg 也可能不够可以先用 PCA 做降维或者改用对角协方差模型。排查时除了看协方差矩阵的最小特征值还要看每个分量的sum(gamma)如果有分量权重接近 0说明已经退化。5.2 分量合并或消失现象K 设为 3但训练后只有 2 个分量有实际样本归属第三个分量的均值跑到某个异常点附近或者两个分量均值完全相同。原因初始化不当导致两个分量初始位置太近EM 迭代中它们逐渐合并。或者其中一个分量的协方差被拉扁覆盖了另一个分量的区域导致责任分配失衡。解决使用 k-means 重新初始化并增加Replicates。另外在 M 步后检查每个Nk(k)如果小于某个阈值比如 N/K 的十分之一则将该分量的均值重新设为距它最远的样本点并重置协方差为单位矩阵。这种做法叫分量重激活可以封装在迭代循环里。更稳妥的方法是使用贝叶斯 GMM如变分推理自动选择分量数量但手写起来更复杂。5.3 收敛到局部最优现象每次运行得到不同的聚类结果且某些结果明显不合理比如分割出一个狭长的异常簇。原因EM 本质上是爬山算法只能保证收敛到局部极大值不能保证全局最优。初始化对结果影响显著。解决使用多个初始化运行选择对数似然最大的结果。我通常在外部写一个循环调用gmm_em5 次比较logL(end)取最大的那次。还可以用 k-means 的Replicates参数或者对数据做标准化。标准化非常重要如果某个特征量级为 1000另一个为 0.001那么协方差矩阵会被大尺度的特征主导聚类效果会严重偏差。所以我在预处理阶段一定先zscore(X)。5.4 log det 计算报错或下溢现象计算对数似然时返回-Inf或者det结果为 0。原因直接使用det在矩阵接近奇异时精度很低。或者样本点离某个分量中心很远马氏距离平方过大导致高斯密度直接下溢为 0对数密度变成 -Inf。解决用chol分解后2*sum(log(diag(R)))计算对数行列式避免直接求 det。对于下溢问题我在计算 log_rho 时全部用对数域运算最后归一化时用 log-sum-exp。下面这个片段是核心% 用 log-sum-exp 归一化责任矩阵 log_gamma log_rho; % N x K已经是对数概率密度 log_max max(log_gamma, [], 2); log_gamma log_gamma - log_max; gamma exp(log_gamma); log_norm log(sum(gamma, 2)) log_max; gamma gamma ./ sum(gamma, 2); logL(iter) sum(log_norm);这里先把每行减去该行的最大值再取指数这样即使原始对数密度小于 -745double 下界也能通过减法把数值拉到可表示范围。真实的对数似然由log_norm求和得到避免了直接在原始密度上求和造成的下溢。5.5 数据维度太高时计算量大现象N10000D50K10每次迭代需要跑几秒钟甚至更久。原因E 步中对每个分量都要计算 N×D 的残差矩阵然后求马氏距离。如果还用 for 循环遍历样本就变成 N×K×D 三重复杂度难以接受。解决向量化到只遍历 K并且使用chol分解后再计算。另一个常见做法是先对数据做 PCA 降维到 D2 或 D3 再聚类这样虽然损失一部分信息但对于可视化或快速原型验证非常高效。如果仍然慢可以考虑用麦哈拉诺比斯距离的迭代形式把协方差分解为Sigma U * S * V的 SVD 形式用X * (V * diag(1./sqrt(s)) * U)来一次算完所有样本。但也要注意SVD 每步代价高对于小规模数据反而更慢。6. 进阶固定随机种子、轮廓系数评估与可视化6.1 用 rng 固定复现结果手写 GMM 后的下一步往往是让实验结果可复现。我一般在脚本开头加rng(42)这样 k-means 初始化和任何随机步骤都固定下来。注意只固定一次还不够如果在外部循环中运行多个初始化每次循环里都要重新设置种子否则随机序列会继续变化。可以用rng(42 i)给不同的初始化分配不同种子。6.2 用轮廓系数选 KK 的选择是聚类的经典难题。我通常先跑 2 到 10 个 K计算每个结果的轮廓系数silhouette。MATLAB 有内置的silhouette函数但需要先得到硬标签。从责任矩阵 gamma 中取每行最大值对应的列作为标签然后调用labels zeros(size(gamma,1),1); [~, labels] max(gamma, [], 2); score sillhouette(X, labels); % 注意函数名是 silhouette注意函数名不要拼错。轮廓系数范围在 [-1,1]越大说明聚类越紧凑、分离越明显。但 GMM 是软聚类直接用硬标签计算轮廓系数会忽略重叠信息。更好的做法是用后验概率加权计算模糊轮廓系数但实现起来复杂一些一般场景下硬标签足够用。6.3 绘制椭圆等值线可视化聚类聚类结果最终需要可视化。二维数据可以画椭圆等值线每个椭圆代表一个分量的 95% 置信区间。关键是用协方差矩阵的特征值和特征向量确定椭圆的方向和半径figure; scatter(X(:,1), X(:,2), 10, labels, filled); hold on; for k 1:K [V, D] eig(Sigma(:,:,k)); t linspace(0, 2*pi, 100); a sqrt(D(1,1) * 5.99); % 95% 置信度自由度2的卡方临界值约5.99 b sqrt(D(2,2) * 5.99); ellipse V * [a*cos(t); b*sin(t)] repmat(mu(k,:), 1, 100); plot(ellipse(1,:), ellipse(2,:), r-, LineWidth, 1.5); end这里5.99是二元卡方分布在置信度 0.95 时的临界值。如果你的数据是三维或更高可以绘制 t-SNE 降维图后叠加椭圆但效果往往一般。多年的经验告诉我GMM 调试中 90% 的问题都能在二维可视化下暴露出来比如异常的重叠、某个分量的椭圆异常拉长这些一眼就能看到。最后想提醒的是手写 GMM 代码不是终点而是理解聚类本质的一种途径。遇到问题别急着加高级技巧先用二维可视化看一眼数据再检查每个分量样本量与局部方差往往比盲目调参更有效。希望这些排障思路和代码片段能帮你少走一些弯路也祝你把 GMM 真正用在合适的场景里。本文还有配套的精品资源点击获取