潜在狄利克雷分布LDA

潜在狄利克雷分布(latent Dirichlet allocation,LDA)，作为基于贝叶斯学习的话题模型，是潜在语义分析、概率潜在语义分析的扩展。

LDA模型是文本集合的生成概率模型。假设每个文本由话题的一个多项分布表示，每个话题由单词的一个多项分布表示，特别假设文本的话题分布的先验分布是狄利克雷分布，话题的单词分布的先验分布也是狄利克雷分布。

先验分布的导入使LDA能够更好的应对话题模型学习中的过拟合现象。

LDA的文本集合的生成过程如下：首先，随机生成一个文本的话题分布，之后在该文本的每个位置，依据该文本的话题分布随机生成一个话题，然后在该位置依据该话题的单词分布生成一个单词，直至文本的最后一个位置，生成整个文本。重复以上过程生成所有文本。

LDA模型是含有隐变量的概率图模型。每个话题的单词分布，每个文本的话题分布，文本的每个位置的话题是隐变量；文本的每个位置的单词是观测变量。

LDA模型的学习与推理无法直接求解，通常使用吉布斯抽样和变分EM算法(variational EM algorithm)，前者是蒙特卡罗法，而后者是近似算法。

狄利克雷分布

多项分布是一种多元离散随机变量的概率分布，是二项分布的扩展。

假设重复进行n次独立随机试验，每次试验可能出现的结果有k种，第i种结果出现的概率为pi，第i种结果出现的次数为ni。如果用随机变量X=(X1,X2,…,Xk)表示试验所有可能结果的次数，其中Xi表示第i种结果出现的次数，那么随机变量X服从多项分布。

定义20.1(多项分布) 若多元离散变量X=(X1,X2,…,Xk)的概率质量函数为:

P(X1=n1,X2=n2,…,Xk=nk)=n!n1!n2!…nk!pn11pn22…pnkk=n!∏ki=1ni!k∏i=1pnIi

其中p=(p1,p2,…,pk),pi≥0 ,i=1,2,…,k ,∑ki=1pi=1 ,∑ki=1ni=n，则称随机变量X服从参数为(n,p)的多项分布，记作X∼Mult(n,p)。

当试验的次数n为1时，多项分布变成类别分布(categorical distribution)。类别分布表示试验可能出现的k种结果的概率。

狄利克雷分布

狄利克雷分布是一种多元连续随机变量概率分布，是贝塔分布(beta distribution)的扩展。在贝叶斯学习中，狄利克雷分布常作为多项式分布的先验分布使用。

定义20.2(狄利克雷分布) 若多元连续随机变量θ(θ1,θ2,…,θk)的概率密度函数为：

p(θ|α)=Γ(∑ki=1αi)∏ki=1Γ(αi)k∏i=1θαi−1i

其中∑ki=1θi=1 ,θi≥0 ,α=(α1,α2,…,αk) ,αi>0 ,i=1,2,…,k，则称随机变量θ服从参数为α的狄利克雷分布，记作 θ∼Dir(α)。

式中Γ(s)是伽马函数，定义为：

Γ(s)=∫∞0xs−1e−xdx ,s>0 Γ(s+1)=sΓ(s)

当s是自然数时，有

Γ(s+1)=s!

由于满足条件

θi≥0 ,k∑i=1θi=1

所以狄利克雷分布θ存在于(k-1)维单纯形上。

B(α)=∏ki=1Γ(αi)Γ(∑ki=1αi)

则狄利克雷分布的密度函数可以写成：

p(θ|α)=1B(α)k∏i=1θαi−1i

其中，B(α)是规范化因子，称为多元贝塔函数(或扩展的贝塔函数)。由密度函数的性质：

∫Γ(∑ki=1αi)∏ki=1Γ(αi)k∏i=1θαi−1idθ=Γ(∑ki=1αi)∏ki=1Γ(αi)∫k∏i=1θαi−1idθ=1 B(α)=∫k∏i=1θαi−1idθ

所以式(20-5)是多元贝塔函数的积分表示。

二项分布和贝塔分布

二项分布是多项分布的特殊情况，贝塔分布是狄利克雷分布的特殊情况。

二项分布是指如下概率分布。X为离散随机变量，取值为m，其概率质量函数为:

P(X=m)=(nm)pm(1−p)n−m ,m=0,1,2,…,n

其中n和p(0≤p≤1)是参数。

贝塔分布是指如下概率分布，X为连续随机变量，取值范围为[0,1]，其概率密度函数为:

p(x)={1B(s,t)xs−1(1−x)t−1,0≤x≤10,其它

其中，s>0和t>0是参数，B(s,t)=Γ(s)Γ(t)Γ(s+t)是贝塔函数，定义为:

B(s,t)=∫10xs−1(1−x)t−1dx

当s,t是自然数时，

B(s,t)=(s−1)!(t−1)!(s+t−1)!

当n=1时，二项分布变成伯努利分布或0-1分布。伯努利分布表示试验可能出现的2种结果的概率。

下图给出了几种概率分布的关系：

狄利克雷分布的一些重要性质：

1. 狄利克雷分布属于指数分布族；
2. 狄利克雷分布是多项分布的共轭先验

贝叶斯学习中常使用共轭分布。如果后验分布与先验分布属于同类，则先验分布与后验分布称为共轭分布(conjugate distribution)，先验分布称为共轭先验(conjugate prior)。

作为先验分布的狄利克雷分布的参数又称为超参数。使用共轭分布的好处是便于从先验分布计算后验分布。

设W={w1,w2,…,wk}是由k个元素组成的集合。随机变量X服从W上的多项分布，X∼Mult(n,θ)，其中n=(n1,n2,…,nk)和θ=(θ1,θ2,…,θk)是参数。参数n为从W中重复独立抽取样本的次数，ni为样本中wi出现的次数(i=1,2,…,k)；参数θi为wi出现的概率(i=1,2,…,k)。

将样本数据表示为D，目标是计算在样本数据D给定条件下参数θ的后验概率p(θ|D)。对于给定的样本数据D，似然函数是:

p(D|θ)=θn11θn22…θnkk=k∏i=1θnii

假设随机变量θ服从狄利克雷分布p(θ|α)，其中α=(α1,α2,…,αk)为参数。则θ的先验分布为：

p(θ|α)=Γ(∑ki=1αi)∏ki=1Γ(αi)k∏i=1θαi−1i=1B(α)k∏i=1θαi−1i=Dir(θ|α) ,αi>0

根据贝叶斯规则，在给定样本数据D和参数α条件下，θ 的后验概率分布是：

p(θ|D,α)=p(D|θ)p(θ|α)p(D|α)=∏ki=1θnii1B(α)∏ki=1θαi−1i∫∏ki=1θnii1B(α)∏ki=1θαi−1idθ=1B(α+n)k∏i=1θαi+ni−1i=Dir(θ|α+n)

狄利克雷后验分布的参数等于狄利克雷先验分布参数α=(α1,α2,…,αk)加上多项分布的观测计数n=(n1,n2,…,nk)。好像试验之前就已经观察到计数α=(α1,α2,…,αk)，因此也把α叫做先验伪计数(prior pseudo-counts)。

潜在狄利克雷分配模型

潜在狄利克雷分配是文本集合的生成概率模型。

LDA模型表示文本集合的自动生成过程：

1. 首先，基于单词分布的先验分布(狄利克雷分布)生成多个单词分布，即决定多个主题内容；
2. 之后，基于话题分布的先验分布(狄利克雷分布)生成多个话题分布，即决定多个文本内容；
3. 然后，基于每一话题分布生成话题序列，针对每个话题，基于话题的单词分布生成单词，整体构成一个单词序列，即生成文本。重复这个过程生成所有文本。

下图显示LDA的文本生成过程：

LDA模型是概率图模型，其特点是以狄利克雷分布为多项分布的先验分布，学习就是给定文本集合，通过后验概率分布的估计，推断模型的所有参数。利用LDA进行话题分析，就是对给定文本集合，学习到每个文本的话题分布，以及每个话题的单词分布。

LDA与PLSA的区别：

不同点：

1. LDA使用狄利克雷分布作为先验分布，而PLSA不使用先验分布(或者说假设先验分布是均匀分布)，两者对文本生成过程有不同假设。
2. 学习过程LDA基于贝叶斯学习，而PLSA基于极大似然估计

相同点：两者都假设话题是单词的多项分布，文本是话题的多项分布。

LDA的优点是，使用先验概率分布，可以防止学习过程中产生的过拟合。

潜在狄利克雷分配使用三种集合:

1. 单词集合W={w1,…,wv,…,wV}，其中wv是第v个单词，v=1,2,…,V，V是单词的个数。
2. 文本集合D={w1,…,wm,…,wM}，其中wm是第m个文本,m=1,2,…,M，M是文本的个数。文本wm是一个单词序列wm=(wm1,…,wmn,…,wmNm)，其中wmn是文本wm的第n个单词，n=1,2,…,Nm，Nm是文本wm中单词的个数。
3. 话题集合Z={z1,…,zk,…,zK}，其中zk是第k个话题，k=1,2,…,K，K是话题的个数。

每个话题zk由一个单词的条件概率分布p(w|zk)决定，w∈W。分布p(w|zk)服从多项分布(严格意义上类别分布)，其参数为φk。

参数φk服从狄利克雷分布(先验分布)，其超参数为β。参数φk是一个V维向量φk=(φk1,φk2,…,φkV)，其中φkv表示话题zk生成单词wv的概率。

所有话题的参数向量构成一个KxV矩阵φ={φk}Kk=1。超参数β也是一个V维向量β=(β1,β2,…,βV)。

每个文本wm由一个话题的条件概率分布p(z|wm)决定的，z∈Z。分布p(z|wm)服从多项分布(严格意义上类别分布)，其参数为θm。

参数θm服从狄利克雷分布(先验分布)，其超参数为α。参数θm是一个K维向量θm=(θm1,θm2,…,θmK)，其中θmk表示文本wm生成话题zk的概率。

所有文本的参数向量构成一个MxK矩阵θ={θm}Mm=1。超参数α也是一个K维向量α=(α1,α2,…,αK)。

每个文本wm中的每个单词wmn由该文本的话题分布p(z|wm)以及所有话题的单词分布p(w|zk)决定。

LDA文本集合的生成过程如下：

给定单词集合W，文本集合D，话题集合Z，狄利克雷分布的超参数α和β。

1. 生成话题的单词分布

   随机生成K个话题的单词分布。具体过程如下：

   按照狄利克雷分布Dir(β)随机生成一个参数向量φk，φ∼Dir(β)，作为话题zk的单词分布p(w|zk),w∈W ，k=1,2,…,K

2. 生成文本的话题分布

   随机生成M个文本的话题分布。具体过程如下：

   按照狄利克雷分布Dir(α)随机生成一个参数向量θm，θm∼Dir(α)，作为文本wm的话题分布p(z|wm)，m=1,2,…,M。

3. 生成文本的单词序列

   随机生成M个文本的Nm个单词。文本wm(m=1,2,…,M)的单词wmn(n=1,2,…,Nm)的生成过程如下：

   1. 首先，按照多项分布Mult(θm)随机生成一个话题zmn ,zmn∼Mult(θm)。
   2. 然后，按照多项分布Mult(φzmn)随机生成一个单词wmn ,wmn∼Mult(φzmn)。

   文本wm本身是单词序列wm=(wm1,…,wmn,…,wmNm)，对应着隐式的话题序列:zm=(zm1,zm2,…,zmNm)。

算法20.1(LDA的文本生成算法)

1. 对于话题zk(k=1,2,…,K)：

   生成多项分布参数φk∼Dir(β)，作为话题的单词分布p(w|zk)

2. 对于文本wm(m=1,2,…,M)：

   生成多项分布参数θm∼Dir(α)，作为文本的话题分布p(z|wm)。

3. 对于文本wm的单词wmn(m=1,2,…,M,n=1,2,…,Nm)

   1. 生成话题zmn∼Mult(θm)，作为单词对应的话题。
   2. 生成单词wmn∼Mult(φzmn)

LDA的文本生成过程中，假定话题个数K给定(通常通过试验给定)，狄利克雷分布的超参数α和β通常也是事先给定的(在没有其它先验知识情况下,可以假设向量α和β的所有分量均为1)。这时的文本的话题分布θm是对称的，话题的单词分布φk也是对称的。

概率图模型

如下图所示，结点表示随机变量，双边圆是观测变量，单边圆是隐变量；有向边表示概率依存关系；矩形表示重复，矩形内的数字表示重复次数。

上图的LDA板块表示，结点α和β是模型的超参数，结点φk表示话题zk的单词分布的参数，结点θm表示文本wm的话题分布的参数，结点zmn表示话题，结点wmn表示单词。

结点β指向结点φk，重复K次，表示根据超参数β生成K个话题的单词分布的参数φk;

结点α指向结点θm，重复M次，表示根据超参数α生成M个文本的话题分布的参数θm；

结点θm指向结点zmn，重复Nm次，表示根据文本的话题分布θm生成Nm个话题zmn；

结点zmn指向结点wmn，同时K个结点φk也指向结点wmn，表示根据话题zmn以及K个话题的单词分布φk生成单词wmn。

把板块图展开得到以下的有向图表示形式：

随机变量序列的可交换性

一个有限的随机变量序列是可交换的，是指随机变量的联合概率分布对随机变量的排列不变：

P(x1,x2,…,xN)=P(xπ(1),xπ(2),…,xπ(N))

这里π(1),π(2),…,π(N)表示自然数1,2,…,N的任意一个排列。一个无限的随机变量序列是无限可交换的，是指它的任意一个有限子序列都是可交换的。

如果一个随机变量序列X1,X2,…,XN,…是独立同分布的，那么它们是无限可交换的。反之不然。

根据De Finetti定理，任意一个无限可交换的随机变量序列对一个随机参数是条件独立同分布的。即任意一个无限可交换的随机变量序列X1,X2,…,Xi,…的i基于一个随机参数Y的条件概率，等于基于这个随机参数Y的各个随机变量X1,X2,…,Xi,…的条件概率的乘积：

P(X1,X2,…,Xi,…|Y)=P(X1|Y)P(X2|Y)…P(Xi|Y)…

LDA的假设前提为：文本中的话题对一个随机参数是条件独立同分布的。因此，在参数给定的情况下，文本中的话题的顺序可以忽略。

LDA模型整体是由观测变量和隐变量组成的联合概率分布，可以表示为：

p(w,z,θ,φ|α,β)=K∏k=1p(φk|β)M∏m=1p(θm|α)Nm∏n=1p(zmn|θm)p(wmn|zmn,φ)

其中观测变量w表示所有文本中的单词序列，隐变量z表示所有文本中的话题序列，隐变量θ表示所有文本的话题分布的参数，隐变量φ表示所有话题的单词分布的参数，α和β是超参数。

• p(φk|β)表示超参数β给定的情况下，第k个话题的单词分布的参数φk的生成概率；
• p(θm|α)表示超参数α给定的情况下，第m个文本的话题分布的参数θm的生成概率；
• p(zmn|θm)表示第m个文本的话题分布θm给定的情况下，文本的第n个位置的话题zmn的生成概率；
• p(wmn|zmn,φ)，表示在第m个文本的第n个位置的话题zmn及所有话题的单词分布的参数φ给定的情况下，第m个文本的第n个位置的单词wmn的生成概率。

第m个文本的联合概率分布可以表示为：

p(wm,zm,θm,φ|α,β)=K∏k=1p(φk|β)p(θm|α)Nm∏n=1p(zmn|θm)p(wmn|zmn,φ)

其中，wm表示该文本中的单词序列，zm表示该文本的话题序列，θm表示该文本的话题分布参数。

LDA模型的联合分布含有隐变量，对隐变量进行积分得到边缘分布。

参数θm,φ给定的情况下，第m个文本的生成概率是：

p(wm|θm,φ)=Nm∏n=1[K∑k=1p(zmn=k|θm)p(wmn|φk)]

超参数α,β给定的条件下，第m个文本的生成概率是：

p(wm|α,β)=K∏k=1∫p(φk|β){∫p(θm|α)Nm∏n=1[K∑l=1p(zmn=l|θm)p(wmn|φl)]dθm}dφk

超参数α,β给定的条件下，所有文本的生成概率是

p(w|α,β)=K∏k=1∫p(φk|β){M∏m=1∫p(θm|α)Nm∏n=1[K∑l=1p(zmn=l|θm)p(wmn|φl)]dθm}dφk

LDA的吉布斯抽样算法

潜在狄利克雷分布的学习(参数估计)是一个复杂的最优化问题，很难精确求解，只能近似求解。

常用的近似求解方法有吉布斯抽样和变分推理。

吉布斯抽样优点是实现简单，缺点是收敛速度慢。

给定文本(单词序列)的集合D={w1,…,wm,…,wM}，其中wm是第m个文本,m=1,2,…,M，M是文本的个数。文本wm是一个单词序列wm=(wm1,…,wmn,…,wmNm)，其中wmn是文本wm的第n个单词，n=1,2,…,Nm，Nm是文本wm中单词的个数。

超参数α,β已知。目标是要推断：

1. 话题序列的集合z={z1,…,zm,…,zM}的后验概率分布，其中zm是第m个文本的话题序列，zm=(zm1,zm2,…,zmNm)；
2. 参数θ={θ1,…,θm,…,θM}，其中θm是文本wm的话题分布的参数；
3. 参数φ={φ1,…,φk,…,φK}，其中φk是话题zk的单词分布的参数。

也就是说，要对联合概率分布p(w,z,θ,φ|α,β)进行估计，其中w是观测变量，而z,θ,φ是隐变量。

使用吉布斯抽样对多元随机变量x的联合分布进行估计：

为了估计多元随机变量x的联合分布p(x)，吉布斯抽样法选择x的一个分量，固定其它分量，按照其条件概率分布进行随机抽样，依次循环对每一分量执行这个操作，得到联合分布p(x)的一个随机样本，重复这个过程，在燃烧期之后，得到联合概率分p(x)的样本集合。

LDA模型的学习通常采用收缩的吉布斯抽样(collapsed Gibbs sampling)方法。基本想法是，

1. 通过对隐变量θ和φ积分，得到边缘概率分布p(w,z|α,β)，其中w是可观测变量，变量z是e不可观测的；
2. 对后验概率分布p(z|w,α,β)进行吉布斯抽样，得到分布p(z|w,α,β)的样本集合；
3. 利用这个样本集合对参数θ和φ进行估计，最终得到LDA模型p(w,z,θ,φ|α,β)的所有参数估计。

算法 的主要部分

根据前面的分析，问题转化为对后验概率分布p(z|w,α,β)的吉布斯抽样，该分布表示在所有文本的单词序列给定的情况下，所有可能话题序列的条件概率。

抽样分布的表达式

p(z|w,α,β)=p(w,z|α,β)p(w|α,β)∝p(w,z|α,β)

这里变量w,α,β已知。联合分布p(w,z|α,β)的表达式可进一步分解为：

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

两个因子可以分开处理。

推导第一个因子p(w|z,β)的表达式。首先：

p(w|z,φ)=K∏k=1V∏v=1φnkvkv

其中，φkv是第k个话题生成单词集合第v个单词的概率，nkv是数据中第k个话题生成第v个单词的次数。于是

p(w|z,β)=∫p(w|z,φ)p(φ|β)dφ=∫K∏k=11B(β)V∏v=1φnkv+βv−1kvdφ=K∏k=1B(nk+β)B(β)

其中，nk={nk1,nk2,…,nkV}

第二个因子p(z|α)的表达式可以类似推导。首先：

p(z|θ)=M∏m=1K∏k=1θnmkmk

其中，θmk是第m个文本生成第k个话题的概率，nmk是数据中第m个文本生成第k个话题的次数。于是：

p(z|α)=∫p(z|θ)p(θ|α)dθ=∫M∏m=11B(α)K∏k=1θnmk+αk−1mkdθ=M∏m=11B(α)∫K∏k=1θnmk+αk−1mkdθ=M∏m=1B(nm+α)B(α)

其中nm={nm1,nm2,…,nmK}。由式(20-23)和式(20-25)得

p(w,z|α,β)=K∏k=1B(nk+β)B(β)∙M∏m=1B(nm+α)B(α)

故由式(20-20)和式(20-26)，得收缩的吉布斯抽样分布的公式：

p(z|w,α,β)∝K∏k=1B(nk+β)B(β)∙M∏m=1B(nm+α)B(α)

满条件分布 的表达式

分布p(z|w,α,β)的满条件分布可以写成：

p(zi|z−i,w,α,β)=1Zzip(z|w,α,β)

这里wi表示所有文本的单词序列的第i个位置的单词，zi表示单词wi对应的话题，i=(m,n)，i=1,2,…,I,z−i={zj:j≠i}，Zzi表示分布p(z|w,α,β)对变量zi的边缘化因子。

式(20-28)是在所有文本单词序列、其它位置话题序列给定条件下，第i个位置的话题的条件概率分布。由式(20-27)和式(20-28)可以推出：

p(zi|z−i,w,α,β)∝nkv+βv∑Vv=1(nkv+βv)∙nmk+αk∑Kk=1(nmk+αk)

其中第m个文本的第n个位置的单词wi是单词集合的第v个单词，其话题zi是话题 集合的第k个话题，nkv表示第k个话题中第v个单词的计数，但减去当前单词的计数，nmk表示第m个文本中第k个话题的计数，但减去当前单词的话题的计数。

算法的后处理

通过吉布斯抽样得到的分布p(z|w,α,β)的a样本，可以得到变量z的分配值，也可以估计变量θ和φ。

• 参数θ={θm}

  根据LDA模型的定义，后验概率满足：

  p(θm|zm,α)=1ZθmNm∏n=1p(zmn|θm)p(θm|α)=Dir(θm|nm+α)

这里nm={nm1,nm2,…,nmK}是第m个文本的话题计数，Zθm表示 分布p(θm，zm|α)对变量θm的边缘化因子。于是得到参数θ={θm}的估计式：

θmk=nmk+αk∑Kk=1(nmk+αk) ,m=1,2,…,M;k=1,2,…,K

• 参数φ={φk}

  后验概率满足

  p(φk|w,z,β)=1ZφkI∏i=1p(wi|φk)p(φk|β)=Dir(φk|nk+β)

这里nk={nk1,nk2,…,nkV}是第k个话题的单词计数,Zφk表示分布p(φk,w|z,β)对变量φk的边缘化因子，I是文本集合单词序列w的单词总数。于是，得到参数的估计式：

φkv=nkv+βv∑Vv=1(nkv+βv) ,k=1,2,…,K;v=1,2,…,V

算法20.2(LDA吉布斯抽样算法)

输入：文本的单词序列w={w1,…,wm,…,wM}，wm=(wm1,…,wmn,…,wmNm)；

输出：文本的话题序列z={z1,…,zm,…,zM}，zm=(zm1,zm2,…,zmNm)的后验概率分布p(z|w,α,β)的样本计数，模型的参数φ和θ的估计值

参数：超参数α,β 和话题个数K

1. 设所有计数矩阵的元素nmk,nkv，计数向量的元素nm,nk初值为0；

2. 对所有文本wm，m=1,2,…,M

   对第m个文本中的所有单词wmn，n=1,2,…,Nm

   1. 抽样话题zmn=zk∼Mult(1K)

      增加文本-话题计数nmk=nmk+1,

      增加文本-话题和计数nm=nm+1，

      增加话题-单词计数nkv=nkv+1,

      增加话题-单词和计数nk=nk+1

3. 循环执行以下操作，直到进入燃烧期

   对所有文本wm,m=1,2,…,M

   对第m个文本中的所有单词wmn，n=1,2,…,Nm

   1. 当前的单词wmn是第v个单词，话题指派zmn是第k个话题

      减少计数nmk=nmk−1,nm=nm−1,nkv=nkv−1,nk=nk−1

   2. 按照满条件分布进行抽样

      p(zi|z−i,w,α,β)∝nkv+βv∑Vv=1(nkv+βv)∙nmk+αk∑Kk=1(nmk+αk)

      得到新的第k′个话题，分配给zmn

   3. 增加计数nmk′=nmk′+1,nm=nm+1,nk′v=nk′v+1,nk=nk+1

   4. 得到更新的两个计数矩阵NK×V=[nkv]和NM×K=[nmk]，表示后验概率分布p(z|w,α,β)的样本计数

4. 利用得到的样本计数，计算模型参数：

   θmk=nmk+αk∑Kk=1(nmk+αk)φkv=nkv+βv∑Vv=1(nkv+βv)

LDA的变分EM算法

变分推理(variational inference)是贝叶斯学习中常用的、含有隐变量模型的学习和推理方法。

MCMC通过随机抽样的方法近似的计算模型的后验概率；

变分推理则通过解析的方法计算模型的后验概率的近似值。

变分推理的基本想法：假设模型是联合概率分布p(x,z)，其中x是观测变量，z是隐变量，包括参数。目标是学习模型的后验概率分布p(z|x)，用模型进行概率推理。由于目标分布很复杂，因此通过使用概率分布q(z)近似条件概率分布p(z|x)，用KL散度D(q(z)‖p(z|x))计算两者的相似度，q(z)称为变分分布(variational distribution)。

如果能找到与p(z|x)在KL散度意义下最近的分布q∗(z)，则可以用这个分布近似p(z|x)。

p(z|x)≈q∗(z)

下图给出了q∗(z)与p(z|x)的关系：

KL散度可以写成以下形式：

D(q(z)‖p(z|x))=Eq[logq(z)]−Eq[logp(z|x)]=Eq[logq(z)]−Eq[logp(x,z)]+logp(x)=logp(x)−{Eq[logp(x,z)]−Eq[logq(z)]}

KL散度大或等于0，当且仅当两个分布一致时为0，由式(20-35)可推出：

logp(x)≥Eq[logp(x,z)]−Eq[logq(z)]

不等式右端是左端的下界，左端称作证据(evidence)，右端称为证据下界(evidence lower bound,ELBO)，证据下界记作：

L(q)=Eq[logp(x,z)]−Eq[logq(z)]

KL散度(20-35)的最小化可以通过证据下界(20-37)最大化实现，因为目标是求q(z)使KL散度最小化，这时logp(x)是常量。

因此，变分推理变成求解证据下界最大化的问题。

假设q(z)对z的所有分量都是相互独立的(实际是条件独立于参数)，即满足：

q(z)=q(z1)q(z2)…q(zn)

这时，变分分布称为平均场(mean field)。

KL散度的最小化或证据下界最大化实际是在平均场的集合，即满足独立假设的分布集合Q={q(z)|q(z)=∏ni=1q(zi)}之中进行的。

变分推理有以下几个步骤：

1. 定义变分分布q(z)
2. 推导其证据下界表达式
3. 用最优化方法对证据下界进行优化，如坐标上升，得到最有分布q∗(z)，作为后验概率分布p(z|x)的近似

变分EM算法

变分推理中，可以通过迭代的方法最大化证据下界，此算法是EM算法的推广，称为变分EM算法。

假设模型是联合概率分布p(x,z|θ)，其中x是观测变量，z是隐变量，θ是参数。目标是通过观测数据的概率(证据)logp(x|θ) 的最大化，估计模型的参数θ。

使用变分推理，导入平均场q(z)=∏ni=1q(zi)，定义证据下界：

L(q,θ)=Eq[logp(x,z|θ)]−Eq[logq(z)]

通过迭代，分别以q和θ为变量，对证据下界进行最大化。

算法20.3(变分EM算法)

循环执行以下E步和M步，直到收敛。

1. E步：固定θ，求L(q,θ)对q的最大化
2. M步：固定q，求L(q,θ)对θ的最大化

给出模型参数θ的估计值。

根据变分推理原理，观测数据的概率和证据下界满足：

logp(x|θ)−L(q,θ)=D(q(z)‖p(z|x,θ))≥0

变分EM算法的迭代过程中，以下关系成立：

logp(x|θ(t−1))=L(q(t),θ(t−1))≤L(q(t),θ(t))≤logp(x|θ(t))

其中上角标t-1和t表示迭代次数。

左边的等式基于E步计算和变分推理原理，中间的不等式基于M步计算，右边的不等式基于变分推理原理。说明了变分EM算法的收敛性。

对照《统计学习方法》9.4节，EM算法的推广，EM算法的推广是求F函数的极大-极大算法，其中的F函数就是证据下界。

EM算法假设q(z) = p(z|x)且p(z|x)容易计算；而变分EM算法则考虑一般情况，使用容易计算的平均场q(z)=∏ni=1q(zi)。

当模型复杂时，EM算法未必可用，但变分EM算法仍然可以使用。

将变分EM算法应用到下图的LDA模型的学习上。

首先定义具体的变分分布，推导证据下界的表达式，接着，推导变分分布的参数和LDA模型的参数的估计式，最后给出LDA模型的变分EM算法。

证据下界的定义

一次只考虑一个文本。记作w

文本的单词序列w=(w1,…,wn,…,wN)，对应的话题序列z=(z1,…,zn,…,zN)，以及话题分布θ，随机变量w,z和θ的联合分布是：

p(θ,z,w|α,φ)=p(θ|α)N∏n=1p(zn|θ)p(wn|zn,φ)

其中，w是观测变量，θ和z是隐变量，α,φ是参数。

定义基于平均场的变分分布：

q(θ,z|γ,η)=q(θ|γ)N∏n=1q(zn|ηn)

其中γ是狄利克雷分布参数，η=(η1,η2,…,ηn)是多项分布的参数，变量θ和z的各个分量都是条件独立的。目标是求KL散度意义下最相近的变分分布q(θ,z|γ,η)，以近似LDA模型的后验分布p(θ,z|w,α,φ)。

下图为变分分布的板块表示：

LDA模型中的隐变量θ,z之间存在依存关系，变分分布中这些依存关系被去掉，变量θ,z条件独立。由此得到一个文本的证据下界：

L(γ,η,α,φ)=Eq[logp(θ,z,w|α,φ)]−Eq[logq(θ,z|γ,η)]

其中，数学期望是对分布q(θ,z|γ,η)定义的，为了方便写作Eq[∙]，γ,η是变分分布的参数，α,φ是LDA模型参数。

所有文本的证据下界为：

Lw(γ,η,α,φ)=M∑m=1{Eqm[logp(θm,zm,wm|α,φ)]−Eqm[logq(θm,zm|γm,ηm)]}

为求解证据下界L(γ,η,α,φ)的最大化，首先写出证据下界的表达式。为此展开证据下界(20-44)：

L(γ,η,α,φ)=Eq[logp(θ|α)]+Eq[logp(z|θ)]+Eq[logp(w|z,φ)]−Eq[logq(θ|γ)]−Eq[logq[z|η]]

根据变分参数γ和η，模型参数α和φ继续展开，并将展开式的每一项写成一行：

L(γ,η,α,φ)=logΓ(K∑l=1αl)−K∑k=1logΓ(αk)+K∑k=1(αk−1)[Ψ(γk)−Ψ(K∑l=1γl)]+N∑n=1K∑k=1ηnk[Ψ(γk)−Ψ(K∑l=1γl)]+N∑n=1K∑k=1V∑v=1ηnkwvnlogφkv−logΓ(K∑l=1γl)+K∑k=1logΓ(γk)−K∑k=1(γk−1)[Ψ(γk)−Ψ(K∑l=1γl)]−N∑n=1K∑k=1ηnklogηnk

式中Ψ(αk)是对数伽马函数的导数，即：

Ψ(αk)=ddαklogΓ(αk)

第一项推导，求Eq[logp(θ|α)]，是关于分布q(θ,z|γ,η)的数学期望。

Eq[logp(θ|α)]=K∑k=1(αk−1)Eq[logθk]+logΓ(K∑l=1αl)−K∑k=1logΓ(αk)

其中θ∼Dir(θ|γ)，根据这里的狄利克雷分布的相关性质，可得：

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

式中αk,γk表示第k个话题的狄利克雷分布参数。

第二项推导，求Eq[logp(z|θ)]，是关于分布q(θ,z|γ,η)的数学期望。

Eq[logp(z|θ)]=N∑n=1Eq[logp(zn|θ)]=N∑n=1Eq(θ,zn|γ,η)[log(zn|θ)]=N∑n=1K∑k=1q(znk|η)Eq(θ|γ)[logθk]=N∑n=1K∑k=1ηnk[Ψ(γk)−Ψ(K∑l=1γl)]

式中ηnk表示文档第n个位置的单词由第k个话题产生的概率，γk表示第k个话题的狄利克雷分布参数。最后一步请参考这里的狄利克雷分布的相关性质。

第三项推导，求Eq[logp(w|z,φ)]，是关于分布q(θ,z|γ,η)的数学期望。

Eq[logp(w|z,φ)]=N∑n=1Eq[logp(wn|zn,φ)]=N∑n=1Eq(zn|η)[logp(wn|zn,φ)]=N∑n=1K∑k=1q(znk|η)logp(wn|zn,φ)=N∑n=1K∑k=1V∑v=1ηnkwvnlogφkv

式中ηnk表示文档第n个位置的单词由第k个话题产生的概率，wvn在第n个位置的单词是单词集合的第v个单词时为1，否则取值为0，φkv表示第k个话题生成单词集合中第v个单词的概率。

第四项推导，求Eq[logq(θ|γ)]，是关于分布q(θ,z|γ,η)的数学期望。

由于θ∼Dir(γ)，类似于(20-50)，可以得到：

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

式中γk表示第k个话题的狄利克雷分布参数。

第五项推导，求Eq[logq[z|η]]，是关于分布q(θ,z|γ,η)的数学期望。

Eq[logq[z|η]]=N∑n=1Eq[logq[zn|η]]=N∑n=1Eq(zn|η)[logq[zn|η]]=N∑n=1K∑k=1Eq(znk|η)[logq[znk|η]]=N∑n=1K∑k=1ηnklogηnk

式中ηnk表示文档第n个位置的单词由第k个话题产生的概率，γk表示第k个话题的狄利克雷分布参数。

变分参数估计

• 首先，通过证据下界最优化估计参数η。ηnk表示第n个位置的单词是由第k个话题生成的概率。

  考虑式(20-47)关于ηnk的最大化，ηnk满足约束条件∑Kl=1ηnl=1。使用拉格朗日乘子法，得：

  L[ηnk]=ηnk[Ψ(γk)−Ψ(K∑k=1γl)]+ηnklogφkv−ηnklogηnk+λn(K∑l=1ηnl−1)

这里φkv是(在第n个位置)由第k个话题生成第v个单词的概率。

对ηnk求偏导数得:

∂L∂ηnk=Ψ(γk)−Ψ(K∑l=1γl)+logφkv−logηnk−1+λn

令偏导数为0，得到参数ηnk的估计值：

ηnk∝φkvexp[φ(γk)−φ(K∑l=1γl)]

• 接着，通过证据下界最优化估计参数γ。γk是第k个话题的狄利克雷分布参数。

  考虑式(20-47)关于γk得最大化：

  L[γk]=K∑k=1(αk−1)[Ψ(γk)−Ψ(K∑l=1γl)]+N∑n=1K∑k=1ηnk[Ψ(γk)−Ψ(K∑l=1γl)]−logΓ(K∑l=1γl)+logΓ(γk)−K∑k=1(γk−1)[Ψ(γk)−Ψ(K∑l=1γl)]

L[γk]=K∑k=1[Ψ(γk)−Ψ(K∑l=1γl)](αk+N∑n=1ηnk−γk)−logΓ(K∑l=1γl)+logΓ(γk)

对γk求偏导数得：

∂L∂γk=[Ψ′(γk)−Ψ′(K∑l=1γl)](αk+N∑n=1ηnk−γk)

令偏导数为0，求解得到参数γk的估计值：

γk=αk+N∑n=1ηnk

据此 ，得到由坐标上升法估计变分参数的方法，具体算法如下：

算法20.4(LDA的变分参数估计算法)

1. 初始化：对所有k和n，η(0)nk=1/K

2. 初始化：对所有k,γk=αk+N/K

3. 1. 对n=1到N

      1. 对k=1到K

         η(t+1)nk=φkvexp[φ(γ(t)k)−φ(K∑l=1γ(t)l)]

      2. 规范化η(t+1)nk使其和为1

   2. γ(t+1)=α+N∑n=1η(t+1)n

模型参数的估计

• 给定一个文本集合D={w1,…,wm,…,wM}，模型参数估计对所有文本同时进行。

  首先，通过证据下界的最大化估计φ。φkv表示第gek话题生成单词集合第v个单词的概率。

  将式(20-47)扩展到所有文本，并考虑关于φ的最大化。满足K个约束条件:

  V∑v=1φkv=1,k=1,2,…,K

从而得到对应的拉格朗日函数：

L[β]=M∑m=1Nm∑n=1K∑k=1V∑v=1ηmnkwvmnlogφkv+K∑k=1λk(V∑v=1φkv−1)

求φkv求偏导数并令其为0，归一化求解，得到参数φkv的估计值为：

φkv=m∑m=1Nm∑n=1ηmnkwvmn

其中，ηmnk为第m个文本的第n个单词属于第k个话题的概率，wvmn 在第m个文本的第n个单词是单词集合的第v个单词时取值为1，否则为0.

• 接着，通过证据下界的最大化估计参数α。αk表示第k个话题的狄利克雷分布参数。

  将式(20-47)扩展到所有文本，并考虑关于α的最大化：

  L[α]=M∑m=1{logΓ(K∑l=1αl)−K∑k=1logΓ(αk)+K∑k=1(αk−1)[Ψ(γmk)−Ψ(K∑l=1γml)]}

对αk求偏导数得：

∂L∂αk=M[Ψ(K∑l=1αl)−Ψ(αk)]+M∑m=1[Ψ(γmk)−Ψ(K∑l=1γml)]

再对αl求偏导数，得：

∂2L∂αk∂αl=M[Ψ′(K∑l=1αl)−δ(k,l)Ψ′(αk)]

这里δ(k,l)是delta函数。

式2(0-66)和式(20-67)分别是函数(20-65)对变量α的梯度g(α)和Hessian矩阵H(α)。

应用牛顿法求该函数的最大化。用以下公式迭代，得到参数α的估计值。

αnew=αold−H(αold)−1g(αold)

​ 据此，得到估计参数α的算法。

根据上面的推导，给出LDA的变分EM算法：

算法20.5(LDA的变分EM算法)

输入：给定文本集合D={w1,…,wm,…,wM}；

输出：变分参数γ,η，模型参数α,φ。

交替迭代E步和M步，直到收敛。

1. 固定模型参数α,φ，通过关于变分参数γ,η的证据下界的最大化，估计变分参数γ,η。具体见算法20.4

2. 固定变分参数γ,η，通过关于模型参数α,φ的证据下界的最大化，估计模型参数α,φ。具体见式(20-64)和式(20-68)。

   根据变分参数(γ,η)可以估计模型参数θ=(θ1,…,θm,…,θM)，z=(z1,…,zm,…,zM)。