NLP --- > LDA

共轭先验分布

在贝叶斯概率理论中，如果后验概率P(θ|x)和先验概率p(θ)满足同样的分布律，那么，先验分布和后验分布被叫做共轭分布，同时，先验分布叫做似然函数的共轭先验分布
Beta分布是二项式分布的共轭先验分布，而狄利克雷(Dirichlet)分布是多项式分布的共轭分布。
共轭的意思是，以Beta分布和二项式分布为例，数据符合二项分布的时候，参数的先验分布和后验分布都能保持Beta分布的形式，这种形式不变的好处是，我们能够在先验分布中赋予参数很明确的物理意义，这个物理意义可以延续到后续分布中进行解释，同时从先验变换到后验过程中从数据中补充的知识也容易有物理解释。

似然函数
在数理统计学中，似然函数是一种关于统计模型中的参数的函数，表示模型参数中的似然性。似然函数在统计推断中有重大作用，如在最大似然估计和费雪信息之中的应用等等。“似然性”与“或然性”或“概率”意思相近，都是指某种事件发生的可能性，但是在统计学中，“似然性”和“概率”（或然性）又有明确的区分。概率用于在已知一些参数的情况下，预测接下来的观测所得到的结果，而似然性则是用于在已知某些观测所得到的结果时，对有关事物的性质的参数进行估计。
在这种意义上，似然函数可以理解为条件概率的逆反

符号含义

N：所有词的数量f
V:所有词的种类数量
K:主题数量
D:文档数
$w^j=1and{w}^{-j}=0$

LDA 建模过程

1. LDA首先在主题的词分布的分布中抽取K次，得到K个主题的分布。并以此对所有文档采用该先验分布。
2. LDA中然后从主题分布的分布抽样1次得到一篇文档的主题分布，然后该主题分布不变的情况下每个词抽样N次得到N的词的主题。
3. 对于所有文档中，每个文档的主题分布的系数，依然未知；对于K个主题中，每个主题的词分布的系数，依然未知；
4. 目的就是已知$\alpha ,\beta$的情况下，求3中参数，使得$p(w|\alpha ,\beta )$最大，最大likelihood。
5. 方法有两种，一种是collapsed gibbs sampling，一种是variational inference。

Gibbs Sampling

对LDA 准确推断很难,一般使用近似推断，其中gibbs sampling更加简单易懂；
gibbs sampling是mcmc算法的一个特例，运行方式是每次选取概率向量的一个维度，给定其他维度的变量值sampling出当前维度的值，直到收敛。

采样算法流程总结

1. 选择合适的主题数K， 以及合适的超参数$\vec{\alpha },\vec{\eta }$

2. 初始时根据先验分布随机地给每个单词分配主题编号z，然后更新5个数据，分别为
   - nd: M*K（M表示文档数量，K表示主题数量）每个doc中各个topic的词数量 （规则矩阵）
   - nw: V*K （V表示语料库中所有词的数量，相同词合并为一个）词word在主题topic中的数量（规则矩阵）
   - Z: M*n（n表示对于文档m中所有词的数量，相同词不合并为一个）每个doc中主题列表（非规则list集合）
   - nwsum K1 （每个topic词的总数）
   - ndsum M1 （每个doc中词的总数）

3. 对每个文档的每个词（按文档顺序逐个词，一篇文档中词可能有重复），利用gibbs为每个词采样主题及词:$z_i=k,w_i=v$
   $$p(z_i=k|\vec{w},{\vec{z}}_{\neg i})=\frac{n_{d,\neg i}^k+{\alpha }_k}{\sum _{s=1}^K{n}_{d,\neg i}^s+{\alpha }_s}\frac{n_{k,\neg i}^v+{\eta }_v}{\sum _{f=1}^V{n}_{k,\neg i}^f+{\eta }_f}$$

4. 重复第3步过程，直至Gibbs采样收敛

5. 统计预料库中各个文档各个词的主题，得到文档主题分布${\theta }_d$，统计语料库中各个主题词的分布，得到LDA的主题与词的分布。

注意

• 第三步中的采样公式其实就是$p(topic|doc)\cdot p(word|topic)$,这个概率是一个$doc\to topic\to word$的路径，由于top 有K个，所有Gibbs 采样总就是在2K条路径中采样，doc-topic-word路径概率：
  
  在矩阵上，因为我们已知文档编号和词编号，所以也就是根据两向量的各自的积生产的向量，然后据此采样出一个主题

• 然后使用更新的词的主题重新统计， 重新计算$n_k^m+=1,n_m+=1,n_k^t+=1,n_k+=1$

关键代码部分：

def sampling(self,i,j):
    topic = self.Z[i][j]
    word = self.dpre.docs[i].words[j]
    self.nw[word][topic] -= 1
    self.nd[i][topic] -= 1
    self.nwsum[topic] -= 1
    self.ndsum[i] -= 1
    
    Vbeta = self.dpre.words_count * self.beta
    Kalpha = self.K * self.alpha
    self.p = (self.nw[word] + self.beta)/(self.nwsum + Vbeta) * \
             (self.nd[i] + self.alpha) / (self.ndsum[i] + Kalpha)
    for k in xrange(1,self.K):
        self.p[k] += self.p[k-1] 
    
    u = random.uniform(0,self.p[self.K-1])
    for topic in xrange(self.K):
        if self.p[topic]>u:
            break
    
    self.nw[word][topic] +=1
    self.nwsum[topic] +=1
    self.nd[i][topic] +=1
    self.ndsum[i] +=1

def est(self):
    # Consolelogger.info(u"迭代次数为%s 次" % self.iter_times)
    for x in xrange(self.iter_times):
        for i in xrange(self.dpre.docs_count):
            for j in xrange(self.dpre.docs[i].length):
                topic = self.sampling(i,j)
                self.Z[i][j] = topic

---

变分推断

再看一下LDA的整体过程

1. 对于任一篇文档d， 其主题分布${\theta }_d$由Dirichlet分布得到，其中$\alpha$为分布的超参数， 是一个K维向量，我们采样得到的$\theta$也是一个K维的向量：
   $${\theta }_d=Dirichlet(\vec{\alpha })$$
2. 然后对于文档d中的每一词，我们从上述多项式分布$\theta$中再采样出一个主题$z_{dn}$，其为一个值:
   $$z_{dn}=multi({\theta }_d)$$
3. 对于任一主题，其对应的词分布的先验分布也是Dirichlet分布，我们从该先验分布中为一主题采样词分布，其中$\vec{\eta }$和${\beta }_k$的维度均为$1\ast K$：
   $${\beta }_k=Dirichlet(\vec{\eta })$$
4. 现在我们已经为文档d中第n个词采样一个主题$z_{dn}$，并且该主题的词分布为${\beta }_{z_{dn}}$，然后我们为该位置采样出一个词：
   $$w_{dn}=multi({\beta }_{z_{dn}})$$

变化推断过程

上述过程是模拟LDA生成文档的过程，而我们在求解的过程中，文档已经生成，我们要求的是上述过程中涉及到的参数，显然，此时利用到极大似然估计的思想。
我们要最大化的直观上就是$${\int }_d{\int }_{n}p(w_{dn}|\alpha ,\eta )$$
对于其中一个词，可知：
$$p(w|\alpha ,\eta )==\frac{p(\theta ,\beta ,z,w|\alpha ,\eta )}{p(\theta ,\beta ,z|w,\alpha ,\eta )}$$
由于$\theta ,\beta ,z$之间的耦合，这些条件概率没法之间求；为此，我们引入变分推断，具体是引入 mean field assumption的变分推断，这个推断假设所有的隐藏变量都是通过各自的独立分布形成的，如下所示：

这里假设隐藏变量$\theta$是由独立Dirichlet分布$\gamma$形；隐藏变量$z$是由独立多项式分布$\varphi$形成的，隐藏变量$\beta$是由独立Dirichlet分布$\lambda$形成的。这样我们得到三个隐藏变量联合变分分布q为：

$$\begin{array}{rl}q(\beta ,z,\theta |\lambda ,\varphi ,\gamma )& =\prod _{k=1}^{K}q({\beta }_k|{\lambda }_k)\prod _{d=1}^{M}q({\theta }_d,z_d|{\gamma }_d,{\varphi }_d)\\ & =\prod _{k=1}^{K}q({\beta }_k|{\lambda }_k)\prod _{d=1}^M(q({\theta }_d|{\gamma }_d)\prod _{n=1}^{N_d}q(z_{dn}|{\varphi }_{dn}))\end{array}$$

我们的目标就是用$q(\beta ,z,\theta |\lambda ,\varphi ,\gamma )$来近似估计$p(\theta ,\beta ,z|w,\alpha ,\eta )$，也就是让两个分布尽可能相似，这里用到KL散度：
$$({\lambda }^{\ast },{\varphi }^{\ast },{\gamma }^{\ast })=\underset{\lambda ,\varphi ,\gamma }{\underbrace{argmin}}D(q(\beta ,z,\theta |\lambda ,\varphi ,\gamma )||p(\theta ,\beta ,z|w,\alpha ,\eta ))$$ —
我们的目的找到合适的${\lambda }^{\ast },{\varphi }^{\ast },{\gamma }^{\ast }$，然后用$q(\beta ,z,\theta |{\lambda }^{\ast },{\varphi }^{\ast },{\gamma }^{\ast })$来近似隐藏变量的条件分布$p(\theta ,\beta ,z|w,\alpha ,\eta )$，进而使用EM算法求解几个近似分布的参数。

怎么求合适的${\lambda }^{\ast },{\varphi }^{\ast },{\gamma }^{\ast }$参数呢，有
$$\begin{array}{rl}D(q(\beta ,z,\theta |\lambda ,\varphi ,\gamma )||p(\theta ,\beta ,z|w,\alpha ,\eta ))& =E_{q}logq(\beta ,z,\theta |\lambda ,\varphi ,\gamma )-E_{q}logp(\theta ,\beta ,z|w,\alpha ,\eta )\\ & =E_{q}logq(\beta ,z,\theta |\lambda ,\varphi ,\gamma )-E_{q}log\frac{p(\theta ,\beta ,z,w|\alpha ,\eta )}{p(w|\alpha ,\eta )}\\ & =-L(\lambda ,\varphi ,\gamma ;\alpha ,\eta )+log(w|\alpha ,\eta )\end{array}$$
由于对数似然部分$logp(w|\alpha ,\eta )$和我们的KL散度无关，可以看做常量，因此我们希望最小化KL散度等价于最大化ELBO，即转化为求ELBO的最大值。

$$\begin{array}{rl}logp(w|\alpha ,\eta )& =log\int \int \sum _{z}p(\theta ,\beta ,z,w|\alpha ,\eta )d\theta d\beta \\ & =log\int \int \sum _z\frac{p(\theta ,\beta ,z,w|\alpha ,\eta )q(\beta ,z,\theta |\lambda ,\varphi ,\gamma )}{q(\beta ,z,\theta |\lambda ,\varphi ,\gamma )}d\theta d\beta \\ & =log{E}_q\frac{p(\theta ,\beta ,z,w|\alpha ,\eta )}{q(\beta ,z,\theta |\lambda ,\varphi ,\gamma )}\\ & \ge E_{q}log\frac{p(\theta ,\beta ,z,w|\alpha ,\eta )}{q(\beta ,z,\theta |\lambda ,\varphi ,\gamma )}\\ & =E_{q}logp(\theta ,\beta ,z,w|\alpha ,\eta )-E_{q}logq(\beta ,z,\theta |\lambda ,\varphi ,\gamma )\end{array}$$
其中不等号部分用到Jensen不等式：
$$f(E(x))\ge E(f(x))f(x)为凹函数$$
我们把上述推导最后的式子记做：
$$L(\lambda ,\varphi ,\gamma ;\alpha ,\eta )=E_{q}logp(\theta ,\beta ,z,w|\alpha ,\eta )-E_{q}logq(\beta ,z,\theta |\lambda ,\varphi ,\gamma )$$
可以看出$$L(\lambda ,\varphi ,\gamma ;\alpha ,\eta )$$是我们对数似然log $p(w|\alpha ,\eta )$的一个下界，所以这个L一般称为ELBO（(Evidence Lower BOund)。

极大化ELBO求解变分参数

对其进行整理
$$\begin{array}{rl}L(\lambda ,\varphi ,\gamma ;\alpha ,\eta )& =E_q[logp(\beta |\eta )]+E_q[logp(z|\theta )]+E_q[logp(\theta |\alpha )]\\ & +E_q[logp(w|z,\beta )]-E_q[logq(\beta |\lambda )]\\ & -E_q[logq(z|\varphi )]-E_q[logq(\theta |\gamma )]\end{array}$$
可见展开后有7项，我们只对第一项展开：
$$\begin{array}{rl}{E}_q[logp(\beta |\eta )]& =E_q[log\prod _{k=1}^K(\frac{\Gamma (\sum _{i=1}^V{\eta }_i)}{{\prod }_{i=1}^V\Gamma ({\eta }_i)}\prod _{i=1}^V{\beta }_{ki}^{{\eta }_i-1})]\\ & =Klog\Gamma (\sum _{i=1}^V{\eta }_i)-K\sum _{i=1}^{V}log\Gamma ({\eta }_i)+\sum _{k=1}^K{E}_q[\sum _{i=1}^V({\eta }_i-1)log{\beta }_{ki}]\end{array}$$

利用指数分布族的性质
$$E_q[\sum _{i=1}^{V}log{\beta }_{ki}]=(log\Gamma ({\lambda }_{ki})-log\Gamma (\sum _{i^{{}^{\prime }}=1}^V{\lambda }_{k{i}^{{}^{\prime }}}){)}^{{}^{\prime }}=\Psi ({\lambda }_{ki})-\Psi (\sum _{i^{{}^{\prime }}=1}^V{\lambda }_{k{i}^{{}^{\prime }}})$$
$$其中，\Psi (x)=\frac{d}{dx}log\Gamma (x)=\frac{{\Gamma }^{{}^{\prime }}(x)}{\Gamma (x)}$$

上式继续化简
$$\begin{array}{rl}原式& =\\ & Klog\Gamma (\sum _{i=1}^V{\eta }_i)-K\sum _{i=1}^{V}log\Gamma ({\eta }_i)+\sum _{k=1}^K\sum _{i=1}^V({\eta }_i-1)(\Psi ({\lambda }_{ki})-\Psi (\sum _{i^{{}^{\prime }}=1}^V{\lambda }_{k{i}^{{}^{\prime }}}))\end{array}$$
同样，我们可以求出其他几项的表达，得到关于变分参数$\lambda ,\varphi ,\gamma$的表达式。

EM算法E步

再回顾一下我们求解的思路，对于这种存在隐变量的场景，我们需要能够写出含隐变量的分布，但由于原始隐变量之间存在关系，没法直接写成相应的形式，因为我们使用若干分布来近似，并且我们需要保证近似分布与真实分布尽可能接近，即是上述提到的KL散度，让其尽可能变小；然后利用得到的近似分布，执行EM算法。

与普通EM算法不同的是，这里的E步需要在包含期望的EBLO计算最佳变分参数。即通过EBLO对各个变分参数$\lambda ,\varphi ,\gamma$分布求导并令偏导数为0，可以得到迭代表达式，多次收敛后即为最佳变分参数。

直接给出变分参数的表达式
$$\begin{array}{rl}{\varphi }_{nk}& \propto exp(\sum _{i=1}^V{w}_n^i(\Psi ({\lambda }_{ki})-\Psi (\sum _{i^{{}^{\prime }}=1}^V{\lambda }_{k{i}^{{}^{\prime }}}))+\Psi ({\gamma }_k)-\Psi (\sum _{k^{{}^{\prime }}=1}^K{\gamma }_{k^{{}^{\prime }}}))\\ {\gamma }_k& ={\alpha }_k+\sum _{n=1}^N{\varphi }_{nk}\\ {\lambda }_{ki}& ={\eta }_i+\sum _{d=1}^M\sum _{n=1}^{N_d}{\varphi }_{dnk}{w}_{dn}^i\end{array}$$

其中，对于第一个式子
$$\begin{array}{rl}{\varphi }_{nk}& \propto exp(\sum _{i=1}^V{w}_n^i(\Psi ({\lambda }_{ki})-\Psi (\sum _{i^{{}^{\prime }}=1}^V{\lambda }_{k{i}^{{}^{\prime }}}))+\Psi ({\gamma }_k)-\Psi (\sum _{k^{{}^{\prime }}=1}^K{\gamma }_{k^{{}^{\prime }}}))\end{array}$$, $w_n^i=1$当且仅当文档中第n个词为词汇表中第i个词。

最终的E步就是用上述三个式子来更新参数，并且不断循环迭代更新，直到三个变分参数收敛。当三个都收敛后，下一步就是M步，固定变分参数，更新模型参数$\alpha ,\eta$

EM算法M步

我们在该步固定变分参数$\lambda ,\varphi ,\gamma$，极大化ELBO得到最优的模型参数$\alpha ,\eta$，求解方法有很多，LDA这里一般使用牛顿法，即通过求出ELBO对于$\alpha ,\eta$的一阶和二阶导数的表达式，然后迭代求解$\alpha ,\eta$在M步的最优解。

对于$\alpha ,它的一阶和二阶导数表达式为$：

$$\begin{gathered} {\nabla }_{{\alpha }_k}L=M(\Psi (\sum _{k^{{}^{\prime }}=1}^K{\alpha }_{k^{{}^{\prime }}})-\Psi ({\alpha }_k))+\sum _{d=1}^M(\Psi ({\gamma }_{dk})-\Psi (\sum _{k^{{}^{\prime }}=1}^K{\gamma }_{d{k}^{{}^{\prime }}})) \\ {\nabla }_{{\alpha }_k{\alpha }_j}L=M({\Psi }^{{}^{\prime }}(\sum _{k^{{}^{\prime }}=1}^K{\alpha }_{k^{{}^{\prime }}})-\delta (k,j){\Psi }^{{}^{\prime }}({\alpha }_k)) \end{gathered}$$
其中看，当且仅当k=j时，$\delta (k,j)=1$，否则$\delta (k,j)=0$

对于$\eta$，它的一阶和二阶导数表达式为：
$$\begin{gathered} {\nabla }_{{\eta }_i}L=K(\Psi (\sum _{i^{{}^{\prime }}=1}^V{\eta }_{i^{{}^{\prime }}})-\Psi ({\eta }_i))+\sum _{k=1}^K(\Psi ({\lambda }_{ki})-\Psi (\sum _{i^{{}^{\prime }}=1}^V{\lambda }_{k{i}^{{}^{\prime }}})) \\ {\nabla }_{{\eta }_i{\eta }_j}L=K({\Psi }^{{}^{\prime }}(\sum _{i^{{}^{\prime }}=1}^V{\eta }_{i^{{}^{\prime }}})-\delta (i,j){\Psi }^{{}^{\prime }}({\eta }_i)) \end{gathered}$$

最终牛顿法迭代公式为：
$$\begin{array}{r}{\alpha }_k={\alpha }_k+\frac{{\nabla }_{{\alpha }_k}L}{{\nabla }_{{\alpha }_k{\alpha }_j}L}\\ {\eta }_i={\eta }_i+\frac{{\nabla }_{{\eta }_i}L}{{\nabla }_{{\eta }_i{\eta }_j}L}\end{array}$$

summary

变分推断

1. 求解过程实际上求解$logp(w|\alpha ,\eta )$,得到其下界函数$L(\lambda ,\varphi ,\gamma ;\alpha ,\eta )=E_{q}logp(\theta ,\beta ,z,w|\alpha ,\eta )-E_{q}logq(\beta ,z,\theta |\lambda ,\varphi ,\gamma )$
2. 可以通过求解该下届，得到逼近该数值
3. 通过计算化简得到该L函数关于$\lambda ,\varphi ,\gamma$的表达式。
4. 通过对ELBO函数求导并令其为0，得到$\lambda ,\varphi ,\gamma$的表达式的值，三者之间存在联系，迭代更新直到收敛。得到固定$\alpha ,\eta$情况下的最大值。
5. 由ELBO函数对$\alpha ,\eta$求导，另其为0，得到其对应的表达式。由4步得到的值对两者进行更新
6. 重复上述过程，即求解。

变分推断算法流程

算法

输入：主题数K，M个文档与对应的词
1. 初始化$\alpha ,\eta$向量
2. 开始EM算法迭代直至收敛
a）. 初始化所有的$\varphi ,\gamma ,\lambda$，进行E步迭代循环，直至$\varphi ,\gamma ,\lambda$收敛
    (i) for d from 1 to M:
            for n from 1 to $N_d$:
                for k from 1 to K:
                    按照E步中${\varphi }_{nk}$公式更新${\varphi }_{dnk}$
            标准化${\varphi }_{nk}$使得该向量各项的和为1.
        按照E步中${\gamma }_k$公式更新${\gamma }_{dk}$
    (ii) for k from 1 to K:
            for i from 1 to V:
                按照E步更新${\lambda }_{ki}$
    (iii)如果$\varphi ,\gamma ,\lambda$收敛，则跳出a) 步，否则回到(i)步。
b) 进行LDA的M步迭代循环， 直到$\alpha ,\eta$收敛
    按照M步中更新$\alpha ,\eta$参数直到收敛
c) 如果所有的参数均收敛，则算法结束，否则回到第2)步。

算法结束后，我们可以得到模型的后验参数$\alpha ,\eta$，以及我们需要的近似模型主题词分布$\lambda$，以及近似训练文档主题分布$\gamma$。

summary

Q函数也即$logp(w|\alpha ,\eta )$

1. 原则上求Q函数，然后再极大化Q函数
2. 对于单隐变量条件下Q函数为 $E_z[logp(z,o|\lambda )|o,{\lambda }^i]={\sum }_{z}p(z|o,{\lambda }^i)logp(z,o|\lambda )$
3. LDA中隐变量分布$p(\theta ,\beta ,z|w,\alpha ,\eta )$不可直接得到,由$q(\theta ,\beta ,z|\lambda ,\varphi ,\gamma )近似得到$，得到该分布后，代入Q函数，极大化求取$\alpha ,\eta$
4. 在求Q函数时，由KL散度代表两分布的差异性，此时第一个分布的参数$w,\alpha ,\eta$已知，求KL散度时，可转化为求ELBO（evidence lower bound）的极大值
5. Q步时，需迭代求取参数$\lambda ,\varphi ,\gamma$，直到收敛
6. M步时，同样需迭代求取参数$\alpha ,\gamma$，直到收敛

reference
LDA数学八卦
变分推断求解 http://blog.csdn.net/happyer88/article/details/46723941
https://segmentfault.com/a/1190000012215533
https://www.cnblogs.com/pinard/p/6873703.html