爱喝兽奶帝天荒

这个屌丝很懒，什么也没留下！

热门标签

机器学习：手撕 EM 算法_手撕em算法

作者：爱喝兽奶帝天荒 | 2024-07-04 03:26:23

踩

手撕em算法

本文首次发表于知乎，欢迎关注作者。

1. 前言

EM 算法是求解带有隐变量的最大似然估计的方法。常见的需要 EM 求解的模型有 GMM，HMM。甚至 K-mean 的迭代过程，也是 EM 的一个特例。本文尝试说明 EM 算法的基本原理，以及它在几个不同模型上的应用。学习 EM 的过程中比较曲折，经常出现” 翻开书，马冬梅；合上书，马冬啥” 的现象。对 EM 的理解一直不够深入，记得曾经老师对 EM 总结时说：EM 算法就像古代武林中的轻功，高手们左脚踩右脚，右脚踩左脚便不断的拔地而起了。当时理解不了这句话。整理出这篇文章，希望加深对 EM 的理解。文本尝试说明：

隐变量给 MLE 带来哪些困难；
EM 算法的主要原理，包括 EM 算法的动机，以及如何求解 E 步和 M 步；
解释 EM 算法中常见两种 q（z）形式为什么是等价；
用 EM 如何求解 GMM 模型；
解释为什么 k-meams 是 GMM 的特例；
用 EM 如何求解 HMM 模型；

如果不知道如何从 EM 的角度解释封面这张图，那么读完本文后，便知道封面这张图的含义了。

2. 隐变量给 MLE 带来的困难

在机器学习任务中，常常会遇到概率密度估计的任务：即给我们一个数据集 $D$ ，其中数据集中每个元素 x ∈ $R^d$ ，我们需要找到一个模型和模型的一组参数 θ，可以很好的描述数据集的概率密度 P（x|θ）。MLE（MaximumLikelihood Estimation）是解决这类问题的方法之一，MLE 将密度估计问题看做一个优化问题，然后搜索到模型的一组参数，使模型对数据的拟合达到最好：

θ* 表示最优的参数，l(D|θ) 表示数据集 D 的对数似然函数。实际求解时会对式 2-1 取负号变为 NLL（negative log likelihood），从而将最大化转化为最小化的问题，然后采用基梯度（gradient-based）的优化方法，求解出式 2-1 的一个局部极小点。这里为了叙述方便，我们以 MLE 的原始形式说明。

但是 MLE 有一个天然的限制，它适用于“完全可观测”的数据集，即数据集中的元素 x ∈ $R^d$ 包含了与问题相关的所有变量。与“完全可观测”数据集相对应的是“部分可观测”的数据集，即与问题相关的所有变量，只有部分变量是可观测的，存在一部分变量不可观测，其中不可观测的变量叫做隐变量。

当 MLE 遇到带有隐变量 z（ z∈ $R^d$ ）的数据集 D 时，优化问题转化为如下形式：

式 2-2 相当于引入隐变量 z，然后用全概率公式消掉参数 z。相比于式 2-1，在式 2-2 中 log 里面多出了对 z 的求和项（若 z 是连续变量，则是对 z 的积分），从而形成 logsum 的形式。因为 logsum 的出现，使式 2-2 不论在求和还是在=求导上，都会变得不方便计算。若是在计算过程中，不考虑隐变量 z，虽然避免了计算上的不便，但对于明确存在隐变量的问题，在求解时不去考虑隐变量，那么求出的 p（x|θ）效果也不会太好。

面对如式 2-2 这种不容易计算的 MLE 公式时，采用 EM 算法可以很方便的解决这个问题。

3. EM 算法的主要原理

3.1 EM 算法背后的动机

在介绍 EM 算法前，把式 2-2 中的 log-likelihood 提取出来：

为了方便说明，我们从数据集 D 中抽取一条数据 x，计算一条数据的 log-likelihood:

现在，我们假设存在一个关于隐变量 z 的概率质量分布 q（z）（如果是 z 是连续变量，则是概率密度函数），我们在式 3-4 中同时乘一个 q（z）和除一个 q（z），不改变结果：

在式 3-5 中 log 后面的部分可以看做是变量 $\frac{P(x,z|\theta )}{q(z)}$ 关于分布 q（z）的期望。这时，可以借用 Jenson’s inequality 在凹函数下形式"期望的函数大于等于函数的期望"，便可以将 log 放到 ∑ 内部：

式 3-6 中红色的部分有一个专有名字 ELBO（Evidence lower bound），即：

很容易看出 ELBO 是 $l(x|\theta )$ 的一个下界，即：

我们知道求解如式 3-4 这种 logsum 的形式比较困难，但求解式 3-7 这种 sumlog 的形式相对容易些。所以我们不再直接求解 log p（x|θ），而是通过最大化 $l(q|\theta )$ ，近似最大化 log p（x|θ）的值。这样优化问题也由原来：

转化为如下形式：

式 3-10 将最大化 ELBO 转化为两个子问题。第一个子问题，在 θ 固定的情况下，寻找一个 q（z）使 $L(q|\theta )$ 最大，这步这称作 Expectation Step，即 E 步：

第二个子问题固定 q（z），寻找一个 θ 使 $L(q|\theta )$ 最大，这步称作 Maximization Step，即 M 步：

这两个子问题不断循环迭代，直到 ELBO 收敛。

3.2 如何求解 q（z）

在求解 q（z）前，我们先进一步对 $L(q|\theta )$ 进行展开。由概率相关知识可知：

于是我们可以对 ELBO 进一步推导：

通过式 3-14 我们可以得到 $L(q|\theta )$ 和 log p（x|θ）的等式关系：

通过散度的定义我们知道 KL ≥ 0，结合式 3-15 我们知道当 KL = 0 时：

此时 $L(q|\theta )$ 取得最大值。由散度的相关知识可知，当

时，有：

这样我们便求解出 q（z）。

3.3 如何求解 θ

接下来在 $\hat{q}(z)$ 已知的情况下，如何求解 θ，为了方便区分我们将 $\hat{q}(z)$ 的参数用 $\theta ^{old}$ 表示。把 $L(q|\theta )$ 以另外一种方式展开：

式 3-18 第三行的第二项是 $\hat{q}(z)$ 的熵，它是一个常数量，与参数 θ 无关，所以可以舍去，得到第四行的表达式。第四行是“完全可观测”数据 log-likelihood 的表达式，所以可以用基于梯度的优化方法求解 θ。

3.4 EM 算法思路总结

首先带有隐变量的 log-likelihood 难以计算，EM 的思想是找到 log-likelihood 的下界 ELBO，然后通过最大化 ELBO，间接增大 log-likelihood。然后在增大 ELBO 的过程中，通过两步迭代计算。通过上面的叙述，我们可以看到，在 E 步：

式 3-19 是在 θ， x 已知的条件下，对隐变量的推断，这是推断过程。同理当已知 q（z），在 M 步我们求得：

式 3-20 是求解参数 θ 的过程，而且在求解过程中，必须求解“完全可观测”数据的 log-likelihood，这个过程叫做学习过程。其中我们将式 3-20 中的最大化项定义为 Q 函数，即：

式 3-21 表示“完全可观测”数据的 log-likelihood 相对于分布 q（z）的期望。

在有的文献中，会看到 $\hat{q}(z) = p(z,x| \theta )$ ，其实这种表述与式 3-19 等价。若将 ELBO 直接看成散度：

最大化式 3-22 便可以得到 $\hat{q}(z) = p(z,x|\theta)$ 。相应的 Q 函数变为

式 3-14 的推导过程相比式 3-22，不仅可以求出 q（z），而且还可以得到 ELBO 与 p（x|θ）的定量关系（如式 3-15 所示）。后面我们以 $\hat{q}(z) = p(z,x|\theta)$ 的表达形式说明。

图 3-1 经常出现在 EM 相关文档中，我们也借助图 3-1 来说明下 EM 过程。图中，横坐标是 θ 值，纵坐标是对数似然值。

4. EM 求解 GMM 模型

4.1 GMM 介绍

GMM（Gaussian Mixture Model）是由多个高斯分布组成的混合分布，我们需要求它的概率密度。它的参数有三类：

每个高斯模型占整体的比重，这是一个标量： $\pi=\{\pi_1,\pi_2,\pi_3,...,\pi_K\}$ ；
每个高斯模型的均值，它的维度和 x 的维度相同： $\mu=\{\mu_1,\mu_2,\mu_3,...,\mu_K\}$ ；
每个高斯模型的协方差矩阵： $\Sigma =\{\Sigma _1,\Sigma _2,\Sigma _3,...,\Sigma _K\}$ ；

其中 K 是高斯混合模型的个数。在 GMM 中，可观测数据是 x ∈ $R^d$ ，其中数据 $x_i$ 对应的隐变量的分布 q（z）记作 $r_i$ ，表示数据 $x_i$ 来自每个高斯分布的概率分布：

根据 $r_i^k$ 的定义，对于每一个数据点，我们有：

我们需要学习出一组参数 π， µ， Σ，估计出整体概率密度：

有可观测数据集 $D=\{x_1,x_2,x_3,...,x_K\}$ ，我们先计算下 GMM 的 log-likelihood:

可以看到 GMM 的 log-likelihood 式 4-27 化简出 log − sum 的形式，不方便直接计算。

4.2 求解 GMM

4.2.1 E 步

我们采用 EM 算法求解 GMM。根据 EM 算法，先初始化一组参数 $\theta ^0=\{\pi ^0,\mu ^0,\Sigma ^0\}$ ，先求隐变量的分布 q（z），E 步：

对应到 GMM 模型中，是数据 $x_i$ 来自各个高斯分布的概率分布 $r_i$ 。比如求来自第 k 个高斯分布的概率 $r_i^k$ :

在已知 $\theta ^{t-1}$ 时，由式 4-29 可以很容易的求出 $r_i$ 。

4.2.2 M 步

接下来我们计算 M 步，即已知 q（z）时，求解参数 θ：

对应到 GMM 中是已知 $r_i$ ，求解 {π， µ， Σ}，我们先计算 GMM 的完全数据的 log-likelihood 的期望 Q，为了更好的理解 GMM 解的含义，我们记作 Q 为整个数据集 D 上的完全数据的 log-likelihood 的期望：

因为数据集 D 中，任意 2 条数据相互独立，他们在求 log-likelihood 时没有耦合，所以对他们的“对数似然的和求期望”，等价于“分别求期望后的和”，于是第一行的 $\Sigma _i$ 可以提到最左侧。对于 GMM 模型，当知道完全数据（x， z）时，完全数据的概率：

所以有了式 4-31 中第三行到第四行的变换。又因为 $p(z_i =k|x_i,\theta _{t-1})$ 已经在 E 步计算求出，直接带入，便得到式 4-31 的第五行，然后按照 log 的运算法则，得到式 4-31 最后一行。

求参数 π 通过式 4-31 可以看到参数 π 和参数 µ， Σ 相互独立，所以我们可以分别求解。由式 4-30，我们可以将 $\pi _k$ 转化下面这个优化问题的解：

式 4-33 的拉格朗日函数：

其中 λ 是拉格朗日乘子，然后分别对 $\pi _k$ ，求导，并令导数为 0:

将式 4-35 变形得：

将式 4-36 前 K 个方程等号左右两边分别相加：

然后结合式 4-25，可以求出：

其中 N 为数据总数。将式 4-38 代入 4-36，求出 $\pi _k$ :

求参数 $\mu _k$ ， $\Sigma _k$ ：因为参数 $\mu _k$ ， $\Sigma _k$ 没有额外的约束，比较容易求解，由式 4-31，可以将 $\mu _k$ ， $\Sigma _k$ 转化为下面这个优化问题的解（忽略常数项）：

分别对 $\mu _k$ ， $\Sigma _k$ 求偏导数，并令偏导数为 0：

在求解 $\Sigma _k$ 的偏导数时，用到了矩阵代数的相关知识，具体可以参考 [1]。求解式 4-41 可得：

总结一下，在第 t 轮迭代中，GMM 的参数求解公式如下：

在式 4-43 中，包含和 E 步和 M 步，每次新求的参数用于下次迭代。

4.3 K-means 是 GMM 的特例

本小节主要解释为什么人们常说“K-means 是 GMM 的特例”。我们先回顾下 K-means 的计算过程：

现在我们给 GMM 模型填加一个约束：每个高斯的形状为“圆形”，即在每个维度上数据分散程度一致，且具有相同的协方差矩阵，即 $\Sigma =\delta ^2I$ 。然后我们让 $\delta ^2$ → 0，直观上每个高斯分布的概率质量更集中，如图 4-2。当 $\delta ^2$ 小到一定程度，会出现每个点在某个特定的高斯分布上概率特别大，在其他高斯分布上概率很小，几乎为 0，此时 4-43 中的 $r_i^k$ 退化为 one-hot 向量：

当 $r_i^k$ 为 one-hot 时，式 4-43 其他参数求解为：

式 4-45 中 $n_k$ 表示落在第 k 个高斯分布中的点的个数，N 表示所有点的个数。 $\mu _k$ 退化为第 k 个高斯分布中点的均值（重心）。这时 EM 算法中的 E 步，即式 4-44 中 $r_i^k$ 的求解过程对应 Kmeans 算法的“计算每个样本所属的类别”，EM 算法的 M 步，即式 4-45 中 $\mu _k$ 的求解过程对应 Kmeans 算法的“更新类别中心”。此时受约束的 GMM 退化为 Kmeans 算法，所以说 Kmeans 算法是 GMM 算法的特例。

5. EM 求解 HMM 模型

5.1 HMM 基本概念

首先我们来回顾下 HMM 的基本概念。HMM 的模型如图 5-3 所示。HMM 是关于随机变量 {O1， O2， O3， ..., OT , Q1, Q2, Q3, ...， QT } 的联合概率密度模型。其中 Qt 表示隐变量，它是离散型数值。Ot 是观测变量，它既可以是离散的数值也可以是连续的数值。HMM 引入 2 个假设，使问题变得简化，容易求解：

在以上两种假设下，HMM 可以用三组参数表示：

转移概率矩阵 $A=\{a_{i,j}\}=p(Q_t=j|Q_{t-1}=i)$ ，表示 2 个隐变量的转移概率。若 Q 的取值有 N 种，则 A 是一个 N × N 的矩阵；
初始状态概率向量 $\pi =\{p(Q_1=i)\}$ ，在 t =0 时刻，隐变量每个状态的概率。其中 π 的维度为 N；
观测概率矩阵 $B=\{b_j(o_t)\}=p(O_t=o_t|Q_t=j)$ 。表示在 t 时刻，隐变量为 j 的条件下，观测变量为 $o_t$ 的概率。其中 $O_t$ 既可以是离散值也可以是连续值，为了方便说明问题，这里我们规定 $O_t$ 为离散值，且维度为 M，则 B 是 N × M 的矩阵；

这样 HMM 表示的概率密度可以写成：

因为隐变量序列 q 是 T 维度的向量，所以式中第二行是每个维度都展开后的表示，方便进一步理解。围绕 HMM 模型有三个基本的问题：

概率计算。给定模型参数 $\lambda=\{A,B,\pi \}$ 和观测数据 $O=\{o_1,o_2,...o_T\}$ ，求观测数据的概率 $p(O|\lambda )$ ；
学习模型参数。已知观测数据 $O=\{o_1,o_2,...o_T\}$ ，学习一组模型参数 $\lambda =\{A,B,\pi \}$ ，使观测数据 O 出现的概率 $p(O|\lambda )$ 最大；
预测隐变量。给定模型参数 λ 和观测序列 O，求出最有可能的隐变量序列 $q=\{q_1,q_2,q_3,...q_T\}$ ；

本文我们主要讲 EM 算法，所以我们只关注问题 2，看看通过 EM 算法，如何学习 HMM 的参数。

5.2 求解 HMM

5.2.1 E 步

我们先初始化一组参数 $\lambda ^0=\{A^0,B^0,\pi ^0\}$ 。然后根据式 3-19 求解隐变量序列分布 q（z）。因为在 HMM 中我们已经用 q 表示隐变量序列了，所以我们改用 p（q）表示隐变量序列 q 的分布（注意这里的 q 是 T 维的向量）：

当知道前一时刻的参数 λ′ 和观测数据 O，对于指定的任意一个隐变量序列 q，很容易求出概率 p（q|λ′， O）：

这样我们完成 E 步的求解，接下来求解 M 步。

5.2.2 M 步

由式 3-20 和式 3-21。我们先计算 HMM 的 Q 函数：

在式 5-51 中：

第一行为 Q 函数的定义，可由式 3-21 推出；
将式 5-48 带入第一行，便得到第二行的形式；
根据 log 的运算法则，便得到第三行；
将第三行的括号打开，便得到第四行的表达形式；

通过式 5-51，我们看到 Q 函数可以表示为 3 个独立项的和，对 Q 函数求最大，可以分别表示为对每一项单独求最大。

求初始状态 π : 因为 π 表示 $q_0$ 的取值概率，所以 π 的维度和 $q_0$ 的取值数量相同为 N。其中式 5-51 第一项可以化简为：

又因为 $\Sigma ^N_{k=1}\pi _k=1$ ，所以我们将求解 $\pi _k$ 转化为下面这个优化问题：

式 5-53 的拉格朗日函数：

其中式 $\alpha$ 是拉格朗日乘子，然后分别对式 5-54 中的 $\pi _i$ ，求导，得：

变形可得：

将式 5-56 的前 N 个方程左右两边分别相加，结合 $\Sigma ^N_{k=1}\pi _k=1$ 可得：

将式 5-57 代入式 5-56，求得：

由式 5-58 和式 5-59 便可以计算出所有的 π 值：

求转移矩阵 $A=\{a_{i,j}\}$ 我们化简式 5-51 中第二项（关于 A 的项）:

在式 5-59 中 $p(q_{t-1}=i,q_t=j|\lambda ',O)$ ，可以通过全概率公式消掉隐变量序列 $q_{1:t-2},q_{t+1:T}$ ，求解：

另外结合约束 $\Sigma ^N_j a_{i,j}=1$ ，得到关于 A 的优化问题：

式 5-62 的拉格朗日函数：

同理，然后分别对式 5-63 中的 $a_{i,j}$ ， $\alpha _i$ 求导，得：

求解方法与式 5-56 相同，求解式 5-64 可得：

通过式 5-65 可以求出状态转移矩阵 A。

求观测矩阵 $B=\{b{q_t}(o_t)\}$ 我们化简式 5-51 中第三项（关于 B 的项）:

根据 B 的定义，我们易知 $\Sigma ^M_j b_i(j)=1$ ，所以关于矩阵 B 的优化问题是：

同理式 5-67 的拉格朗日函数：

然后分别对式 5-68 中的 $b_i(j)$ ， $\alpha _i$ 求导，其中注意式中的红色字体部分 $o_t$ 是观测数据 O 在 t 时刻的状态，所以在求导时，若 $b_i(o_t)$ 与被求导变量 $b_i(j)$ 一致，则导数为 1，否则为 0，用指示函数 $\Pi \{b_i(j)=b_i(o_t)\}$ 表示。综上得：

求解式 5-69 可得（解方法与式 5-56 相同，这里省了）：

在式 5-70 中，在求解 $\alpha _i$ 时，如红色字体所示，对于每一个时刻 t，都存在一个 j 使 $b_i(j)=b_i(o_t)$ ，所以指示函数 $\Pi$ 可以去掉，得到 $\alpha _i$ 的最终表达式。至此我们求解出参数 B。然后我们通过 EM 算法不断的迭代，直到算法收敛，HMM 的所有参数 λ = {A， B， π} 将全部求解出来。

总结在第 k 轮 EM 迭代中，HMM 的参数求解公式如下：

6. 结束

本文主要探讨了 EM 算法的背景和原理。并应用 EM 算法求解机器学习中常用的 GMM 和 HMM 模型。本文侧重 EM 算法在 GMM 和 HMM 上的应用，希望大家对 EM 的原理有更进一步的理解，对 GMM 和 HMM 没有过细的分析。作者能力有限，文中难免存在理解不正确，或者描述不清的地方。欢迎留言讨论。

参考文献

[1] Bilmes: A Gentle Tutorial of the EM Algorithm and its Application to

Parameter Estimation for Gaussian Mixture and Hidden Markov Models

http://imaging.mrc-cbu.cam.ac.uk/methods/BayesianStuff?action=AttachFile&do=get&target=bilmes-em-algorithm.pdf

[2] Yida.Xu: Expectation Maximization

roboticcam/machine-learning-notes

[3] LongMingsheng <deep learning> lecture

Mingsheng Long - Tsinghua University

[4] David Rosenberg: Expectation Maximization Algorithm

https://davidrosenberg.github.io/mlcourse/Archive/2016/Lectures/14a.EM-algorithm.pdf

[5] What is the difference between K-means and the mixture model of

Gaussian?

https://www.quora.com/What-is-the-difference-between-K-means-and-the-mixture-model-of-Gaussian

[6] Murphy: Machine Learning A Probabilistic Perspective

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/爱喝兽奶帝天荒/article/detail/785487