主题模型之LDA

$\newcommand{\X}{\mathcal{X}} \newcommand{\a}{\vec{\alpha}} \newcommand{\b}{\vec{\beta}} \newcommand{\p}{\vec{p}} \newcommand{\w}{\vec{w}} \newcommand{\x}{\vec{x}} \newcommand{\z}{\vec{z}} \newcommand{\vt}{\vec{\vartheta}} \newcommand{\vp}{\vec{\varphi}}$

[TOC]

在原始的pLSA模型中,我们求解出两个参数:“主题-词项”矩阵 $\Phi$ 和“文档-主题”矩阵 $\Theta$ ,但是我们并未考虑参数的先验知识;而LDA的改进之处,是对这俩参数之上分别增加了先验分布,相应参数称作超参数(hyperparamter)。概率图表示如下:

lda pgm

其中,单圆圈表示隐变量;双圆圈表示观察到的变量;把节点用方框(plate)圈起来,表示其中的节点有多种选择。所以这种表示方法也叫做plate notation,具体可参考PRML 8.0 Graphical Models。

对应到上图,只有 $w_{m,n}$ 是观察到的变量,其他都是隐变量或者参数,其中 $\a$ 和 $\b$ 是超参数;方框中, 表示有 $K$ 种“主题-词项”分布; 有 $M$ 种“文档-主题”分布,即对每篇文档都会产生一个 $\vt_m$ 分布;每篇文档 $m$ 中有 $n$ 个词,每个词 $w_{m,n}$ 都有一个主题 $z_{m,n}$ ,该词实际是由 $\vp_{z_{m,n}}$ 产生。具体生成过程下面再说。

贝叶斯估计

在pLSA原参数之上增加先验分布,其实就是用贝叶斯估计取代最大似然估计,具体要了解各种参数估计方法可以参考Heinrich论文的第二部分。简单说,最大似然估计(MLE)和最大后验估计(MAP)都是把待估计的参数看作一个拥有固定值的变量,只是取值未知。通常估计的方法都是找使得相应的函数最大时的参数;由于MAP相比于MLE会考虑先验分布的影响,所以MAP也会有超参数,它的超参数代表的是一种信念(belief),会影响推断(inference)的结果。比如说抛硬币,如果我先假设是公平的硬币,这也是一种归纳偏置(bias),那么最终推断的结果会受我们预先假设的影响。贝叶斯估计是对MAP的扩展,但它不再对参数做直接的估计,而是把待估计的参数看作服从某种分布的随机变量。根据贝叶斯法则:

在MLE和MAP中,由于是要求函数最大值时的参数,所以都不会考虑evidence。但在贝叶斯估计中,不再直接取极值,所以还会考虑evidence,下面的这个积分也是通常贝叶斯估计中最难处理的部分:

evidence相当于对所有的似然概率积分或求和(离散时),所以也称作边界似然

共轭先验

由于有积分的存在,贝叶斯估计常常会很难推算,这里我们就需要利用一种共轭先验(Conjugate Prior)的数学知识。在贝叶斯统计理论中,如果某个随机变量 $\vartheta$ 的先验分布 $p(\vartheta)$ 和后验分布 $p(\vartheta\vert\X)$ 属于统一分布簇(也就是说有同样的函数形式),则称先验分布 $p(\vartheta)$ 和后验分布 $p(\vartheta\vert\X)$ 为共轭分布,先验分布 $p(\vartheta)$ 是似然函数 $p(\X\vert\vartheta)$ 的共轭先验。

由于pLSA中参数 $\Phi$ 和 $\Theta$ 都服从多项分布,因此选择Dirichlet分布作为他们的共轭先验。在pLSA中,假设这俩多项分布是确定的,我们在已经假设确定的分布下,选择具体的主题和词项;但在LDA中,这俩多项分布的参数是被当作随机变量,服从另一层先验分布,所以,Dirichlet分布也称作分布的分布,其定义如下:

其中, $\p$ 是要猜测的随机向量, $\a$ 是超参数, $\Delta(\a)$ 称作Delta函数,可以看作Beta函数的多项式扩展,是Dirichlet分布的归一化系数,定义如下:

相应的多项分布定义:

其中, $\p$ 和 $\vec{n}$ 服从约束 $\sum_k p_k=1$ 和 $\sum_k n_k = N$ 。

由于 $\Gamma(x+1) = x!$ 就是阶乘在实数集上的扩展,显然,公式(1)和(3)有相同的形式,所以,这俩分布也称作Dirichlet-Multinomail共轭。

如果 $\p \sim \mathrm{Dir}(\p\vert\a)$ ,则 $\p$ 中的任一元素 $p_i$ 的期望是:

可以看出,超参数 $\alpha_k$ 的直观意义就是事件先验的伪计数(prior pseudo-count)。

LDA生成模型

在LDA中,“文档-主题”向量 $\vt_m$ 由超参数为 $\a$ 的Dirichlet分布生成,“主题-词项”向量 $\vp_k$ 由超参数为 $\b$ 的Dirichlet分布生成,根据概率图,整个样本集合的生成过程如下:

  • 对所有的主题 $k \in [1,K]$ (生成 $\Phi$ , $K \times N$ 矩阵):
    • 生成“主题-词项”分布 $\vp_k \sim \mathrm{Dir}(\b)$ ( $N$ 维矩阵,对应词表 $\mathcal{V}$ 中的每个词项的概率)
  • 对所有的文档 $m \in [1,M]$ :
    • 生成当前文档 $m$ 相应的“文档-主题”分布 $\vt_m \sim \mathrm{Dir}(\a)$ ( $K$ 维向量,即第 $m$ 篇文档对应的每个主题的概率)
    • 生成当前文档 $m$ 的长度 $N_m \sim \mathrm{Poiss(\xi)}$
    • 对当前文档 $m$ 中的所有词 $n \in [1, N_m]$ :
      • 生成当前位置的词的所属主题 $z_{m,n} \sim \mathrm{Mult}(\vt_m)$
      • 根据之前生成的主题分布 $\Phi$ ,生成当前位置的词的相应词项 $w_{m,n} \sim \mathrm{Mult}(\vp_{z_{m,n}})$

由该生成过程可知,第 $m$ 篇文档中第 $n$ 个词 $t$ 的生成概率:

其中, $\Phi = \brace{\vp_k}_{k=1}^K$ 。

根据所有已知信息和带超参数的隐变量,我们可以写出联合分布:

通过对 $\vt_m$ 和 $\Phi$ 积分以及 $z_{m,n}$ 求和,可以求得 $\w_{m,n}$ 的分布:

整个样本的分布:

符号解释:

  • $M$ 文档数(固定值)
  • $K$ 主题(component)数(固定值)
  • $V$ 词项数(固定值)
  • $\a$ “文档-主题”分布的超参数( $K$ 维向量,如果对称(symmetric)参数,则是一个标量)
  • $\b$ “主题-词项”分布的超参数( $V$ 维向量,如果对称参数,则是一个标量)
  • $\vt_m$ 对于文档 $m$ 的主题分布参数标记 $p(z\vert d=m)$ ,每篇文档均不同, $\Theta=\brace{\vt_m}_{m=1}^M$ 是一个 $M \times K$ 矩阵
  • $\vp_m$ 对于主题 $k$ 的词项分布参数标记 $p(t\vert z=k)$ ,每个主题均不同, $\Phi=\brace{\vp_k}_{k=1}^K$ 是一个 $K \times V$ 矩阵
  • $N_m$ 文档 $m$ 的长度,这里由Poisson分布决定
  • $z_{m,n}$ 文档 $m$ 中第 $n$ 个词所属的主题
  • $w_{m,n}$ 文档 $m$ 中第 $n$ 个词的词项

Gibbs Sampling

Blei的原始论文使用变分法(Variational inference)和EM算法进行贝叶斯估计的近似推断,但不太好理解,并且EM算法可能推导出局部最优解。Heinrich使用了Gibbs抽样法,这也是目前LDA的主流算法。

通常均匀分布 $\mathrm{Uniform}(0,1)$ 的样本,即我们熟悉的类rand()函数,可以由线性同余发生器生成;而其他的随机分布都可以在均匀分布的基础上,通过某种函数变换得到,比如,正态分布可以通过Box-Muller变换得到。然而,这种变换依赖于计算目标分布的积分的反函数,当目标分布的形式很复杂,或者是高维分布时,很难简单变换得到。

当一个问题无法用分析的方法来求精确解,此时通常只能去推断该问题的近似解,而随机模拟(MCMC)就是求解近似解的一种强有力的方法。随机模拟的核心就是对一个分布进行抽样(Sampling)。随机模拟也可用于类pLSA算法,但现在很少有人这么做。

MCMC的基础:Markov链通过转移概率矩阵可以收敛到稳定的概率分布。这意味着MCMC可以借助Markov链的平稳分布特性模拟高维概率分布 $p(\x)$ ;当Markov链经过burn-in阶段,消除初始参数的影响,到达平稳状态后,每一次状态转移都可以生成待模拟分布的一个样本。Gibbs抽样是MCMC的一个特例,它交替的固定某一维度 $x_i$ ,然后通过其他维度 $\x_{\neg{i}}$ 的值来抽样该维度的值。它的基本算法如下:

  1. 选择一个维度 $i$ ,可以随机选择;
  2. 根据分布 $p(x_i\vert\x_{\neg{i}})$ 抽样 $x_i$ 。

所以,如果要完成Gibbs抽样,需要知道如下条件概率:

如果模型包含隐变量 $\z$ ,通常需要知道后验概率分布 $p(\z\vert\x)$ ,所以,包含隐变量的Gibbs抽样器公式如下:

LDA Gibbs Sampler

为了构造LDA Gibbs抽样器,我们需要使用隐变量的Gibbs抽样器公式。在LDA模型中,隐变量为 $z_{m,n}$ ,即样本中每个词 $w_{m,n}$ 所属的主题,而参数 $\Theta$ 和 $\Phi$ 等可以通过观察到的 $w_{m,n}$ 和相应的 $z_{m,n}$ 积分求得,这种处理方法称作collapsed,在Gibbs sampling中经常使用。

下面给出要推断的目标分布,它和联合分布成正比:

这里省略了超参数,这个分布涉及很多离散随机变量,并且分母是 $K_W$ 个项的求和,很难求解。此时,就需要Gibbs sampling发挥用场了,我们期望Gibbs抽样器可以通过Markov链利用全部的条件分布 $p(z_i\vert\z_{\neg{i}},\w)$ 来模拟 $p(\z\vert\w)$ 。根据公式(4),我们需要写出联合概率分布 $p(\w,\z)$ :

由于此公式第一部分独立于 $\a$ ,第二部分独立于 $\b$ ,所以可以分别处理。

第一部分,可以由观察到的词数以及相应主题的多项分布产生:

由于样本中的 $W$ 个词服从参数为主题 $z_i$ 的独立多项分布,这意味着,我们可以把上面的对词的乘积分解成对主题和对词项的两层乘积:

其中, $n_k^{(t)}$ 是词项 $t$ 在主题 $k$ 中出现的次数。目标分布 $p(\w\vert\z,\b)$ 需要对 $\Phi$ 积分,根据 $\Delta(\a)$ 函数公式(2)可得:

这个结果可以看作 $K$ 个Dirichlet-multinomial模型的乘积。

第二部分,类似于 $p(\w\vert\z,\b)$ 的步骤,先写出条件分布,然后分解成两部分的乘积:

其中, $d_i$ 是单词 $i$ 所属的文档, $n_m^{(k)}$ 是主题 $k$ 在文章 $m$ 中出现的次数。对 $\Theta$ 积分可得:

由如上两部分得到联合分布:

根据联合分布,求解下标为 $i=(m,n)$ 的词,即第 $m$ 篇文档中的第 $n$ 个词,的全部的条件概率。令 $\w = \brace{w_i=t, \w_{\neg{i}}}$ , $\z = \brace{z_i=k, \z_{\neg{i}}}$ ,有:

最终,我们需要根据Markov链的状态 $z_i$ 获取多项分布的参数 $\Theta$ 和 $\Phi$ 。根据贝叶斯法则和Dirichlet先验,以及公式(6)和(7):

其中, $\vec{n}_m$ 是构成文档 $m$ 的主题数向量, $\vec{n}_k$ 是构成主题 $k$ 的词项数向量。求解Dirichlet分布的期望,即公式(4)可得:

梳理一下Gibbs sampling中所用到的数据结构:统计量 $n_m^{(z)}$ 和 $n_z^{(t)}$ 分别是 $M \times K$ 和 $K \times V$ 矩阵,它们每行的和分别是 $M$ 维向量 $n_m$ (文档长度)和 $K$ 维向量 $n_z$ 。Gibbs sampling算法有三个阶段:初始化、burn-in和sampling。具体算法如下:

  • 算法: $\mathrm{LdaGibbs}(\brace{\w,\alpha,\beta,K})$
  • 输入:单词向量 $\w$ ,超参数 $\alpha$ 和 $\beta$ ,主题数 $K$
  • 全局变量:统计量 $\brace{n_m^{(k)}}$ 、 $\brace{n_k^{(t)}}$ ,以及它们的总数 $\brace{n_m}$ 、 $\brace{n_k}$ ,全部条件概率数组 $p(z_i\vert\cdot)$
  • 输出:主题向量 $\brace{\z}$ ,多项分布参数 $\Phi$ 和 $\Theta$ ,超参数估计量 $\alpha$ 和 $\beta$

[初始化]

  • 设置全局变量 $n_m^{(k)}$ 、 $n_k^{(t)}$ 、 $n_m$ 、 $n_k$ 为零
  • 对所有文档 $m \in [1, M]$ :
    • 对文档 $m$ 中的所有单词 $n \in [1, N_m]$ :
      • 初始化每个单词对应的主题 $z_{m,n}=k \sim \mathrm{Mult}(1/K)$
      • 增加“文档-主题”计数: $n_m^{(k)}+=1$
      • 增加“文档-主题”总数: $n_m +=1$
      • 增加“主题-词项”计数: $n_k^{(t)}+=1$
      • 增加“主题-词项”总数: $n_k +=1$

[迭代下面的步骤,直到Markov链收敛]

  • 对所有文档 $m \in [1, M]$ :
    • 对文档 $m$ 中的所有单词 $n \in [1, N_m]$ :
      • 删除该单词的主题计数: $n_m^{(k)}-=1; n_m-=1; n_k^{(t)}-=1; n_k-=1;$
      • 根据公式(8)采样出该单词的新主题: $\tilde{k} \sim p(z_i\vert\z_{\neg{i}}, \w)$
      • 增加该单词的新主题计数: $n_m^{(\tilde{k})}+=1; n_m+=1; n_\tilde{k}^{(t)}+=1; n_\tilde{k}+=1;$
  • 如果Markov链收敛:
    • 根据公式(9)生成主题-词项分布 $\Phi$
    • 根据公式(10)生成文档-主题分布 $\Theta$

参考

  • Gregor Heinrich, Parameter estimation for text analysis
  • David M.Blei, Andrew Y.Ng, Michael I.Jordan, Latent Dirichlet Allocation
  • Philip Resnik, Eric Hardisty, Gibbs Sampling for the Uninitiated