文本主题模型之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算法在这里的使用，我们的模型里面有隐藏变量θ,β,zθ,β,z，模型的参数是α,ηα,η。为了求出模型参数和对应的隐藏变量分布，EM算法需要在E步先求出隐藏变量θ,β,zθ,β,z的基于条件概率分布的期望，接着在M步极大化这个期望，得到更新的后验模型参数α,ηα,η。

　　　　问题是在EM算法的E步，由于θ,β,zθ,β,z的耦合，我们难以求出隐藏变量θ,β,zθ,β,z的条件概率分布，也难以求出对应的期望，需要“变分推断“来帮忙，这里所谓的变分推断，也就是在隐藏变量存在耦合的情况下，我们通过变分假设，即假设所有的隐藏变量都是通过各自的独立分布形成的，这样就去掉了隐藏变量之间的耦合关系。我们用各个独立分布形成的变分分布来模拟近似隐藏变量的条件分布，这样就可以顺利的使用EM算法了。

　　　　当进行若干轮的E步和M步的迭代更新之后，我们可以得到合适的近似隐藏变量分布θ,β,zθ,β,z和模型后验参数α,ηα,η，进而就得到了我们需要的LDA文档主题分布和主题词分布。

　　　　可见要完全理解LDA的变分推断EM算法，需要搞清楚它在E步变分推断的过程和推断完毕后EM算法的过程。

 2. LDA的变分推断思路

　　　　要使用EM算法，我们需要求出隐藏变量的条件概率分布如下：

p(θ,β,z|w,α,η)=p(θ,β,z,w|α,η)p(w|α,η)p(θ,β,z|w,α,η)=p(θ,β,z,w|α,η)p(w|α,η)

 

　　　　前面讲到由于θ,β,zθ,β,z之间的耦合，这个条件概率是没法直接求的，但是如果不求它就不能用EM算法了。怎么办呢，我们引入变分推断，具体是引入基于mean field assumption的变分推断，这个推断假设所有的隐藏变量都是通过各自的独立分布形成的，如下图所示：

　　　　我们假设隐藏变量θθ是由独立分布γγ形成的，隐藏变量zz是由独立分布ϕϕ形成的，隐藏变量ββ是由独立分布λλ形成的。这样我们得到了三个隐藏变量联合的变分分布qq为：

q(β,z,θ|λ,ϕ,γ)=∏k=1Kq(βk|λk)∏d=1Mq(θd,zd|γd,ϕd)=∏k=1Kq(βk|λk)∏d=1M(q(θd|γd)∏n=1Ndq(zdn|ϕdn))(1)(2)(1)q(β,z,θ|λ,ϕ,γ)=∏k=1Kq(βk|λk)∏d=1Mq(θd,zd|γd,ϕd)(2)=∏k=1Kq(βk|λk)∏d=1M(q(θd|γd)∏n=1Ndq(zdn|ϕdn))

 

　　　　我们的目标是用q(β,z,θ|λ,ϕ,γ)q(β,z,θ|λ,ϕ,γ)来近似的估计p(θ,β,z|w,α,η)p(θ,β,z|w,α,η)，也就是说需要这两个分布尽可能的相似，用数学语言来描述就是希望这两个概率分布之间有尽可能小的KL距离，即：

(λ∗,ϕ∗,γ∗)=argminλ,ϕ,γD(q(β,z,θ|λ,ϕ,γ)||p(θ,β,z|w,α,η))(λ∗,ϕ∗,γ∗)=argmin⏟λ,ϕ,γD(q(β,z,θ|λ,ϕ,γ)||p(θ,β,z|w,α,η))

 

　　　　其中D(q||p)D(q||p)即为KL散度或KL距离，对应分布qq和pp的交叉熵。即：

D(q||p)=∑xq(x)logq(x)p(x)=Eq(x)(logq(x)−logp(x))D(q||p)=∑xq(x)logq(x)p(x)=Eq(x)(logq(x)−logp(x))

 

　　　　我们的目的就是找到合适的λ∗,ϕ∗,γ∗λ∗,ϕ∗,γ∗,然后用q(β,z,θ|λ∗,ϕ∗,γ∗)q(β,z,θ|λ∗,ϕ∗,γ∗)来近似隐藏变量的条件分布p(θ,β,z|w,α,η)p(θ,β,z|w,α,η)，进而使用EM算法迭代。

　　　　这个合适的λ∗,ϕ∗,γ∗λ∗,ϕ∗,γ∗,也不是那么好求的，怎么办呢？我们先看看我能文档数据的对数似然函数log(w|α,η)log(w|α,η)如下,为了简化表示，我们用Eq(x)Eq(x)代替Eq(β,z,θ|λ,ϕ,γ)(x)Eq(β,z,θ|λ,ϕ,γ)(x)，用来表示xx对于变分分布q(β,z,θ|λ,ϕ,γ)q(β,z,θ|λ,ϕ,γ) 的期望。

 

log(w|α,η)=log∫∫∑zp(θ,β,z,w|α,η)dθdβ=log∫∫∑zp(θ,β,z,w|α,η)q(β,z,θ|λ,ϕ,γ)q(β,z,θ|λ,ϕ,γ)dθdβ=logEqp(θ,β,z,w|α,η)q(β,z,θ|λ,ϕ,γ)≥Eqlogp(θ,β,z,w|α,η)q(β,z,θ|λ,ϕ,γ)=Eqlogp(θ,β,z,w|α,η)−Eqlogq(β,z,θ|λ,ϕ,γ)(3)(4)(5)(6)(7)(3)log(w|α,η)=log∫∫∑zp(θ,β,z,w|α,η)dθdβ(4)=log∫∫∑zp(θ,β,z,w|α,η)q(β,z,θ|λ,ϕ,γ)q(β,z,θ|λ,ϕ,γ)dθdβ(5)=logEqp(θ,β,z,w|α,η)q(β,z,θ|λ,ϕ,γ)(6)≥Eqlogp(θ,β,z,w|α,η)q(β,z,θ|λ,ϕ,γ)(7)=Eqlogp(θ,β,z,w|α,η)−Eqlogq(β,z,θ|λ,ϕ,γ)

 

　　　　其中，从第(5)式到第(6)式用到了Jensen不等式：

f(E(x))≥E(f(x))f(x)为凹函数f(E(x))≥E(f(x))f(x)为凹函数

 

　　　　我们一般把第(7)式记为：

L(λ,ϕ,γ;α,η)=Eqlogp(θ,β,z,w|α,η)−Eqlogq(β,z,θ|λ,ϕ,γ)L(λ,ϕ,γ;α,η)=Eqlogp(θ,β,z,w|α,η)−Eqlogq(β,z,θ|λ,ϕ,γ)

 

　　　　由于L(λ,ϕ,γ;α,η)L(λ,ϕ,γ;α,η)是我们的对数似然的一个下界（第6式），所以这个LL一般称为ELBO(Evidence Lower BOund)。那么这个ELBO和我们需要优化的的KL散度有什么关系呢？注意到：

D(q(β,z,θ|λ,ϕ,γ)||p(θ,β,z|w,α,η))=Eqlogq(β,z,θ|λ,ϕ,γ)−Eqlogp(θ,β,z|w,α,η)=Eqlogq(β,z,θ|λ,ϕ,γ)−Eqlogp(θ,β,z,w|α,η)p(w|α,η)=−L(λ,ϕ,γ;α,η)+log(w|α,η)(8)(9)(10)(8)D(q(β,z,θ|λ,ϕ,γ)||p(θ,β,z|w,α,η))=Eqlogq(β,z,θ|λ,ϕ,γ)−Eqlogp(θ,β,z|w,α,η)(9)=Eqlogq(β,z,θ|λ,ϕ,γ)−Eqlogp(θ,β,z,w|α,η)p(w|α,η)(10)=−L(λ,ϕ,γ;α,η)+log(w|α,η)

 

　　　　在(10)式中，由于对数似然部分和我们的KL散度无关，可以看做常量，因此我们希望最小化KL散度等价于最大化ELBO。那么我们的变分推断最终等价的转化为要求ELBO的最大值。现在我们开始关注于极大化ELBO并求出极值对应的变分参数λ,ϕ,γλ,ϕ,γ。

3. 极大化ELBO求解变分参数

　　　　为了极大化ELBO，我们首先对ELBO函数做一个整理如下：

L(λ,ϕ,γ;α,η)=Eq[logp(β|η)]+Eq[logp(z|θ)]+Eq[logp(θ|α)]+Eq[logp(w|z,β)]−Eq[logq(β|λ)]−Eq[logq(z|ϕ)]−Eq[logq(θ|γ)](11)(12)(13)(11)L(λ,ϕ,γ;α,η)=Eq[logp(β|η)]+Eq[logp(z|θ)]+Eq[logp(θ|α)](12)+Eq[logp(w|z,β)]−Eq[logq(β|λ)](13)−Eq[logq(z|ϕ)]−Eq[logq(θ|γ)]

 

　　　　可见展开后有7项，现在我们需要对这7项分别做一个展开。为了简化篇幅，这里只对第一项的展开做详细介绍。在介绍第一项的展开前，我们需要了解指数分布族的性质。指数分布族是指下面这样的概率分布：

p(x|θ)=h(x)exp(η(θ)∗T(x)−A(θ))p(x|θ)=h(x)exp(η(θ)∗T(x)−A(θ))

 

　　　　其中，A(x)A(x)为归一化因子，主要是保证概率分布累积求和后为1，引入指数分布族主要是它有下面这样的性质：

ddη(θ)A(θ)=Ep(x|θ)[T(x)]ddη(θ)A(θ)=Ep(x|θ)[T(x)]

 

　　　　这个证明并不复杂，这里不累述。我们的常见分布比如Gamma分布，Beta分布，Dirichlet分布都是指数分布族。有了这个性质，意味着我们在ELBO里面一大推的期望表达式可以转化为求导来完成，这个技巧大大简化了计算量。

　　　　回到我们ELBO第一项的展开如下：

Eq[logp(β|η)]=Eq[log∏k=1K(Γ(∑i=1Vηi)∏Vi=1Γ(ηi)∏i=1Vβηi−1ki)]=KlogΓ(∑i=1Vηi)−K∑i=1VlogΓ(ηi)+∑k=1KEq[∑i=1V(ηi−1)logβki](14)(15)(14)Eq[logp(β|η)]=Eq[log∏k=1K(Γ(∑i=1Vηi)∏i=1VΓ(ηi)∏i=1Vβkiηi−1)](15)=KlogΓ(∑i=1Vηi)−K∑i=1VlogΓ(ηi)+∑k=1KEq[∑i=1V(ηi−1)logβki]

 

　　　　第(15)式的第三项的期望部分，可以用上面讲到的指数分布族的性质，转化为一个求导过程。即：

Eq[∑i=1Vlogβki]=(logΓ(λki)−logΓ(∑i′=1Vλki′))′=Ψ(λki)−Ψ(∑i′=1Vλki′)Eq[∑i=1Vlogβki]=(logΓ(λki)−logΓ(∑i′=1Vλki′))′=Ψ(λki)−Ψ(∑i′=1Vλki′)

 

　　　　其中：

Ψ(x)=ddxlogΓ(x)=Γ′(x)Γ(x)Ψ(x)=ddxlogΓ(x)=Γ′(x)Γ(x)

 

　　　　最终，我们得到EBLO第一项的展开式为：

Eq[logp(β|η)]=KlogΓ(∑i=1Vηi)−K∑i=1VlogΓ(ηi)+∑k=1K∑i=1V(ηi−1)(Ψ(λki)−Ψ(∑i′=1Vλki′))(16)(16)Eq[logp(β|η)]=KlogΓ(∑i=1Vηi)−K∑i=1VlogΓ(ηi)+∑k=1K∑i=1V(ηi−1)(Ψ(λki)−Ψ(∑i′=1Vλki′))

 

　　　　类似的方法求解其他6项，可以得到ELBO的最终关于变分参数λ,ϕ,γλ,ϕ,γ的表达式。其他6项的表达式为：

Eq[logp(z|θ)]=∑n=1N∑k=1KϕnkΨ(γk)−Ψ(∑k′=1Kγk′)(17)(17)Eq[logp(z|θ)]=∑n=1N∑k=1KϕnkΨ(γk)−Ψ(∑k′=1Kγk′)

 

 

Eq[logp(θ|α)]=logΓ(∑k=1Kαk)−∑k=1KlogΓ(αk)+∑k=1K(αk−1)(Ψ(γk)−Ψ(∑k′=1Kγk′))(18)(18)Eq[logp(θ|α)]=logΓ(∑k=1Kαk)−∑k=1KlogΓ(αk)+∑k=1K(αk−1)(Ψ(γk)−Ψ(∑k′=1Kγk′))

 

 

Eq[logp(w|z,β)]=∑n=1N∑k=1K∑i=1Vϕnkwin(Ψ(λki)−Ψ(∑i′=1Vλki′))(19)(19)Eq[logp(w|z,β)]=∑n=1N∑k=1K∑i=1Vϕnkwni(Ψ(λki)−Ψ(∑i′=1Vλki′))

 

 

Eq[logq(β|λ)]=∑k=1K(logΓ(∑i=1Vλki)−∑i=1VlogΓ(λki))+∑k=1K∑i=1V(λki−1)(Ψ(λki)−Ψ(∑i′=1Vλki′))(20)(20)Eq[logq(β|λ)]=∑k=1K(logΓ(∑i=1Vλki)−∑i=1VlogΓ(λki))+∑k=1K∑i=1V(λki−1)(Ψ(λki)−Ψ(∑i′=1Vλki′))

 

 

Eq[logq(z|ϕ)]=∑n=1N∑k=1Kϕnklogϕnk(21)(21)Eq[logq(z|ϕ)]=∑n=1N∑k=1Kϕnklogϕnk

 

 

Eq[logq(θ|γ)]=logΓ(∑k=1Kγk)−∑k=1KlogΓ(γk)+∑k=1K(γk−1)(Ψ(γk)−Ψ(∑k′=1Kγk′))(22)(22)Eq[logq(θ|γ)]=logΓ(∑k=1Kγk)−∑k=1KlogΓ(γk)+∑k=1K(γk−1)(Ψ(γk)−Ψ(∑k′=1Kγk′))

 

　　　　有了ELBO的具体的关于变分参数λ,ϕ,γλ,ϕ,γ的表达式，我们就可以用EM算法来迭代更新变分参数和模型参数了。

4. EM算法之E步：获取最优变分参数

　　　　有了前面变分推断得到的ELBO函数为基础，我们就可以进行EM算法了。但是和EM算法不同的是这里的E步需要在包含期望的EBLO计算最佳的变分参数。如何求解最佳的变分参数呢？通过对ELBO函数对各个变分参数λ,ϕ,γλ,ϕ,γ分别求导并令偏导数为0，可以得到迭代表达式，多次迭代收敛后即为最佳变分参数。

　　　　这里就不详细推导了，直接给出各个变分参数的表达式如下：

 

ϕnk∝exp(∑i=1Vwin(Ψ(λki)−Ψ(∑i′=1Vλki′))+Ψ(γk)−Ψ(∑k′=1Kγk′))(23)(23)ϕnk∝exp(∑i=1Vwni(Ψ(λki)−Ψ(∑i′=1Vλki′))+Ψ(γk)−Ψ(∑k′=1Kγk′))

 

 　　　　其中，win=1wni=1当且仅当文档中第nn个词为词汇表中第ii个词。

 

γk=αk+∑n=1Nϕnk(24)(24)γk=αk+∑n=1Nϕnk

 

 

λki=ηi+∑n=1Nϕnkwin(25)(25)λki=ηi+∑n=1Nϕnkwni

 

　　　　由于变分参数λλ决定了ββ的分布，对于整个语料是共有的，因此我们有：

 

λki=ηi+∑d=1M∑n=1Ndϕdnkwidn(26)(26)λki=ηi+∑d=1M∑n=1Ndϕdnkwdni

 

　　　　最终我们的E步就是用（23）（24）（26）式来更新三个变分参数。当我们得到三个变分参数后，不断循环迭代更新，直到这三个变分参数收敛。当变分参数收敛后，下一步就是M步，固定变分参数，更新模型参数α,ηα,η了。

5. EM算法之M步：更新模型参数

　　　　由于我们在E步，已经得到了当前最佳变分参数，现在我们在M步就来固定变分参数，极大化ELBO得到最优的模型参数α,ηα,η。求解最优的模型参数α,ηα,η的方法有很多，梯度下降法，牛顿法都可以。LDA这里一般使用的是牛顿法，即通过求出ELBO对于α,ηα,η的一阶导数和二阶导数的表达式，然后迭代求解α,ηα,η在M步的最优解。

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

∇αkL=M(Ψ(∑k′=1Kαk′)−Ψ(αk))+∑d=1M(Ψ(γdk)−Ψ(∑k′=1Kγdk′))∇αkL=M(Ψ(∑k′=1Kαk′)−Ψ(αk))+∑d=1M(Ψ(γdk)−Ψ(∑k′=1Kγdk′))

 

 

∇αkαjL=M(Ψ′(∑k′=1Kαk′)−δ(k,j)Ψ′(αk))∇αkαjL=M(Ψ′(∑k′=1Kαk′)−δ(k,j)Ψ′(αk))

 

　　　　其中，当且仅当k=jk=j时，δ(k,j)=1δ(k,j)=1,否则δ(k,j)=0δ(k,j)=0。

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

∇ηiL=K(Ψ(∑i′=1Vηi′)−Ψ(ηi))+∑k=1K(Ψ(λki)−Ψ(∑i′=1Vλki′))∇ηiL=K(Ψ(∑i′=1Vηi′)−Ψ(ηi))+∑k=1K(Ψ(λki)−Ψ(∑i′=1Vλki′))

 

 

∇ηiηjL=K(Ψ′(∑i′=1Vηi′)−δ(i,j)Ψ′(ηi))∇ηiηjL=K(Ψ′(∑i′=1Vηi′)−δ(i,j)Ψ′(ηi))

 

　　　　其中，当且仅当i=ji=j时，δ(i,j)=1δ(i,j)=1,否则δ(i,j)=0δ(i,j)=0。

　　　　最终牛顿法法迭代公式为：

αk+1=αk+∇αkL∇αkαjL(27)(27)αk+1=αk+∇αkL∇αkαjL

 

 

ηi+1=ηi+∇ηiL∇ηiηjL(28)(28)ηi+1=ηi+∇ηiL∇ηiηjL

 

6. LDA变分推断EM算法流程总结

　　　　下面总结下LDA变分推断EM的算法的概要流程。

　　　　输入：主题数KK,M个文档与对应的词。

　　　　1） 初始化α,ηα,η向量。

　　　　2）开始EM算法迭代循环直到收敛。

　　　　　　a) 初始化所有的ϕ,γ,λϕ,γ,λ，进行LDA的E步迭代循环,直到λ,ϕ,γλ,ϕ,γ收敛。

　　　　　　　　(i) for d from 1 to M:

　　　　　　　　　  　　for n from 1 to NdNd:

　　　　　　　　　　　  　　for k from 1 to K:

　　　　　　　　　　　　　　　　按照(23)式更新ϕnkϕnk

　　　　　　　　　　　  标准化ϕnkϕnk使该向量各项的和为1.

　　　　　　　　　　按照(24) 式更新γkγk。

　　　　　　　　(ii) for k from 1 to K:

　　　　　　　　　　　　for i from 1 to V:

　　　　　　　　　　按照(26) 式更新λkiλki。

　　　　　　　　(iii)如果ϕ,γ,λϕ,γ,λ均已收敛，则跳出a)步，否则回到(i)步。

　　　　　　b) 进行LDA的M步迭代循环， 直到α,ηα,η收敛

　　　　　　　　(i) 按照(27)(28)式用牛顿法迭代更新α,ηα,η直到收敛

　　　　　　c) 如果所有的参数均收敛，则算法结束，否则回到第2)步。

 

　　　　　　算法结束后，我们可以得到模型的后验参数α,ηα,η,以及我们需要的近似模型主题词分布λλ,以及近似训练文档主题分布γγ。

转载 文本主题模型之LDA(三) LDA求解之变分推断EM算法 - 刘建平Pinard - 博客园  http://www.cnblogs.com/pinard/p/6873703.html