【作者主页】Francek Chen
【专栏介绍】$\lceil$Python机器学习$\rfloor$ 机器学习是一门人工智能的分支学科，通过算法和模型让计算机从数据中学习，进行模型训练和优化，做出预测、分类和决策支持。Python成为机器学习的首选语言，依赖于强大的开源库如Scikit-learn、TensorFlow和PyTorch。本专栏介绍机器学习的相关算法以及基于Python的算法实现。
【GitCode】专栏资源保存在我的GitCode仓库：https://gitcode.com/Morse_Chen/Python_machine_learning。

---

  本文我们继续介绍概率相关模型的算法。在前面的文章中，我们已经讲解了贝叶斯网络与最大似然估计（MLE）。设数据集的样本为${\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_N$，标签为$y_1,y_2,\cdots ,y_N$，我们用$\mathbf{X}$和$\mathbf{y}$分别表示全体样本和全体标签。对于监督学习的任务，我们可以通过最大化对数似然 $l(\mathbf{w})=\log P(\mathbf{y}|\mathbf{X},\mathbf{w})$ 来求解模型的参数$\mathbf{w}$。而对于无监督学习任务，我们要求解样本$\mathbf{X}$的分布。这时，我们通常需要先假设数据服从某种分布，再求解这个分布的参数。例如假设数据呈现高斯分布$\mathcal{N}(\mathbf{\mu },\mathbf{\Sigma })$，然后可以通过MLE的方式来求得最佳的高斯分布参数$\mathbf{\mu }$和$\mathbf{\Sigma }$。

  更进一步思考，真实世界的数据分布往往较为复杂，其背后往往具有一定的结构性，直接使用一个概率分布模型无法有效刻画数据分布。例如，我们可以假设数据服从高斯混合模型（Gaussian mixture model，GMM），即$\mathcal{N}({\mathbf{\mu }}_1,{\mathbf{\Sigma }}_1,\cdots ,{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)$，该模型是由$k$个相互独立的高斯分布 ${\mathcal{N}}_i({\mathbf{\mu }}_i,{\mathbf{\Sigma }}_i)(i=1,\cdots ,k)$ 组合而成的，数据集中的每个样本${\mathbf{x}}_j$都从其中的某个高斯分布采样得到。在现实生活中，符合GMM的数据集有很多。例如，我们统计了某学校中所有学生的身高。通常认为，人的身高是在某个均值附近的高斯分布，然而男生和女生身高的均值是不同的。因此，我们可以认为男生身高和女生身高分别服从不同的高斯分布，而总的数据集就符合GMM。

  在GMM中，我们要求解的参数共有两种，一种是每个高斯分布的参数${\mathbf{\mu }}_i$和${\mathbf{\Sigma }}_i$，另一种是每个高斯分布在GMM中的占比。记 $z\in 1,\cdots ,k$ 是高斯分布的编号，$z$出现的次数越多，从分布$\mathcal{N}({\mathbf{\mu }}_z,{\mathbf{\Sigma }}_z)$采样的数据在数据集中的占比就越大。所以，后者相当于求解$z$的多项分布$p(z)$。

  从贝叶斯网络的角度来看，上面的分析过程建立了如下图1所示的贝叶斯网络。其中，$\mathbf{\varphi }$是多项分布$p(z)$的参数，$z$取 $1,\cdots ,k$ 的概率分别是 ${\varphi }_1,\cdots ,{\varphi }_k$。而对每个样本${\mathbf{x}}_i$，我们先从多项分布中采样${\mathbf{z}}_i$，确定样本属于哪个高斯分布，再从该高斯分布$\mathcal{N}({\mathbf{\mu }}_{z_i},{\mathbf{\Sigma }}_{z_i})$中采样出样本${\mathbf{x}}_i$。于是，我们可以利用中间变量$z$把$\mathbf{x}$的分布拆分为 $p(\mathbf{x})=p(\mathbf{x}|z)p(z)$。

图1 高斯混合模型的贝叶斯网络

  按照MLE的思想，参数的似然就是在此参数条件下出现观测到的数据分布的概率，也即 $L(\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })=P(\mathbf{X}|\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })$。我们将似然取对数，得到
$$\begin{array}{rl}l(\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })& =\log L(\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })=\log P(\mathbf{X}|\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })\\ & =\log \prod _{i=1}^{N}P({\mathbf{x}}_i|\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })\\ & =\sum _{i=1}^N\log P({\mathbf{x}}_i|\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })\\ & =\sum _{i=1}^N\log \sum _{j=1}^{k}P({\mathbf{x}}_i|z=j,{\mathbf{\mu }}_j,{\mathbf{\Sigma }}_j)P(z=j|{\varphi }_j)\\ & =\sum _{i=1}^N\log \sum _{j=1}^k\mathcal{N}({\mathbf{x}}_i|{\mathbf{\mu }}_j,{\mathbf{\Sigma }}_j){\varphi }_j\end{array}$$ 为了求出使对数似然最大化的参数，我们需要令对数似然对3个参数的梯度均为零：$${\nabla }_{\mathbf{\varphi }}l(\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })=0,{\nabla }_{\mathbf{\mu }}l(\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })=0,{\nabla }_{\mathbf{\Sigma }}l(\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })=0$$

  然而，上述方程的求解非常复杂，并且没有解析解。因此，我们需要寻找其他方法求解。本文介绍的EM算法是一种求解复杂分布参数的通用算法。

一、高斯混合模型的EM算法

  我们在最开始为了拆分$\mathbf{X}$的分布引入了中间变量$z$，用来指示每个样本属于哪个高斯分布，而$z$的分布就是混合高斯分布中每个分布的占比。像这样虽然不能直接被观测到、但是可以直接影响最后观测结果的变量，就称为隐变量（latent variable）。通常来说，隐变量比可观测的变量更本质。因此，我们经常会用引入隐变量的方法来把复杂的问题简单化。同时，隐变量也对问题的求解有关键作用。在对数似然的表达式中，如果每个样本${\mathbf{x}}_i$所属的高斯分布编号$z_i$已知，那么推导的倒数第二步中 $P({\mathbf{x}}_i|z=j,{\mathbf{\mu }}_j,{\mathbf{\Sigma }}_j)$ 一项在 $j\ne z_i$ 时就变为0。因此，对数似然可以写为
$$\begin{array}{rl}l(\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })& =\sum _{i=1}^N\log P({\mathbf{x}}_i|z=z_i,{\mathbf{\mu }}_i,{\mathbf{\Sigma }}_i)P(z=z_i|\mathbf{\varphi })\\ & =\sum _{i=1}^N(\log \mathcal{N}({\mathbf{x}}_i|{\mathbf{\mu }}_i,{\mathbf{\Sigma }}_i)+\log {\varphi }_{z_i})\end{array}$$ 求这一函数的最大值就相对容易了。首先，由于$z$已知，其多项分布的参数$\mathbf{\varphi }$可以直接通过统计数据集中属于每个高斯分布的样本比例得到：
$$\begin{array}{rl}{N}_i& =\sum _{j=1}^N\mathbb{I}(z_j=i)\\ {\varphi }_i& =\frac{N_i}{N}\end{array}$$ 其中，$N_i$表示数据集中属于第$i$个高斯分布的样本数目。每个高斯分布之间的参数是相互独立的，并且在$z$已知的条件下，它们也和$\mathbf{\varphi }$独立。这一点从图1的贝叶斯网络示意图中也可以看出来，当中间的节点$z$已知时，左右两端的节点之间不存在依赖关系。高斯分布的参数分别是均值和协方差，也可以从数据集中属于该分布的样本直接计算出来：$${\mathbf{\mu }}_i=\frac{1}{N_i}\sum _{j=1}^N\mathbb{I}(z_j=i){\mathbf{x}}_j$$ $${\mathbf{\Sigma }}_i=\frac{1}{N_i-1}\sum _{j=1}^N\mathbb{I}(z_j=i)({\mathbf{x}}_j-{\mathbf{\mu }}_i{)}^{\mathrm{T}}({\mathbf{x}}_j-{\mathbf{\mu }}_i)$$ 求和中每一项的因子 $\mathbb{I}(z_j=i)$ 是为了筛选出属于第$i$个高斯分布的样本。如果样本${\mathbf{x}}_j$不属于分布$i$，那么这一项就为0，对求和没有贡献。如果我们已知的是隐变量$z$的分布$q(z|\mathbf{x})$，也即第$j$个样本${\mathbf{x}}_j$属于第$i$个高斯分布的概率是 $q(z=i|{\mathbf{x}}_j)$，那么上面3个计算公式也相应的改为
$$\begin{array}{rl}{\varphi }_i& =\frac{1}{N}\sum _{j=1}^{N}q(z=i|{\mathbf{x}}_j)\\ {\mathbf{\mu }}_i& =\frac{1}{N{\varphi }_i}\sum _{j=1}^{N}q(z=i|{\mathbf{x}}_j){\mathbf{x}}_j\\ {\mathbf{\Sigma }}_i& =\frac{1}{(N-1){\varphi }_i}\sum _{j=1}^{N}q(z=i|{\mathbf{x}}_j)({\mathbf{x}}_j-{\mathbf{\mu }}_i{)}^{\mathrm{T}}({\mathbf{x}}_j-{\mathbf{\mu }}_i)\end{array}$$

  反过来考虑，如果我们已经知道了分布的参数$\mathbf{\varphi }$、$\mathbf{\mu }$和$\mathbf{\Sigma }$我们同样也可以推断出每个样本${\mathbf{x}}_i$属于第$j$个高斯分布的概率 $P(z_i=j|{\mathbf{x}}_i,\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })$。根据贝叶斯公式可以得到$z$的后验分布为
$$\begin{array}{rl}P(z_i=j|{\mathbf{x}}_i,\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })& =\frac{P({\mathbf{x}}_i|z_i=j,\mathbf{\mu },\mathbf{\Sigma })P(z_i=j|\mathbf{\varphi })}{P({\mathbf{x}}_i,\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })}\\ & =\frac{P({\mathbf{x}}_i|z_i=j,{\mathbf{\mu }}_j,{\mathbf{\Sigma }}_j)P(z_i=j|\mathbf{\varphi })}{\sum _{l=1}^{k}P({\mathbf{x}}_i|z_i=l,{\mathbf{\mu }}_l,{\mathbf{\Sigma }}_l)P(z_i=l|\mathbf{\varphi })}\\ & =\frac{\mathcal{N}({\mathbf{x}}_i|{\mathbf{\mu }}_j,{\mathbf{\Sigma }}_j){\varphi }_j}{\sum _{l=1}^k\mathcal{N}({\mathbf{x}}_i|{\mathbf{\mu }}_l,{\mathbf{\Sigma }}_l){\varphi }_l}\end{array}$$ 在各个参数已知的情况下，分子分母都可以直接计算出来。到此为止，看起来我们似乎在进行循环论证，在已知$z$的前提下可以计算出参数$\mathbf{\varphi }$、$\mathbf{\mu }$和$\mathbf{\Sigma }$，而在已知参数的前提下可以计算出$z$，但是这两者都是未知的。因此，我们可以采用之前用过许多次的近似方式，先固定其中一方来优化另一方，再反过来固定另一方，如此迭代，这就形成了EM算法的基本思想。具体在本例中，EM算法每一次迭代由两步构成：

• 期望步骤（E-step）：固定各个参数，由数据集中的样本统计计算隐变量$z$的后验分布$p(z|\mathbf{X},\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })$。

• 最大化步骤（M-step）：固定隐变量，最大化参数的对数似然$l(\mathbf{\varphi },\mathbf{\mu },\mathbf{\Sigma })$。

  由于两步分别在计算样本的统计期望和最大化参数的对数似然，这样交替优化隐变量和参数的方法就称为期望最大化算法（Expectation-maximization algorithm），简称EM算法。

二、动手求解GMM拟合数据分布

  虽然GMM是由高斯分布组成的，然而理论上它可以用来拟合任意的数据分布。如图2所示，左上角是由 $y=\sin x$ 加上随机的高斯噪声生成的数据，显然这些数据并不符合高斯分布。如果我们试图用两个高斯分布来拟合数据，如右上角所示，每个椭圆形区域表示拟合出的一个高斯分布，可以看出该结果与实际数据偏差仍然较大。但是，当我们继续增加GMM中高斯分布的数目时，拟合结果也会越来越精确。下图的左下方和右下方分别展示了用5个和10个高斯分布拟合出的结果，已经可以基本覆盖所有的数据点，且几乎没有多出的区域。事实上我们可以证明，任何数据分布都可以用GMM来无限逼近。而在现实场景中，由于高斯分布简单易算、性质良好，GMM也就成为了最常用的模型之一。

图2 用高斯混合分布拟合正弦函数

  从上图的拟合结果中我们还可以发现，如果我们有一组分布未知的数据，并且还希望能按照同样的分布生成一些新的数据，那么也可以先用GMM对数据进行拟合，再由它来继续做数据生成。因此，从数据中拟合GMM就变得十分重要。下面，我们就来用上文介绍的EM算法来求解GMM模型。简单起见，我们就采用与上图左上角中相同的正弦数据集。

import numpy as np
import matplotlib.pyplot as plt

X = np.loadtxt('sindata.csv', delimiter=',')

  在实现EM算法之前，我们先来定义计算高斯分布概率密度$\mathcal{N}(\mathbf{x}|\mathbf{\mu },\mathbf{\Sigma })$的函数。设 $\mathbf{x}\in {\mathbb{R}}^d$，那么概率密度的表达式为$$\mathcal{N}(\mathbf{x}|\mathbf{\mu },\mathbf{\Sigma })=\frac{1}{\sqrt{(2\pi {)}^d|\mathbf{\Sigma }|}}{\mathrm{e}}^{-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^{\mathrm{T}}{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathbf{\mu })}$$ 式中，$|\mathbf{\Sigma }|$表示矩阵$\mathbf{\Sigma }$的行列式，可以用numpy.linalg中的工具直接计算得到。此外，我们通常会计算概率密度的对数而非概率密度本身，这是因为形如 $\log \sum _x{\mathrm{e}}^{f(x)}$ 的式子在计算时，前面的对数-求和-指数部分可以进行优化，其速度相对更快。如果我们计算的是 $f(\mathbf{x})=\log \mathcal{N}(\mathbf{x}|\mathbf{\mu },\mathbf{\Sigma })$，就可以将对数似然计算中的第二个求和转化为对数-求和-指数的形式。

# 计算多元高斯分布的概率密度的对数
def log_gaussian_prob(x, mu, sigma):
    d = x.shape[-1]
    det = np.linalg.det(sigma)
    diff = x - mu
    # 由于x可能包含多个样本
    # x.T @ inv(sigma) @ x 的第二个乘法只需要保留前后样本一致的部分，
    # 所以第二个乘法用点乘再求和代替
    # 此外，由于数据存储维度的问题，这里的转置和上面公式中的转置相反
    log_prob = -d / 2 * np.log(2 * np.pi) - 0.5 * np.log(det) - 0.5 * np.sum(diff @ np.linalg.inv(sigma) * diff, axis=-1)
    return log_prob

  接下来是EM算法的核心部分。为了使算法的调用和参数设置尽可能灵活，我们把EM算法封装在GMM类中。上文中只介绍了EM算法的迭代过程，但没有提及该如何初始化算法的各个参数。为了使算法从比较良好的初始状态出发，我们可以用k-means算法先对数据做聚类，得到每个点的聚类标签。如果将每个聚类看成一个高斯分布，那么这就相当于计算出了每个样本属于哪个分布，也就是隐变量$z$。最后，再从$z$统计计算出每个分布的均值和协方差即可完成初始化。

  下面的代码中，我们分别实现随机初始化和k-means初始化两种方法。需要注意，协方差矩阵必须是半正定的，因此我们先在$[-1,1]$范围内均匀采样矩阵的每个元素，再加上$d$倍的单位矩阵$d{\mathbf{I}}_d$，从而保证矩阵的半正定性。

from sklearn.cluster import KMeans
from scipy.special import logsumexp

class GMM:

    def __init__(self, n_components=2, eps=1e-4, max_iter=100, init='random'):
        # n_components：GMM中高斯分布的数目
        # eps：迭代精度，当对数似然的变化小于eps时迭代终止
        # max_iter：最大迭代次数
        # init：初始化方法，random或kmeans
        self.k = n_components
        self.eps = eps
        self.max_iter = max_iter
        self.init = init
        self.phi = None # 隐变量的先验分布，即每个高斯分布的占比
        self.means = None # 每个高斯分布的均值
        self.covs = None # 每个高斯分布的协方差

    def EM_fit(self, X):
        # 用EM算法求解GMM的参数
        # 参数初始化
        if self.init == 'random': 
            self._random_init_params(X) 
        elif self.init == 'kmeans':
            self._kmeans_init_params(X)
        else:
            raise NotImplementedError
        ll = self._calc_log_likelihood(X) # 当前的对数似然
        n, d = X.shape
        # 开始迭代
        qz = np.zeros((n, self.k)) # z的后验分布
        for t in range(self.max_iter):
            # E步骤，更新后验分布
            for i in range(self.k):
                # 计算样本属于第i类的概率
                log_prob = log_gaussian_prob(X, self.means[i], self.covs[i])
                qz[:, i] = self.phi[i] * np.exp(log_prob)
            # 归一化
            qz = qz / np.sum(qz, axis=1).reshape(-1, 1)

            # M步骤，统计更新参数，最大化对数似然
            self.phi = np.sum(qz, axis=0) / n # 更新隐变量分布
            for i in range(self.k):
                # 更新均值
                self.means[i] = np.sum(qz[:, i, None] * X, axis=0) / n / self.phi[i]
                # 更新协方差
                diff = X - self.means[i]
                self.covs[i] = (qz[:, i, None] * diff).T @ diff / (n - 1) / self.phi[i]
            
            # 判断对数似然是否收敛
            new_ll = self._calc_log_likelihood(X)
            # assert new_ll >= ll, new_ll
            if new_ll - ll <= self.eps:
                break
            ll = new_ll
            
    def _calc_log_likelihood(self, X):
        # 计算当前的对数似然
        ll = 0
        for i in range(self.k):
            log_prob = log_gaussian_prob(X, self.means[i], self.covs[i])
            # 用logsumexp简化计算
            # 该函数底层对对数-求和-指数形式的运算做了优化
            ll += logsumexp(log_prob + np.log(self.phi[i]))
        return ll

    def _random_init_params(self, X):
        self.phi = np.random.uniform(0, 1, self.k) # 随机采样phi
        self.phi /= np.sum(self.phi)
        self.means = np.random.uniform(np.min(X), np.max(X), (self.k, X.shape[1])) # 随机采样均值
        self.covs = np.random.uniform(-1, 1, (self.k, X.shape[1], X.shape[1])) # 随机采样协方差
        self.covs += np.eye(X.shape[1]) * X.shape[1] # 加上维度倍的单位矩阵

    def _kmeans_init_params(self, X):
        # 用Kmeans算法初始化参数
        # 简单起见，我们直接调用sklearn库中的Kmeans方法
        kmeans = KMeans(n_clusters=self.k, init='random', n_init='auto', random_state=0).fit(X)
        # 计算高斯分布占比
        data_in_cls = np.bincount(kmeans.labels_, minlength=self.k)
        self.phi = data_in_cls / len(X)
        # 计算均值和协方差
        self.means = np.zeros((self.k, X.shape[1]))
        self.covs = np.zeros((self.k, X.shape[1], X.shape[1]))
        for i in range(self.k):
            # 取出属于第i类的样本
            X_i = X[kmeans.labels_ == i]
            self.means[i] = np.mean(X_i, axis=0)
            diff = X_i - self.means[i]
            self.covs[i] = diff.T @ diff / (len(X_i) - 1)

  最后，我们设置好超参数，用GMM模型拟合正弦数据集，并绘图观察拟合效果。由于绘图过程要用到从协方差矩阵计算椭圆参数、用matplotlib中的工具绘制椭圆等方法，较为复杂，我们在这里直接给出绘制椭圆的函数和使用方法，不再具体讲解绘制过程。我们分别用2、3、5、10个高斯分布和两种初始化方法进行拟合，并把结果画在一起进行对比。

import matplotlib as mpl

def plot_elipses(gmm, ax):
    # 绘制椭圆
    # gmm：GMM模型
    # ax：matplotlib的画布
    covs = gmm.covs
    for i in range(len(covs)):
        # 计算椭圆参数
        cov = covs[i][:2, :2]
        v, w = np.linalg.eigh(cov)
        u = w[0] / np.linalg.norm(w[0])
        ang = np.arctan2(u[1], u[0])
        ang = ang * 180 / np.pi
        v = 2.0 * np.sqrt(2.0) * np.sqrt(v)
        # 设置椭圆的绘制参数
        # facecolor和edgecolor分别是填充颜色和边缘颜色
        # 可以自由调整
        elp = mpl.patches.Ellipse(gmm.means[i, :2], v[0], v[1], angle=180 + ang, facecolor='orange', edgecolor='none')
        elp.set_clip_box(ax.bbox)
        # 设置透明度
        elp.set_alpha(0.5)
        ax.add_artist(elp)

  从结果图中可以看出，用k-means初始化参数可以让初始的聚类中心比较接近最优的高斯分布均值，而随机初始化的参数则很容易陷入局部的驻点。虽然EM算法理论上是收敛的，但是它并不保证收敛到最优解。并且，如果我们设置了对数似然的变化精度为算法的终止条件，那么算法同样可能在变化较为缓慢的驻点出停止，得到较差的解。因此，EM算法的参数初始化同样重要。除了k-means之外，我们还可以用随机性更小、效果更好的k-means++算法进行初始化。

小故事：
  GMM通常用来拟合数据分布后，再生成相同分布的数据。像这样从数据中学习出分布、再从分布出发完成后续任务的模型称为生成模型（generative model）。与之相对，直接学习样本特征和样本标签直接关联的模型称为判别模型（discriminative model）。我们之前介绍的线性回归等有监督学习模型都属于判别模型。从概率分布的角度出发，假设训练样本为 $\{({\mathbf{x}}_1,y_1),\cdots ,({\mathbf{x}}_n,y_n)\}$，判别模型会将${\mathbf{x}}_i$与$y_i$关联起来，直接学习由样本$\mathbf{x}$给出标签$y$的条件概率分布$p(y|\mathbf{x})$；而生成模型将所有样本关联起来，学习了联合分布$p(\mathbf{x},y)$，然后由此构建对其中标签的条件概率分布$p(y|\mathbf{x})$。在无监督学习下，生成模型只需要学习$p(\mathbf{x})$即可。
  现代深度学习中，生成模型同样有广泛的应用。2013年，迪德里克·金玛（Diederik Kingma）和 马克斯·韦林（Max Welling）提出了变分自编码器；2014年，伊恩·古德费洛（Ian Goodfellow）提出了划时代的生成对抗网络。这两类模型衍生出了许多算法，让生成模型在深度学习时代保有了自己的一席之地。如今，最新的扩散模型在图像任务上取得了令人震惊的效果，并催生了足以以假乱真的人工智能绘画。用于文本任务的生成式预训练变换器可以生成高质量的自然语言，以它为基础的ChatGPT掀起了新一轮人工智能热潮。

三、一般情况下的EM算法

  对于更一般的、样本分布不服从GMM的情况，我们同样可以模仿上面的步骤使用EM算法推进学习。设样本为${\mathbf{x}}_1,\cdots ,{\mathbf{x}}_N$，隐变量为$z$，样本服从参数为$\theta$的某个分布$p(\mathbf{X}|\theta )$，那么参数的对数似然为
$$\begin{array}{rl}l(\theta )& =\log P(\mathbf{X}|\theta )=\log \prod _{i=1}^{N}P({\mathbf{x}}_i|\theta )\\ & =\sum _{i=1}^N\log \sum _{z_i}P({\mathbf{x}}_i,z_i|\theta )\end{array}$$ 上式的第二个求和是对隐变量$z_i$的所有可能取值求和。如果隐变量连续，这里的求和应当改为积分。当每个样本的隐变量$z$给定时，上式对$z$的所有可能性求和就只剩下一项：$$l(\theta )=\sum _{i=1}^N\log P({\mathbf{x}}_i,z_i|\theta )=\sum _{i=1}^N\log P({\mathbf{x}}_i|z_i,\theta )$$ 其中，$z_i$表示样本${\mathbf{x}}_i$对应的隐变量的值。

  然而到此为止，我们并没有论证为什么EM算法是合理的。为什么可以假设隐变量已知？这样的迭代方式为何能给出收敛的结果？为了探究这一问题，我们来从头考虑MLE的求解。我们不妨先对原本的$l(\theta )$进行放缩，尝试求出其下界。设对每个样本${\mathbf{x}}_i$，隐变量$z$的分布是$q_i(z)$，对任意$z$，满足归一性 $\sum _z{q}_i(z)=1$ 和非负性 $q_i(z)\ge 0$。这样对数似然可以写为
$$\begin{array}{rl}l(\theta )& =\sum _{i=1}^N\log \sum _{z_i}P({\mathbf{x}}_i,z_i|\theta )\\ & =\sum _{i=1}^N\log \sum _{z_i}{q}_i(z_i)\frac{P({\mathbf{x}}_i,z_i|\theta )}{q_i(z_i)}\end{array}$$

  为了进行放缩，我们介绍一个重要的常用不等式：延森不等式（Jensen’s inequality）。设函数$f(x)$是凸函数，那么任取自变量$x_1$、$x_2$和$0<\alpha <1$，都有 $f(\alpha x_1+(1-\alpha )x_2)\le \alpha f(x_1)+(1-\alpha )f(x_2)$。反之，如果函数$f(x)$是凹函数，那么任取自变量$x_1$、$x_2$和$0<\alpha <1$，都有 $f(\alpha x_1+(1-\alpha )x_2)\ge \alpha f(x_1)+(1-\alpha )f(x_2)$。如果$f(x)$是严格凸或凹函数，那么上式的不等号在 $x_1\ne x_2$ 时严格成立，当且仅当 $x_1=x_2$ 时，等号成立。

  图3以$\log x$的图像为例展示了延森不等式的含义。从图中可以看出，如果$f(x)$是凹函数，那么在图像上任意取两点连线，必定在函数曲线的下方。因此，自变量的组合 $x^{\prime }=\alpha x_1+(1-\alpha )x_2$ 所对应的函数值$f(x^{\prime })$要高于函数值的组合 $\alpha f(x_1)+(1-\alpha )f(x_2)$。

图3 延森不等式的几何理解(以对数函数为例)

  延森不等式还可以推广到多个点的情况。设$f(x)$是凹函数，任取$x_1,\cdots ,x_n$及其组合系数 $0<{\alpha }_1,\cdots ,{\alpha }_n<1$，且组合系数满足 $\sum _{i=1}^n{\alpha }_i=1$，那么$$f({\alpha }_1{x}_1+\cdots +{\alpha }_n{x}_n)\ge {\alpha }_{1}f(x_1)+\cdots +{\alpha }_{n}f(x_n)$$

  注意到在上面对数似然的表达式中，$\log x$是凹函数，内部求和的系数$q_i(z_i)$之和为1，于是由延森不等式可得$$l(\theta )=\sum _{i=1}^N\log \sum _{z_i}{q}_i(z_i)\frac{P({\mathbf{x}}_i,z_i|\theta )}{q_i(z_i)}\ge \sum _{i=1}^N\sum _{z_i}{q}_i(z_i)\log \frac{P({\mathbf{x}}_i,z_i|\theta )}{q_i(z_i)}$$

  上式给出了对数似然的下界。虽然最大化$l(\theta )$较为困难，但是如果我们能提升其下界，就能为$l(\theta )$的值提供最低保证。如果下界不断提升，那么最后的$l(\theta )$就很可能接近其真正的最大值。因此，我们需要设法选取合适的$q_i(z)$，使其下界尽可能大。从延森不等式的形式中容易看出，如果 $x_1=x_2=\cdots =x_n$，不等式显然可以取到等号。而在下界的表达式中，$\begin{array}{r}\frac{P({\mathbf{x}}_i,z_i|\theta )}{q_i(z_i)}\end{array}$相当于自变量，为了使等号成立，我们令$$q_i(z_i)=\frac{1}{C}P({\mathbf{x}}_i,z_i|\theta )$$ 其中，$C$是常数。这样，所有的“自变量”都等于$C$，不等式等号成立。而常数$C$可以由归一化条件 $\sum _z{q}_i(z)=1$ 得到。所以，$q_i(z_i)$的表达式为$$q_i(z_i)=\frac{P({\mathbf{x}}_i,z_i|\theta )}{\sum _{z}P({\mathbf{x}}_i,z|\theta )}=\frac{P({\mathbf{x}}_i,z_i|\theta )}{P({\mathbf{x}}_i|\theta )}=\frac{P({\mathbf{x}}_i|\theta )P(z_i|{\mathbf{x}}_i,\theta )}{P({\mathbf{x}}_i|\theta )}=P(z_i|{\mathbf{x}}_i,\theta )$$ 这恰好是隐变量$z$的后验分布。这也就解释了为什么EM算法的E步骤要计算$z$的后验分布，再用后验分布得出的$z$来优化$l(\theta )$。本质上，E步骤通过计算后验得出了$l(\theta )$的一个最好的下界，M步骤通过调整参数$\theta$来优化这一下界。而当M步一旦调整了参数$\theta$，那么当前的下界 $\begin{array}{r}\sum _{i=1}^N\sum _{z_i}{q}_i(z_i)\log \frac{P({\mathbf{x}}_i,z_i|\theta )}{q_i(z_i)}\end{array}$ 就不再贴合而是严格小于$l(\theta )$，则又需要再次做E步骤，计算$z$的后验，如此往复。那么这样的迭代是否保证能收敛呢？答案是肯定的，我们在下一节具体讨论EM算法的收敛性。

四、EM算法的收敛性

  上一节中，我们通过放缩解释了EM算法的含义以及合理性，本节来继续证明算法的收敛性。由于我们最终要求解的是参数$\theta$，记迭代第$t$步的M步骤优化前参数为$\theta (t)$，那么我们需要证明优化后的参数${\theta }^{(t+1)}$会使对数似然函数增大，即 $l({\theta }^{(t)})\le l({\theta }^{(t+1)})$。由于对数似然函数是一些不超过1的概率相乘再取对数，其上界为0，因此证明其单调递增就可以保证其收敛。同时，记第$t$步隐变量$z$的分布为 $q_i^{(t)}(z_i)=P(z_i|{\mathbf{x}}_i,{\theta }^{(t)})$，上面已经证明过，这一选择可以使延森不等式取到等号。

  证明的关键在于，$q_i^{(t)}(z_i)$并不是参数${\theta }^{(t+1)}$对应的最优选择，从而
$$\begin{array}{rl}l\left({\theta }^{(t+1)}\right)& =\sum _{i=1}^N\log \sum _{z_i}{q}_i^{(t)}(z_i)\frac{P({\mathbf{x}}_i,z_i|{\theta }^{(t+1)})}{q_i^{(t)}(z_i)}\\ & \ge \sum _{i=1}^N\sum _{z_i}{q}_i^{(t)}(z_i)\log \frac{P({\mathbf{x}}_i,z_i|{\theta }^{(t+1)})}{q_i^{(t)}(z_i)}\\ & \ge \sum _{i=1}^N\sum _{z_i}{q}_i^{(t)}(z_i)\log \frac{P({\mathbf{x}}_i,z_i|{\theta }^{(t)})}{q_i^{(t)}(z_i)}\\ & =l\left({\theta }^{(t)}\right)\end{array}$$ 上式的第一个不等号是延森不等式，且在参数为${\theta }^{(t+1)}$时$q_i^{(t)}(z_i)$不一定能使等号成立。第二个不等号是由于M步骤对参数进行了优化，得到的${\theta }^{(t+1)}$对应的值一定不小于${\theta }^{(t)}$对应的值。最后一个等号是由于$q_i^{(t)}(z_i)$可以使参数为${\theta }^{(t)}$的延森不等式取到等号。于是，数列 $l({\theta }^{(1)}),l({\theta }^{(2)}),\cdots ,l({\theta }^{(t)}),l({\theta }^{(t+1)}),\cdots$ 是单调递增的，并且是有上界的（元素都为非正数），根据单调收敛原理，该数列一定收敛。这样，我们就证明了EM算法的迭代过程必定是收敛的。

  如果记$$J(\theta ,q)=\sum _{i=1}^N\sum _{z_i}{q}_i(z_i)\log \frac{P({\mathbf{x}}_i,z_i|\theta )}{q_i(z_i)}$$ 那么，EM算法的E步骤就可以看作固定$\theta$、优化$q$，而M步骤可以看作固定$q$、优化$\theta$。这样交替优化的方式又称作坐标上升（coordinate ascent）。在支持向量机一文中，我们求解用到的SMO算法其实就是一种特殊的坐标上升算法。当然，坐标上升法并非对所有二元函数都适用。简单来说，如果目标函数是凹且光滑的，那么坐标上升法就能收敛到全局最大值。对于更复杂的情况和收敛性讨论，感兴趣的话可以自行查阅相关数学资料。

附：以上文中的数据集及相关资源下载地址：
链接：https://pan.quark.cn/s/c8dac204cab4
提取码：sviY