EM算法详解：从隐变量到GMM聚类的参数估计实战

1. EM算法到底在解决什么问题EM算法全称Expectation-Maximization中文叫期望最大化算法。初次接触这玩意儿的人十有八九都会被那一堆公式劝退什么Q函数、琴生不等式、隐变量……看着就像天书。但我想先说句大实话EM算法的思想内核极其简单简单到可以用一句话概括——先猜一个结果根据结果反推原因再用反推出来的原因去修正结果然后循环。是不是很像生活中我们干的事举个例子。你新接手了一个团队手下有几个员工你不太清楚每个人的真实水平。第一周你观察他们的产出摸个大概算是初始猜测。然后根据产出反推谁强谁弱再给不同的人分配不同难度的任务。过两周再看产出修正你的判断。循环几次你对每个人的能力评估就越来越准任务分配也越来越合理。EM算法干的事情跟这个如出一辙。它解决的问题是在存在隐变量的情况下估计概率模型的参数。所谓隐变量就是那些你观测不到、但切实影响着观测结果的变量。我在实际工作中接触到的场景典型的像用户聚类你能看到用户的购买行为但看不到用户属于哪个潜在人群。这个人群归属就是隐变量。又比如做图像分割你能看到像素点的颜色但看不到这个像素属于前景还是背景这也是隐变量。再比如语音识别里的HMM模型、文本主题挖掘里的LDA模型核心都有隐变量的影子。EM算法就是为了这类问题而生的。它能在你完全不知道隐变量的取值时仅凭观测数据硬生生地把模型参数估计出来。那没有隐变量的情况呢直接用极大似然估计就好了跟EM没什么关系。这一点必须先分清楚很多人学EM学得云里雾里就是因为没搞明白隐变量这三个字的分量。这篇文章的定位是通俗但不过分注水的详细解读。我会先抛一个经典到不能再经典的例子帮你把直觉建立起来然后从数学上把E步和M步一步步拆开讲清楚每一步在干什么、为什么这么干接着用高斯混合模型GMM做完整推导给出一个跑通的全过程最后整理一些工程实战中躲不开的坑和调试经验。适合谁看正在学机器学习、统计模型的同学以及工作中要用到聚类、缺失数据填补、混合分布的建模工程师。看这篇文章之前你只要知道极大似然估计的基本概念就行公式部分我会尽量解释到位。2. 从三硬币问题建立直觉2.1 三硬币问题是什么先看统计学习里那个著名的例子。假设有A、B、C三枚硬币正面朝上的概率分别为p、q、r。实验过程分两步先掷硬币A如果正面朝上就掷硬币B并记录结果如果反面朝上就掷硬币C并记录结果。重复做n次实验你只能看到第二步的结果正或反看不到A掷出了哪一面。这个过程中A的结果就是隐变量——你观测不到但它决定了后续你用的是B还是C。如果隐变量可见问题好办得多。你能统计出使用B时得到正面的次数和使用B的总次数比值就是q的估计p和r同理。这种分组统计的思路就是极大似然估计的直接解法。但问题恰恰在于你根本看不到哪些数据来自B、哪些来自C。数据混在一起正面朝上的概率是qr(1-q)p这样一个混合表达式。你要在一个混合表达式里同时估计p、q、r三个参数直接求导会发现方程很难解而且没有解析解。这就是隐变量带来的核心困境。2.2 EM的应对策略EM的思路很巧妙既然不知道数据来自B还是C那就先假定一组参数比如随便猜p0.5q0.5r0.5然后根据这组参数去计算每一个样本来自B的概率和来自C的概率。注意这个词概率。不是硬性地把样本归给B或C而是软性地分配。比如某个样本是正面在当前的参数下它有0.7的可能是B掷出来的有0.3的可能是C掷出来的。有了这个软分配接下来就好办了。我们拿着这些带权重的样本去重新统计B和C各自出现正面的加权频率得到一组新的参数估计。然后是新一轮循环用新参数重新计算每个样本的归属概率再重新估计参数……就这样一轮接一轮参数会逐渐逼近一个稳定值。最终的p、q、r就是EM算法给出的答案。这里面有个非常关键的心理转变EM并没有试图一次性解决所有问题而是把一个难啃的骨头分解成猜归属和估参数两个交替进行的简单环节。每个环节单独拿出来都不难难的是原来那个混合在一起的联合估计问题。这种分而治之的思路在工程里太常见了。2.3 为什么这种交替能奏效我当年第一次接触EM时最大的疑问是这轮流猜来猜去凭什么保证最后能收敛会不会越猜越偏答案在于一个保证不下滑的机制。EM的每一轮迭代都会保证观测数据的似然函数值不下降或者说只会上升或持平。这就好比爬山每一步都确保自己不比刚才的位置低虽然不保证一步跨到山顶但它保证永远在往上走或原地踏步。有这条性质兜底反复迭代下去自然能到达一个局部最高点。要理解这个性质的数学证明得引入一些符号和不等式这部分放在下一节详细展开。这里先记住结论EM是单调上升的最后会收敛到一个局部最优解。注意我说的是局部最优。EM算法本身不保证找到全局最优。初始参数选得不好很可能收敛到某个不太好的局部极值点。这个特性在工程上坑过很多人后面我会专门讲应对方法。3. 数学拆解E步和M步到底在算什么3.1 极大似然为什么失效先回顾一下标准极大似然估计。假设观测数据是X参数是θ我们要最大化的是对数似然函数log P(X|θ)如果没有隐变量这个函数可以直接对θ求导令导数为零解方程完事。但有了隐变量Z之后事情变了。我们不知道Z的具体值只能把Z所有可能的情况加总起来log P(X|θ) log [Σ_Z P(X, Z|θ)]这里的Σ_Z是个求和或积分作用在P(X,Z|θ)上然后再取对数。问题就出在这个先求和再取对数上。对数里面有个求和导致我们对θ求导时分母会非常复杂而且耦合在一起很难分离出各个参数的解析更新式。如果能把取对数挪到求和里面就好了。比如log P(X,Z|θ)这种形式对θ求导就容易得多。EM的核心操作本质上就是想办法构造这样一个好处理的下界。3.2 琴生不等式与ELBO证据下界这里需要引入一个工具琴生不等式。它的含义用大白话说就是对于一个凹函数f有f( E[Y] ) ≥ E[ f(Y) ]对数函数是凹函数所以如果我们对一个正随机变量取对数这个不等式成立。现在做一个技巧性操作。我们引入一个关于隐变量Z的分布q(Z)然后把log P(X|θ)改写成log P(X|θ) log [Σ_Z q(Z) * P(X,Z|θ) / q(Z)]这一手看起来只是乘了个1但它巧妙地把P(X,Z|θ)/q(Z)变成了一个关于Z的随机变量权重就是q(Z)。根据琴生不等式log P(X|θ) ≥ Σ_Z q(Z) log [P(X,Z|θ) / q(Z)]右边这一坨就是传说中的ELBO证据下界。它是对数似然的一个下界而且注意右边不再有和的对数全是对数的和好处理太多了。到了这一步EM的数学策略已经呼之欲出直接最大化左边难搞那就去最大化右边的下界下界抬高了左边自然也被推高。3.3 E步让下界变紧下界右边的表达式中除了θ还有一个自由变量q(Z)。我们怎么选q琴生不等式成立是不分q的随便什么q都能给一个下界。但不同q给出的下界松紧程度不同。为了做有效迭代我们希望在一个固定的θ下下界和真实值越接近越好——最好是相等。什么时候等号成立琴生不等式的等号条件是随机变量为常数也就是P(X,Z|θ)/q(Z)不随Z变化。稍微整理就能发现这要求q(Z) P(Z|X, θ)也就是说q要取给定当前参数θ后隐变量的后验分布。这个选择会让ELBO恰好等于log P(X|θ)。这一整步就是E步的数学本质固定当前参数θ_old让q(Z)等于后验分布P(Z|X, θ_old)从而使得当前点上下界与真实值相等。实践中我们甚至不显式算q只需要算后验分布下的某个期望就够用后面讲GMM时会看到具体操作。3.4 M步最大化Q函数选定q之后ELBO里跟θ相关的部分可以单独拎出来Σ_Z q(Z) log P(X,Z|θ)这坨东西就是EM教科书里的Q函数。注意ELBO还有另一项-Σ_Z q(Z) log q(Z)这一项是对q的熵不依赖θ最大化时不用管。M步就是在q固定的情况下把Q函数当成θ的函数求它的最大值得到新的θ_new。因为Q函数的形式通常是对数线性组合求导往往能得到闭式解比原始问题友好得多。3.5 为什么每轮必然上升现在把两件事情串起来看。在E步我们用θ_old求出后验q后这个点上下界和真实值相等log P(X|θ_old) ELBO(θ_old, q)在M步我们最大化Q函数等价于最大化ELBO得到了θ_new。因为θ_new是使Q最大的点所以ELBO(θ_new, q) ≥ ELBO(θ_old, q) log P(X|θ_old)另一方面ELBO永远是对数似然的下界所以log P(X|θ_new) ≥ ELBO(θ_new, q)把两条线连起来log P(X|θ_new) ≥ log P(X|θ_old)完美。每一轮迭代对数似然只升不降。这就解释了EM的单调性。你只要记住这个链条的逻辑比背多少遍公式都管用。这个下界收紧-最大化下界的套路后来也被变分推断直接继承所以理解了EM你等于提前打了个底子。4. 实操案例高斯混合模型GMM手把手推到底4.1 GMM 模型长什么样GMM是EM最广为人知的应用场景甚至可以说提到EM就绕不开GMM。假设数据由K个高斯分布混合产生。每个样本先按概率π_k选一个簇k再从高斯分布N(μ_k, Σ_k)里采样一个点。我们观测到的只有样本点不知道它来自哪个簇这个簇指示变量就是隐变量。模型参数是混合系数π (π_1,...,π_K)均值μ (μ_1,...,μ_K)协方差Σ (Σ_1,...,Σ_K)。约束条件是所有π_k加起来等于1。在工程上GMM常被用来做无监督聚类但跟K-means硬聚类不同的是GMM给的是软聚类——每个样本属于每个簇的概率这在处理有重叠的数据分布时优势非常明显。比如用户行为有交叉的场景硬切一刀往往不符合真实情况GMM的软分配就自然得多。4.2 E步算责任值E步要求出后验分布P(z_ik | x_i, θ_old)。这是熟悉的条件概率问题。一个样本x_i属于第k个簇的后验概率由贝叶斯公式给出γ_ik π_k * N(x_i | μ_k, Σ_k) / Σ_j [π_j * N(x_i | μ_j, Σ_j)]这个γ_ik被称为责任值直观理解就是第i个样本对第k个簇的归属程度。分母是把所有簇的可能性加起来做归一化保证所有γ加起来等于1。如果你让我说一个GMM中最重要的量就是这个γ。整个E步的工程计算基本都在算它。这里有一个隐藏问题N(x_i|μ_k,Σ_k)是高斯密度值当维度高、数据离均值远的时候这个值可能小到超出浮点数精度变成0。分母一变成0γ就变成NaN。后面我会讲应对这种数值问题的log-sum-exp技巧现场写代码时经常被这个坑打乱节奏。4.3 M步更新三个参数有了γ这个软归属度M步的参数更新就非常符合直觉了。新均值用责任值加权的中心。每个样本对均值的贡献不是均等的而是按γ_ik加权μ_k_new Σ_i γ_ik x_i / Σ_i γ_ik这个式子的直觉太直接了——如果你把γ当作样本i对簇k的隶属权重那么新均值就是权重下的加权平均。新协方差也是加权的外积和Σ_k_new Σ_i γ_ik (x_i - μ_k_new)(x_i - μ_k_new)^T / Σ_i γ_ik注意这里用的是刚更新完的新均值。这一步要放在均值更新之后。新混合系数簇k在所有样本中的平均责任值π_k_new Σ_i γ_ik / N它反映了平均来看有多少比例的样本属于簇k。这是在约束Σπ_k1下的自然结果可以用拉格朗日乘子法严格推出来但直觉上也不难接受。至于为什么最大化Q函数会得到这样的式子简单推导一下Q函数展开后是一个关于μ_k的二项式对μ_k求导并令其为零整理后可解出μ_k的闭式解就是上面那个加权平均形式。协方差的推导也是对称的套路。4.4 完整跑一遍两轮迭代的数字演示用个极端简化的小例子让你直观看到数值在怎么变。假设只有5个二维数据点K2x1(1.0, 1.0), x2(1.5, 0.5), x3(5.0, 5.0), x4(5.5, 5.5), x5(5.0, 4.5)一眼能看出两个簇一个在(1,1)附近一个在(5,5)附近。但EM算法没有眼睛它只有数学。初始化假设π(0.5,0.5)μ1(0,0)μ2(3,3)协方差矩阵都用单位阵I。第一轮E步计算x1的γ。高斯密度N(x1|μ1,I)在(1,1)处正比于exp(-16/2)也就是e^{-8}N(x1|μ2,I)在(1,1)到(3,3)的距离平方是8正比于e^{-4}。归一化后x1对簇1的γ约为e^{-8}/(e^{-8}e^{-4}) ≈ 0.018。也就是说当前参数下x1几乎被分给簇2。其他点类似算。x3(5,5)到μ2(3,3)的距离平方是8到μ1(0,0)的距离平方是50所以它强烈归簇2。这样第一轮E步算下来所有点几乎都归簇2簇1的权重总和接近于0。M步一看这个局势μ1会往空转的方向走——分子是各个点乘上极小的权重数值上非常不稳定事实上此时的μ1计算出来仍然接近0但μ2会被所有点拉向大约(3.6, 3.3)的位置。第二轮E步以μ2(3.6,3.3)为基准重新算γ此时左边的点x1、x2到μ2的距离仍然比到μ1近归属依旧倾向簇2。但注意随着μ2向右移动x1、x2相对距离也在变化γ_11会开始慢慢变大。再迭代几轮μ1会被拉到左边两个点的中心(1.25,0.75)μ2被拉到右边三个点的中心(5.17,5.0)。最终收敛就是这两组均值π相应接近(0.4,0.6)。这个数值示例虽然粗糙但你能清楚看到EM的两大动态一个是责任值决定参数走向另一个是参数重新调整责任值。两者共同作用最后自动收敛到直观上正确的答案。5. 常见问题与调试技巧实录5.1 对数似然下降了怎么办理论保证EM单调上升那你代码里为什么会出现对数似然下降我遇到过的情况主要是这几种log计算出了问题概率值下溢成0然后取log得到-inf或者出现了NaN。检查你的事件有没有用log-sum-exp技巧。E步和M步用的参数不是一套某些实现里一个不小心上一轮更新的参数没有传进E步导致Q函数和实际计算对不上。责任值没有归一化归一化是求后验的基本要求漏掉会导致步长乱套。排查技巧在每轮迭代前后都打印对数似然值一旦出现下降直接用二分查找定位到是哪一步、哪一组参数导致的。这个习惯帮我省了无数时间。5.2 协方差矩阵奇异了一个常见场景某个簇被分配到的样本数太少或者某些样本几乎线性相关导致协方差矩阵变成奇异矩阵求逆的时候直接炸掉。应对思路有三条在M步更新协方差时加一个小的正则项比如Σ_k λI数值上更稳健。这相当于给协方差施加一个微弱的高斯先验。监控每个簇的有效样本量N_k Σ_i γ_ik如果某个N_k小于某个阈值比如特征维度d的10倍说明这个簇快空了可以停止迭代或合并这个簇。如果协方差全用全矩阵确实不稳可以先试对角协方差很多业务场景下结果差别不大但数值稳定性好得多。5.3 初始化选不好收敛到垃圾解EM对初值极其敏感。我做的第一个GMM就吃过这个亏随机初始化结果两个簇的中心被初始化成同一个点整个模型最后退化成了一个簇对数似然毫无起伏。现在我的习惯是三步走先跑一遍K-means用K-means的聚类中心作为GMM的初始均值混合系数用各簇样本占比。这个初始化的组合拳在绝大多数场景下都很稳。然后用多个随机种子比如5到10个分别跑EM比较最后收敛的对数似然选最大的那个。虽然有点暴力但实验证明这个策略基本不会让你错过好解。如果数据量极大还有一个办法用一小部分数据先完成初始化再对全量数据跑EM收敛速度通常也不错。5.4 K值怎么定GMM里的K通常事先不知道怎么设。一种比较靠谱的做法是画BIC贝叶斯信息准则曲线对K从1到K_max分别跑EM计算BIC画出来后找拐点。BIC的公式是BIC -2 * log P(X|θ) m * log N其中m是模型参数个数。K增大时第一项似然惩罚通常会变小第二项参数数量惩罚会变大最优K出现在两者平衡的位置。不要只依赖似然值选K。似然值永远是K越大越高因为模型越灵活但过拟合风险增加。BIC、AIC这类带惩罚项的指标才是工程上的决定依据。5.5 死循环与收敛判定的坑早期我犯过一个很经典的错误用参数变化量作为收敛判据阈值设了1e-6。结果有一次在接近极值点的地方参数变化量一直在1e-5量级徘徊程序跑了几个小时没停。正确做法是同时监控对数似然的变化量。如果两轮迭代对数似然的变化绝对值小于一个阈值比如1e-3说明收益已经微乎其微可以停了。参数变化量作为辅助判据。另外尽管理论上EM单调递增实际上受浮点数精度影响最后几轮可能出现微小波动。所以停止条件里建议加上连续若干轮变化都小于阈值而不是某一轮变化小于阈值。6. 变体与应用场景延伸说句实话纯EM在你实际工作里大概率不会以教科书里的样子出现但你理解了它很多高级且花哨的东西就有了地基。先说变体变分EM如果后验P(Z|X,θ)本身不可解E步无法精确完成就用变分分布去近似它这就是变分推断的雏形。MCEM蒙特卡洛EME步的期望没有闭式解时用采样的方法近似计算期望比如用MCMC采样抽取隐变量样本。在线EM / SGD-EM数据量太大不能全部进内存时把E步和M步放到小批量式更新框架里逐块更新参数。这些变体虽然各显神通但骨架跟我前面讲的完全一致。再看应用HMM隐马尔可夫模型的参数学习大名鼎鼎的Baum-Welch算法本质就是EM在HMM上的展开E步对应前向后向算法算状态后验M步更新转移概率和发射概率。LDA主题模型的求解虽然具体用的是变分EM或吉布斯采样但底层思路仍是隐变量-观测文本-参数估计的三件套。缺失数据填补问卷、传感器数据里常见的缺失值可以建模成隐变量用EM迭代估计出合理填补值。推荐系统的隐反馈把用户未点击但也未必讨厌的物品建模成隐变量借助EM思路做更精细的建模。你会发现只要问题里藏着看不见的状态变量EM就有用武之地。它的意义不只是单一算法更是一种处理不完整数据的通用范式。关于实际工程应用我想多说两句。我在做用户消费行为分析时用GMM做人群分层的频率比K-means高得多。原因很简单用K-means做硬聚类边界上的用户分类结果很不稳定模型一次微调就可能让大批用户跳跃到别的群组运营侧的人根本没法接受。GMM给出的是属于群A的概率是0.6、群B的概率是0.4这样的描述业务可以设置一个阈值来决定怎么处理边界用户。这种信息量硬聚类给不了。另一个心得体会是跑EM之前强烈建议先用K-means预热。K-means本身在GMM初始化场景下几乎成了标准的配置你说它是EM的一部分也没错。用K-means给出的中心做初始均值再用样本点到中心的最小距离做初始协方差估计比如用类内样本的协方差后面的EM收敛速度和稳定性都会好很多。还有一条写EM代码时要极力避免显式存储所有N×K的责任值矩阵。那玩意在数据量大的时候就是个内存炸弹。更优雅的做法是在M步需要统计量的时候一遍过数据、边算边累加把中间统计量直接累积起来。这种做法不需要保留全量γ工程上可伸缩性强得多。当初我接手一个千万级用户的数据建模项目第一版代码就是因为这个矩阵直接内存溢出了。改成流式更新后整个训练流程跑完只用了几百兆内存。对比一下EM和K-means的差异用一个单口相声式的总结K-means是让每个人站好队队长往队的中心挪EM是让大家按概率站在各队之间队长按加权中心挪站队概率也顺势更新。前者刚烈后者柔和。但在一个高度混合、簇之间有重叠的数据集上柔和的那套往往才是真正贴近现实的那套。最后如果你想自我检查是否真正理解了EM一个特别有效的测试是你能不能向一个不懂数学的同事用两三句话讲清楚EM在干什么。我自己的版本是先把数据按当前猜测的模型软分配到各个簇然后重新估计每个簇的参数反复做直到参数不再明显变化。只要这句话你能发自内心地说出来EM算法就已经真正属于你了。