文字主題模型之LDA(二) LDA求解之Gibbs取樣演算法
文字主題模型之LDA(三) LDA求解之變分推斷EM演算法
本文是LDA主題模型的第三篇，讀這一篇之前建議先讀文字主題模型之LDA(一) LDA基礎，同時由於使用了EM演算法，如果你對EM演算法不熟悉，建議先熟悉EM演算法的主要思想。LDA的變分推斷EM演算法求解，應用於Spark MLlib和Scikit-learn的LDA演算法實現，因此值得好好理解。
1. 變分推斷EM演算法求解LDA的思路
首先，回顧LDA的模型圖如下：
變分推斷EM演算法希望通過“變分推斷(Variational Inference)”和EM演算法來得到LDA模型的文件主題分佈和主題詞分佈。首先來看EM演算法在這裡的使用，我們的模型裡面有隱藏變數$\theta,\beta, z$，模型的引數是$\alpha,\eta$。為了求出模型引數和對應的隱藏變數分佈，EM演算法需要在E步先求出隱藏變數$\theta,\beta, z$的基於條件概率分佈的期望，接著在M步極大化這個期望，得到更新的後驗模型引數$\alpha,\eta$。
問題是在EM演算法的E步，由於$\theta,\beta, z$的耦合，我們難以求出隱藏變數$\theta,\beta, z$的條件概率分佈，也難以求出對應的期望，需要“變分推斷“來幫忙，這裡所謂的變分推斷，也就是在隱藏變數存在耦合的情況下，我們通過變分假設，即假設所有的隱藏變數都是通過各自的獨立分佈形成的，這樣就去掉了隱藏變數之間的耦合關係。我們用各個獨立分佈形成的變分分佈來模擬近似隱藏變數的條件分佈，這樣就可以順利的使用EM演算法了。
當進行若干輪的E步和M步的迭代更新之後，我們可以得到合適的近似隱藏變數分佈$\theta,\beta, z$和模型後驗引數$\alpha,\eta$，進而就得到了我們需要的LDA文件主題分佈和主題詞分佈。
可見要完全理解LDA的變分推斷EM演算法，需要搞清楚它在E步變分推斷的過程和推斷完畢後EM演算法的過程。
2. LDA的變分推斷思路
要使用EM演算法，我們需要求出隱藏變數的條件概率分佈如下：$$p(\theta,\beta, z | w, \alpha, \eta) = \frac{p(\theta,\beta, z, w| \alpha, \eta)}{p(w|\alpha, \eta)}$$
前面講到由於$\theta,\beta, z$之間的耦合，這個條件概率是沒法直接求的，但是如果不求它就不能用EM演算法了。怎麼辦呢，我們引入變分推斷，具體是引入基於mean field assumption的變分推斷，這個推斷假設所有的隱藏變數都是通過各自的獨立分佈形成的，如下圖所示：
我們假設隱藏變數$\theta$是由獨立分佈$\gamma$形成的，隱藏變數$z$是由獨立分佈$\phi$形成的，隱藏變數$\beta$是由獨立分佈$\lambda$形成的。這樣我們得到了三個隱藏變數聯合的變分分佈$q$為：$$\begin{align} q(\beta, z, \theta|\lambda,\phi, \gamma) & = \prod_{k=1}^Kq(\beta_k|\lambda_k)\prod_{d=1}^Mq(\theta_d, z_d|\gamma_d,\phi_d) \\ & = \prod_{k=1}^Kq(\beta_k|\lambda_k)\prod_{d=1}^M(q(\theta_d|\gamma_d)\prod_{n=1}^{N_d}q(z_{dn}| \phi_{dn})) \end{align}$$
我們的目標是用$q(\beta, z, \theta|\lambda,\phi, \gamma) $來近似的估計$p(\theta,\beta, z | w, \alpha, \eta)$，也就是說需要這兩個分佈儘可能的相似，用數學語言來描述就是希望這兩個概率分佈之間有儘可能小的KL距離，即：$$(\lambda^*,\phi^*, \gamma^*) = \underbrace{arg \;min}_{\lambda,\phi, \gamma} D(q(\beta, z, \theta|\lambda,\phi, \gamma) || p(\theta,\beta, z | w, \alpha, \eta))$$
其中$D(q||p)$即為KL散度或KL距離，對應分佈$q$和$p$的交叉熵。即：$$D(q||p) = \sum\limits_{x}q(x)log\frac{q(x)}{p(x)} = E_{q(x)}(log\;q(x) - log\;p(x))$$
我們的目的就是找到合適的$\lambda^*,\phi^*, \gamma^*$,然後用$q(\beta, z, \theta|\lambda^*,\phi^*, \gamma^*)$來近似隱藏變數的條件分佈$p(\theta,\beta, z | w, \alpha, \eta)$，進而使用EM演算法迭代。
這個合適的$\lambda^*,\phi^*, \gamma^*$,也不是那麼好求的，怎麼辦呢？我們先看看我能文件資料的對數似然函式$log(w|\alpha,\eta)$如下,為了簡化表示，我們用$E_q(x)$代替$E_{q(\beta, z, \theta|\lambda,\phi, \gamma) }(x)$，用來表示$x$對於變分分佈$q(\beta, z, \theta|\lambda,\phi, \gamma)$ 的期望。
$$\begin{align} log(w|\alpha,\eta) & = log \int\int \sum\limits_z p(\theta,\beta, z, w| \alpha, \eta) d\theta d\beta \\ & = log \int\int \sum\limits_z \frac{p(\theta,\beta, z, w| \alpha, \eta) q(\beta, z, \theta|\lambda,\phi, \gamma)}{q(\beta, z, \theta|\lambda,\phi, \gamma)}d\theta d\beta \\ & = log\;E_q \frac{p(\theta,\beta, z, w| \alpha, \eta) }{q(\beta, z, \theta|\lambda,\phi, \gamma)} \\ & \geq E_q\; log\frac{p(\theta,\beta, z, w| \alpha, \eta) }{q(\beta, z, \theta|\lambda,\phi, \gamma)} \\ & = E_q\; log{p(\theta,\beta, z, w| \alpha, \eta) } - E_q\; log{q(\beta, z, \theta|\lambda,\phi, \gamma)} \end{align} $$
其中，從第(5)式到第(6)式用到了Jensen不等式：$$f(E(x)) \geq E(f(x)) \;\; f(x)為凹函式$$
我們一般把第(7)式記為：$$L(\lambda,\phi, \gamma; \alpha, \eta) = E_q\; log{p(\theta,\beta, z, w| \alpha, \eta) } - E_q\; log{q(\beta, z, \theta|\lambda,\phi, \gamma)}$$
由於$ L(\lambda,\phi, \gamma; \alpha, \eta)$是我們的對數似然的一個下界（第6式），所以這個$L$一般稱為ELBO(Evidence Lower BOund)。那麼這個ELBO和我們需要優化的的KL散度有什麼關係呢？注意到：$$\begin{align} D(q(\beta, z, \theta|\lambda,\phi, \gamma) || p(\theta,\beta, z | w, \alpha, \eta)) & = E_q logq(\beta, z, \theta|\lambda,\phi, \gamma) - E_q log p(\theta,\beta, z | w, \alpha, \eta) \\& =E_q logq(\beta, z, \theta|\lambda,\phi, \gamma) - E_q log \frac{p(\theta,\beta, z, w| \alpha, \eta)}{p(w|\alpha, \eta)} \\& = - L(\lambda,\phi, \gamma; \alpha, \eta) + log(w|\alpha,\eta) \end{align} $$
在(10)式中，由於對數似然部分和我們的KL散度無關，可以看做常量，因此我們希望最小化KL散度等價於最大化ELBO。那麼我們的變分推斷最終等價的轉化為要求ELBO的最大值。現在我們開始關注於極大化ELBO並求出極值對應的變分引數$\lambda,\phi, \gamma$。
3. 極大化ELBO求解變分引數
為了極大化ELBO，我們首先對ELBO函式做一個整理如下：$$\begin{align} L(\lambda,\phi, \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|\phi)] - E_q[logq(\theta|\gamma)] \end{align} $$
可見展開後有7項，現在我們需要對這7項分別做一個展開。為了簡化篇幅，這裡只對第一項的展開做詳細介紹。在介紹第一項的展開前，我們需要了解指數分佈族的性質。指數分佈族是指下面這樣的概率分佈：$$p(x|\theta) = h(x) exp(\eta(\theta)*T(x) -A(\theta))$$
其中，$A(x)$為歸一化因子，主要是保證概率分佈累積求和後為1，引入指數分佈族主要是它有下面這樣的性質：$$\frac{d}{d \eta(\theta)} A(\theta) = E_{p(x|\theta)}[T(x)]$$
這個證明並不複雜，這裡不累述。我們的常見分佈比如Gamma分佈，Beta分佈，Dirichlet分佈都是指數分佈族。有了這個性質，意味著我們在ELBO裡面一大推的期望表示式可以轉化為求導來完成，這個技巧大大簡化了計算量。
回到我們ELBO第一項的展開如下：$$\begin{align} E_q[logp(\beta|\eta)] & = E_q[log\prod_{k=1}^K(\frac{\Gamma(\sum\limits_{i=1}^V\eta_i)}{\prod_{i=1}^V\Gamma(\eta_i)}\prod_{i=1}^V\beta_{ki}^{\eta_i-1})] \\ & = Klog\Gamma(\sum\limits_{i=1}^V\eta_i) - K\sum\limits_{i=1}^Vlog\Gamma(\eta_i) + \sum\limits_{k=1}^KE_q[\sum\limits_{i=1}^V(\eta_i-1) log\beta_{ki}] \end{align} $$
第(15)式的第三項的期望部分，可以用上面講到的指數分佈族的性質，轉化為一個求導過程。即：$$E_q[\sum\limits_{i=1}^Vlog\beta_{ki}] = (log\Gamma(\lambda_{ki} ) - log\Gamma(\sum\limits_{i^{'}=1}^V\lambda_{ki^{'}}))^{'} = \Psi(\lambda_{ki}) - \Psi(\sum\limits_{i^{'}=1}^V\lambda_{ki^{'}})$$
其中：$$\Psi(x) = \frac{d}{d x}log\Gamma(x) = \frac{\Gamma^{'}(x)}{\Gamma(x)}$$
最終，我們得到EBLO第一項的展開式為：$$\begin{align} E_q[logp(\beta|\eta)] & = Klog\Gamma(\sum\limits_{i=1}^V\eta_i) - K\sum\limits_{i=1}^Vlog\Gamma(\eta_i) + \sum\limits_{k=1}^K\sum\limits_{i=1}^V(\eta_i-1)(\Psi(\lambda_{ki}) - \Psi(\sum\limits_{i^{'}=1}^V\lambda_{ki^{'}}) ) \end{align} $$
類似的方法求解其他6項，可以得到ELBO的最終關於變分引數$\lambda,\phi, \gamma$的表示式。其他6項的表示式為：$$\begin{align} E_q[logp(z|\theta)] = \sum\limits_{n=1}^N\sum\limits_{k=1}^K\phi_{nk}\Psi(\gamma_{k}) - \Psi(\sum\limits_{k^{'}=1}^K\gamma_{k^{'}}) \end{align} $$
$$\begin{align} E_q[logp(\theta|\alpha)] & = log\Gamma(\sum\limits_{k=1}^K\alpha_k) - \sum\limits_{k=1}^Klog\Gamma(\alpha_k) + \sum\limits_{k=1}^K(\alpha_k-1)(\Psi(\gamma_{k}) - \Psi(\sum\limits_{k^{'}=1}^K\gamma_{k^{'}})) \end{align} $$
$$\begin{align} E_q[logp(w|z, \beta)] & = \sum\limits_{n=1}^N\sum\limits_{k=1}^K\sum\limits_{i=1}^V\phi_{nk}w_n^i(\Psi(\lambda_{ki}) - \Psi(\sum\limits_{i^{'}=1}^V\lambda_{ki^{'}}) ) \end{align} $$
$$\begin{align} E_q[logq(\beta|\lambda)] = \sum\limits_{k=1}^K(log\Gamma(\sum\limits_{i=1}^V\lambda_{ki}) - \sum\limits_{i=1}^Vlog\Gamma(\lambda_{ki})) + \sum\limits_{k=1}^K\sum\limits_{i=1}^V (\lambda_{ki}-1)(\Psi(\lambda_{ki}) - \Psi(\sum\limits_{i^{'}=1}^V\lambda_{ki^{'}}) )\end{align} $$
$$\begin{align} E_q[logq(z|\phi)] & = \sum\limits_{n=1}^N\sum\limits_{k=1}^K\phi_{nk}log\phi_{nk} \end{align} $$
$$\begin{align} E_q[logq(\theta|\gamma)] & = log\Gamma(\sum\limits_{k=1}^K\gamma_k) - \sum\limits_{k=1}^Klog\Gamma(\gamma_k) + \sum\limits_{k=1}^K(\gamma_k-1)(\Psi(\gamma_{k}) - \Psi(\sum\limits_{k^{'}=1}^K\gamma_{k^{'}}))\end{align} $$
有了ELBO的具體的關於變分引數$\lambda,\phi, \gamma$的表示式，我們就可以用EM演算法來迭代更新變分引數和模型引數了。
4. EM演算法之E步：獲取最優變分引數
有了前面變分推斷得到的ELBO函式為基礎，我們就可以進行EM演算法了。但是和EM演算法不同的是這裡的E步需要在包含期望的EBLO計算最佳的變分引數。如何求解最佳的變分引數呢？通過對ELBO函式對各個變分引數$\lambda,\phi, \gamma$分別求導並令偏導數為0，可以得到迭代表示式，多次迭代收斂後即為最佳變分引數。
這裡就不詳細推導了，直接給出各個變分引數的表示式如下：
$$\begin{align} \phi_{nk} & \propto exp(\sum\limits_{i=1}^Vw_n^i(\Psi(\lambda_{ki}) - \Psi(\sum\limits_{i^{'}=1}^V\lambda_{ki^{'}}) ) + \Psi(\gamma_{k}) - \Psi(\sum\limits_{k^{'}=1}^K\gamma_{k^{'}}))\end{align} $$
其中，$w_n^i =1$當且僅當文件中第$n$個詞為詞彙表中第$i$個詞。
$$\begin{align} \gamma_k & = \alpha_k + \sum\limits_{n=1}^N\phi_{nk} \end{align} $$
$$\begin{align} \lambda_{ki} & = \eta_i + \sum\limits_{n=1}^N\phi_{nk}w_n^i \end{align} $$
由於變分引數$\lambda$決定了$\beta$的分佈，對於整個語料是共有的，因此我們有：
$$\begin{align} \lambda_{ki} & = \eta_i + \sum\limits_{d=1}^M\sum\limits_{n=1}^{N_d}\phi_{dnk}w_{dn}^i \end{align} $$
最終我們的E步就是用（23）（24）（26）式來更新三個變分引數。當我們得到三個變分引數後，不斷迴圈迭代更新，直到這三個變分引數收斂。當變分引數收斂後，下一步就是M步，固定變分引數，更新模型引數$\alpha,\eta$了。
5. EM演算法之M步：更新模型引數
由於我們在E步，已經得到了當前最佳變分引數，現在我們在M步就來固定變分引數，極大化ELBO得到最優的模型引數$\alpha,\eta$。求解最優的模型引數$\alpha,\eta$的方法有很多，梯度下降法，牛頓法都可以。LDA這裡一般使用的是牛頓法，即通過求出ELBO對於$\alpha,\eta$的一階導數和二階導數的表示式，然後迭代求解$\alpha,\eta$在M步的最優解。
對於$\alpha$,它的一階導數和二階導數的表示式為：$$\nabla_{\alpha_k}L = M(\Psi(\sum\limits_{k^{'}=1}^K\alpha_{k^{'}}) - \Psi(\alpha_{k}) ) + \sum\limits_{d=1}^M(\Psi(\gamma_{dk}) - \Psi(\sum\limits_{k^{'}=1}^K\gamma_{dk^{'}}))$$
$$\nabla_{\alpha_k\alpha_j}L = M(\Psi^{'}(\sum\limits_{k^{'}=1}^K\alpha_{k^{'}})- \delta(k,j)\Psi^{'}(\alpha_{k}) )$$
其中，當且僅當$k=j$時，$\delta(k,j)=1$,否則$\delta(k,j)=0$。
對於$\eta$,它的一階導數和二階導數的表示式為：$$\nabla_{\eta_i}L = K(\Psi(\sum\limits_{i^{'}=1}^V\eta_{i^{'}}) - \Psi(\eta_{i}) ) + \sum\limits_{k=1}^K(\Psi(\lambda_{ki}) - \Psi(\sum\limits_{i^{'}=1}^V\lambda_{ki^{'}}))$$
$$\nabla_{\eta_i\eta_j}L = K(\Psi^{'}(\sum\limits_{i^{'}=1}^V\eta_{i^{'}}) - \delta(i,j)\Psi^{'}(\eta_{i}) )$$
其中，當且僅當$i=j$時，$\delta(i,j)=1$,否則$\delta(i,j)=0$。
最終牛頓法法迭代公式為：$$\begin{align} \alpha_{k+1} = \alpha_k + \frac{\nabla_{\alpha_k}L}{\nabla_{\alpha_k\alpha_j}L} \end{align}$$
$$\begin{align} \eta_{i+1} = \eta_i+ \frac{\nabla_{\eta_i}L}{\nabla_{\eta_i\eta_j}L} \end{align}$$
6. LDA變分推斷EM演算法流程總結
下面總結下LDA變分推斷EM的演算法的概要流程。
輸入：主題數$K$,M個文件與對應的詞。
1） 初始化$\alpha,\eta$向量。
2）開始EM演算法迭代迴圈直到收斂。
a) 初始化所有的$\phi, \gamma, \lambda$，進行LDA的E步迭代迴圈,直到$\lambda,\phi, \gamma$收斂。
(i) for d from 1 to M:
for n from 1 to $N_d$:
for k from 1 to K:
按照(23)式更新$\phi_{nk}$
標準化$\phi_{nk}$使該向量各項的和為1.
按照(24) 式更新$\gamma_{k}$。
(ii) for k from 1 to K:
for i from 1 to V:
按照(26) 式更新$\lambda{ki}$。
(iii)如果$\phi, \gamma, \lambda$均已收斂，則跳出a)步，否則回到(i)步。
b) 進行LDA的M步迭代迴圈， 直到$\alpha,\eta$收斂
(i) 按照(27)(28)式用牛頓法迭代更新$\alpha,\eta$直到收斂
c) 如果所有的引數均收斂，則演算法結束，否則回到第2)步。
演算法結束後，我們可以得到模型的後驗引數$\alpha,\eta$,以及我們需要的近似模型主題詞分佈$\lambda$,以及近似訓練文件主題分佈$\gamma$。
（歡迎轉載，轉載請註明出處。歡迎溝通交流： liujianping-ok@163.com）