机器学习笔记十一

第十二讲：$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)}:=\arg \underset{j}{\min }{\Vert x^{(i)}-{\mu }_j\Vert }^2$
• 对每一个簇$j$，令 μ j : = ∑ i = 1 m 1 { c ( i ) = j } x ( i ) ∑ i = 1 m 1 { c ( i ) = j } \mu_j:=\frac{\sum_{i=1}^m1\left\{c^{(i)}=j\right\}x^{(i)}}{\sum_{i=1}^m1\left\{c^{(i)}=j\right\}} μj​:=∑i=1m​1{c(i)=j}∑i=1m​1{c(i)=j}x(i)​

}

在这个算法中，参数$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{\Vert x^{(i)}-{\mu }_{c^{(i)}}\Vert }^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}(\varphi )$（即${\varphi }_j>0,{\sum }_{j=1}^k{\varphi }_j=1,{\varphi }_j=p\left(z^{(i)}=j\right)$，当只有两个高斯分布时，$z^{(i)}$只有两种取值的可能，此时它是一个伯努利分布。），并且有$x^{(i)}\mid z^{(i)}=j\sim \mathcal{N}\left({\mu }_j,{\Sigma }_j\right)$。令$k$表示事件$z^{(i)}$可能取的值的总数，因此，我们的模型假设每个事件$x^{(i)}$均由事件“$z^{(i)}$随机的在 { 1 , ⋯   , k } \{1,\cdots,k\} {1,⋯,k}做出选择”生成，而$x^{(i)}$可以表示为$k$个“由$z^{(i)}$确定的高斯分布”中的一个。这就是混合高斯模型（mixture of Gaussians model）。应注意$z^{(i)}$是**潜在（latent）**随机变量，这意味着它们是隐藏的、未知的，也正是这一点加大了我们估值难度。（变量$z^{(i)}$类似于样本$x^{(i)}$的标记，只是我们并不知道这个变量是多少，也就是说我们不知道每个样本属于哪个标记类，在后面的算法中，我们会尝试猜测这个标记。）

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

$$\begin{array}{rl}l\left(\varphi ,\mu ,\Sigma \right)& =\sum _{i=1}^m\log p\left(x^{(i)};\varphi ,\mu ,\Sigma \right)\\ & =\sum _{i=1}^m\log \sum _{z^{(i)}=1}^{k}p\left(x^{(i)}\mid z^{(i)};\mu ,\Sigma \right)p\left(z^{(i)};\varphi \right)\end{array}$$

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

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

$$l\left(\varphi ,\mu ,\Sigma \right)=\sum _{i=1}^m\log p\left(x^{(i)}\mid z^{(i)};\mu ,\Sigma \right)+\log p\left(z^{(i)};\varphi \right)$$

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

$$\begin{array}{rl}{\varphi }_j& =\frac{1}{m}\sum _{i=1}^{m}1\left\{z^{(i)}=j\right\}\\ {\mu }_j& =\frac{{\sum }_{i=1}^{m}1\left\{z^{(i)}=j\right\}{x}^{(i)}}{{\sum }_{i=1}^{m}1\left\{z^{(i)}=j\right\}}\\ {\Sigma }_j& =\frac{{\sum }_{i=1}^{m}1\left\{z^{(i)}=j\right\}\left(x^{(i)}-{\mu }_j\right){\left(x^{(i)}-{\mu }_j\right)}^T}{{\sum }_{i=1}^{m}1\left\{z^{(i)}=j\right\}}\end{array}$$

如果能确定$z^{(i)}$，则这里的参数估计与前面第五讲高斯判别分析模型中的参数估计几乎是一样的，只不过$z^{(i)}$代替了原来的类标记$y^{(i)}$。（式子里还有其它细节不同，在的高斯判别分析中我们一了解到：第一，不同与问题集中 y ( i ) ∈ { − 1 , 1 } y^{(i)}\in\{-1,1\} y(i)∈{−1,1}服从伯努利分布，我们已经将$z^{(i)}$泛化为多项分布；第二，不同于问题集中不同的$y$共用一个协方差矩阵，我们已经给每个高斯分布赋上了不同的协方差矩阵${\Sigma }_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)};\varphi ,\mu ,\Sigma \right)$（样本$x^{(i)}$属于第$j$个高斯分布的概率）
• (M步骤)更新模型中的参数：
  $$\begin{array}{rl}{\varphi }_j& =\frac{1}{m}\sum _{i=1}^m{w}_j^{(i)}\\ {\mu }_j& =\frac{{\sum }_{i=1}^m{w}_j^{(i)}{x}^{(i)}}{{\sum }_{i=1}^m{w}_j^{(i)}}\\ {\Sigma }_j& =\frac{{\sum }_{i=1}^m{w}_j^{(i)}\left(x^{(i)}-{\mu }_j\right){\left(x^{(i)}-{\mu }_j\right)}^T}{{\sum }_{i=1}^m{w}_j^{(i)}}\end{array}$$

}

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

$$p\left(z^{(i)}=j\mid x^{(i)};\varphi ,\mu ,\Sigma \right)=\frac{p\left(x^{(i)}\mid z^{(i)}=j;\mu ,\Sigma \right)p\left(z^{(i)}=j;\varphi \right)}{{\sum }_{l=1}^{k}p\left(x^{(i)}\mid z^{(i)}=l;\mu ,\Sigma \right)p\left(z^{(i)}=l;\varphi \right)}$$

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

E步骤中求出的$w_j^{(i)}$表示我们对$z^{(i)}$的“软”猜测。（“软”意味着一种概率上的猜测，从$[0,1]$区间取值；相对的“硬”猜测意味着最佳的一次性猜测，是一种对于值的猜测，比如从集合 { 0 , 1 } \{0,1\} {0,1}或 { 1 , ⋯   , k } \{1,\cdots,k\} {1,⋯,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^{\prime \prime }(x)\ge 0$则$f$为凸函数（中文，英文）。在这里，一般化的$f$的输入为向量，而函数的海森矩阵$H$应为半正定矩阵（$H\ge 0$）。如果如果$\forall x\in \mathbb{R},f^{\prime \prime }(x)>0$，则称其为**严格（strictly）**凸函数（再输入为向量的情形下，函数对应的$H$应为为正定矩阵，记为$H>0$）。于是，延森不等式为：

**定理：**令$f$为凸函数，$X$为随机变量，则有$\mathrm{E}\left[f(x)\right]\ge 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{E}X)=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)]\ge f(\mathrm{E}X)$。（有些同学可能会把这个不等式的符号方向记反，记住上面的图就可以快速推导出符号的方向。）

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

2. EM算法

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

$$\begin{array}{rl}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{array}$$

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

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

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

$$\begin{array}{rl}\sum _i\log p\left(x^{(i)};\theta \right)& =\sum _i\log \sum _{z^{(i)}}p\left(x^{(i)},z^{(i)};\theta \right)\\ & =\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)}\\ & \ge \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)}\end{array}\text{\hspace{1em}(1)(2)(3)}$$

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

$$f\left(\underset{z^{(i)}\sim Q_i}{E}\left[\frac{p\left(x^{(i)},z^{(i)};\theta \right)}{Q_i\left(z^{(i)}\right)}\right]\right)\ge \underset{z^{(i)}\sim Q_i}{E}\left[f\left(\frac{p\left(x^{(i)},z^{(i)};\theta \right)}{Q_i\left(z^{(i)}\right)}\right)\right]$$

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

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

对于特定的$\theta$猜测，要使下界紧贴着目标函数，可以再观察引入延森不等式的那一步推导，我们的目标是使不等式取到等号。而上面已经介绍了，如果要取等号，只能是对一个“常量随机变量”求期望，也就是需要令$\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$是一个概率分布，有$\sum _z{Q}_i\left(z^{(i)}\right)=1=\sum _z\frac{p\left(x^{(i)},z^{(i)};\theta \right)}{c}=\frac{{\sum }_{z}p\left(x^{(i)},z^{(i)};\theta \right)}{c}⟹c=\sum _{z}p\left(x^{(i)},z^{(i)};\theta \right)$，则有：

$$\begin{array}{rl}{Q}_i\left(z^{(i)}\right)& =\frac{p\left(x^{(i)},z^{(i)};\theta \right)}{{\sum }_{z}p\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{array}$$

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

对于用上面的后验概率确定的$Q_i$，$(3)$式就能给出一个紧贴目标函数$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步骤)更新参数，令$\theta :=\arg \underset{\theta }{\max }\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算法的介绍。