機器學習主題模型之LDA引數求解——Gibbs取樣

阿新 • • 發佈：2019-01-04

LDA引數推導的Gibbs取樣方法基於馬爾科夫鏈蒙特卡洛方法，因此首先學習MCMC方法。

一、馬爾科夫鏈蒙特卡洛方法

MCMC（Markov Chain Monte Carlo）方法是構造適合的馬爾科夫鏈，使其平穩分佈為待估引數的後驗分佈，抽樣並使用蒙特卡洛方法進行積分計算，實現了抽樣分佈隨模擬的進行而改變的動態模擬，彌補了傳統蒙特卡洛積分只能靜態模擬的缺陷。

1、蒙特卡洛方法

蒙特卡洛方法是以概率統計的理論方法為基礎的一種數值計算方法，將所要求解的問題同一定的概率模型相聯絡，用計算機實現統計模擬或抽樣，以獲得問題的近似解，故又稱隨機抽樣法或統計試驗法。這種方法求得的是近似解，對於一些複雜問題往往是唯一可行的方法，模擬的樣本數越大越接近真實值，但同時計算量也大幅上升。

蒙特卡洛方法的三個步驟：

（1）構造/描述概率過程

對於本身就具有隨機性質的問題，如粒子運輸問題，只要正確描述和模擬這個概率過程即可；對於本來不是隨機性質的確定性問題，如定積分計算，就需要構造一個人為的概率過程，將不具有隨機性質的問題轉化為隨機性質的問題。

例如要求f(x)在[a,b]上的積分，但原函式F(x)難以求解；如果可以得到x在[a,b]上的概率分佈函式p(x)，則定積分可以近似求解為：

$\int_{a}^{b}f(x)dx=\int_{a}^{b}\frac{f(x)}{p(x)}p(x)dx\approx \frac{1}{n}\sum_{i=0}^{n-1}\frac{f(x_{i})}{p(x_{i})}$

（2）實現從已知概率分佈中抽樣

構造概率模型後，就可以基於概率分佈去取樣所需的n個x的樣本集，帶入求和式子即可求解。最基本的一個概率分佈是(0,1)上的均勻分佈，可以使用線性同餘發生器LCG生成偽隨機數樣本：

x:=(a*x+c)%m, 0<a<m, 0≤c<m, m>0

對於一些常見的連續分佈，如正態分佈、t分佈、F分佈、Beta分佈、 $\Gamma$ 分佈等都可以通過uniform(0,1)的樣本轉換而得。

對於不是常見的分佈，可以使用接受-拒絕取樣來得到該分佈的樣本，即使用一個可以取樣的分佈q(x)與一個常數c相乘後，使得真實分佈總在cq(x)曲線的下方。

重複以下步驟直到獲得n個被接受的取樣點：

①從cq(x)中獲得一個隨機取樣點 $x_{0}$ ，值為 $cq(x_{0})$ ；

②計算接受概率acceptance probability： $\alpha =\frac{p(x_{0})}{cq(x_{0})}$ ；

③從uniform(0,1)中隨機生成一個值u，如果ɑ≥u，則接受取樣 $x_{0}$

，否則拒絕 $x_{0}$ 並回到第一步。

從圖上可以看出，cq(x)與p(x)曲線接近的地方接受概率較高，所以這樣的地方採到的點就會比較多；而在接受概率較低的、差距較大的地方採到的點就會比較少，這也就保證了這個方法的有效性。

（3）建立各種估計量。

在構造了概率模型並從中抽樣後，就可以確定一個隨機變數，作為所要求問題的解——無偏估計。建立各種估計量，相當於對模擬實驗的結果進行考察和登記，從中得到問題的解。

2、馬爾科夫鏈

當後驗分佈的樣本相互獨立時，可以使用大數定律求解；當樣本不獨立時，就要使用馬爾科夫鏈的平穩分佈進行抽樣。

對於一個非週期的馬爾科夫鏈，任何兩個狀態i,j∈E是相通的（不可約），那麼它的狀態轉移概率的極限存在且與i無關，則稱此馬爾科夫鏈具有遍歷性： $\lim_{n\rightarrow \infty }p^{(n)}(i,j)=p_{j}$ ，此時若滿足 $p_{j}\geqslant 0,\; \sum p_{j}=1$ ，則稱 $p_{j}$ 為轉移概率的極限分佈。

平穩分佈：若有限或無限數列 $\left \{ q_{j},j=1,2,... \right \}$ 滿足 $q_{j}\geq 0,\; \sum q_{j}=1$ ，則稱它是概率分佈；如果此概率分佈滿足 $q_{j}=\sum_{i}^{ }q_{i}p(i,j)$ ，則稱它是平穩分佈。

有限馬爾科夫鏈轉移概率的極限分佈一定是平穩分佈；無限馬爾科夫鏈轉移概率的極限分佈不一定是平穩分佈。

在得到平穩分佈的馬爾科夫鏈狀態轉移矩陣後，就可以從這個平穩分佈中取樣蒙特卡洛方法所需的樣本集了。

3、Metropolis-Hastings演算法

對不滿足細緻平穩條件的馬爾科夫鏈狀態轉移矩陣p，引入ɑ(i,j)使其可以取等號，即得到滿足條件的目標矩陣P：

$q_{i}p(i,j)\neq q_{j}p(j,i)\\\\q_{i}p(i,j)\alpha (i,j)= q_{j}p(j,i)\alpha (j,i)\\\\\Rightarrow \left\{\begin{matrix} \alpha (i,j)=q_{j}p(j,i)\\ \alpha (j,i)=q_{i}p(i,j) \end{matrix}\right.\\\\\therefore P(i,j)=\alpha (i,j)\cdot p(i,j)$

ɑ(i,j)∈[0,1]稱為接受率，P可以通過任意一個馬爾科夫鏈狀態轉移矩陣p乘以ɑ(i,j)得到，注意每次轉移後需要歸一化矩陣使得每一行之和為1。但由於ɑ(i,j)可能非常小，導致大部分取樣都被拒絕轉移，馬爾科夫鏈需要取樣很多次才能收斂，取樣效率很低。

MH取樣將接受率做了如下改進，解決了接受率過低的問題：

$\alpha (i,j)=\min \left \{ \frac{q_{j}p(j,i)}{q_{i}p(i,j)},1 \right \}$

MH演算法在特徵很多的情況下，接受率的計算量很大；且在特徵維度很高，難以求出特徵的聯合概率分佈，只能得到特徵之間的條件概率分佈時，需要採取Gibbs取樣的方法解決這類問題。

4、Gibbs Sampling

馬爾科夫鏈通過轉移概率矩陣可以收斂到穩定的概率分佈，這意味著MCMC可以藉助馬爾科夫鏈的平穩分佈特性模擬高維概率分佈p(x)；當馬爾科夫鏈經過burn-in階段，消除初始引數的影響、達到平穩狀態後，每一次狀態轉移都可以生成待模擬分佈的一個樣本。

Gibbs取樣是MCMC演算法的一個特例，只對2維以上的情況有效，它交替隨機固定某一維度，通過其他維度的值來抽樣該維度的值。設p(x,y)是一個二維的聯合概率分佈，由下圖可知：

$p(x_{1},y_{1})p(y_{2}|x_{1})=p(x_{1})p(y_{1}|x_{1})p(y_{2}|x_{1})=p(x_{1},y_{2})p(y_{1}|x_{1})$

$p(x_{1},y_{1})p(x_{2}|y_{1})=p(y_{1})p(x_{1}|y_{1})p(x_{2}|y_{1})=p(x_{2},y_{1})p(x_{1}|y_{1})$

$\Rightarrow \left\{\begin{matrix} p(A)p(y_{2}|x_{1})=p(B)p(y_{1}|x_{1})\\ p(A)p(x_{2}|y_{1})=p(C)p(x_{1}|y_{1}) \end{matrix}\right.$

所以在 $x=x_{1}$ 這條直線上，如果用條件概率分佈 $p(y|x_{1})$ 作為馬爾科夫鏈的狀態轉移概率，則任意兩個點之間的轉移滿足細緻平穩條件；在 $y=y_{1}$ 這條直線上也是如此。Gibbs取樣和MH的不同之處在於，MH只是單獨考慮某兩個狀態的取樣（相當於直接從A到D），而Gibbs取樣是沿著不同維度進行輪換取樣的，且接受率ɑ=1 是MH演算法的特例情況。最終的近似平穩分佈可以通過取樣結果 $\left \{(x_{n},y_{n}),(x_{n+1},y_{n+1}),... \right \}$ 近似表達。

二、Gibbs取樣演算法求解LDA思路

在Gibbs取樣演算法求解LDA的方法中，假設α、η是已知的先驗輸入，目的是推斷隱變數的分佈，即估計p(θ,z,β,w|ɑ,η)，可以通過對θ、β積分將它們邊緣化去掉，這種方法稱為collapse Gibbs sampling。

$p(\beta ,z,\theta ,w|\alpha ,\eta )=\prod_{k=1}^{K}p(\beta _{k}|\eta )\prod_{m=1}^{M}p(\theta _{m}|\alpha )\prod_{n=1}^{N}p(z_{mn}|\theta _{m})p(w_{mn}|z_{mn},\beta _{k})$

$p(z,w|\alpha ,\eta )\\\\=\int \int p(\theta ,z,w,\beta |\alpha ,\eta )d\theta d\beta \\\\=\int \prod_{m=1}^{M}p(\theta _{m}|\alpha )\prod_{n=1}^{N}p(z_{mn}|\theta_{m})d\mathbf{\theta}\cdot \int \prod_{k=1}^{K}p(\beta _{k}|\eta)\prod_{m=1}^{M}\prod_{n=1}^{N}p(w_{mn}|z_{mn},\beta_{k})d\beta$

假設每篇文件長度N相同，分別考慮θ、β這兩部分的積分，令 $n_{n}^{k}$ 代表第m篇文件中分配給第k個主題的詞的個數，對應的多項分佈的計數為： $n_{m}=\left ( n_{m}^{1}, n_{m}^{2},..., n_{m}^{K} \right )$

$\sum_{n=1}^{N}p(z_{mn}|\theta _{m})=\prod_{k=1}^{K}\theta _{mk}^{n_{m}^{k}}$

先對θ部分進行計算，根據α→θ→z的Dir-Mult共軛分佈關係，省略M的連乘可得：

$\int p(\theta _{m}|\alpha )\prod_{n=1}^{N}p(z_{mn}|\theta _{m})d\theta _{m}\\\\=\int \frac{\Gamma (\sum_{k=1}^{K}\alpha _{k})}{\prod_{k=1}^{K}\Gamma (\alpha _{k})}\prod_{k=1}^{K}\theta _{mk}^{\alpha _{k}-1}\prod_{k=1}^{K}\theta _{mk}^{n_{m}^{k}}d\theta _{m}\\\\=\int \frac{\Gamma (\sum_{k=1}^{K}\alpha _{k})}{\prod_{k=1}^{K}\Gamma (\alpha _{k})}\prod_{k=1}^{K}\theta _{mk}^{n_{m}^{k}+\alpha _{k}-1}d\theta _{m}$

根據Dirichlet分佈的積分性質可知：

$\int \frac{\Gamma (\sum_{k=1}^{K}n_{m}^{k}+\alpha _{k})}{\prod_{k=1}^{K}\Gamma (n_{m}^{k}+\alpha _{k})}\prod_{k=1}^{K}\theta _{mk}^{n_{m}^{k}+\alpha _{k}-1}d\theta _{m}=1$

因此θ部分積分可化為：

$\int p(\theta _{m}|\alpha )\prod_{n=1}^{N}p(z_{mn}|\theta _{m})d\theta _{m}= \frac{\Gamma (\sum_{k=1}^{K}\alpha _{k})}{\prod_{k=1}^{K}\Gamma (\alpha _{k})}\cdot \frac{\prod_{k=1}^{K}\Gamma (n_{m}^{k}+\alpha _{k})}{\Gamma (\sum_{k=1}^{K}n_{m}^{k}+\alpha _{k})}$

同樣地，令 $n_{k}^{v}$ 代表第k個主題對應詞典中第v個詞的詞個數，對應的多項分佈的計數為： $n_{k}=\left ( n_{k}^{1}, n_{k}^{2},..., n_{k}^{V} \right )$ ，可將β部分的積分寫為：

$\int \prod_{k=1}^{K}p(\beta _{k}|\eta)\prod_{m=1}^{M}\prod_{n=1}^{N}p(w_{mn}|z_{mn},\beta_{k})d\beta \\\\=\prod_{k=1}^{K}\int p(\beta _{k}|\eta)\prod_{m=1}^{M}\prod_{n=1}^{N}p(w_{mn}|z_{mn},\beta_{k})d\beta_{k}\\\\=\prod_{k=1}^{K}\int \frac{\Gamma (\sum_{v=1}^{V}\eta_{v})}{\prod_{v=1}^{V}\Gamma (\eta_{v})}\prod_{v=1}^{V}\beta _{mk}^{n_{k}^{v}+\eta _{v}-1}d\beta _{k}\\\\=\prod_{k=1}^{K}\frac{\Gamma (\sum_{v=1}^{V}\eta_{v})}{\prod_{v=1}^{V}\Gamma (\eta_{v})}\cdot \frac{\prod_{v=1}^{V}\Gamma (n_{k}^{v}+\eta _{v})}{\Gamma (\sum_{v=1}^{V}n_{k}^{v}+\eta _{v})}$

最終得到主題和詞的聯合分佈如下：

$p(z,w|\alpha ,\eta )=\prod_{m=1}^{M}\frac{\Gamma (\sum_{k=1}^{K}\alpha _{k})}{\prod_{k=1}^{K}\Gamma (\alpha _{k})}\cdot \frac{\prod_{k=1}^{K}\Gamma (n_{m}^{k}+\alpha _{k})}{\Gamma (\sum_{k=1}^{K}n_{m}^{k}+\alpha _{k})} \times \prod_{k=1}^{K}\frac{\Gamma (\sum_{v=1}^{V}\eta_{v})}{\prod_{v=1}^{V}\Gamma (\eta_{v})}\cdot \frac{\prod_{v=1}^{V}\Gamma (n_{k}^{v}+\eta _{v})}{\Gamma (\sum_{v=1}^{V}n_{k}^{v}+\eta _{v})}$

Gibbs取樣的目標是接近分佈p(z|w,ɑ,η)，由於w,ɑ,η已知，令 $z_{mn}$ 表示第m篇文件的第n個詞對應的主題， $z_{-mn}$ 代表z去掉 $z_{mn}$ 後的主題分佈。

$p(z_{mn}|z_{-mn},w,\alpha ,\eta )=\frac{p(z_{mn},z_{-mn},w|\alpha ,\eta)}{p(z_{-mn},w|\alpha ,\eta )}$

注意到Gibbs取樣只需要為 $z_{mn}$ 採一個樣，因此不用得到上式準確的值，只需得到 $z_{mn}$ 取值的比值即可，將上式簡化為：

$p(z_{mn}=k|z_{-mn},w,\alpha ,\eta )\\\\\propto p(z_{mn}=k,z_{-mn},w|\alpha ,\eta )\\\\=\left ( \frac{\Gamma (\sum_{i=1}^{K}\alpha _{i})}{\prod_{i=1}^{K}\Gamma (\alpha _{i})} \right )^{M}\prod_{j\neq m}^{ }\frac{\prod_{i=1}^{K}\Gamma (n_{j}^{i}+\alpha _{i})}{\Gamma (\sum_{i=1}^{K}n_{j}^{i}+\alpha _{i})}\left ( \frac{\Gamma (\sum_{r=1}^{V}\eta _{r})}{\prod_{r=1}^{V}\Gamma (\eta _{r})} \right )^{K}\\\\\cdot \prod_{i=1}^{K}\prod_{r\neq v}^{ }\Gamma (n_{i}^{r}+\eta _{r})\frac{\prod_{i=1}^{K}\Gamma (n_{m}^{i}+\alpha _{i})}{\Gamma (\sum_{i=1}^{K}n_{m}^{i}+\alpha _{i})}\prod_{i=1}^{K}\frac{\Gamma (n_{i}^{v}+\eta _{v})}{\Gamma (\sum_{r=1}^{V}n_{i}^{r}+\eta _{r})}$

$\propto \prod_{i=1}^{K}\Gamma (n_{m}^{i}+\alpha _{i})\frac{\Gamma (n_{i}^{v}+\eta _{v})}{\Gamma (\sum_{r=1}^{V}n_{i}^{r}+\eta _{r})}$

$\propto \prod_{i\neq k}^{ }\Gamma (n_{m}^{i,-mn}+\alpha _{i})\frac{\Gamma (n_{i,-mn}^{v}+\eta _{v})}{\Gamma (\sum_{r=1}^{V}n_{i,-mn}^{r}+\eta _{r})}\cdot \Gamma (n_{m}^{k,-mn}+\alpha _{k}+1)\frac{\Gamma (n_{k,-mn}^{v}+\eta _{v}+1)}{\Gamma (\sum_{r=1}^{V}n_{k,-mn}^{r}+\eta _{r}+1)}$

$= \prod_{i}^{ }\Gamma (n_{m}^{i,-mn}+\alpha _{i})\frac{\Gamma (n_{i,-mn}^{v}+\eta _{v})}{\Gamma (\sum_{r=1}^{V}n_{i,-mn}^{r}+\eta _{r})}\cdot (n_{m}^{k,-mn}+\alpha _{k})\frac{ (n_{k,-mn}^{v}+\eta _{v})}{(\sum_{r=1}^{V}n_{k,-mn}^{r}+\eta _{r})}$

$\propto (n_{m}^{k,-mn}+\alpha _{k})\frac{ (n_{k,-mn}^{v}+\eta _{v})}{(\sum_{r=1}^{V}n_{k,-mn}^{r}+\eta _{r})}$

Gibbs取樣公式的物理意義就是在K條doc→topic→word的路徑中進行取樣，當取樣收斂後可得所有詞最終的主題z，統計每篇文件中的主題概率p(topic|doc)可得文件主題分佈 $\theta _{m}$ ，統計每個主題中所有的詞概率p(word|topic)可得主題詞分佈 $\beta _{k}$ 。

三、LDA Gibbs取樣演算法流程總結

（1）選擇適當的主題總數K，及LDA模型引數ɑ、η；

（2）對應文集中的每一篇文件的每個詞w，隨機賦予一個主題編號z；

（3）重新掃描文集，對每個詞w使用Gibbs取樣公式更新它的主題編號，並更新文集中的詞；

（4）重複取樣過程直至收斂；

（5）統計文集中topic-word共現頻率矩陣，即可得到LDA模型。

通常訓練LDA模型時，取Gibbs取樣收斂之後的n個迭代的結果進行平均來做引數估計，這樣模型質量更高。

當LDA模型訓練好後，要確定一篇新文件的主題時，由於各個主題的詞分佈 $\beta _{k}$ 已經確定，只需要對 $\theta _{m}$ 進行取樣計算即可。

選擇適合的超引數K非常關鍵，對於簡單的語義區分可選擇較小的K，對於複雜的語義區分則選擇較大的K，且需要足夠多的語料。

Gibbs取樣容易並行化，可以很方便的使用大資料平臺來分散式地訓練海量文件的LDA模型。

參考資料

機器學習主題模型之LDA引數求解——Gibbs取樣

機器學習主題模型之LDA引數求解——Gibbs取樣

文字主題模型之LDA(二) LDA求解之Gibbs取樣演算法

文字主題模型之LDA(一) LDA基礎

機器學習精簡教程之九——用scikit-learn的網格搜尋快速找到最優模型引數

機器學習三要素之資料、模型、演算法

【機器學習】最小二乘法求解線性迴歸引數

文字主題模型之潛在語義分析（LDA:Latent Dirichlet Allocation）

李巨集毅機器學習2016 第十五講無監督學習生成模型之 VAE

Java機器學習庫ML之六關於模型迭代訓練的思考

機器學習面試題之——簡單解釋正則化為什麼能減小模型複雜度

機器學習Python實現之線性模型

機器學習：模型評估和選擇

機器學習：模型性能度量(performance measure)(待補充)

人工智能與機器學習的不同之處

機器學習：模型性能評估與參數調優

機器學習數學基礎之矩陣理論（三）

機器學習：模型泛化（L1、L2 和彈性網絡）

機器學習基礎——模型參數評估與選擇

機器學習--boosting家族之GBDT

python 機器學習中模型評估和調參

機器學習主題模型之LDA引數求解——Gibbs取樣

相關推薦