NO.33.tip： 边界优化

背景

我们将讨论一类称为边界优化（bound optimization）或者 MM 的算法。在最小化的上下文中，MM 代表多数化最小化（majorize–minimize）。在最大化的上下文中，MM 代表少数化最大化（minorize–maximize）。MM 算法有很多例子，例如期望最大化算法、近端梯度法、聚类的均值偏移算法等。

通用算法

我们假设目标是最大化某个函数的相对于其参数 $\mathbf{\theta }$ 的 $\ell (\mathbf{\theta })$。MM 算法的基本方法是构造一个代理函数（surrogate function）$Q(\mathbf{\theta },{\mathbf{\theta }}^t)$，代理函数是 $\ell (\mathbf{\theta })$ 的紧下界。

使得 $Q(\mathbf{\theta },{\mathbf{\theta }}^t)\le \ell (\mathbf{\theta })$ 并且 $Q({\mathbf{\theta }}^t,{\mathbf{\theta }}^t)=\ell ({\mathbf{\theta }}^t)$。如果满足这些条件，我们称 $Q$ 最小化 $\ell$。然后，我们在每个步骤执行以下更新：

$${\mathbf{\theta }}^{t+1}=\arg \underset{\mathbf{\theta }}{\max }Q(\mathbf{\theta },{\mathbf{\theta }}^t)\text{\hspace{1em}(1)}$$

这保证了我们在原始目标中单调增加：

$$\ell ({\mathbf{\theta }}^{t+1})\ge Q({\mathbf{\theta }}^{t+1},{\mathbf{\theta }}^t)\ge Q({\mathbf{\theta }}^t,{\mathbf{\theta }}^t)=\ell ({\mathbf{\theta }}^t)\text{\hspace{1em}(2)}$$

其中，第一个不等式成立，因为 $Q({\mathbf{\theta }}^t,{\mathbf{\theta }}^t)$ 是任何 ${\mathbf{\theta }}^t$ 的 $\ell ({\mathbf{\theta }}^t)$ 的下界；第二个不等式由式 (1) 得出；最后的等式遵循紧密性性质。基于这个结果，如果没有观察到目标的单调增加，那么所采用的数学和/或代码一定存在错误。这是一个功能强大的调试工具。

示例：逻辑回归

如果 $\ell (\mathbf{\theta })$ 是我们想要最大化的凹函数，那么获得有效下界的一种方法是在其黑塞矩阵上使用一个界，即找到一个负定矩阵 $\mathbf{B}$，使得 $\mathbf{H}(\mathbf{\theta })\succ \mathbf{B}$。在这种情况下，可以证明：

$$\ell (\mathbf{\theta })\ge \ell ({\mathbf{\theta }}^t)+(\mathbf{\theta }-{\mathbf{\theta }}^t{)}^{\mathsf{T}}\mathbf{g}({\mathbf{\theta }}^t)+\frac{1}{2}(\mathbf{\theta }-{\mathbf{\theta }}^t{)}^{\mathsf{T}}\mathbf{B}(\mathbf{\theta }-{\mathbf{\theta }}^t)\text{\hspace{1em}(3)}$$

其中，$\mathbf{g}({\mathbf{\theta }}^t)=\nabla \ell ({\mathbf{\theta }}^t)$。因此，以下函数是一个有效的下界：

$$Q(\mathbf{\theta },{\mathbf{\theta }}^t)={\mathbf{\theta }}^{\mathsf{T}}(\mathbf{g}({\mathbf{\theta }}^t)-\mathbf{B}{\mathbf{\theta }}^t)+\frac{1}{2}{\mathbf{\theta }}^{\mathsf{T}}\mathbf{B}\mathbf{\theta }\text{\hspace{1em}(4)}$$

对应的更新变为：

$${\mathbf{\theta }}^{t+1}={\mathbf{\theta }}^t-{\mathbf{B}}^{-1}\mathbf{g}({\mathbf{\theta }}^t)\text{\hspace{1em}(5)}$$

这类似于牛顿更新，区别在于我们使用 $\mathbf{B}$ 而不是 $\mathbf{H}({\mathbf{\theta }}^t)$，$\mathbf{B}$ 是一个固定矩阵，这个固定矩阵在每次迭代时都会变化。这种方法可以在较低的计算成本下为我们提供二阶方法的一些优势。

例如，让我们使用 MM 拟合一个多类别逻辑回归模型。样例 $n$ 属于类别 $c\in \{1,\dots ,C\}$ 的概率为：

$$p(y_n=c\mid {\mathbf{x}}_n,\mathbf{w})=\frac{\exp ({\mathbf{w}}_c^{\mathsf{T}}{\mathbf{x}}_n)}{{\sum }_{i=1}^C\exp ({\mathbf{w}}_i^{\mathsf{T}}{\mathbf{x}}_n)}\text{\hspace{1em}(6)}$$

由于归一化条件 ${\sum }_{c=1}^{C}p(y_n=c\mid {\mathbf{x}}_n,\mathbf{w})=1$，我们可以设置 ${\mathbf{w}}_C=0$。（例如，在二元逻辑回归中，其中 $C=2$，我们只学习单个权重向量。）因此，参数 $\mathbf{\theta }$ 对应于大小为 $D(C-1)$ 的权重矩阵 $\mathbf{w}$，其中 $\mathbf{x}\in {\mathbb{R}}^D$。

如果我们设 ${\mathbf{p}}_n(\mathbf{w})=[p(y_n=1\mid {\mathbf{x}}_n,\mathbf{w}),\dots ,p(y_n=C-1\mid {\mathbf{x}}_n,\mathbf{w})]$ 以及 ${\mathbf{y}}_n=[\mathbb{I}(y_n=1),\dots ,\mathbb{I}(y_n=C-1)]$，那么可以将似然记作：

$$\ell (\mathbf{w})=\sum _{n=1}^N\left[\sum _{c=1}^{C-1}{y}_{nc}{\mathbf{w}}_c^{\mathsf{T}}{\mathbf{x}}_n-\log \sum _{c=1}^C\exp ({\mathbf{w}}_c^{\mathsf{T}}{\mathbf{x}}_n)\right]\text{\hspace{1em}(7)}$$

梯度由以下公式给出：

$$\mathbf{g}(\mathbf{w})=\sum _{n=1}^N({\mathbf{y}}_n-{\mathbf{p}}_n(\mathbf{w}))\otimes {\mathbf{x}}_n\text{\hspace{1em}(8)}$$

其中，$\otimes$ 表示 Kronecker 乘积（在这种情况下，该乘积只是两个向量的外积）。黑塞矩阵由以下公式给出：

$$\mathbf{H}(\mathbf{w})=-\sum _{n=1}^N(\mathrm{diag}({\mathbf{p}}_n(\mathbf{w}))-{\mathbf{p}}_n(\mathbf{w}){\mathbf{p}}_n(\mathbf{w}{)}^{\mathsf{T}})\otimes ({\mathbf{x}}_n{\mathbf{x}}_n^{\mathsf{T}})\text{\hspace{1em}(9)}$$

我们可以在黑塞矩阵上构造一个下界：

$$\mathbf{H}(\mathbf{w})\succ -\frac{1}{2}\left[\mathbf{I}-11^{\mathsf{T}}/C\right]\otimes \left(\sum _{n=1}^N{\mathbf{x}}_n{\mathbf{x}}_n^{\mathsf{T}}\right)\triangleq \mathbf{B}\text{\hspace{1em}(10)}$$

其中，$\mathbf{I}$ 是 $(C-1)$ 维单位矩阵，$1$ 是取值全为 $1$ 的 $(C-1)$ 维向量。在二元的情况下，这变成：

$$\mathbf{H}(\mathbf{w})\succ -\frac{1}{2}\left(1-\frac{1}{2}\right)\left(\sum _{n=1}^N{\mathbf{x}}_n^{\mathsf{T}}{\mathbf{x}}_n\right)=-\frac{1}{4}{\mathbf{X}}^{\mathsf{T}}\mathbf{X}\text{\hspace{1em}(11)}$$

这也遵循以下结论：由于 $p_n\le 0.5$，因此 $-(p_n-p_n^2)\ge -0.25$。

我们可以使用这个下界来构造一个 MM 算法，以查找最大似然估计。公式更新为：

$${\mathbf{w}}^{t+1}={\mathbf{w}}^t-{\mathbf{B}}^{-1}\mathbf{g}({\mathbf{w}}^t)\text{\hspace{1em}(12)}$$

例如，让我们考虑二元的情况，因此 ${\mathbf{g}}^t=\nabla \ell ({\mathbf{w}}^t)={\mathbf{X}}^{\mathsf{T}}(\mathbf{y}-{\mathbf{\mu }}^t)$，其中 ${\mathbf{\mu }}^t=[p_n({\mathbf{w}}^t),(1-p_n({\mathbf{w}}^t)){]}_{n=1}^N$。公式更新为：

$${\mathbf{w}}^{t+1}={\mathbf{w}}^t-4({\mathbf{X}}^{\mathsf{T}}\mathbf{X}{)}^{-1}{\mathbf{g}}^t\text{\hspace{1em}(13)}$$

上述算法（每一步）比 IRLS（Iteratively Reweighted Least Squares，迭代重加权最小二乘法）算法（即牛顿法）更快，后者是拟合广义线性模型的标准方法。为了理解这一点，请注意牛顿更新具有以下形式：

$${\mathbf{w}}^{t+1}={\mathbf{w}}^t-{\mathbf{H}}^{-1}\mathbf{g}({\mathbf{w}}^t)={\mathbf{w}}^t-({\mathbf{X}}^{\mathsf{T}}{\mathbf{S}}^t\mathbf{X}{)}^{-1}{\mathbf{g}}^t\text{\hspace{1em}(14)}$$

其中，${\mathbf{S}}^t=\mathrm{diag}({\mathbf{\mu }}^t\odot (1-{\mathbf{\mu }}^t))$。可以看出，式 (13) 的计算速度更快，因为我们可以预先计算常数矩阵 $({\mathbf{X}}^{\mathsf{T}}\mathbf{X}{)}^{-1}$。

源与流

期望最大化算法

我们将讨论期望最大化（Expectation Maximization, EM）算法。该算法针对具有缺失数据和/或隐藏变量的概率模型，用于计算最大似然估计或最大后验参数估计，是 MM 算法的一个特例。

期望最大化算法背后的基本思想是在期望步骤（E 步骤）期间估计隐藏变量（或缺失值），然后在最大化步骤（M 步骤）期间使用完全观察到的数据来计算最大似然估计。当然，我们需要迭代这个过程，因为期望值取决于参数，但参数取决于期望值。

我们将证明期望最大化算法是一个边界优化算法，这意味着该迭代过程将收敛到对数似然的局部最大值。收敛速度取决于缺失数据的数量，这会影响边界的紧密性。

接下来，我们将描述一个通用模型的期望最大化算法。设 $y_n$ 是样例 $n$ 的可见数据，而 $z_n$ 是隐藏数据。

下界

期望最大化算法的目标是最大化所观测数据的对数似然：

$$\ell (\mathbf{\theta })=\sum _{n=1}^N\log p(y_n\mid \mathbf{\theta })=\sum _{n=1}^N\log \left[\sum _{z_n}p(y_n,z_n\mid \mathbf{\theta })\right]\text{\hspace{1em}(15)}$$

其中，$y_n$ 是可见变量，$z_n$ 是隐藏变量。遗憾的是，这很难优化，因为无法将 $\log$（对数）推入求和中。

期望最大化算法通过以下方式解决了这个问题。首先，考虑每个隐藏变量 $z_n$ 上的一组任意分布 $q_n(z_n)$。观测数据对数似然可以记作：

$$\ell (\mathbf{\theta })=\sum _{n=1}^N\log \left[\sum _{z_n}{q}_n(z_n)\frac{p(y_n,z_n\mid \mathbf{\theta })}{q_n(z_n)}\right]\text{\hspace{1em}(16)}$$

使用詹森不等式，我们可以将 $\log$（这是一个凹函数）推入期望内，以获得以下对数似然的下界：

$$\begin{array}{rl}\ell (\mathbf{\theta })& \ge \sum _n\sum _{z_n}{q}_n(z_n)\log \frac{p(y_n,z_n\mid \mathbf{\theta })}{q_n(z_n)}\\ & =\sum _n{\mathbb{E}}_{q_n}\left[\log p(y_n,z_n\mid \mathbf{\theta })\right]+\mathbb{H}(q_n)\\ & =\sum _n\mathcal{L}(\mathbf{\theta },q_n\mid y_n)\triangleq \mathcal{L}(\mathbf{\theta },\{q_n\}\mid \mathcal{D})\end{array}\text{\hspace{1em}(17)(18)(19)}$$

其中，$\mathbb{H}(q)$ 是概率分布 $q$ 的熵，$\mathcal{L}(\mathbf{\theta },\{q_n\}\mid \mathcal{D})$，被称为证据下界，因为该值是对数边缘似然 $\log p(y_{1:N}\mid \mathbf{\theta })$（也称为证据）的下界。优化这个下界是变分推理的基础。

期望步骤

我们看到下界是 $N$ 个项的和，每个项都具有以下形式：

$$\begin{array}{rl}\mathcal{L}(\mathbf{\theta },q_n\mid y_n)& =\sum _{z_n}{q}_n(z_n)\log \frac{p(y_n,z_n\mid \mathbf{\theta })}{q_n(z_n)}\\ & =\sum _{z_n}{q}_n(z_n)\log \frac{p(z_n\mid y_n,\mathbf{\theta })p(y_n\mid \mathbf{\theta })}{q_n(z_n)}\\ & =\sum _{z_n}{q}_n(z_n)\log \frac{p(z_n\mid y_n,\mathbf{\theta })}{q_n(z_n)}+\sum _{z_n}{q}_n(z_n)\log p(y_n\mid \mathbf{\theta })\\ & =-D_{\mathrm{KL}}(q_n(z_n)\Vert p(z_n\mid y_n,\mathbf{\theta }))+\log p(y_n\mid \mathbf{\theta })\end{array}\text{\hspace{1em}(20)(21)(22)(23)}$$

其中，$D_{\mathrm{KL}}(q\Vert p)\triangleq {\sum }_{z}q(z)\log \frac{q(z)}{p(z)}$ 是概率分布 $q$ 和 $p$ 之间的 KL 散度。这里我们需要的关键性质是 $D_{\mathrm{KL}}(q\Vert p)\ge 0$，当且仅当 $q=p$ 时，$D_{\mathrm{KL}}(q\Vert p)=0$。因此，我们可以通过将每个值设置为 $q_n^{\ast }=p(z_n\mid y_n,\mathbf{\theta })$ 来最大化相对于 $\{q_n\}$ 的下界 $\mathcal{L}(\mathbf{\theta },\{q_n\}\mid \mathcal{D})$。这被称为期望步骤。这确保了证据下界是一个紧下界：

$$\mathcal{L}(\mathbf{\theta },\{q_n^{\ast }\}\mid \mathcal{D})=\sum _n\log p(y_n\mid \mathbf{\theta })=\ell (\mathbf{\theta }\mid \mathcal{D})\text{\hspace{1em}(24)}$$

为了了解这与绑定优化之间的联系，让我们定义：

$$Q(\mathbf{\theta },{\mathbf{\theta }}^t)=\mathcal{L}(\mathbf{\theta },\{p(z_n\mid y_n;{\mathbf{\theta }}^t)\})$$

于是，我们有 $Q(\mathbf{\theta },{\mathbf{\theta }}^t)\le \ell (\mathbf{\theta })$ 和 $Q({\mathbf{\theta }}^t,{\mathbf{\theta }}^t)=\ell ({\mathbf{\theta }}^t)$，这符合要求。

然而，如果不能精确地计算后验 $p(z_n\mid y_n;{\mathbf{\theta }}^t)$，我们仍然可以使用近似分布 $q(z_n\mid y_n,{\mathbf{\theta }}^t)$。这将产生对数似然性的非紧下界。期望最大化算法的这种广义形式被称为变分期望最大化。

最大化步骤

在最大化步骤中，我们需要最大化相对于 $\mathbf{\theta }$ 的 $\mathcal{L}(\mathbf{\theta },\{q_n^t\})$，其中 $q_n^t$ 是在迭代 $t$ 的期望步骤中计算的分布。由于熵项 $\mathbb{H}(q_n)$ 是相对于 $\mathbf{\theta }$ 的常数，我们可以在最大化步骤中省略这些常数。因此，只剩下：

$${\ell }^t(\mathbf{\theta })=\sum _n{\mathbb{E}}_{q_n^t(z_n)}\left[\log p(y_n,z_n\mid \mathbf{\theta })\right]\text{\hspace{1em}(25)}$$

这被称为期望的完全数据对数似然（expected complete data log likelihood）。如果联合概率在指数族中，我们可以将其重写为：

$${\ell }^t(\mathbf{\theta })=\sum _n(\mathbb{E}[\mathcal{T}(y_n,z_n){]}^{\mathsf{T}}\mathbf{\theta }-A(\mathbf{\theta }))=\sum _n(\mathbb{E}[\mathcal{T}(y_n,z_n){]}^{\mathsf{T}}\mathbf{\theta }-A(\mathbf{\theta }))\text{\hspace{1em}(26)}$$

其中，$\mathbb{E}[\mathcal{T}(y_n,z_n)]$ 被称为期望的充分统计量。

在最大化步骤中，我们最大化期望的完全数据对数似然，以获得：

$${\mathbf{\theta }}^{t+1}=\arg \underset{\mathbf{\theta }}{\max }\sum _n{\mathbb{E}}_{q_n^t}\left[\log p(y_n,z_n\mid \mathbf{\theta })\right]\text{\hspace{1em}(27)}$$

在指数族的情况下，最大化可以通过匹配期望充分统计量的矩以闭合形式求解。

从上面的讨论可以看出，期望步骤实际上不需要返回后验分布 $\{q(z_n)\}$ 的完整集合，而是可以只返回期望的充分统计量的总和 ${\sum }_n{\mathbb{E}}_{q(z_n)}[\mathcal{T}(y_n,z_n)]$。

期望最大化算法的一个常见应用是对混合物模型进行拟合。

示例：缺失数据多元正态分布的期望最大化算法

当我们有一个完整观察到的数据矩阵时，很容易计算多元正态分布的最大似然估计：我们只计算样本的均值和协方差。在本节中，我们将考虑缺失数据或部分观测（partially observed data）到的数据情况。例如，我们可以将 $\mathbf{Y}$ 的条目视为调查问卷的答案，其中一些答案可能是未知的。存在许多类型的缺失数据。在本节中，为了简单起见，我们做出了随机缺失（Missing At Random, MAR）假设。在随机缺失的假设下，可见数据的对数似然具有以下形式：

$$\log p(\mathbf{X}\mid \mathbf{\theta })=\sum _n\log p({\mathbf{x}}_n\mid \mathbf{\theta })=\sum _n\log \left[\int p({\mathbf{x}}_n,{\mathbf{z}}_n\mid \mathbf{\theta })d{\mathbf{z}}_n\right]\text{\hspace{1em}(28)}$$

其中，${\mathbf{x}}_n$ 是样例 $n$ 中的可见变量，${\mathbf{z}}_n$ 是隐藏变量，${\mathbf{y}}_n=({\mathbf{z}}_n,{\mathbf{x}}_n)$ 是所有变量。遗憾的是，这个目标很难最大化。因为我们无法将对数运算移动到期望值内。幸运的是，我们可以很容易地应用期望最大化运算，接下来将展开阐述。

期望步骤

假设我们有来自上一次迭代的参数 ${\mathbf{\theta }}^{t-1}$，那么，我们可以计算第 $t$ 次迭代的期望完全数据对数似然，如下所示：

$$\begin{array}{rl}Q(\mathbf{\theta },{\mathbf{\theta }}^{t-1})& =\mathbb{E}\left[\sum _{n=1}^N\log \mathcal{N}({\mathbf{y}}_n\mid \mathbf{\mu },\mathbf{\Sigma })|\mathcal{D},{\mathbf{\theta }}^{t-1}\right]\\ & =-\frac{N}{2}\log |2\pi \mathbf{\Sigma }|-\frac{1}{2}\sum _{n=1}^N\mathbb{E}\left[({\mathbf{y}}_n-\mathbf{\mu }{)}^{\mathsf{T}}{\mathbf{\Sigma }}^{-1}({\mathbf{y}}_n-\mathbf{\mu })\right]\end{array}\text{\hspace{1em}(29)(30)}$$

$$=-\frac{N}{2}\log |2\pi \mathbf{\Sigma }|-\frac{1}{2}\mathrm{tr}\left({\mathbf{\Sigma }}^{-1}\sum _n\mathbb{E}[({\mathbf{y}}_n-\mathbf{\mu })({\mathbf{y}}_n-\mathbf{\mu }{)}^{\mathsf{T}}]\right)\text{\hspace{1em}(31)}$$

$$=-\frac{N}{2}\log |\mathbf{\Sigma }|-\frac{ND}{2}\log (2\pi )-\frac{1}{2}\mathrm{tr}\left({\mathbf{\Sigma }}^{-1}\mathbb{E}[\mathbf{S}(\mathbf{\mu })]\right)\text{\hspace{1em}(32)}$$

其中：

$$\mathbb{E}[\mathbf{S}(\mathbf{\mu })]\triangleq \sum _n(\mathbb{E}[{\mathbf{y}}_n{\mathbf{y}}_n^{\mathsf{T}}]+\mathbf{\mu }{\mathbf{\mu }}^{\mathsf{T}}-2\mathbf{\mu }\mathbb{E}[{\mathbf{y}}_n{]}^{\mathsf{T}})\text{\hspace{1em}(33)}$$

（为了简洁起见，我们省略对 $\mathcal{D}$ 和 ${\mathbf{\theta }}^{t-1}$ 的期望条件。）我们发现，需要计算 ${\sum }_n\mathbb{E}[{\mathbf{y}}_n]$ 和 ${\sum }_n\mathbb{E}[{\mathbf{y}}_n{\mathbf{y}}_n^{\mathsf{T}}]$；这些都是期望的充分统计量。

于是：

$$p({\mathbf{z}}_n\mid {\mathbf{x}}_n,\mathbf{\theta })=\mathcal{N}({\mathbf{z}}_n\mid {\mathbf{m}}_n,{\mathbf{V}}_n)\text{\hspace{1em}(34)}$$

$${\mathbf{m}}_n\triangleq {\mathbf{\mu }}_h+{\mathbf{\Sigma }}_{hv}{\mathbf{\Sigma }}_{vv}^{-1}({\mathbf{x}}_n-{\mathbf{\mu }}_v)\text{\hspace{1em}(35)}$$

$${\mathbf{V}}_n\triangleq {\mathbf{\Sigma }}_{hh}-{\mathbf{\Sigma }}_{hv}{\mathbf{\Sigma }}_{vv}^{-1}{\mathbf{\Sigma }}_{vh}\text{\hspace{1em}(36)}$$

其中，我们基于隐藏和可见索引 $h$ 和 $v$ 将 $\mathbf{\mu }$ 和 $\mathbf{\Sigma }$ 划分为块。因此，期望的充分统计量为：

$$\mathbb{E}[{\mathbf{y}}_n]=(\mathbb{E}[{\mathbf{z}}_n];{\mathbf{x}}_n)=({\mathbf{m}}_n;{\mathbf{x}}_n)\text{\hspace{1em}(37)}$$

为了计算 $\mathbb{E}[{\mathbf{y}}_n{\mathbf{y}}_n^{\mathsf{T}}]$，我们使用结论 $\mathrm{Cov}[\mathbf{y}]=\mathbb{E}[\mathbf{y}{\mathbf{y}}^{\mathsf{T}}]-\mathbb{E}[\mathbf{y}]\mathbb{E}[{\mathbf{y}}^{\mathsf{T}}]$。因此：

$$\mathbb{E}[{\mathbf{y}}_n{\mathbf{y}}_n^{\mathsf{T}}]=\mathbb{E}\left[\left(\begin{array}{c}{\mathbf{z}}_n\\ {\mathbf{x}}_n\end{array}\right)\left(\begin{array}{cc}{\mathbf{z}}_n^{\mathsf{T}}& {\mathbf{x}}_n^{\mathsf{T}}\end{array}\right)\right]=\left(\begin{array}{cc}\mathbb{E}[{\mathbf{z}}_n{\mathbf{z}}_n^{\mathsf{T}}]& \mathbb{E}[{\mathbf{z}}_n]{\mathbf{x}}_n^{\mathsf{T}}\\ {\mathbf{x}}_n\mathbb{E}[{\mathbf{z}}_n{]}^{\mathsf{T}}& {\mathbf{x}}_n{\mathbf{x}}_n^{\mathsf{T}}\end{array}\right)\text{\hspace{1em}(38)}$$

$$\mathbb{E}[{\mathbf{z}}_n{\mathbf{z}}_n^{\mathsf{T}}]=\mathbb{E}[{\mathbf{z}}_n]\mathbb{E}[{\mathbf{z}}_n{]}^{\mathsf{T}}+{\mathbf{V}}_n\text{\hspace{1em}(39)}$$

最大化步骤

通过求解 $\nabla Q(\mathbf{\theta },{\mathbf{\theta }}^{(t-1)})=0$，我们可以证明最大化步骤等效于将这些期望的充分统计量插入通常的最大似然估计方程中，从而得到：

$${\mathbf{\mu }}^t=\frac{1}{N}\sum _n\mathbb{E}[{\mathbf{y}}_n]\text{\hspace{1em}(40)}$$

$${\mathbf{\Sigma }}^t=\frac{1}{N}\sum _n\mathbb{E}[{\mathbf{y}}_n{\mathbf{y}}_n^{\mathsf{T}}]-{\mathbf{\mu }}^t({\mathbf{\mu }}^t{)}^{\mathsf{T}}\text{\hspace{1em}(41)}$$

因此，我们看到期望最大化并不等同于简单地使用变量的期望值替换变量并应用标准最大似然估计公式，这将忽略后验方差并且将导致不正确的估计。相反，我们必须计算充分统计量的期望值，并将其插入最大似然估计的一般公式中。

初始化

为了启动算法，我们可以根据数据矩阵中完整观察到的行来计算最大似然估计。如果不存在这样的行，我们可以使用观察到的边缘统计量来估计 $\mathbf{\Sigma }$ 的对角项。然后我们准备启动期望最大化算法。

作为该程序的一个例子，让我们考虑一个缺失值插值问题，其中有 $N=100$ 个 $10$ 维数据样例，假设这些样例来自高斯分布。在我们生成的合成数据中，$50\%$ 的观测值是随机缺失的。首先，我们使用期望最大化算法拟合参数。将得到的参数称为 $\hat{\mathbf{\theta }}$。现在可以使用我们的模型通过计算 $\mathbb{E}[{\mathbf{z}}_n\mid {\mathbf{x}}_n,\hat{\mathbf{\theta }}]$ 来进行预测。使用学习参数获得的结果几乎与使用真实参数一样好。理所当然地，性能会随着数据的增加而提高，或者随着缺失数据的减少而提高。

示例：使用学生似然的稳健线性回归

我们将讨论如何使用期望最大化算法来拟合线性回归模型，为了实现鲁棒性，该模型使用学生分布来表示其似然，而不是更常见的高斯分布。更准确地说，似然为：

$$p(y\mid \mathbf{x},\mathbf{w},{\sigma }^2,\nu )=\mathcal{T}(y\mid {\mathbf{w}}^{\mathsf{T}}\mathbf{x},{\sigma }^2,\nu )\text{\hspace{1em}(42)}$$

乍一看，实现方法可能并不直观，因为没有缺失数据，也没有隐藏变量。然而，事实证明，我们可以引入“人为的”隐藏变量，使问题更容易解决。这是一个常见的处理技巧。其中的关键是，我们可以将学生分布表示为高斯缩放的混合模型。

可以将学生分布的高斯缩放混合模型应用于我们的问题，方法是将潜在缩放 $z_n\in {\mathbb{R}}_{+}$ 与每个样例联系起来。因此，完全的数据对数似然由下式给出：

$$\log p(\mathbf{y},\mathbf{z}\mid \mathbf{X},\mathbf{w},{\sigma }^2,\nu )=\sum _n-\frac{1}{2}\log (2\pi z_n{\sigma }^2)-\frac{1}{2z_n{\sigma }^2}(y_i-{\mathbf{w}}^{\mathsf{T}}{\mathbf{x}}_i{)}^2+\text{\hspace{1em}(43)}$$

$$\left(\frac{\nu }{2}-1\right)\log (z_n)-z_n\frac{\nu }{2}+\text{常数}\text{\hspace{1em}(44)}$$

忽略不涉及 $\mathbf{w}$ 的项，并考虑期望值，我们有：

$$Q(\mathbf{\theta },{\mathbf{\theta }}^t)=-\sum _n\frac{{\lambda }_n^t}{2{\sigma }^2}(y_n-{\mathbf{w}}^{\mathsf{T}}{\mathbf{x}}_n{)}^2\text{\hspace{1em}(45)}$$

其中，${\lambda }_n^t\triangleq \mathbb{E}[1/z_n\mid y_n,{\mathbf{x}}_n,{\mathbf{w}}^t]$。我们发现，这是一个加权最小二乘目标，每个数据点的权重为 ${\lambda }_n^t$。

我们现在讨论如何计算这些权重。可以证明：

$$p(z_n\mid y_n,{\mathbf{x}}_n,\mathbf{\theta })=\mathrm{IG}\left(\frac{\nu +1}{2},\frac{\nu +{\delta }_n}{2}\right)\text{\hspace{1em}(46)}$$

其中，${\delta }_n=\frac{(y_n-{\mathbf{w}}^{\mathsf{T}}{\mathbf{x}}_n{)}^2}{{\sigma }^2}$ 是标准化残差。因此：

$${\lambda }_n=\mathbb{E}[1/z_n]=\frac{{\nu }^{\prime }+1}{{\nu }^{\prime }+{\delta }_n^{\prime }}\text{\hspace{1em}(47)}$$

所以，如果残差 ${\delta }_n^{\prime }$ 很大，则该点将被赋予低权重 ${\lambda }_n^{\prime }$，这是显而易见的，因为这个点可能是一个异常值。

期望最大化算法的扩展

期望最大化算法有许多变体和扩展。接下来，我们将讨论其中的一些扩展。

变分期望最大化算法

假设在期望步骤中，我们选择 $q_n^{\ast }=\arg {\min }_{q_n\in \mathcal{Q}}{D}_{\mathrm{KL}}(q_n\Vert p({\mathbf{z}}_n\mid {\mathbf{x}}_n,\mathbf{\theta }))$。因为我们在函数空间上进行优化，所以这被称为变分推理（见 10.1 节）。如果分布族 $\mathcal{Q}$ 足够丰富，足以包含真正的后验，$q_n=p({\mathbf{z}}_n\mid {\mathbf{x}}_n,\mathbf{\theta })$，那么我们可以使 KL 散度为零。但一般来说，出于计算原因，我们可能会选择一个限制性更强的类别。例如，即使真正的后验是相关的，我们也可以使用 $q_n({\mathbf{z}}_n)=\mathcal{N}({\mathbf{z}}_n\mid {\mathbf{\mu }}_n,\mathrm{diag}({\mathbf{\sigma }}_n))$。

在期望最大化算法的期望步骤中，使用受限后验分布族 $\mathcal{Q}$ 称为变分期望最大化（variational EM）算法。与常规期望最大化算法不同，变分期望最大化算法本身不能保证增加实际对数似然，但算法确实单调增加了变分下界。我们可以通过改变变分分布族 $\mathcal{Q}$ 来控制这个下界的紧密性。在 $q_n=p_n$ 的极限中，对应于精确推理，我们化简为与正则期望最大化算法相同的行为。

硬期望最大化算法

假设我们在变分期望最大化算法的背景下使用退化后验分布近似，对应于点估计，$q(\mathbf{z}\mid {\mathbf{x}}_n)={\delta }_{{\hat{\mathbf{z}}}_n}(\mathbf{z})$，其中 ${\hat{\mathbf{z}}}_n=\arg {\max }_{\mathbf{z}}p(\mathbf{z}\mid {\mathbf{x}}_n)$。这相当于硬期望最大化（hard EM）算法，我们在期望步骤中忽略了 $z_n$ 的不确定性。

这种退化方法的问题是非常容易过拟合，因为潜在变量的数量与数据样例的数量成比例。

蒙特卡罗期望最大化算法

处理棘手的期望步骤的另一种方法是对期望的充分统计量使用蒙特卡罗近似。也就是说，我们从后验分布 ${\mathbf{z}}_n^s\sim p({\mathbf{z}}_n\mid {\mathbf{x}}_n,{\mathbf{\theta }}^t)$ 抽取样本，然后计算每个完整向量 $({\mathbf{x}}_n,{\mathbf{z}}_n^s)$ 的充分统计量，最后对结果进行平均。这被称为蒙特卡罗期望最大化（Monte Carlo EM, MCEM）算法 [WT90; Nea12]。

抽取样本的一种方法是使用马尔可夫链蒙特卡罗（Markov Chain Monte Carlo, MCMC）近似方法。然而，如果我们必须等待马尔可夫链蒙特卡罗在每个期望步骤内收敛，则该方法会变得非常缓慢。另一种可替代的方法是使用随机近似，只在期望步骤中执行“简洁”采样，然后进行部分参数更新。这被称为随机近似期望最大化（stochastic approximation EM）算法，并且往往比蒙特卡罗期望最大化算法工作得更好。

广义期望最大化算法

有时我们可以精确地执行期望步骤，但不能精确地执行最大化步骤。然而，我们仍然可以通过执行“部分”最大化步骤来单调增加对数似然，在该步骤中，我们只增加期望的完全数据对数似然而不是最大化该似然。例如，我们可以遵循几个梯度步骤。这被称为广义期望最大化（Generalized EM, GEM）算法。我们可以执行一个牛顿–拉弗森步骤：

$${\mathbf{\theta }}_{t+1}={\mathbf{\theta }}_t-{\eta }_t{\mathbf{H}}_t^{-1}{\mathbf{g}}_t\text{\hspace{1em}(48)}$$

其中，$0 < \eta_t \le 1$ 是步长，并且

$${\mathbf{g}}_t={\frac{\partial }{\partial \mathbf{\theta }}Q(\mathbf{\theta },{\mathbf{\theta }}_t)|}_{\mathbf{\theta }={\mathbf{\theta }}_t}\text{\hspace{1em}(49)}$$

$${\mathbf{H}}_t={\frac{{\partial }^2}{\partial \mathbf{\theta }\partial {\mathbf{\theta }}^{\mathsf{T}}}Q(\mathbf{\theta },{\mathbf{\theta }}_t)|}_{\mathbf{\theta }={\mathbf{\theta }}_t}\text{\hspace{1em}(50)}$$

如果 ${\eta }_t=1$，我们称之为梯度期望最大化（gradient EM algorithm）算法。然而，可以使用更大的步长来加快算法，我们提出的拟牛顿期望最大化（quasi-Newton EM algorithm）算法。该方法还用 BFGS 近似代替了式 (50) 中的黑塞矩阵，该项可能不是负定的（对于非指数族模型）。这确保了整个算法是上升算法。然而，请注意，当最大化步长不能以闭合形式计算时，期望最大化算法在使用基于梯度的求解器直接优化边缘似然方面失去了一些吸引力。

ECM算法

ECM 算法代表“期望条件最大化”（Expectation Conditional Maximization）指的是在最大化步骤中，如果各个参数是相关的，则按顺序对这些参数进行优化。ECME 算法代表“期望条件最大化选择其一”（ECM Either，是 ECM 的一种变体，在一个或多个条件最大化步骤中，我们像通常一样最大化期望的完全数据对数似然（Q 函数），或者最大化观测的数据对数似然。对观测的数据对数似然进行最大化可能运算速度更快，因为这种方法忽略了期望步骤的结果，并直接优化了感兴趣的目标。这方面的一个标准示例是拟合学生分布。对于固定的 $\nu$，我们可以像往常一样更新 $\mathbf{\Sigma }$，但为了更新 $\nu$，我们将形式为 ${\nu }^{t+1}=\arg {\max }_{\nu }Q(({\mathbf{\mu }}^{t+1},{\mathbf{\Sigma }}^{t+1},\nu ),{\mathbf{\theta }}^t)$ 的标准更新替换为 ${\nu }^{t+1}=\arg {\max }_{\nu }\log p(\mathcal{D}\mid {\mathbf{\mu }}^{t+1},{\mathbf{\Sigma }}^{t+1},\nu )$。

在线期望最大化算法

当处理大型或者流式数据集时，能够在线学习很重要，有两种主要的在线期望最大化（online EM）算法方法。第一种方法称为增量期望最大化（incremental EM）算法，优化下界 $Q(\mathbf{\theta },q_1,\dots ,q_N)$，每次一个 $q_n$。然而，这需要为每个数据样例存储期望的充分统计量。

第二种方法称为逐步期望最大化（stepwise EM）算法，基于随机梯度下降。这种方法优化了每个步骤中 ${\ell }_n(\mathbf{\theta })=\log p({\mathbf{x}}_n\mid \mathbf{\theta })$ 的局部上界。