轉載請注明出處，該文章的官方來源：

LDA | Teaching ML

前言

LDA是一種概率主題模型：隱式狄利克雷分布（Latent Dirichlet Allocation，簡稱LDA）。LDA是2003年提出的一種主題模型，它可以將文檔集中每篇文檔的主題以概率分布的形式給出。 通過分析一些文檔，我們可以抽取出它們的主題（分布），根據主題（分布）進行主題聚類或文本分類。同時，它是一種典型的詞袋模型，即一篇文檔是由一組詞構成，詞與詞之間沒有先后順序的關系。一篇文檔可以包含多個主題，文檔中每一個詞都由其中的一個主題生成。

舉一個簡單的例子，比如假設事先給定了這幾個主題：Arts、Budgets、Children、Education，然后通過學習的方式，獲取每個主題Topic對應的詞語，如下圖所示：

然后以一定的概率選取上述某個主題，再以一定的概率選取那個主題下的某個單詞，不斷的重復這兩步，最終生成如下圖所示的一篇文章（不同顏色的詞語分別表示不同主題）。

我們看到一篇文章后，往往會推測這篇文章是如何生成的，我們通常認為作者會先確定幾個主題，然后圍繞這幾個主題遣詞造句寫成全文。LDA要干的事情就是根據給定的文檔，判斷它的主題分布。在LDA模型中，生成文檔的過程有如下幾步：

• 從狄利克雷分布 中生成文檔i的主題分布 ；
• 從主題的多項式分布 中取樣生成文檔i第j個詞的主題 ；
• 從狄利克雷分布 中取樣生成主題 對應的詞語分布 ；
• 從詞語的多項式分布 中采樣最終生成詞語 .

LDA的圖模型結構如下圖所示：

LDA會涉及很多數學知識，后面的章節我會首先介紹LDA涉及的數學知識，然后在這些數學知識的基礎上詳細講解LDA的原理。

1 數學預備

1.1 Gamma函數

在高等數學中，有一個長相奇特的Gamma函數

通過分部積分，可以推導gamma函數有如下遞歸性質

通過該遞歸性質，我們可以很容易證明，gamma函數可以被當成階乘在實數集上的延拓，具有如下性質

1.2 Digamma函數

如下函數被稱為Digamma函數，它是Gamma函數對數的一階導數

這是一個很重要的函數，在涉及Dirichlet分布相關的參數的極大似然估計時，往往需要用到這個函數。Digamma函數具有如下一個漂亮的性質

1.3 二項分布（Binomial distribution）

二項分布是由伯努利分布推出的。伯努利分布，又稱兩點分布或0-1分布，是一個離散型的隨機分布，其中的隨機變量只有兩類取值，即0或者1。二項分布是重復n次的伯努利試驗。簡言之，只做一次實驗，是伯努利分布，重復做了n次，是二項分布。二項分布的概率密度函數為：

對于k=1,2，...,n，其中C(n,k)是二項式系數（這就是二項分布的名稱的由來）

1.4 多項分布

多項分布是二項分布擴展到多維的情況。多項分布是指單次試驗中的隨機變量的取值不再是0-1，而是有多種離散值可能（1,2,3...,k）。比如投擲6個面的骰子實驗，N次實驗結果服從K=6的多項分布。其中：

多項分布的概率密度函數為：

1.5 Beta分布

1.5.1 Beta分布

首先看下面的問題1（問題1到問題4都取自于文獻【1】）。

問題1：

為解決這個問題，可以嘗試計算

落在區間[x,x+delta x]的概率。首先，把[0,1]區間分成三段[0,x),[x,x+delta x]，(x+delta x,1]，然后考慮下簡單的情形：即假設n個數中只有1個落在了區間[x,x+delta x]內，由于這個區間內的數X(k)是第k大的，所以[0,x)中應該有k?1個數，(x+delta x,1]這個區間中應該有n?k個數。 如下圖所示：

上述問題可以轉換為下述事件E：

對于上述事件E，有：

其中，o(delta x)表示delta x的高階無窮小。顯然，由于不同的排列組合，即n個數中有一個落在[x,x+delta x]區間的有n種取法，余下n?1個數中有k?1個落在[0,x)的有C(n-1,k-1)種組合。所以和事件E等價的事件一共有nC(n-1,k-1)個。

文獻【1】中證明，只要落在[x,x+delta x]內的數字超過一個，則對應的事件的概率就是o(delta x)。所以

的概率密度函數為：

利用Gamma函數，我們可以將f(x)表示成如下形式：

在上式中，我們用alpha=k，beta=n-k+1替換，可以得到beta分布的概率密度函數

1.5.2 共軛先驗分布

什么是共軛呢？軛的意思是束縛、控制。共軛從字面上理解，則是共同約束，或互相約束。在貝葉斯概率理論中，如果后驗概率P(z|x)和先驗概率p(z)滿足同樣的分布，那么，先驗分布和后驗分布被叫做共軛分布，同時，先驗分布叫做似然函數的共軛先驗分布。

1.5.3 Beta-Binomial 共軛

我們在問題1的基礎上增加一些觀測數據，變成問題2：

第2步的條件可以用另外一句話來表述，即“Yi中有m1個比X(k)小，m2個比X(k)大”，所以X(k)是

中k+m1大的數。

根據1.5.1的介紹，我們知道事件p服從beta分布,它的概率密度函數為：

按照貝葉斯推理的邏輯，把以上過程整理如下：

• 1、p是我們要猜測的參數，我們推導出p的分布為f(p)=Beta(p|k,n-k+1),稱為p的先驗分布
• 2、根據Yi中有m1個比p小，有m2個比p大，Yi相當是做了m次伯努利實驗，所以m1服從二項分布B(m,p)
• 3、在給定了來自數據提供(m1,m2)知識后，p的后驗分布變為f(p|m1,m2)=Beta(p|k+m1,n-k+1+m2)

貝葉斯估計的基本過程是：

先驗分布 + 數據的知識 = 后驗分布

以上貝葉斯分析過程的簡單直觀的表示就是：

Beta(p|k,n-k+1) + BinomCount(m1,m2) = Beta(p|k+m1,n-k+1+m2)

更一般的，對于非負實數alpha和beta，我們有如下關系

Beta(p|alpha,beta) + BinomCount(m1,m2) = Beta(p|alpha+m1,beta+m2)

針對于這種觀測到的數據符合二項分布，參數的先驗分布和后驗分布都是Beta分布的情況，就是Beta-Binomial共軛。換言之，Beta分布是二項式分布的共軛先驗概率分布。二項分布和Beta分布是共軛分布意味著，如果我們為二項分布的參數p選取的先驗分布是Beta分布，那么以p為參數的二項分布用貝葉斯估計得到的后驗分布仍然服從Beta分布。

1.6 Dirichlet 分布

1.6.1 Dirichlet 分布

Dirichlet分布，是beta分布在高維度上的推廣。Dirichlet分布的的密度函數形式跟beta分布的密度函數類似：

其中

至此，我們可以看到二項分布和多項分布很相似，Beta分布和Dirichlet分布很相似。并且Beta分布是二項式分布的共軛先驗概率分布。那么Dirichlet分布呢？Dirichlet分布是多項式分布的共軛先驗概率分布。下文來論證這點。

1.6.2 Dirichlet-Multinomial 共軛

在1.5.3章問題2的基礎上，我們更進一步引入問題3：

類似于問題1的推導，我們可以容易推導聯合分布。為了簡化計算，我們取x3滿足x1+x2+x3=1,x1和x2是變量。如下圖所示。

概率計算如下：

于是我們得到聯合分布為：

觀察上述式子的最終結果，可以看出上面這個分布其實就是3維形式的Dirichlet分布。令alpha1=k1,alpha2=k2,alpha3=n-k1-k2+1，分布密度函數可以寫為：

為了論證Dirichlet分布是多項式分布的共軛先驗概率分布，在上述問題3的基礎上再進一步，提出問題4。

為了方便計算，我們記

根據問題中的信息，我們可以推理得到p1,p2在X;Y這m+n個數中分別成為了第k1+m1,k1+k2+m1+m2大的數。后驗分布p應該為

同樣的，按照貝葉斯推理的邏輯，可將上述過程整理如下：

• 1 我們要猜測參數P=(p1,p2,p3)，其先驗分布為Dir(p|k);
• 2 數據Yi落到三個區間[0,p1),[p1,p2],(p2,1]的個數分別是m1,m2,m3,所以m=(m1,m2,m3)服從多項分布Mult(m|p);
• 3 在給定了來自數據提供的知識m后，p的后驗分布變為Dir(P|k+m)

上述貝葉斯分析過程的直觀表述為：

Dir(p|k) + Multcount(m) = Dir(p|k+m)

針對于這種觀測到的數據符合多項分布，參數的先驗分布和后驗分布都是Dirichlet分布的情況，就是Dirichlet-Multinomial共軛。這意味著，如果我們為多項分布的參數p選取的先驗分布是Dirichlet分布，那么以p為參數的多項分布用貝葉斯估計得到的后驗分布仍然服從Dirichlet分布。

1.7 Beta和Dirichlet分布的一個性質

如果p=Beta(t|alpha,beta)，那么

上式右邊的積分對應到概率分布Beta(t|alpha+1,beta)，對于這個分布，我們有

把上式帶人E(p)的計算式，可以得到：

這說明，對于Beta分布的隨機變量，其期望可以用上式來估計。Dirichlet分布也有類似的結論。對于p=Dir(t|alpha)，有

這個結論在后文的推導中會用到。

1.8 總結

LDA涉及的數學知識較多，需要認真體會，以上大部分的知識來源于文獻【1,2,3】,如有不清楚的地方，參見這些文獻以了解更多。

2 主題模型LDA

在介紹LDA之前，我們先介紹幾個基礎模型：Unigram model、mixture of unigrams model、pLSA model。為了方便描述，首先定義一些變量：

• 1 w表示詞，V表示所有詞的個數
• 2 z表示主題，k表示主題的個數
• 3 表示語料庫，M表示語料庫中的文檔數。
• 4 表示文檔，N表示文檔中詞的個數。

2.1 一元模型(Unigram model)

對于文檔

,用 表示 的先驗概率，生成文檔W的概率為：

其圖模型為（圖中被涂色的w表示可觀測變量，N表示一篇文檔中總共N個單詞，M表示M篇文檔）：

2.2 混合一元模型(Mixture of unigrams model)

該模型的生成過程是：給某個文檔先選擇一個主題Z，再根據該主題生成文檔，該文檔中的所有詞都來自一個主題。生成文檔的概率為：

其圖模型為（圖中被涂色的w表示可觀測變量，未被涂色的z表示未知的隱變量，N表示一篇文檔中總共N個單詞，M表示M篇文檔）：

2.3 pLSA模型

在混合一元模型中，假定一篇文檔只由一個主題生成，可實際中，一篇文章往往有多個主題，只是這多個主題各自在文檔中出現的概率大小不一樣。在pLSA中，假設文檔由多個主題生成。下面通過一個投色子的游戲（取自文獻【2】的例子）說明pLSA生成文檔的過程。

首先，假定你一共有K個可選的主題，有V個可選的詞。假設你每寫一篇文檔會制作一顆K面的“文檔-主題”骰子（扔此骰子能得到K個主題中的任意一個），和K個V面的“主題-詞項”骰子（每個骰子對應一個主題，K個骰子對應之前的K個主題，且骰子的每一面對應要選擇的詞項，V個面對應著V個可選的詞）。 比如可令K=3，即制作1個含有3個主題的“文檔-主題”骰子，這3個主題可以是：教育、經濟、交通。然后令V = 3，制作3個有著3面的“主題-詞項”骰子，其中，教育主題骰子的3個面上的詞可以是：大學、老師、課程，經濟主題骰子的3個面上的詞可以是：市場、企業、金融，交通主題骰子的3個面上的詞可以是：高鐵、汽車、飛機。

其次，每寫一個詞，先扔該“文檔-主題”骰子選擇主題，得到主題的結果后，使用和主題結果對應的那顆“主題-詞項”骰子，扔該骰子選擇要寫的詞。先扔“文檔-主題”的骰子，假設以一定的概率得到的主題是：教育，所以下一步便是扔教育主題篩子，以一定的概率得到教育主題篩子對應的某個詞大學。

• 上面這個投骰子產生詞的過程簡化一下便是：“先以一定的概率選取主題，再以一定的概率選取詞”。事實上，一開始可供選擇的主題有3個：教育、經濟、交通，那為何偏偏選取教育這個主題呢？其實是隨機選取的，只是這個隨機遵循一定的概率分布。比如可能選取教育主題的概率是0.5，選取經濟主題的概率是0.3，選取交通主題的概率是0.2，那么這3個主題的概率分布便是{教育：0.5，經濟：0.3，交通：0.2}，我們把各個主題z在文檔d中出現的概率分布稱之為主題分布，且是一個多項分布。
• 同樣的，從主題分布中隨機抽取出教育主題后，依然面對著3個詞：大學、老師、課程，這3個詞都可能被選中，但它們被選中的概率也是不一樣的。比如大學這個詞被選中的概率是0.5，老師這個詞被選中的概率是0.3，課程被選中的概率是0.2，那么這3個詞的概率分布便是{大學：0.5，老師：0.3，課程：0.2}，我們把各個詞語w在主題z下出現的概率分布稱之為詞分布，這個詞分布也是一個多項分布。
• 所以，選主題和選詞都是兩個隨機的過程，先從主題分布{教育：0.5，經濟：0.3，交通：0.2}中抽取出主題：教育，然后從該主題對應的詞分布{大學：0.5，老師：0.3，課程：0.2}中抽取出詞：大學。

最后，你不停的重復扔“文檔-主題”骰子和”主題-詞項“骰子，重復N次（產生N個詞），完成一篇文檔，重復這產生一篇文檔的方法M次，則完成M篇文檔。

上述過程抽象出來即是pLSA的文檔生成模型。在這個過程中，我們并未關注詞和詞之間的出現順序，所以pLSA是一種詞袋模型。定義如下變量:

• 表示隱藏的主題；
• 表示海量文檔中某篇文檔被選中的概率；
• 表示詞 在文檔 中出現的概率；針對海量文檔，對所有文檔進行分詞后，得到一個詞匯列表，這樣每篇文檔就是一個詞語的集合。對于每個詞語，用它在文檔中出現的次數除以文檔中詞語總的數目便是它在文檔中出現的概率；
• 表示主題 在文檔 中出現的概率；
• 表示詞 在主題 中出現的概率。與主題關系越密切的詞其條件概率越大。

我們可以按照如下的步驟得到“文檔-詞項”的生成模型：

• 1 按照 選擇一篇文檔 ；
• 2 選定文檔 之后，從主題分布中按照概率 選擇主題；
• 3 選定主題后，從詞分布中按照概率 P(w_{j}|z_{k})選擇一個詞。

利用看到的文檔推斷其隱藏的主題（分布）的過程，就是主題建模的目的：自動地發現文檔集中的主題（分布）。文檔d和單詞w是可被觀察到的，但主題z卻是隱藏的。如下圖所示（圖中被涂色的d、w表示可觀測變量，未被涂色的z表示未知的隱變量，N表示一篇文檔中總共N個單詞，M表示M篇文檔）。

上圖中，文檔d和詞w是我們得到的樣本，可觀測得到，所以對于任意一篇文檔，其

是已知的。根據這個概率可以訓練得到文檔-主題概率以及主題-詞項概率。即：

故得到文檔中每個詞的生成概率為：

可以直接得出，而 和 未知，所以 就是我們要估計的參數,我們要最大化這個參數。因為該待估計的參數中含有隱變量z，所以我們可以用EM算法來估計這個參數。

2.4 LDA模型

LDA的不同之處在于，pLSA的主題的概率分布P(c|d)是一個確定的概率分布，也就是雖然主題c不確定，但是c符合的概率分布是確定的，比如符合某個多項分布，這個多項分布的各參數是確定的。 但是在LDA中，這個多項分布都是不確定的，高斯分布又服從一個狄利克雷先驗分布(Dirichlet prior)。即LDA就是pLSA的貝葉斯版本,正因為LDA被貝葉斯化了，所以才會加的兩個先驗參數。

LDA模型中一篇文檔生成的方式如下所示:

• 1 按照 選擇一篇文檔 ；
• 2 從狄利克雷分布 中生成文檔i的主題分布 ；
• 3 從主題的多項式分布 中取樣生成文檔i第j個詞的主題 ；
• 4 從狄利克雷分布 中取樣生成主題 對應的詞語分布 ；
• 5 從詞語的多項式分布 中采樣最終生成詞語

從上面的過程可以看出，LDA在pLSA的基礎上，為主題分布和詞分布分別加了兩個Dirichlet先驗。

拿之前講解pLSA的例子進行具體說明。如前所述，在pLSA中，選主題和選詞都是兩個隨機的過程，先從主題分布{教育：0.5，經濟：0.3，交通：0.2}中抽取出主題：教育，然后從該主題對應的詞分布{大學：0.5，老師：0.3，課程：0.2}中抽取出詞：大學。 在LDA中，選主題和選詞依然都是兩個隨機的過程。但在LDA中，主題分布和詞分布不再唯一確定不變，即無法確切給出。例如主題分布可能是{教育：0.5，經濟：0.3，交通：0.2}，也可能是{教育：0.6，經濟：0.2，交通：0.2}，到底是哪個我們不能確定，因為它是隨機的可變化的。 但再怎么變化，也依然服從一定的分布，主題分布和詞分布由Dirichlet先驗確定。

舉個文檔d產生主題z的例子。

在pLSA中，給定一篇文檔d，主題分布是一定的，比如{ P(zi|d), i = 1,2,3 }可能就是{0.4,0.5,0.1}，表示z1、z2、z3這3個主題被文檔d選中的概率都是個固定的值：P(z1|d) = 0.4、P(z2|d) = 0.5、P(z3|d) = 0.1，如下圖所示:

在LDA中，主題分布（各個主題在文檔中出現的概率分布）和詞分布（各個詞語在某個主題下出現的概率分布）是唯一確定的。LDA為提供了兩個Dirichlet先驗參數，Dirichlet先驗為某篇文檔隨機抽取出主題分布和詞分布。

給定一篇文檔d，現在有多個主題z1、z2、z3，它們的主題分布{ P(zi|d), i = 1,2,3 }可能是{0.4,0.5,0.1}，也可能是{0.2,0.2,0.6}，即這些主題被d選中的概率都不再是確定的值，可能是P(z1|d) = 0.4、P(z2|d) = 0.5、P(z3|d) = 0.1，也有可能是P(z1|d) = 0.2、P(z2|d) = 0.2、P(z3|d) = 0.6，而主題分布到底是哪個取值集合我們不確定，但其先驗分布是dirichlet分布，所以可以從無窮多個主題分布中按照dirichlet先驗隨機抽取出某個主題分布出來。如下圖所示

LDA在pLSA的基礎上給兩參數

和 加了兩個先驗分布的參數。這兩個分布都是Dirichlet分布。 下面是LDA的圖模型結構：

3 LDA 參數估計

在spark中，提供了兩種方法來估計參數，分別是變分EM（期望最大）算法（見文獻【3】【4】）和在線學習算法（見文獻【5】）。下面將分別介紹這兩種算法以及其源碼實現。

3.1 變分EM算法

變分貝葉斯算法的詳細信息可以參考文獻【9】。

在上文中，我們知道LDA將變量theta和phi（為了方便起見，我們將上文LDA圖模型中的beta改為了phi）看做隨機變量，并且為theta添加一個超參數為alpha的Dirichlet先驗，為phi添加一個超參數為eta的Dirichlet先驗來估計theta和beta的最大后驗（MAP）。 可以通過最優化最大后驗估計來估計參數。我們首先來定義幾個變量：

• 下式的gamma表示詞為w，文檔為j時，主題為k的概率，如公式**（3.1.1）**

• 表示詞w在文檔j中出現的次數；
• 表示詞w在主題k中出現的次數，如公式**（3.1.2）**

• 表示主題k在文檔j中出現的次數，如公式**（3.1.3）**

• 表示主題k中包含的詞出現的總次數，如公式**（3.1.4）**

• 表示文檔j中包含的主題出現的總次數,如公式**（3.1.5）**

根據文獻【4】中2.2章節的介紹，我們可以推導出如下更新公式**(3.1.6)**，其中alpha和eta均大于1：

收斂之后，最大后驗估計可以得到公式**(3.1.7)**：

變分EM算法的流程如下：

• 1 初始化狀態，即隨機初始化 和
• 2 E-步，對每一個（文檔，詞匯）對i，計算 ，更新gamma值
• 3 M-步，計算隱藏變量phi和theta。即計算 和
• 4 重復以上2、3兩步，直到滿足最大迭代數

第4.2章會從代碼層面說明該算法的實現流程。

3.2 在線學習算法

3.2.1 批量變分貝葉斯

在變分貝葉斯推導(VB)中，根據文獻【3】，使用一種更簡單的分布q(z,theta,beta)來估計真正的后驗分布，這個簡單的分布使用一組自由變量(free parameters)來定義。 通過最大化對數似然的一個下界（Evidence Lower Bound (ELBO)）來最優化這些參數，如下公式**(3.2.1)**

最大化ELBO就是最小化q(z,theta,beta)和p(z,theta,beta|w,alpha,eta)的KL距離。根據文獻【3】，我們將q因式分解為如下**（3.2.2）**的形式：

后驗z通過phi來參數化，后驗theta通過gamma來參數化，后驗beta通過lambda來參數化。為了簡單描述，我們把lambda當作“主題”來看待。公式**(3.2.2)分解為如下(3.2.3)**形式：

我們現在將上面的期望擴展為變分參數的函數形式。這反映了變分目標只依賴于

，即詞w出現在文檔d中的次數。當使用VB算法時，文檔可以通過它們的詞頻來匯總（summarized），如公式**(3.2.4)**

上面的公式中，W表示詞的數量，D表示文檔的數量。l表示文檔d對ELBO所做的貢獻。L可以通過坐標上升法來最優化，它的更新公式如**(3.2.5)**:

log(theta)和log(beta)的期望通過下面的公式**(3.2.6)**計算：

通過EM算法，我們可以將這些更新分解成E-步和M-步。E-步固定lambda來更新gamma和phi；M-步通過給定phi來更新lambda。批VB算法的過程如下**(算法1)**所示：

3.2.2 在線變分貝葉斯

批量變分貝葉斯算法需要固定的內存，并且比吉布斯采樣更快。但是它仍然需要在每次迭代時處理所有的文檔，這在處理大規模文檔時，速度會很慢，并且也不適合流式數據的處理。 文獻【5】提出了一種在線變分推導算法。設定gamma(n_d,lambda)和phi(n_d,lambda)分別表示gamma_d和phi_d的值，我們的目的就是設定phi來最大化下面的公式**(3.2.7)**

我們在算法2中介紹了在線VB算法。因為詞頻的第t個向量

是可觀察的，我們在E-步通過固定lambda來找到gamma_t和phi_t的局部最優解。 然后，我們計算lambda_cap。如果整個語料庫由單個文檔重復D次組成，那么這樣的lambda_cap設置是最優的。之后，我們通過lambda之前的值以及lambda_cap來更新lambda。我們給lambda_cap設置的權重如公式**(3.2.8)**所示：

在線VB算法的實現流程如下算法2所示

那么在在線VB算法中，alpha和eta是如何更新的呢？參考文獻【8】提供了計算方法。給定數據集，dirichlet參數的可以通過最大化下面的對數似然來估計

其中，

有多種方法可以最大化這個目標函數，如梯度上升，Newton-Raphson等。Spark使用Newton-Raphson方法估計參數，更新alpha。Newton-Raphson提供了一種參數二次收斂的方法， 它一般的更新規則如下公式**(3.3.3)**:

其中，H表示海森矩陣。對于這個特別的對數似然函數，可以應用Newton-Raphson去解決高維數據，因為它可以在線性時間求出海森矩陣的逆矩陣。一般情況下，海森矩陣可以用一個對角矩陣和一個元素都一樣的矩陣的和來表示。 如下公式**(3.3.4)**，Q是對角矩陣，C11是元素相同的一個矩陣。

為了計算海森矩陣的逆矩陣，我們觀察到，對任意的可逆矩陣Q和非負標量c，有下列式子**(3.3.5)**:

因為Q是對角矩陣，所以Q的逆矩陣可以很容易的計算出來。所以Newton-Raphson的更新規則可以重寫為如下**(3.3.6)**的形式

其中b如下公式**(3.3.7)**，

4 LDA代碼實現

4.1 LDA使用實例

我們從官方文檔【6】給出的使用代碼為起始點來詳細分析LDA的實現。

import org.apache.spark.mllib.clustering.{LDA, DistributedLDAModel} import org.apache.spark.mllib.linalg.Vectors // 加載和處理數據 val data = sc.textFile("data/mllib/sample_lda_data.txt") val parsedData = data.map(s => Vectors.dense(s.trim.split(' ').map(_.toDouble))) // 為文檔編號，編號唯一。List（（id，vector）....） val corpus = parsedData.zipWithIndex.map(_.swap).cache() // 將文檔聚類為3類 val ldaModel = new LDA().setK(3).run(corpus) val topics = ldaModel.topicsMatrix for (topic <- Range(0, 3)) {print("Topic " + topic + ":")for (word <- Range(0, ldaModel.vocabSize)) { print(" " + topics(word, topic)); }println() }

以上代碼主要做了兩件事：加載和切分數據、訓練模型。在樣本數據中，每一行代表一篇文檔，經過處理后，corpus的類型為List((id,vector)*)，一個(id,vector)代表一篇文檔。將處理后的數據傳給org.apache.spark.mllib.clustering.LDA類的run方法， 就可以開始訓練模型。run方法的代碼如下所示：

def run(documents: RDD[(Long, Vector)]): LDAModel = {val state = ldaOptimizer.initialize(documents, this)var iter = 0val iterationTimes = Array.fill[Double](maxIterations)(0)while (iter < maxIterations) {val start = System.nanoTime()state.next()val elapsedSeconds = (System.nanoTime() - start) / 1e9iterationTimes(iter) = elapsedSecondsiter += 1}state.getLDAModel(iterationTimes)}

這段代碼首先調用initialize方法初始化狀態信息，然后循環迭代調用next方法直到滿足最大的迭代次數。在我們沒有指定的情況下，迭代次數默認為20。需要注意的是， ldaOptimizer有兩個具體的實現類EMLDAOptimizer和OnlineLDAOptimizer，它們分別表示使用EM算法和在線學習算法實現參數估計。在未指定的情況下，默認使用EMLDAOptimizer。

4.2 變分EM算法的實現

在spark中，使用GraphX來實現EMLDAOptimizer，這個圖是有兩種類型的頂點的二分圖。這兩類頂點分別是文檔頂點（Document vertices）和詞頂點（Term vertices）。

• 文檔頂點使用大于0的唯一的指標來索引，保存長度為k（主題個數）的向量
• 詞頂點使用{-1, -2, ..., -vocabSize}來索引，保存長度為k（主題個數）的向量
• 邊（edges）對應詞出現在文檔中的情況。邊的方向是document -> term，并且根據document進行分區

我們可以根據3.1節中介紹的算法流程來解析源代碼。

4.2.1 初始化狀態

spark在EMLDAOptimizer的initialize方法中實現初始化功能。包括初始化Dirichlet參數alpha和eta、初始化邊、初始化頂點以及初始化圖。

//對應超參數alphaval docConcentration = lda.getDocConcentration//對應超參數etaval topicConcentration = lda.getTopicConcentrationthis.docConcentration = if (docConcentration == -1) (50.0 / k) + 1.0 else docConcentrationthis.topicConcentration = if (topicConcentration == -1) 1.1 else topicConcentration

上面的代碼初始化了超參數alpha和eta，根據文獻【4】，當alpha未指定時，初始化其為(50.0 / k) + 1.0，其中k表示主題個數。當eta未指定時，初始化其為1.1。

//對于每個文檔，為每一個唯一的Term創建一個(Document->Term)的邊 val edges: RDD[Edge[TokenCount]] = docs.flatMap { case (docID: Long, termCounts: Vector) =>// Add edges for terms with non-zero counts.termCounts.toBreeze.activeIterator.filter(_._2 != 0.0).map { case (term, cnt) =>//文檔id，termindex，詞頻Edge(docID, term2index(term), cnt)} } //term2index將term轉為{-1, -2, ..., -vocabSize}索引private[clustering] def term2index(term: Int): Long = -(1 + term.toLong)

上面的這段代碼處理每個文檔，對文檔中每個唯一的Term（詞）創建一個邊，邊的格式為(文檔id，詞索引，詞頻)。詞索引為{-1, -2, ..., -vocabSize}。

//創建頂點val docTermVertices: RDD[(VertexId, TopicCounts)] = {val verticesTMP: RDD[(VertexId, TopicCounts)] =edges.mapPartitionsWithIndex { case (partIndex, partEdges) =>val random = new Random(partIndex + randomSeed)partEdges.flatMap { edge =>val gamma = normalize(BDV.fill[Double](k)(random.nextDouble()), 1.0)//此處的sum是DenseVector，gamma*N_wjval sum = gamma * edge.attr//srcId表示文獻id，dstId表是詞索引Seq((edge.srcId, sum), (edge.dstId, sum))}}verticesTMP.reduceByKey(_ + _)}

上面的代碼創建頂點。我們為每個主題隨機初始化一個值，即gamma是隨機的。sum為gamma * edge.attr，這里的edge.attr即N_wj,所以sum用gamma * N_wj作為頂點的初始值。

this.graph = Graph(docTermVertices, edges).partitionBy(PartitionStrategy.EdgePartition1D)

上面的代碼初始化Graph并通過文檔分區。

4.2.2 E-步：更新gamma

val eta = topicConcentrationval W = vocabSizeval alpha = docConcentrationval N_k = globalTopicTotalsval sendMsg: EdgeContext[TopicCounts, TokenCount, (Boolean, TopicCounts)] => Unit =(edgeContext) => {// 計算 N_{wj} gamma_{wjk}val N_wj = edgeContext.attr// E-STEP: 計算 gamma_{wjk} 通過N_{wj}來計算//此處的edgeContext.srcAttr為當前迭代的N_kj , edgeContext.dstAttr為當前迭代的N_wk,//后面通過M-步，會更新這兩個值,作為下一次迭代的當前值val scaledTopicDistribution: TopicCounts = computePTopic(edgeContext.srcAttr, edgeContext.dstAttr, N_k, W, eta, alpha) *= N_wjedgeContext.sendToDst((false, scaledTopicDistribution))edgeContext.sendToSrc((false, scaledTopicDistribution))}

上述代碼中，W表示詞數，N_k表示所有文檔中，出現在主題k中的詞的詞頻總數，后續的實現會使用方法computeGlobalTopicTotals來更新這個值。N_wj表示詞w出現在文檔j中的詞頻數，為已知數。E-步就是利用公式**(3.1.6)**去更新gamma。 代碼中使用computePTopic方法來實現這個更新。edgeContext通過方法sendToDst將scaledTopicDistribution發送到目標頂點， 通過方法sendToSrc發送到源頂點以便于后續的M-步更新的N_kj和N_wk。下面我們看看computePTopic方法。

private[clustering] def computePTopic(docTopicCounts: TopicCounts,termTopicCounts: TopicCounts,totalTopicCounts: TopicCounts,vocabSize: Int,eta: Double,alpha: Double): TopicCounts = {val K = docTopicCounts.lengthval N_j = docTopicCounts.dataval N_w = termTopicCounts.dataval N = totalTopicCounts.dataval eta1 = eta - 1.0val alpha1 = alpha - 1.0val Weta1 = vocabSize * eta1var sum = 0.0val gamma_wj = new Array[Double](K)var k = 0while (k < K) {val gamma_wjk = (N_w(k) + eta1) * (N_j(k) + alpha1) / (N(k) + Weta1)gamma_wj(k) = gamma_wjksum += gamma_wjkk += 1}// normalizeBDV(gamma_wj) /= sum}

這段代碼比較簡單，完全按照公式**(3.1.6)**表示的樣子來實現。val gamma_wjk = (N_w(k) + eta1) * (N_j(k) + alpha1) / (N(k) + Weta1)就是實現的更新邏輯。

4.2.3 M-步：更新phi和theta

// M-STEP: 聚合計算新的 N_{kj}, N_{wk} counts. val docTopicDistributions: VertexRDD[TopicCounts] =graph.aggregateMessages[(Boolean, TopicCounts)](sendMsg, mergeMsg).mapValues(_._2)

我們由公式**(3.1.7)**可知，更新隱藏變量phi和theta就是更新相應的N_kj和N_wk。聚合更新使用aggregateMessages方法來實現。請參考文獻【7】來了解該方法的作用。

4.3 在線變分算法的代碼實現

4.3.1 初始化狀態

在線學習算法首先使用方法initialize方法初始化參數值

override private[clustering] def initialize(docs: RDD[(Long, Vector)],lda: LDA): OnlineLDAOptimizer = {this.k = lda.getKthis.corpusSize = docs.count()this.vocabSize = docs.first()._2.sizethis.alpha = if (lda.getAsymmetricDocConcentration.size == 1) {if (lda.getAsymmetricDocConcentration(0) == -1) Vectors.dense(Array.fill(k)(1.0 / k))else {require(lda.getAsymmetricDocConcentration(0) >= 0,s"all entries in alpha must be >=0, got: $alpha")Vectors.dense(Array.fill(k)(lda.getAsymmetricDocConcentration(0)))}} else {require(lda.getAsymmetricDocConcentration.size == k,s"alpha must have length k, got: $alpha")lda.getAsymmetricDocConcentration.foreachActive { case (_, x) =>require(x >= 0, s"all entries in alpha must be >= 0, got: $alpha")}lda.getAsymmetricDocConcentration}this.eta = if (lda.getTopicConcentration == -1) 1.0 / k else lda.getTopicConcentrationthis.randomGenerator = new Random(lda.getSeed)this.docs = docs// 初始化變分分布 q(beta|lambda)this.lambda = getGammaMatrix(k, vocabSize)this.iteration = 0this}

根據文獻【5】，alpha和eta的值大于等于0，并且默認為1.0/k。上文使用getGammaMatrix方法來初始化變分分布q(beta|lambda)。

private def getGammaMatrix(row: Int, col: Int): BDM[Double] = {val randBasis = new RandBasis(new org.apache.commons.math3.random.MersenneTwister(randomGenerator.nextLong()))//初始化一個gamma分布val gammaRandomGenerator = new Gamma(gammaShape, 1.0 / gammaShape)(randBasis)val temp = gammaRandomGenerator.sample(row * col).toArraynew BDM[Double](col, row, temp).t}

getGammaMatrix方法使用gamma分布初始化一個隨機矩陣。

4.3.2 更新參數

override private[clustering] def next(): OnlineLDAOptimizer = {//返回文檔集中采樣的子集//默認情況下，文檔可以被采樣多次，且采樣比例是0.05val batch = docs.sample(withReplacement = sampleWithReplacement, miniBatchFraction,randomGenerator.nextLong())if (batch.isEmpty()) return thissubmitMiniBatch(batch)}

以上的next方法首先對文檔進行采樣，然后調用submitMiniBatch對采樣的文檔子集進行處理。下面我們詳細分解submitMiniBatch方法。

• 1 計算log(beta)的期望，并將其作為廣播變量廣播到集群中

val expElogbeta = exp(LDAUtils.dirichletExpectation(lambda)).t //廣播變量 val expElogbetaBc = batch.sparkContext.broadcast(expElogbeta) //參數alpha是dirichlet參數 private[clustering] def dirichletExpectation(alpha: BDM[Double]): BDM[Double] = {val rowSum = sum(alpha(breeze.linalg.*, ::))val digAlpha = digamma(alpha)val digRowSum = digamma(rowSum)val result = digAlpha(::, breeze.linalg.*) - digRowSumresult}

上述代碼調用exp(LDAUtils.dirichletExpectation(lambda))方法實現參數為lambda的log beta的期望。實現原理參見公式**(3.2.6)**。

• 2 計算phi以及gamma，即算法2中的E-步

//對采樣文檔進行分區處理 val stats: RDD[(BDM[Double], List[BDV[Double]])] = batch.mapPartitions { docs =>//val nonEmptyDocs = docs.filter(_._2.numNonzeros > 0)val stat = BDM.zeros[Double](k, vocabSize)var gammaPart = List[BDV[Double]]()nonEmptyDocs.foreach { case (_, termCounts: Vector) =>val ids: List[Int] = termCounts match {case v: DenseVector => (0 until v.size).toListcase v: SparseVector => v.indices.toList}val (gammad, sstats) = OnlineLDAOptimizer.variationalTopicInference(termCounts, expElogbetaBc.value, alpha, gammaShape, k)stat(::, ids) := stat(::, ids).toDenseMatrix + sstatsgammaPart = gammad :: gammaPart}Iterator((stat, gammaPart))}

上面的代碼調用OnlineLDAOptimizer.variationalTopicInference實現算法2中的E-步,迭代計算phi和gamma。

private[clustering] def variationalTopicInference(termCounts: Vector,expElogbeta: BDM[Double],alpha: breeze.linalg.Vector[Double],gammaShape: Double,k: Int): (BDV[Double], BDM[Double]) = {val (ids: List[Int], cts: Array[Double]) = termCounts match {case v: DenseVector => ((0 until v.size).toList, v.values)case v: SparseVector => (v.indices.toList, v.values)}// 初始化變分分布 q(theta|gamma) val gammad: BDV[Double] = new Gamma(gammaShape, 1.0 / gammaShape).samplesVector(k) // K//根據公式（3.2.6）計算 E(log theta)val expElogthetad: BDV[Double] = exp(LDAUtils.dirichletExpectation(gammad)) // Kval expElogbetad = expElogbeta(ids, ::).toDenseMatrix // ids * K//根據公式（3.2.5）計算phi，這里加1e-100表示并非嚴格等于val phiNorm: BDV[Double] = expElogbetad * expElogthetad :+ 1e-100 // idsvar meanGammaChange = 1Dval ctsVector = new BDV[Double](cts) // ids// 迭代直至收斂while (meanGammaChange > 1e-3) {val lastgamma = gammad.copy//依據公式(3.2.5)計算gammagammad := (expElogthetad :* (expElogbetad.t * (ctsVector :/ phiNorm))) :+ alpha//根據更新的gamma，計算E(log theta)expElogthetad := exp(LDAUtils.dirichletExpectation(gammad))// 更新phiphiNorm := expElogbetad * expElogthetad :+ 1e-100//計算兩次gamma的差值meanGammaChange = sum(abs(gammad - lastgamma)) / k}val sstatsd = expElogthetad.asDenseMatrix.t * (ctsVector :/ phiNorm).asDenseMatrix(gammad, sstatsd)}

3 更新lambda

val statsSum: BDM[Double] = stats.map(_._1).reduce(_ += _)val gammat: BDM[Double] = breeze.linalg.DenseMatrix.vertcat(stats.map(_._2).reduce(_ ++ _).map(_.toDenseMatrix): _*)val batchResult = statsSum :* expElogbeta.t// 更新lambda和alphaupdateLambda(batchResult, (miniBatchFraction * corpusSize).ceil.toInt)

updateLambda方法實現算法2中的M-步,更新lambda。實現代碼如下：

private def updateLambda(stat: BDM[Double], batchSize: Int): Unit = {// 根據公式（3.2.8）計算權重val weight = rho()// 更新lambda，其中stat * (corpusSize.toDouble / batchSize.toDouble)+eta表示rho_caplambda := (1 - weight) * lambda +weight * (stat * (corpusSize.toDouble / batchSize.toDouble) + eta)} // 根據公式（3.2.8）計算rho private def rho(): Double = {math.pow(getTau0 + this.iteration, -getKappa)}

4 更新alpha

private def updateAlpha(gammat: BDM[Double]): Unit = {//計算rhoval weight = rho()val N = gammat.rows.toDoubleval alpha = this.alpha.toBreeze.toDenseVector//計算log p_hatval logphat: BDM[Double] = sum(LDAUtils.dirichletExpectation(gammat)(::, breeze.linalg.*)) / N//計算梯度為N（-phi(alpha)+log p_hat）val gradf = N * (-LDAUtils.dirichletExpectation(alpha) + logphat.toDenseVector)//計算公式（3.3.4）中的c，trigamma表示gamma函數的二階導數val c = N * trigamma(sum(alpha))//計算公式（3.3.4）中的qval q = -N * trigamma(alpha)//根據公式(3.3.7)計算bval b = sum(gradf / q) / (1D / c + sum(1D / q))val dalpha = -(gradf - b) / qif (all((weight * dalpha + alpha) :> 0D)) {alpha :+= weight * dalphathis.alpha = Vectors.dense(alpha.toArray)}}

5 參考文獻

【1】LDA數學八卦

【2】通俗理解LDA主題模型

【3】[Latent Dirichlet Allocation](docs/Latent Dirichlet Allocation.pdf)

【4】[On Smoothing and Inference for Topic Models](docs/On Smoothing and Inference for Topic Models.pdf)

【5】[Online Learning for Latent Dirichlet Allocation](docs/Online Learning for Latent Dirichlet Allocation.pdf)

【6】Spark官方文檔

【7】Spark GraphX介紹

【8】Maximum Likelihood Estimation of Dirichlet Distribution Parameters

【9】Variational Bayes

星環科技：機器學習算法2020線上訓練營 | 首營即將開班，限額報名！?zhuanlan.zhihu.com

總結

以上是生活随笔為你收集整理的狄利克雷分布公式_深入机器学习系列11－隐式狄利克雷分布的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。