现在我们开始进入机器学习课程新的一章——无监督学习算法的应用。
在聚类问题(clustering problem)中,对于给定的训练集$\left\{x^{(1)},\cdots,x^{(m)}\right\}$,并且希望对样本分组,并找到能够使样本“紧密聚集”在一起的簇(clusters)。这里的$x^{(i)}\in\mathbb{R}^n$像原来一样,但是没有$y^{(i)}$。因此,这是一个无监督学习问题。
$k$-means聚类算法的过程如下:
{
}
在这个算法中,参数$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)$最小的一组聚类结果即可。
在这一节我们将讨论对聚类样本密度估计(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步骤中,我们使用$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)} $$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算法更广泛的应用,并展示如何将其应用在含有潜在变量的估值问题族中。我们从延森不等式(Jensen's inequality)这一通途广泛的的结论开始本节的讨论。
令函数$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)$等)。
对于给定的包含$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算法了:
重复直到收敛:{
}
下一讲将继续EM算法的介绍。