主题模型（3）——PLSA模型及其EM算法求解

之前整理过两篇关于主题模型的博客《文本建模之Unigram Model，PLSA与LDA》和《再看LDA主题模型》，主要是整理了主题模型的由来和推导过程，关于模型参数怎么计算没有过多涉及，因此接下来将分两篇博客，分别整理PLSA模型和EM算法求解，LDA模型和Gibbs Sample求解。

PLSA

首先回顾下PLSA，作为生成模型，其在文本生成过程中，引入主题的概念，即先从$K$个主题中选定一个主题，然后在该主题下生成对应的单词。这里，我们将文档$d_i$下每个主题$z_k$的生成概率表示为$p(z_k|d_i)$，每个主题$z_k$下单词（这里的单词是只词汇表$V$所对应的每个不同单词）$w_j$的生成概率表示为$p(w_j|z_k)$。因此文档$d_i$中单词$w_j$的生成概率可以表示为
$$p(d_i,w_j)=p(d_i)p(w_j|d_i)$$

根据概率论公式：边缘分布等于联合分布的和有：
$$\begin{array}{rl}p(d_i,w_j)& =p(d_i)\sum _{k=1}^{K}p(w_j,z_k|d_i)\\ & =p(d_i)\sum _{k=1}^{K}p(w_j|z_k)p(z_k|d_i)\end{array}$$

单词之间相互独立，因此文档$d_i$中所有单词的生成概率为
$$p(d_i,\vec{w})=\prod _{j=1}^{N}p(d_i,w_j{)}^{n(d_i,w_j)}$$

其中，$\vec{w}$是一个$N$维向量$(n(d_i,w_1),n(d_i,w_2),\cdots ,n(d_i,w_N))$，每个分量是单词$w_j$在文档$d_i$中的出现次数，$N$是词汇表$V$的大小。
文档之间同样相互独立，因此语料的生成概率为
$$\begin{array}{rl}p(D)& =\prod _{i=1}^{M}p(d_i,\vec{w})\\ & =\prod _{i=1}^M\prod _{j=1}^{N}p(d_i,w_j{)}^{n(d_i,w_j)}\end{array}$$

采用最大似然估计，最大化$\log p(D)$，对应的$p(z_k|d_i)$和$p(w_j|z_k)$就是我们要找的最优参数。
$$\begin{array}{rl}\mathcal{L}=\log p(D)& =\log \prod _{i=1}^M\prod _{j=1}^{N}p(d_i,w_j{)}^{n(d_i,w_j)}\\ & =\sum _{i=1}^M\sum _{j=1}^{N}n(d_i,w_j)\log p(d_i,w_j)\\ & =\sum _{i=1}^M\sum _{j=1}^{N}n(d_i,w_j)\log [p(d_i)\sum _{k=1}^{K}p(w_j|z_k)p(z_k|d_i)]\\ & =\sum _{i=1}^M\sum _{j=1}^{N}n(d_i,w_j)[\log p(d_i)+\log \sum _{k=1}^{K}p(w_j|z_k)p(z_k|d_i)]\\ & =\sum _{i=1}^M(\sum _{j=1}^{N}n(d_i,w_j)\log p(d_i)+\sum _{j=1}^{N}n(d_i,w_j)\log \sum _{k=1}^{K}p(w_j|z_k)p(z_k|d_i))\\ & =\sum _{i=1}^M(n(d_i)\log p(d_i)+\sum _{j=1}^{N}n(d_i,w_j)\log \sum _{k=1}^{K}p(w_j|z_k)p(z_k|d_i))\\ & =\sum _{i=1}^{M}n(d_i)[\log p(d_i)+\sum _{j=1}^N\frac{n(d_i,w_j)}{n(d_i)}\log \sum _{k=1}^{K}p(w_j|z_k)p(z_k|d_i)]\end{array}$$

因为文档长度$n(d_i)$和文档概率p(d_i)可以单独计算，将其去掉不会影响似然函数的最优化。所以在后面部分，我们统一将$\mathcal{L}$写成如下形式
$$\mathcal{L}=\sum _{i=1}^M\sum _{j=1}^{N}n(d_i,w_j)\log \sum _{k=1}^{K}p(w_j|z_k)p(z_k|d_i)$$

可以看到上式对数函数中包含参数的和，进行求导时将非常的复杂，即直接找到对数似然函数的最大值将不那么容易。庆幸的是，EM（Expectation-Maximization，简称EM）算法能够帮助我们解决这个问题。

EM算法

EM算法是一种启发式迭代算法，每次迭代分为两步：E步，求期望（expectation）；M步，极大化对数似然（maximization），这也是算法名称的由来。回到我们的问题，我们将对数函数包含参数和的情况进行一般化，首先目标函数可以写成下式：
$$\theta =\arg \underset{\theta }{\max }\sum _{i=1}^m\log p(x^{(i)};\theta )$$

$p(x^{(i)};\theta )$是在给定参数$\theta$下第$i$个样本的产生概率，有些情况下这个概率值受到一些隐含因素的影响，如主题模型中文档生成过程中的主题，将其用$z^{(i)}$表示，此时$p(x^{(i)};\theta )$需要累加所有的$p(x^{(i)},z^{(i)};\theta )$，而优化目标也变成：
$$\theta =\arg \underset{\theta }{\max }\sum _{i=1}^m\log p(x^{(i)};\theta )=\arg \underset{\theta }{\max }\sum _{i=1}^m\log \sum _{z^{(i)}}p(x^{(i)},z^{(i)};\theta )$$

上式的难办之处就在于对数函数中出现了和的形式，那有没有一个办法将求和符号提到$\log$外面吗？答案是有的。

Jensen不等式

设$f$是定义在实数域上的函数，如果对于任意的实数$x$，都有
$$f^{\prime \prime }\ge 0$$

那么$f$是凸函数，此时对于随机变量$x$，有
$$E[f(x)]\ge f(E(x))$$

等号成立的条件是随机变量$x$是常量。Jensen不等式用语言描述或许方便记忆：对于凸函数，函数的期望大于等于期望的函数。 证明如下：

如上图所示，对于$(x,f(x))$，$(y,f(y))$线段内的任意一点，其横纵坐标可以分别表示为$tx+(1-t)y$，$tf(x)+(1-t)f(y)$，而函数$f(x)$在$tx+(1-t)y$对应的纵坐标为$f(tf(x)+(1-t)f(y))$，因为函数$f(x)$是凸函数，有
$$tf(x)+(1-t)f(y)\ge f(tx+(1-t)y)$$

因为$t+(1-t)=1$，所以我们可以将$t$和$1-t$看作随机变量$x,y$对应的概率，即上式表达的是
$$E[f(x)]\ge f(E(x))$$

若函数$f$是凹函数，上述不等式符号相反。

EM算法推导

有了Jensen不等式，我们可以对上面的对数似然函数进行缩放如下：
$$\begin{array}{rl}\sum _{i=1}^m\log \sum _{z^{(i)}}p(x^{(i)},z^{(i)};\theta )& =\sum _{i=1}^m\log \sum _{z^{(i)}}Q(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta )}{Q(z^{(i)})}\end{array}\text{\hspace{1em}(1)}$$$$\ge \sum _{i=1}^m\sum _{z^{(i)}}Q(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};\theta )}{Q(z^{(i)})}\text{\hspace{1em}(2)}$$

其中${\sum }_{z^{(i)}}Q(z^{(i)})=1$，即隐变量$z^{(i)}$的概率分布。在上式中，将$Q(z^{(i)})$看作随机变量$\frac{p(x^{(i)},z^{(i)};\theta )}{Q(z^{(i)})}$的概率，函数$f$对应$\log$函数，根据Jensen不等式就能得到上面不等式。
OK，到这里我们就成功的将求和符号提到了对数函数外面，因而式$(2)$将容易求导。但是式$(2)$和式$(1)$是不等号关系，式$(2)$的最大值不一定是式$(1)$的最大值，那怎么才能得到式$(1)$的最大值呢？我们知道，当给定参数$\theta$时，对数似然函数$\mathcal{L}$的值就确定了，也就是式$(2)$的上界就确定了，同时$p(x^{(i)},z^{(i)})$的值就确定了，这时我们修改$Q(z^{(i)})$的值，使下界等于上界。

如果要满足Jensen不等式的等号，需要：
$$\begin{array}{rl}& \frac{p(x^{(i)},z^{(i)};\theta )}{Q(z^{(i)})}=c,c为常量\\ & \Rightarrow p(x^{(i)},z^{(i)};\theta )=cQ(z^{(i)})\\ & \Rightarrow \sum _{z^{(i)}}p(x^{(i)},z^{(i)};\theta )=c\sum _{z^{(i)}}Q(z^{(i)})\\ & \Rightarrow \sum _{z^{(i)}}p(x^{(i)},z^{(i)};\theta )=c\end{array}$$

所以，我们可以得到，等号成立需要满足的条件是
$$Q(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta )}{{\sum }_{z^{(i)}}p(x^{(i)},z^{(i)};\theta )}=p(z^{(i)}|x^{(i)};\theta )$$

即隐变量的后验概率。知道了$Q(z^{(i)})$，将其带入式$(2)$，然后找出使其最大的$\theta$，这时我们就确定了在该$\theta$值处的一个下界，同样我们调整$Q(z^{(i)})$使这个下界最大化，同时也在尝试提升我们的对数似然函数，即我们需要最大化下式：
$$\begin{array}{rl}& \arg \underset{\theta }{\max }\sum _{i=1}^m\sum _{z^{(i)}}Q(z^{(i)})\log \frac{p(x^{(i)},z^{(i)};\theta )}{Q(z^{(i)})}\\ & \Rightarrow \arg \underset{\theta }{\max }\sum _{i=1}^m\sum _{z^{(i)}}Q(z^{(i)})\log p(x^{(i)},z^{(i)};\theta )\\ & =\arg \underset{\theta }{\max }\sum _{i=1}^m\sum _{z^{(i)}}p(z^{(i)}|x^{(i)};\theta )\log p(x^{(i)},z^{(i)};\theta )\end{array}$$

上式可以理解为对数似然函数$\log p(x^{(i)},z^{(i)};\theta )$基于条件概率$p(z^{(i)}|x^{(i)};\theta )$的期望，正好对应EM算法的E步，而最大化过程，正好对应M步。重复上述过程，直到收敛到对数似然函数$\mathcal{L}$的最大处${\theta }^{\ast }$，至于为什么会收敛？证明可以参考李航《统计学习方法》等。注意，与所有迭代算法一样，EM算法不能保证得到全局最优解。

EM算法的形象过程如下图：

总结EM算法的流程如下：
输入：观测数据$x=(x^{(1)},x^{(2)},\cdots ,x^{(m)})$，联合分布$p(x^{(i)},z^{(i)};\theta )$，条件概率$p(z^{(i)}|x^{(i)};\theta )$，最大迭代次数$J$。

1. 随机初始化模型参数${\theta }^0$;
2. E步：计算联合分布的条件概率期望：
   $$Q_j(z^{(i)})=p(z^{(i)}|x^{(i)};{\theta }^j)$$$$E=\sum _{i=1}^m\sum _{z^{(i)}}{Q}_j(z^{(i)})\log p(x^{(i)},z^{(i)};\theta )$$
3. M步：最大化期望，得到${\theta }^{j+1}$
   $${\theta }^{j+1}=\arg \underset{\theta }{\max }\sum _{i=1}^m\sum _{z^{(i)}}{Q}_j(z^{(i)})\log p(x^{(i)},z^{(i)};\theta )$$
4. 重复第（2）步和第（3）步，直至收敛或达到最大迭代次数$J$。

PLSA模型的EM算法步骤

$$\begin{array}{rl}\mathcal{L}& =\log p(D)=\sum _{i=1}^M\sum _{j=1}^{N}n(d_i,w_j)\log \sum _{k=1}^{K}p(w_j|z_k)p(z_k|d_i)\\ & \ge \sum _{i=1}^M\sum _{j=1}^{N}n(d_i,w_j)\sum _{k=1}^{K}Q(z_k)\log \frac{p(w_j|z_k)p(z_k|d_i)}{Q(z_k)}\end{array}$$

根据上面的推导得
$$Q(z_k)=p(z_k|d_i,w_j)$$

然后开始EM算法迭代过程。
E步：计算联合分布的条件概率期望
（1）计算隐变量的后验概率
$$\begin{array}{rl}p(z_k|d_i,w_j)& =\frac{p(z_k,d_i,w_j)}{{\sum }_{k=1}^{K}p(z_k,d_i,w_j)}\\ & =\frac{p(w_j|z_k)p(z_k|d_i)}{{\sum }_{k=1}^{K}p(w_j|z_k)p(z_k|d_i)}\end{array}$$

（2）将隐变量的后验概率带入，得到期望函数
$$\begin{array}{rl}{\mathcal{L}}^c& =\sum _{i=1}^M\sum _{j=1}^{N}n(d_i,w_j)\sum _{k=1}^{K}p(z_k|d_i,w_j)\log p(z_k,d_i,w_j)\\ & =\sum _{i=1}^M\sum _{j=1}^{N}n(d_i,w_j)\sum _{k=1}^{K}p(z_k|d_i,w_j)\log p(w_j|z_k)p(z_k|d_i)\end{array}$$

M步：最大化期望函数，这是一个带约束的最优化问题：
$$\begin{gathered} \sum _{j=1}^{N}p(w_j|z_k)=1 \\ \sum _{k=1}^{K}p(z_k|d_i)=1 \end{gathered}$$

可以采用拉格朗日乘子法求解，拉格朗日函数如下：
$$\mathcal{H}={\mathcal{L}}^c+\sum _{k=1}^K{\tau }_k(1-\sum _{j=1}^{N}p(w_j|z_k))+\sum _{i=1}^M\rho (1-\sum _{k=1}^{K}p(z_k|d_i))$$

分别求偏导，令其结果等于0，最终可以估计出参数$p(w_j|z_k)$和$p(z_k|d_i)$
$$p(w_j|z_k)=\frac{{\sum }_{i=1}^{M}n(d_i,w_j)p(z_k|d_i,w_j)}{{\sum }_{j=1}^N{\sum }_{i=1}^{M}n(d_i,w_j)p(z_k|d_i,w_j)}$$$$p(z_k|d_i)=\frac{{\sum }_{j=1}^{N}n(d_i,w_j)p(z_k|d_i,w_j)}{n(d_i)}$$

参考文献

EM算法原理总结
从最大似然到EM算法浅解