第十二讲：$k$-means算法、高斯混合模型及最大期望算法¶

现在我们开始进入机器学习课程新的一章——无监督学习算法的应用。

$k$-means聚类算法（$k$-means clustering algorithm）¶

在聚类问题（clustering problem）中，对于给定的训练集$\left\{x^{(1)},\cdots,x^{(m)}\right\}$，并且希望对样本分组，并找到能够使样本“紧密聚集”在一起的簇（clusters）。这里的$x^{(i)}\in\mathbb{R}^n$像原来一样，但是没有$y^{(i)}$。因此，这是一个无监督学习问题。

$k$-means聚类算法的过程如下：

1. 随机的初始化簇质心（cluster centroids）$\mu_1,\mu_2,\cdots,\mu_k\in\mathbb{R}^n$；
2. 重复直到收敛

{

• 对每一个样本$i$，令$c^{(i)}:=\mathrm{arg}\displaystyle\operatorname*{min}_j\left\lVert x^{(i)}-\mu_j\right\rVert^2$
• 对每一个簇$j$，令$\mu_j:=\frac{\sum_{i=1}^m1\left\{c^{(i)}=j\right\}x^{(i)}}{\sum_{i=1}^m1\left\{c^{(i)}=j\right\}}$

}

在这个算法中，参数$k$表示我们想要聚类的簇的数量，$\mu_j$表示我们当前对于簇质心位置的猜测。为了初始化簇质心（即算法第一步），我们通常随机选择$k$个训练样本，然后令$k$个质心等于这$k$个训练样本。（方法不唯一，其他初始化方法也是可以的。）

算法的内层循环重复这两步：(i)将每个训练样本$x^{(i)}$分配给距离它们最近的簇质心$\mu_j$，(ii)移动每个簇质心$\mu_j$到其下辖样本集合的中心。下图为$k$-means算法的简单演示：

图中圆点为训练样本，叉号为簇质心，演示的过程为：(a)原数据集；(b)随机初始化的簇质心（在本例中并不是由某两个随机样本作为簇质心）；(c-f)演示了$k$-means算法的两个迭代，在每次迭代中，我们都将训练样本分配给距离它们最近的簇质心（图中将样本与对应的簇质心标记为同一颜色）；然后将每个簇质心移动到其下辖样本的中心。

那么，$k$-means算法一定收敛吗？在特定的场景下，答案是肯定的。我们先定义失真函数（distortion function）：

$$ J(c,\mu)=\sum_{i=1}^m\left\lVert x^{(i)}-\mu_{c^{(i)}}\right\rVert^2 $$

$J$为每个样本$x^{(i)}$到其簇质心$\mu_{c^{(i)}}$的距离的平方之和，可以证明$k$-means正是在$J$上的坐标下降过程（参考第八讲）。尤其是$k$-means算法的内循环实际上是在重复：保持$\mu$不变最小化$J$关于$c$的函数，然后再保持$c$不变最小化$J$关于$\mu$的函数。因此，$J$是单调递减的，它的函数值一定收敛。（通常这意味着$c$和$\mu$也会收敛。但在理论上，$k$-means算法的最终结果可能在少数簇之间振荡，比如少数几组$c$或/且$\mu$具有完全一样的$J$值。不过在实际操作中，这种现象几乎不会发生。）

失真函数$J$是一个非凸函数，所以对于$J$的坐标下降并不能保证其收敛于全局最小值，即$k$-means会被局部最小值影响。尽管如此，$k$-means算法通常都都能很好的完成聚类任务。不过，如果真的怕算法落入“不好的”局部最优中，一个可行的解决方案是为簇质心$\mu_j$赋上不同的初始化值，多运行几次，然后选出失真函数$J(c,\mu)$最小的一组聚类结果即可。

混合高斯模型（Mixtures of Gaussians）与最大期望算法（EM algorithm）¶

在这一节我们将讨论对聚类样本密度估计（density estimation）的最大期望（EM: Expectation-Maximization）。

假设有给定训练集$\left\{x^{(1)},\cdots,x^{(m)}\right\}$，因为现在研究无监督学习算法，所以训练样本没有对应的$y$标记。

我们希望使用一个联合分布$p\left(x^{(i)},z^{(i)}\right)=p\left(x^{(i)}\mid z^{(i)}\right)p\left(z^{(i)}\right)$（链式法则）对数据建模。这里的$z^{(i)}\sim\mathrm{Multinomial}(\phi)$（即$\phi_j\gt0,\ \sum_{j=1}^k\phi_j=1,\ \phi_j=p\left(z^{(i)}=j\right)$，当只有两个高斯分布时，$z^{(i)}$只有两种取值的可能，此时它是一个伯努利分布。），并且有$x^{(i)}\mid z^{(i)}=j\sim\mathcal{N}\left(\mu_j,\varSigma_j\right)$。令$k$表示事件$z^{(i)}$可能取的值的总数，因此，我们的模型假设每个事件$x^{(i)}$均由事件“$z^{(i)}$随机的在$\{1,\cdots,k\}$做出选择”生成，而$x^{(i)}$可以表示为$k$个“由$z^{(i)}$确定的高斯分布”中的一个。这就是混合高斯模型（mixture of Gaussians model）。应注意$z^{(i)}$是潜在（latent）随机变量，这意味着它们是隐藏的、未知的，也正是这一点加大了我们估值难度。（变量$z^{(i)}$类似于样本$x^{(i)}$的标记，只是我们并不知道这个变量是多少，也就是说我们不知道每个样本属于哪个标记类，在后面的算法中，我们会尝试猜测这个标记。）

模型的参数为$\phi,\mu,\varSigma$，为了估计这些参数，可以写出参数在给定数据上的对数似然估计：

$$ \begin{align} \mathscr{l}\left(\phi,\mu,\varSigma\right)&=\sum_{i=1}^m\log p\left(x^{(i)};\phi,\mu,\varSigma\right)\\ &=\sum_{i=1}^m\log\sum_{z^{(i)}=1}^kp\left(x^{(i)}\mid z^{(i)};\mu,\varSigma\right)p\left(z^{(i)};\phi\right) \end{align} $$

然而，如果对式中每个参数求偏导并置为零，并尝试解出极值时会发现，我们无法求出这些参数的最大似然估计的解析解。

随机变量$z^{(i)}$指出了$x^{(i)}$服从$k$个高斯分布中的哪一个。可以看出，如果我们能够确定$z^{(i)}$，那么这个最大似然问题就会好解很多，我们可以将似然函数写成：

$$ \mathscr{l}\left(\phi,\mu,\varSigma\right)=\sum_{i=1}^m\log p\left(x^{(i)}\mid z^{(i)};\mu,\varSigma\right)+\log p\left(z^{(i)};\phi\right) $$

最大化上面这个关于$\phi,\mu,\varSigma$的式子就能求得各参数：

$$ \begin{align} \phi_j&=\frac{1}{m}\sum_{i=1}^m1\left\{z^{(i)}=j\right\}\\ \mu_j&=\frac{\sum_{i=1}^m1\left\{z^{(i)}=j\right\}x^{(i)}}{\sum_{i=1}^m1\left\{z^{(i)}=j\right\}}\\ \varSigma_j&=\frac{\sum_{i=1}^m1\left\{z^{(i)}=j\right\}\left(x^{(i)}-\mu_j\right)\left(x^{(i)}-\mu_j\right)^T}{\sum_{i=1}^m1\left\{z^{(i)}=j\right\}} \end{align} $$

如果能确定$z^{(i)}$，则这里的参数估计与前面（第五讲）高斯判别分析模型中的参数估计几乎是一样的，只不过$z^{(i)}$代替了原来的类标记$y^{(i)}$。（式子里还有其它细节不同，在问题集1的高斯判别分析中我们一了解到：第一，不同与问题集中$y^{(i)}\in\{-1,1\}$服从伯努利分布，我们已经将$z^{(i)}$泛化为多项分布；第二，不同于问题集中不同的$y$共用一个协方差矩阵，我们已经给每个高斯分布赋上了不同的协方差矩阵$\varSigma_j$。）

然而，在我们的密度估计问题中，$z^{(i)}$是未知的，我们该怎么办呢？

EM算法（中文，英文）是一种迭代算法，有两个主要步骤。对于上面的问题，E步骤中，算法尝试猜测每个$z^{(i)}$的值，在M步骤中算法会根据猜测更新模型。因为M步骤会假设E步骤中的猜测是正确的，那么，有了缺失的$z^{(i)}$，最大化就很简单了。算法的流程：

重复直到收敛：{

• (E步骤)对于每个$i,j$，令$w_j^{(i)}:=p\left(z^{(i)}=j\mid x^{(i)};\phi,\mu,\varSigma\right)$（样本$x^{(i)}$属于第$j$个高斯分布的概率）
• (M步骤)更新模型中的参数： $$ \begin{align} \phi_j&=\frac{1}{m}\sum_{i=1}^mw_j^{(i)}\\ \mu_j&=\frac{\sum_{i=1}^mw_j^{(i)}x^{(i)}}{\sum_{i=1}^mw_j^{(i)}}\\ \varSigma_j&=\frac{\sum_{i=1}^mw_j^{(i)}\left(x^{(i)}-\mu_j\right)\left(x^{(i)}-\mu_j\right)^T}{\sum_{i=1}^mw_j^{(i)}} \end{align} $$

}

在E步骤中，我们使用$x^{(i)}$以及当前的参数值计算所有的$z^{(i)}$，即使用贝叶斯定理：

$$ p\left(z^{(i)}=j\mid x^{(i)};\phi,\mu,\varSigma\right)=\frac{p\left(x^{(i)}\mid z^{(i)}=j;\mu,\varSigma\right)p\left(z^{(i)}=j;\phi\right)}{\sum_{l=1}^kp\left(x^{(i)}\mid z^{(i)}=l;\mu,\varSigma\right)p\left(z^{(i)}=l;\phi\right)} $$

• 式中的$p\left(x^{(i)}\mid z^{(i)}=j;\mu,\varSigma\right)$为某个以$\mu_j$为期望、$\varSigma_j$为协方差的高斯分布在$x^{(i)}$处的密度估计，即$\frac{1}{\sqrt{(2\pi)^n\left\lvert\varSigma_j\right\rvert}}\exp\left(-\frac{1}{2}\left(x^{(i)}-\mu_j\right)^T\varSigma_j^{-1}\left(x^{(i)}-\mu_j\right)\right)$；
• $p\left(z^{(i)}=j;\phi\right)$为$\phi_j$。

E步骤中求出的$w_j^{(i)}$表示我们对$z^{(i)}$的“软”猜测。（“软”意味着一种概率上的猜测，从$[0,1]$区间取值；相对的“硬”猜测意味着最佳的一次性猜测，是一种对于值的猜测，比如从集合$\{0,1\}$或$\{1,\cdots,k\}$中取值。）

我们也应该对比M步骤中$z^{(i)}$已知时参数更新的式子，能够发现，只有式子里的由高斯分布中每个样本点求出的指示函数$1\left\{z^{(i)}=j\right\}$变为了$w_j^{(i)}$，其余部分同原来一样。

同样也可以用前面的$k$-means聚类算法做对比，不同之处在于$k$-means算法将样本“硬”分配给不同的簇，EM算法用$w_j^{(i)}$做“软”分配。相似之处在于，EM算法也存在落入局部最优的可能。所以，使用不同的初始参数多运行几次的方法对EM算法同样适用。

显然，可以很自然的解释EM算法对于$z^{(i)}$的迭代猜测，但是这个点子来自哪里？我们能对算法属性做出证明吗，比如证明其收敛性？在下一节，我们将给出一个更加泛化的EM算法的描述，它会使我们轻易的将EM算法应用在不同的具有潜在随机变量的估值问题中，而且也能够证明算法的收敛性。

第九部分：EM算法¶

在前面我们介绍了如何用EM算法应用在高斯混合模型的拟合上，在本节，我们将介绍EM算法更广泛的应用，并展示如何将其应用在含有潜在变量的估值问题族中。我们从延森不等式（Jensen's inequality）这一通途广泛的的结论开始本节的讨论。

1. 延森不等式¶

令函数$f$的值域为实数集，回忆以前的知识，如果$\forall x\in\mathbb{R},f''(x)\geq0$则$f$为凸函数（中文，英文）。在这里，一般化的$f$的输入为向量，而函数的海森矩阵$H$应为半正定矩阵（$H\geq0$）。如果如果$\forall x\in\mathbb{R},f''(x)\gt0$，则称其为严格（strictly）凸函数（再输入为向量的情形下，函数对应的$H$应为为正定矩阵，记为$H\gt0$）。于是，延森不等式为：

定理：令$f$为凸函数，$X$为随机变量，则有$\mathrm{E}\left[f(x)\right]\geq f(\mathrm{E}X)$。

此外，如果$f$严格凸，则当且仅当$X=\mathrm{E}[X]$的概率为$1$时（即$p(X=\mathrm{E}[X])=1$，也就是说$X$是常量），有$\mathrm{E}\left[f(x)\right]=f(\mathrm{E}X)$。

根据惯例，我们可以在写期望的时候省略方括号，$f(\mathrm EX)=f(\mathrm E[X])$。

可以用下图解释这个定理：

图中的实现就是凸函数$f$。X是一个随机变量，有$50%$的概率取到$a$、$50%$的概率取到$b$（沿$x$轴），则$X$的期望应在$a,b$的中点。

再看$y$轴上的$f(a),f(b),f(\mathrm{E}[X])$和$\mathrm E[f(x)]$，其中$\mathrm E[f(x)]$在$f(a),f(b)$的中点，从图中可以看出，因为$f$是凸函数，必然有$\mathrm E[f(x)]\geq f(\mathrm EX)$。（有些同学可能会把这个不等式的符号方向记反，记住上面的图就可以快速推导出符号的方向。）

注意：回忆以前的知识，当且仅当$-f$是[严格]凸函数时，$f$是[严格]凹函数（即$f''(x)\leq0$或$H\leq0$）。延森不等式同样适用于凹函数，只不过所有不等式的不等号的方向相反（如$\mathrm E[f(x)]\leq f(\mathrm EX)$等）。

2. EM算法¶

对于给定的包含$m$个相互独立的样本的训练集$\left\{x^{(1)},\cdots,x^{(m)}\right\}$，假设有估值问题：我们希望使用模型$p(x,z)$拟合样本数据，对数似然函数为：

$$ \begin{align} \mathscr l(\theta)&=\sum_{i=1}^m\log p\left(x^{(i)};\theta\right)\\ &=\sum_{i=1}^m\log\sum_{z^{(i)}}p\left(x^{(i)},z^{(i)};\theta\right) \end{align} $$

但是直接计算出参数$\theta$的最大似然估计可能会很难，因为这里的$z^{(i)}$是潜在随机变量。不过，如果能够确定$z^{(i)}$的值，参数的最大似然估计通常很容易计算。

在这种情形下，EM算法给出了一个能够有效的求出参数最大似然估计的方法。直接最大化$\mathscr l(\theta)$可能很难，EM算法的思路是不停的构造$\mathscr l$函数的下界（E步骤），然后最大化这些下界（M步骤）。（也就是说，算法先初始化一套参数$\theta^{(0)}$，然后构造一个下界$\mathscr l$并求这个$\mathscr l$的最大值，而取最大值时的的这套参数就是$\theta^{(1)}$；算法再从$\theta^{(1)}$开始构造下一个下界并最大化得到$\theta^{(2)}$……直到收敛于一个局部最优解）。下面介绍EM算法的形式描述：

对每一个$i$，令$Q_i$为某个关于$z$的分布（$\displaystyle\sum_{z^{(i)}}Q_i\left(z^{(i)}\right)=1,Q_i\left(z^{(i)}\right)\geq0$），考虑下面的式子（若$z$是连续变量，则$Q_i$为概率密度，此时讨论中的所有对$z$的求和应变为对$z$的积分）：

$$ \begin{align} \sum_i\log p\left(x^{(i)};\theta\right)&=\sum_i\log\sum_{z^{(i)}}p\left(x^{(i)},z^{(i)};\theta\right)\tag{1}\\ &=\sum_i\log\sum_{z^{(i)}}Q_i\left(z^{(i)}\right)\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}\tag{2}\\ &\geq\sum_i\sum_{z^{(i)}}Q_i\left(z^{(i)}\right)\log\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}\tag{3} \end{align} $$

$(2)$式到$(3)$式的推导使用了延森不等式，对于$f(x)=\log(x)$有$f''(x)=-\frac{1}{x^2}\lt0$，所以$f(x)$是一个定义在$x\in\mathbb R^+$上的凹函数，则有$f(\mathrm EX)\geq \mathrm E[f(X)]$。再看$(2)$式中的$\displaystyle\sum_{z^{(i)}}\underbrace{Q_i\left(z^{(i)}\right)}_\textrm{probability}\underbrace{\left[\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}\right]}_\textrm{variable}$项，恰好是$z^{(i)}$服从$Q_i$的分布时，变量$\left[\displaystyle\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}\right]$的期望表达式（期望的定义：$\mathrm E[g(x)]=\sum_xp(x)g(x)$）。那么，根据延森不等式，应有：

$$ f\left(\operatorname*E_{z^{(i)}\sim Q_i}\left[\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}\right]\right)\geq\operatorname*E_{z^{(i)}\sim Q_i}\left[f\left(\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}\right)\right] $$

式中的$\displaystyle\operatorname*E_{z^{(i)}\sim Q_i}$指在随机变量$z^{(i)}$服从$Q_i$分布的情况下求期望。最终我们得到$(3)$式。

对于任意的分布$Q_i$，$(3)$式都能给出一个$\mathscr l(\theta)$的下界。而$Q_i$有很多种选择，我们接下来该怎么做？如果我们对当前的参数$\theta$做一次猜测，我们希望这个猜测能够紧贴着目标函数$\mathscr l(\theta)$，因为只有这样，我们才能确定每次迭代都是在逼近目标函数的局部最优解。也就是说我们希望特定的$\theta$猜测能够使这个“下界不等式”取到等号部分。（后面会解释这个等号如何保证$\mathscr l(\theta)$在每一次迭代中都单调递增。）

对于特定的$\theta$猜测，要使下界紧贴着目标函数，可以再观察引入延森不等式的那一步推导，我们的目标是使不等式取到等号。而上面已经介绍了，如果要取等号，只能是对一个“常量随机变量”求期望，也就是需要令$\displaystyle\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}=c$，$c$是某个与$z^{(i)}$无关的常量（对于任意$z^{(i)}$，都取同一个常量$c$）。于是，令$Q_i\left(z^{(i)}\right)\propto p\left(x^{(i)},z^{(i)};\theta\right)$即可（令分母与分子保持同一个比例）。

实际上，因为我们知道$Q_i$是一个概率分布，有$\displaystyle\sum_zQ_i\left(z^{(i)}\right)=1=\sum_z\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{c}=\frac{\sum_zp\left(x^{(i)},z^{(i)};\theta\right)}{c}\implies c=\sum_zp\left(x^{(i)},z^{(i)};\theta\right)$，则有：

$$ \begin{align} Q_i\left(z^{(i)}\right)&=\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{\sum_zp\left(x^{(i)},z;\theta\right)}\\ &=\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{p\left(x^{(i)};\theta\right)}\\ &=p\left(z^{(i)}\mid x^{(i)};\theta\right) \end{align} $$

因此，我们只需要将$Q_i$设置为“对于给定的$x^{(i)}$下$z^{(i)}$的后验概率”即可，然后再此时的$\theta$就行了。

对于用上面的后验概率确定的$Q_i$，$(3)$式就能给出一个紧贴目标函数$\mathscr l(\theta)$的“下界函数”，这就是E步骤，而后面我们会尝试最大化这个“下界函数”，也就是在M步骤中，我们会尝试最大化$(3)$式得出的式子（用前面E步骤的$\theta$确定的式子），并得到一组新的参数$\theta$。不断重复这两部就是EM算法了：

重复直到收敛：{

• (E步骤)对于每个$i$，令$Q_i\left(z^{(i)}\right):=p\left(z^{(i)}\mid x^{(i)};\theta\right)$
• (M步骤)更新参数，令$\displaystyle\theta:=\arg\max_\theta\sum_i\sum_{z^{(i)}}Q_i\left(z^{(i)}\right)\log\frac{p\left(x^{(i)},z^{(i)};\theta\right)}{Q_i\left(z^{(i)}\right)}$

}

下一讲将继续EM算法的介绍。