4.3 高斯混合与期望最大化

动机

位于两个群组之间的数据，未必明显属于其中一组。硬分配只保留选中的群组；混合模型则结合多个概率密度，为每条观测保留属于各分量的概率，再用这些概率加权重新估计分量。

困难在于似然的每一项都包含“总和的对数”。即使单个高斯分量很简单，不同分量的参数仍通过总和相互关联。期望最大化算法，简称 EM，引入责任度，将参数更新拆成加权估计问题。它的两个步骤共同改善一个在当前参数处与观测对数似然相切的下界。

本节先推导权重和更新式，再完整计算一次迭代。软概率与硬标签的区别是关键：最大化选中分量的密度，与最大化混合密度具有不同的目标。加权协方差的分母也必须从所最大化的似然中推导。

几何软权重

软权重也可以由几何定义。设互异站点 $S=\{s_1,\ldots,s_K\}$ 的旧 Voronoi 单元为 $V_j$，距离相同时固定分给其中一个单元。插入查询点 $x\notin S$，用 $V_x$ 表示它的新单元，并假设该单元有界，体积有限且为正。定义

$$0<\operatorname{vol}(V_x)<\infty,\qquad w_j=\frac{\operatorname{vol}(V_x\cap V_j)} {\operatorname{vol}(V_x)},\qquad \sum_{j=1}^K w_j=1.$$

旧单元将新单元分割，所以重叠体积之和等于新单元体积。各权重非负，没有正体积重叠的站点得到零权重；这些比例衡量新区域原来分属各旧站点的份额。

一维体积就是长度。取两个站点和一个新查询点：

$$S=\{0,4\},\quad x=1,\quad V_x=\left[\frac12,\frac52\right],\quad (w_1,w_2)=\left(\frac{3/2}{2},\frac{1/2}{2}\right) =\left(\frac34,\frac14\right).$$

旧边界为 $2$，新区间在其左右的重叠长度分别是 $3/2$ 和 $1/2$，总长为 $2$。这是几何分配权重；下文的高斯责任度是混合模型下的条件概率。两者都跨分量规范，但定义所用的信息不同。

定义

设 $x_1,\ldots,x_n\in\mathbb R^d$ 为独立观测，固定高斯分量数目 $K$。分量 $k$ 具有权重 $\pi_k>0$、均值 $\mu_k\in\mathbb R^d$ 和正定协方差 $\Sigma_k$，且 $\sum_{k=1}^K\pi_k=1$。用 $\theta$ 表示全部参数，以下对数均为自然对数。

定义

高斯混合与观测对数似然

高斯密度和混合密度分别为

$$\phi_d(x;\mu,\Sigma) =\frac{\exp\!\left(-\frac12(x-\mu)^\top\Sigma^{-1}(x-\mu)\right)} {(2\pi)^{d/2}\det(\Sigma)^{1/2}}, \qquad p_\theta(x)=\sum_{k=1}^K\pi_k\phi_d(x;\mu_k,\Sigma_k).$$

观测对数似然为

$$\ell(\theta)=\sum_{i=1}^n\log p_\theta(x_i).$$

潜在分量标签 $z_i$ 以概率 $\pi_k$ 取值 $k$；分量密度描述给定该标签时 $x_i$ 的分布。

混合权重本身属于模型，不是选定分量后才加入的规范。即使两个分量在某条观测处的密度相同，权重仍然影响结果。对全部可能标签求和，才得到观测密度。

如果标签也已观测，对数似然就是 $\log\pi_{z_i}+\log\phi_d(x_i;\mu_{z_i},\Sigma_{z_i})$ 的总和，也可以用各分量的零一指示量选出对应项。未知标签使我们无法直接按这种方式分组；EM 转而取这些指示量的条件期望，得到 $Q$ 中的责任度。因此，分数权重具有明确的概率含义，是当前模型下潜在指示量的期望。

这也区分了推导中的三个函数：观测对数似然在每个对数内对标签求和；完整数据对数似然使用已知标签；$Q$ 则在当前后验分布下平均完整数据表达式。分清观测量、潜在量和暂时固定的量，可以理解加权高斯拟合为什么能够改善更复杂的混合似然。

定义

责任度与分量的有效样本数

在当前参数 $\theta^{\mathrm{old}}$ 下，Bayes 公式给出

$$r_{ik} =\Pr_{\theta^{\mathrm{old}}}(z_i=k\mid x_i) =\frac{\pi_k^{\mathrm{old}} \phi_d(x_i;\mu_k^{\mathrm{old}},\Sigma_k^{\mathrm{old}})} {\sum_{j=1}^K\pi_j^{\mathrm{old}} \phi_d(x_i;\mu_j^{\mathrm{old}},\Sigma_j^{\mathrm{old}})}, \qquad N_k=\sum_{i=1}^n r_{ik}.$$

每一行满足 $\sum_k r_{ik}=1$。混合权重为正、协方差正定时，每个 $r_{ik}$ 都为正，所以 $N_k>0$，且 $\sum_kN_k=n$。

$N_k$ 是在当前后验标签下，预期归于分量 $k$ 的观测数量，可以不是整数。责任度是给定当前模型与观测后的概率，不是随数据提供的已知类别标签。

这里存在两种不同的归一化。E 步针对每条观测跨分量规范，因此每行总和为一；均值更新则针对一个分量跨观测规范，分母是列总和 $N_k$。若均值更新除以 $n$，就会把分给其他分量的责任度也计入本分量总量。两个分母回答不同的条件问题。

一点的高斯密度本身不是成员概率：密度可以大于一，也会随测量单位改变尺度。责任度则是比较同一观测的加权密度后得到的无量纲概率。一般协方差下，比较同时包含指数中的 Mahalanobis 距离和规范中的行列式。更分散的分量具有更大体积和更低峰值，$\pi_k$ 控制其总体混合权重。在等权重、共同各向同性协方差下，分量密度排序可以化为欧几里得距离比较。

定理/命题

定理

高斯 M 步使用有效样本数规范

固定责任度，定义完整数据对数似然的条件期望

$$Q(\theta\mid\theta^{\mathrm{old}}) =\sum_{i=1}^n\sum_{k=1}^K r_{ik}\left[\log\pi_k+\log\phi_d(x_i;\mu_k,\Sigma_k)\right].$$

当下面的加权协方差均为正定时，最大化更新为

$$\pi_k^{\mathrm{new}}=\frac{N_k}{n},\qquad \mu_k^{\mathrm{new}}=\frac1{N_k}\sum_{i=1}^nr_{ik}x_i,$$$$\Sigma_k^{\mathrm{new}} =\frac1{N_k}\sum_{i=1}^nr_{ik} (x_i-\mu_k^{\mathrm{new}}) (x_i-\mu_k^{\mathrm{new}})^\top.$$

若模型预先指定共同协方差为 $\sigma^2I_d$、$\sigma^2>0$，则保持该协方差固定，只使用相同的权重和均值更新。

当且仅当加权中心化观测张成全部特征方向，这个协方差才正定；否则所得矩阵奇异，不属于正定高斯协方差。固定方差模型使我们能够在不估计协方差矩阵的情况下，逐步计算有效模型更新。

定理

有效的 EM 迭代不会降低观测对数似然

先在 $\theta^{\mathrm{old}}$ 处计算精确后验责任度，再选择满足下式的有效新参数：

$$Q(\theta^{\mathrm{new}}\mid\theta^{\mathrm{old}}) \ge Q(\theta^{\mathrm{old}}\mid\theta^{\mathrm{old}}).$$

则

$$\ell(\theta^{\mathrm{new}})\ge\ell(\theta^{\mathrm{old}}).$$

精确 M 步最大化 $Q$；只使它增加的更新称为广义 EM 步。结论针对观测对数似然的值；全局最优与参数向量收敛是另外的数学主张。

证明思路

推导加权更新

权重部分是在 $\sum_k\pi_k=1$ 下最大化 $\sum_kN_k\log\pi_k$。引入乘数 $\lambda$，驻点条件为 $N_k/\pi_k=\lambda$。对 $k$ 求和得到 $\lambda=n$，所以 $\pi_k=N_k/n$。

各均值的相关项是加权 Mahalanobis 距离平方和的负数，求导得到

$$\Sigma_k^{-1}\sum_i r_{ik}(x_i-\mu_k)=0.$$

由于 $\Sigma_k$ 可逆且 $N_k>0$，均值必须是命题中的加权平均。预先指定共同各向同性协方差时，这就是加权最小二乘的最小点。每条观测按责任度贡献给各分量，不只是贡献给责任度最大的分量。

推导协方差时，令

$$S_k=\sum_i r_{ik}(x_i-\mu_k^{\mathrm{new}}) (x_i-\mu_k^{\mathrm{new}})^\top,\qquad \Lambda_k=\Sigma_k^{-1}.$$

$Q$ 的相关部分为

$$\frac{N_k}{2}\log\det\Lambda_k -\frac12\operatorname{tr}(\Lambda_kS_k).$$

驻点方程是 $N_k\Lambda_k^{-1}=S_k$。当 $S_k$ 正定时，由对数行列式的凹性，最大点为 $\Sigma_k=S_k/N_k$。这是似然估计，所以分母为 $N_k$，不减一；省略分母得到的是尺度不同的加权散布矩阵。

与当前似然相切的下界

对任意正概率行 $q_i=(q_{i1},\ldots,q_{iK})$，Jensen 不等式给出

$$\log\sum_k\pi_k\phi_d(x_i;\mu_k,\Sigma_k) \ge\sum_kq_{ik}\log \frac{\pi_k\phi_d(x_i;\mu_k,\Sigma_k)}{q_{ik}}.$$

将右侧对观测求和，定义

$$\mathcal F(q,\theta) =\sum_{i,k}q_{ik} \left[\log\pi_k+\log\phi_d(x_i;\mu_k,\Sigma_k)-\log q_{ik}\right].$$

精确差距为

$$\ell(\theta)-\mathcal F(q,\theta) =\sum_i\operatorname{KL}\!\left(q_i\Vert r_i(\theta)\right)\ge0,$$

其中 $\operatorname{KL}(q\Vert r)=\sum_kq_k\log(q_k/r_k)$，$r_i(\theta)$ 为参数 $\theta$ 下的后验行。当 $q_i=r_i(\theta)$ 时，差距为零。

E 步选择 $q=r(\theta^{\mathrm{old}})$，使下界在当前观测对数似然处取等号。固定 $q$ 后，其熵 $-\sum_{i,k}q_{ik}\log q_{ik}$ 为常数，所以增加 $Q$ 也会增加 $\mathcal F$，得到

$$\ell(\theta^{\mathrm{new}}) \ge\mathcal F(r^{\mathrm{old}},\theta^{\mathrm{new}}) \ge\mathcal F(r^{\mathrm{old}},\theta^{\mathrm{old}}) =\ell(\theta^{\mathrm{old}}).$$

两个步骤各有作用：E 步使下界在当前模型处精确，M 步在固定下界下改善模型。用硬标签代替后验概率，会失去证明所依赖的相切等号。

一次 M 步始终固定责任度。新均值、新协方差和新权重，都使用完整旧参数算出的同一后验矩阵。协方差公式用新加权均值中心化，因为它最大化的是同一个固定 $Q$ 问题。形成新模型之后，下一次 E 步才重新计算后验；不等式链中间的比较，是保留一个概率矩阵而改变参数。

按 $0\log0=0$ 并通过连续延伸解释零权重项，下界与 KL 恒等式也适用于含零元素的概率行。因此可以比较精确后验的相切下界与硬分配产生的较松下界。

例题

例题

1. 相同的分量密度仍可具有不同混合权重

考虑一维混合模型，均值为 $\mu_1=0$、$\mu_2=2$，共同方差为 $\sigma^2=1$，权重为 $(\pi_1,\pi_2)=(3/4,1/4)$。在 $x=1$，到两个均值的距离相同，因此高斯密度相同。代入 Bayes 公式，共同密度抵消：

$$r_1=\frac{(3/4)\phi_1(1;0,1)} {(3/4)\phi_1(1;0,1)+(1/4)\phi_1(1;2,1)} =\frac34,\qquad r_2=\frac14.$$

中点与两个分量形状同样相容，但模型给第一分量三倍的先验权重。省略混合权重、只规范两个密度，会得到 $1/2$ 和 $1/2$；那是另一个混合模型的后验。所以软成员概率同时取决于密度和权重。

例题

2. 固定方差下的完整一次 E 步与 M 步

取观测 $x_1=0$、$x_2=2$，初始均值为 $(0,2)$，混合权重相同，并指定共同方差为 $1$。令

$$a=\frac1{1+e^{-2}}\approx0.880797078.$$

每条观测到另一均值的距离平方，比到对应均值大 $4$，所以密度比为 $e^{-2}$。因此

$$R=\begin{bmatrix}a&1-a\\1-a&a\end{bmatrix}, \qquad N_1=N_2=1.$$

M 步给出

$$\pi^{\mathrm{new}}=(1/2,1/2),\qquad \mu^{\mathrm{new}}=(2(1-a),2a) \approx(0.238405844,1.761594156).$$

每条观测最初虽然恰好位于一个分量均值处，对另一个分量仍有正责任度，从而把两个新均值向内拉。共同方差保持为 $1$。使用完整高斯密度计算，观测对数似然从 $-2.970315405$ 升至 $-2.886562696$，约增加 $0.083752710$。

下面的 NumPy 程序重现同一次迭代。对数规范同时使用全部分量计算后验，不会先挑选某个分量再形成分母。

import numpy as np

X = np.array([0.0, 2.0])
means = np.array([0.0, 2.0])
weights = np.array([0.5, 0.5])
variance = 1.0

def log_components(means, weights):
    return (np.log(weights)[None, :]
            - 0.5 * np.log(2.0 * np.pi * variance)
            - 0.5 * (X[:, None] - means[None, :]) ** 2 / variance)

joint = log_components(means, weights)
normalizer = np.logaddexp.reduce(joint, axis=1)
responsibilities = np.exp(joint - normalizer[:, None])
mass = responsibilities.sum(axis=0)
new_weights = mass / len(X)
new_means = (responsibilities.T @ X) / mass
before = normalizer.sum()
after = np.logaddexp.reduce(
    log_components(new_means, new_weights), axis=1
).sum()
assert np.allclose(responsibilities.sum(axis=1), 1.0)
assert after >= before
print("responsibilities =", np.round(responsibilities, 6).tolist())
print("mass =", np.round(mass, 6).tolist())
print("weights =", np.round(new_weights, 6).tolist())
print("means =", np.round(new_means, 6).tolist())
print(f"log_likelihood: {before:.6f} -> {after:.6f}")

输出为

responsibilities = [[0.880797, 0.119203], [0.119203, 0.880797]]
mass = [1.0, 1.0]
weights = [0.5, 0.5]
means = [0.238406, 1.761594]
log_likelihood: -2.970315 -> -2.886563

例题

3. 协方差塌缩可以使似然无上界

取观测 $0$ 和 $3$。两个分量的权重各为 $1/2$，均值为 $0$ 和 $3$，方差分别为 $\varepsilon^2$ 和 $1$。每个 $\varepsilon>0$ 都给出有效高斯分量，其混合密度满足

$$p_\theta(0)\ge\frac1{2\sqrt{2\pi}\varepsilon}, \qquad p_\theta(3)\ge\frac1{2\sqrt{2\pi}}.$$

因此

$$\ell(\theta)\ge-\log\varepsilon-\log(8\pi) \longrightarrow+\infty \qquad(\varepsilon\downarrow0).$$

将一个分量收缩到某条观测附近，可以使该处密度任意高；另一个分量继续解释其余观测。这个不受限制的似然没有有限最大值。问题在统计目标本身，不是要求计算奇异协方差的高斯密度。像例题 2 那样预先指定正的共同方差，是另一个不会出现这种塌缩的模型。

如何理解 EM 的进度

指定共同各向同性方差时，各分量密度至多为 $(2\pi\sigma^2)^{-d/2}$；它们的加权平均也有相同上界，所以

$$\ell(\theta)\le-\frac{nd}{2}\log(2\pi\sigma^2).$$

因此，这个模型的不下降似然值存在有限极限。这是针对实数目标值序列的结论，不会指定最佳混合模型，也不会建立每个参数坐标的极限。交换分量标签也会描述同一混合密度。

初始模型决定第一次后验权重，这些权重再决定第一次加权拟合问题。不同初始模型可以沿不同的相切下界序列前进。似然用于检查所选模型内部的优化进度，不会自动选择分量数量或比较不同的协方差假设。

等权重、共同各向同性方差下，责任度是负距离平方的指数再作规范。方差越小，责任度越集中于唯一最近的均值。对于每个正方差，后验仍然是软分配；硬最近中心更新则改变优化问题。这解释了加权均值更新与中心式聚类的相似之处，同时保留两种目标的区别。

常见错误

常见错误

用硬 argmax 代替 E 步

从例题 2 的新均值出发，硬分配将 $0$ 归于第一分量、$2$ 归于第二分量，再按各组观测均值更新，就会回到 $(0,2)$。共同方差与等权重不变时，观测对数似然从约 $-2.886563$ 降至 $-2.970315$。硬分配不是使 EM 下界与当前似然相切的后验分布。

常见错误

将散布当成协方差，或省略混合权重

高斯协方差更新将加权散布除以 $N_k$；E 步则跨分量规范 $\pi_k\phi_d(x_i;\mu_k,\Sigma_k)$。两个因子都来自混合似然。改变任何一个规范，都会改变拟合方程，使上述 EM 推导不再成立。

总结

高斯混合将分量密度加权求和。E 步使用规范后的后验责任度，M 步按有效样本数更新权重、均值和需要估计的协方差。Jensen 下界在当前后验处精确，M 步增加 $Q$ 时下界也增加，从而证明观测对数似然不下降。正权重、正定协方差和正有效样本数，使公式保持在指定模型内。估计协方差或预先指定方差，是目标的一部分；不受限制的协方差可以因分量塌缩而使似然无上界。

练习

思考检查

1. 在 $x=1$，两个分量的均值为 $0$ 和 $2$、共同方差为 $1$、权重为 $(3/4,1/4)$。求两个责任度。只保留最大责任度对应的分量，会丢失什么信息？

分清相同分量密度与不同混合权重。

思考检查

2. 对观测 $(0,2)$，固定责任度为 $R=\begin{bmatrix}0.9&0.1\\0.3&0.7\end{bmatrix}$。指定共同方差时，求 M 步的权重和均值。

先计算每个 $N_k$，再用它规范加权总和。

思考检查

3. 对观测 $(0,2,4)$，某分量的责任度为 $(4/5,1/2,1/5)$。求有效样本数、新均值，以及估计协方差时的一维协方差值。

将围绕新均值的加权距离平方和除以有效样本数。

思考检查

4. 对例题 2 更新后的模型作硬分配和均值更新。比较得到的观测对数似然与当前值；EM 证明中的哪个等号不再成立？

当前两个均值的中点仍为 $1$。

解答

解答 · 1. 后验保留两个比例

责任度为 $(3/4,1/4)$，因为共同密度会抵消。只保留第一分量，就只记录一个选定标签，丢失剩余后验质量 $1/4$。

解答 · 2. 分数样本量也会更新混合权重

有效样本数为 $(N_1,N_2)=(6/5,4/5)$。所以 $\pi^{\mathrm{new}}=(3/5,2/5)$，且 $\mu_1^{\mathrm{new}}=(0.9\cdot0+0.3\cdot2)/(6/5)=1/2$， 而 $\mu_2^{\mathrm{new}}=(0.1\cdot0+0.7\cdot2)/(4/5)=7/4$。

解答 · 3. 用有效样本数规范加权散布

有效样本数为 $N=3/2$，均值为 $((4/5)0+(1/2)2+(1/5)4)/(3/2)=6/5$。 加权距离平方和为 $(4/5)(6/5)^2+(1/2)(4/5)^2+(1/5)(14/5)^2=76/25$， 所以协方差估计为 $(76/25)/(3/2)=152/75$。

解答 · 4. 下界不再与当前似然相切

硬分配使每个分量各有一条观测，均值回到 $(0,2)$，权重回到 $(1/2,1/2)$。观测对数似然从 $-2.886562696$ 降至 $-2.970315405$。硬概率行 $q$ 不同于严格为正的后验，因此 $\mathcal F(q,\theta^{\mathrm{old}})=\ell(\theta^{\mathrm{old}})$ 不再成立。改善这个硬分配下界，无法用 EM 的不等式链保证观测似然改善。

可以继续阅读最小二乘与残差， 或查看课程阅读书目。