主题模型

机器学习

---

常见使用场景

在机器学习和自然语言处理中，主题模型是最近比较热门的一个话题，它能够从一个给定文档集合中发现在这些文档中隐藏的“主题”。当给定某个主题的文档，则在这篇文档中，和该主题相关的词语会更加频繁出现。例如，关于体育的新闻中，像“棒球”，“赛点”这些词语相比于“火锅”，“川菜”更加容易出现。而且通常一个文档不仅仅含有一个主题，它更有可能同时包含多个主题。例如，常见的武侠小说中，可能80%的文字在描写打斗，剩下的15%和5%的文字分别描述了爱情或者政治等主题。主题模型会通过训练集学习到组成每篇文章的主题分布和每个词语代表着各个主题的概率，即主题-词项分布和文章-主题分布，这样当遇到一篇新的文档时，我们就可以通过主题-词项分布确定新文档的主题分布。

LSA算法

LSA的目的是要从文本中发现隐含的语义维度，即主题。我们知道，在文档的空间向量模型（VSM）中，文档被表示成由特征词出现概率组成的多维向量，这种方法的好处是可以将查询和文档转化成同一空间下的向量计算相似度，可以对不同词项赋予不同的权重，在文本检索、分类、聚类问题中都得到了广泛应用。然而，向量空间模型没有能力处理一词多义和一义多词问题，因此会造成计算文档间相似度的不准确性。LSA一定程度上解决上述所描述的问题。具体说来就是对一个大型的文档集合使用一个合理的维度建模，并将词项和文档都表示到该空间中。比如有2000个文档，包含7000个索引词，LSA使用一个维度为100的向量空间将文档和词项表示到该空间，进而在该空间进行信息检索。而将文档表示到此空间的过程就是SVD奇异值分解和降维的过程。降维是LSA分析中最重要的一步，通过降维，去除了文档中的“噪音”，也就是无关信息（比如词的误用或不相关的词偶尔出现在一起），语义结构逐渐呈现。相比传统向量空间，潜在语义空间的维度更小，语义关系更明确。因此将文档和词语都在潜在语义空间表示，我们就能计算Document-Document，Document-Term，Term-Term的相似度。

LSA的步骤如下：
1. 分析文档集合，建立Term-Document矩阵。
2. 对Term-Document矩阵进行奇异值分解。
3. 对SVD分解后的矩阵进行降维，保留前$k$个特征值，后面$r-k$个置零，也就是低阶近似。
4. 使用降维后的矩阵构建潜在语义空间，或重建Term-Document矩阵。新得到的Term-Document矩阵就是我们经过LSA模型提取低维隐含语义空间。该空间中，每个奇异值对应的是每个“语义”维度的权重，我们刚才将不太重要的权重置为零，只保留最重要的维度信息，因而可以得到文档的一种更优表达形式。

假设我们有五篇文档和一个查询：
$d_1$ = Romeo and Juliet.
$d_2$ = Juliet: O happy dagger!
$d_3$ = Romeo died by dagger.
$d_4$ = "Live free or die", that’s the New-Hampshire’s motto.
$d_5$ = Did you know, New-Hampshire is in New-England.
$q$ = dies, dagger.
很明显当我们查询dies, dagger的时候，$d_3$应该是最相关的，因为在$d_3$中同时出现dies, dagger，$d_2$和$d_4$各自包含一个单词紧随其后。但是$d_1$和$d_5$该怎么排列？我们都知道$d_1$应该在$d_5$之前，但是机器怎么识别呢？解决方法正是LSA，在这个例子中，$d_1$中的Romeo和Juliet分别在$d_2$和$d_3$和单词dagger有关联，而单词dagger也是同理，对于$d_1$和$d_5$分别通过Romeo和New-Hampshire与dagger相联系。下面我们会通过更正式的方式来表达上述过程。

假设$A$是Term-Document矩阵，那么$A$的每列对应一篇文章，每行是一个单词。则$A$矩阵看上去是如下的样子：

$$\begin{matrix} &d_1&d_2&d_3&d_4&d_5\\ romeo&1&0&1&0&0\\ juliet&1&1&0&0&0\\ happy&0&1&0&0&0\\ dagger&0&1&1&0&0\\ live&0&0&0&1&0\\ die&0&0&1&1&0\\ free&0&0&0&1&0\\ new\ hampshire&0&0&0&1&1\\ \end{matrix}$$
很显然$B=A^TA$是Document-Document矩阵，如果文档i和j有b个相同的单词，则$B[i,j]=b$。另一方面$C=AA^T$是Term-Term矩阵，如果单词i和j同时出现在一篇文档的频数是c，则$C[i,j]=c$。很明显，矩阵B和C都是对称矩阵。我们对A做SVD分解。
$$A=U\Sigma V^T$$
其中U是矩阵B的特征向量，V是矩阵C的特征向量，$\Sigma$是矩阵B的特征值的平法根。
$$\Sigma=\begin{bmatrix} 2.285&0&0&0&0\\ 0&2.201&0&0&0\\ 0&0&1.361&0&0\\ 0&0&0&1.187&0\\ 0&0&0&0&0.797\\ \end{bmatrix}$$
我们只保留$\Sigma$前$K$个特征值，其余替换成0形成新的矩阵$\Sigma_k$。对矩阵$U$和$V$也只保留对于$K$项。然后计算“去噪”后的矩阵$A_k$。
$$A_k=U_k\Sigma_k V_k^T$$
因此，Term-Term与Document-Document在潜在语义空间的相关性矩阵可以表示为：
$$\begin{align} &A_kA_k^T = (U_k\Sigma_kV_k^T)(U_k\Sigma_kV_k^T)^T = U_k\Sigma_kV_k^TV_k\Sigma_k^TU_k^T = U_k\Sigma_k\Sigma_k^TU_k^T\\ &A_k^TA_k = (U_k\Sigma_kV_k^T)^T(U_k\Sigma_kV_k^T) = V_k\Sigma_k^TU_k^TU_k\Sigma_kV_k^T = V_k\Sigma_k^T\Sigma_kV_k^T \end{align}$$
因此我们可以把$U_k\Sigma_k$矩阵的行看作是term在潜在语义空间的坐标。同理我们认为$V_k\Sigma_k^T$的行表示了文档的坐标。

回到之前例子中，我们只取前$k=2$个特征值，则:

$$\Sigma_2=\begin{bmatrix} 2.285&0\\ 0&2.201\\ \end{bmatrix}$$
$$U_2=\begin{bmatrix} 0.396&0.280\\ 0.314&0.450\\ 0.178&0.269\\ 0.438&0.369\\ 0.264&0.346\\ 0.524&0.246\\ 0.264&0.346\\ 0.326&0.460 \end{bmatrix}$$
$$V_2^T=\begin{bmatrix} 0.311&0.407&0.594&0.603&0.143\\ 0.363&0.541&0.200&0.695&0.229 \end{bmatrix}$$
我们知道$U_2\Sigma_2$的行代表单词在潜在语义空间的坐标，同理$V_2\Sigma_2^T$表示文档的坐标。则：
$$\begin{align} &romeo=\begin{bmatrix} 0.905\\ 0.563 \end{bmatrix}, juliet=\begin{bmatrix} 0.717\\ 0.905 \end{bmatrix}, happy=\begin{bmatrix} 0.407\\ 0.541 \end{bmatrix}, dagger=\begin{bmatrix} 1.197\\ 0.494 \end{bmatrix},\\ &live=\begin{bmatrix} 0.603\\ 0.695 \end{bmatrix}, die=\begin{bmatrix} 1.001\\ 0.742 \end{bmatrix}, free=\begin{bmatrix} 0.603\\ 0.695 \end{bmatrix}, new\ hampshire=\begin{bmatrix} 0.745\\ 0.925 \end{bmatrix}, \end{align}$$
和
$$d_1=\begin{bmatrix} 0.711\\ 0.730 \end{bmatrix}, d_2=\begin{bmatrix} 0.930\\ 1.087 \end{bmatrix}, d_3=\begin{bmatrix} 1.357\\ 0.402 \end{bmatrix}, d_4=\begin{bmatrix} 1.378\\ 1.397 \end{bmatrix}, d_5=\begin{bmatrix} 0.327\\ 0.460 \end{bmatrix}$$
因此查询die, dagger也可以表示成两个 质心：
$$q=\frac{ \begin{bmatrix} 1.001\\ 0.742 \end{bmatrix}+ \begin{bmatrix} 1.197\\ 0.494 \end{bmatrix}}{2}= \begin{bmatrix} 1.099\\ 0.124 \end{bmatrix}$$
然后我们可以计算查询词和每片文档的距离：
$$\frac{d_iq}{|d_i||q|}$$
下图就是所有文档、词语和查询在潜在语义空间的坐标：

对于新文档，如何计算它在潜在语义空间的坐标呢？事实上，不需要重新进行SVD分解，由：

$$\begin{align} &A=U\Sigma V^T\\ &\Rightarrow U^TA=U^TU\Sigma V^T\\ &\Rightarrow \Sigma^{-1}U^TA=V^T\\ &\Rightarrow V=A^TU\Sigma^{-1} \tag 1 \end{align}$$
对于新的文档仅需在$A$矩阵加一列再代入(1)式重新计算出新文档的$V$然后再乘以$\Sigma$即可。

pLSA算法

我们了解到通过SVD可以进行LSA，把给定文档投影到语义空间，但其缺乏严谨的数理统计基础并且SVD分解非常耗时。pLSA则是从概率分布的角度建模的一种方法，它假设在词和文档之间有一层主题隐语义，而主题符合多项分布，一个主题中的词项也符合多项分布，由这两层分布的模型，生成各种文档。想象某个人要写$N$篇文档，他需要确定每篇文档里每个位置上的词。假定他一共有$K$个可选的主题，有$V$个可选的词项，所以，他制作了$K$个$V$面的“主题-词项”骰子，每个骰子对应一个主题，骰子每一面对应要选择的词项。然后，每写一篇文档会再制作一颗$K$面的“文档-主题”骰子；每写一个词，先扔该骰子选择主题；得到主题的结果后，使用和主题结果对应的那颗“主题-词项”骰子，扔该骰子选择要写的词。他不停的重复如上两个扔骰子步骤，最终完成了这篇文档。重复该方法$N$次，则写完所有的文档。在这个过程中，我们并未关注词和词之间的出现顺序，所以pLSA也是一种词袋方法。PLSA的概率图模型如下：

实心的节点$d$和$w$表示我们能观察到的文档和单词，空心的$z$表示我们观察不到的隐藏变量，用来表示隐含的主题。图中用了所谓的“盘子记法”，即用方框表示随机变量的重复。这里方框右下角的字母$M$和$N$分别表示有$M$篇文档，第$j$篇文档有$N$个单词。每条有向边表示随机变量间的依赖关系。也就是说，pLSA 假设每对$(d,w)$都是由下面的过程产生的：
1. 以$p(d)$的先验概率选择一篇文档
2. 选定$d$后，以$p(z|d)$的概率选中主题
3. 选中主题$z$后，以$p(w|z)$的概率选中单词

而我们感兴趣的正是其中的$p(z|d)$和$p(w|z)$，利用前者我们可以知道每篇文章中各主题所占的比重，利用后者我们则能知道各单词在各主题中出现的概率，从而进一步找出各主题的“关键词”。记$\theta=(p(z|d),p(w|z))$， 表示我们希望估计的模型参数。当然$\theta$不仅仅代表两个数，而是对于每对$(z^k,d^i)$和$(z^k,w^j)$， 我们都要希望知道$p(z^k|d^i)$和$p(w^j|z^k)$的值。也就是说，模型中共有$|Z|*|D|+|Z|*|W|$个参数。我们还知道：

$$\begin{align} &p(d,w) = p(d)p(w|d)\\ &p(w|d) = \sum_z p(w|z)p(z|d) \end{align}$$
于是，可以很容易写出样本分布的对数似然函数：
$$\begin{align} L(\theta) &= log\prod_{d,w} p(d,w)^{n(d,w)}\\ &= \sum_{d,w}n(d,w)log \sum_zp(d)p(w|z)p(z|d)\\ &= \sum_{d,w}n(d,w)log \sum_zp(w|z)p(z|d) \end{align}$$

因为$p(d)$和参数$\theta$无关，所以第三步略掉，这里出现$log$套$\sum$的形式，导致很难直接拿它做最大似然。但假如能观察到$z$，问题就很简单了。于是对于这种含有隐变量的最大似然估计，我们还是需要使用EM方法。

EM算法的步骤是：
1. E步骤：求隐含变量在给定当前估计的参数条件下的后验概率。
2. M步骤：最大化complete data对数似然函数的期望，此时我们使用E步骤里计算的隐含变量的后验概率，得到新的参数值，两步迭代进行直到收敛。

先解释一下什么是incomplete data和complete data。当原始数据的似然函数很复杂时，我们通过增加一些隐含变量来增强我们的数据，得到complete data，而complete data的似然函数更加简单，方便求极大值。于是，原始的数据就成了incomplete data。我们将会看到，我们可以通过最大化complete data似然函数的期望来最大化incomplete data的似然函数，以便得到求似然函数最大值更为简单的计算途径。

针对我们pLSA参数估计问题，在E步骤中，直接使用贝叶斯公式计算隐含变量在当前参数取值条件下的后验概率，有：

$$P(z^k|d^i,w^j;\theta_t) = \frac{P_t(z^k|d^i)P_t(w^j|z^k)} {\sum_z P_t(z|d^i)P_t(w^j|z)}$$
而$P_t(z^k|d^i)$和$P_t(w^j|z^k)$就是上轮迭代求出的，这样就完成了E-step。接下来M-step就是带入隐变量的后验概率，最大化样本分布的对数似然函数，更新参数$P_t(z^k|d^i)$和$P_t(w^j|z^k)$。
$$\begin{align} Q(\theta,\theta_{old})=E[log\ p(w,z,d)]&=\sum_{d,w}n(d,w)\sum_k p(z^k|d^i,w^j;\theta_t)log\ p(w,z,d)\\ &= \sum_{d,w}n(d,w)\sum_k p(z^k|d^i,w^j;\theta_t)log\ p(z|d)p(w|z) \end{align}$$
注意这里$p(z^k|d^i,w^j;\theta_t)$是已知的，取的是前面E步骤里面的估计值。下面我们来最大化期望，这又是一个带约束条件求极值的问题，可以用拉格朗日乘数法。拉格朗日乘数法可以把条件极值问题转化为无条件极值问题，在pLSA中目标函数就是,约束条件是：
$$\begin{align} \sum_k p(z^k|d^i) &= 1\\ \sum_j p(w^j|z^k) &= 1 \end{align}$$
这里我们引入两个乘子$\lambda$和$u$，可以写出拉格朗日函数，如下：
$$\mathcal{H}=E[log\ p(w,z,d)]+\lambda\sum_i(1-\sum_k p(z^k|d^i))+u\sum_k(1-\sum_j p(w^j|z^k))$$
这是一个关于$\theta=(p(z|d),p(w|z))$的函数，分别对其求偏导数，我们可以得到：
$$\begin{align*} P_{t+1}(w^{(j)}|z^{(k)}) &= \frac {\sum_d n(d,w^{(j)})P(z^{(k)}|d,w^{(j)};\theta_t)} {\sum_{d,w} n(d,w)P(z^{(k)}|d,w;\theta_t)} \\ P_{t+1}(z^{(k)}|d^{(i)}) &= \frac {\sum_w n(d^{(i)},w)P(z^{(k)}|d^{(i)},w;\theta_t)} {\sum_{w,z} n(d,w)P(z|d^{(i)},w;\theta_t)} \end{align*}$$
然后使用更新后的参数值$\theta=(p_{t+1}(z|d),p_{t+1}(w|z))$，我们又进入E步骤，计算隐含变量给出当前估计的参数条件下的后验概率。如此不断迭代，直到满足终止条件。

LDA算法

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

和pLSA模型类似，其中，单圆圈表示隐变量，双圆圈表示观察到的变量。对应到上图，只有$w_{m,n}$是观察到的变量，其他都是隐变量或者参数，其中$\alpha$和$\beta$是超参数。方框中$\Phi=\{\vec\varphi_k\}_{k=1}^K$表示有$K$种“主题-词项”分布，$\Theta=\{\vec\vartheta_m\}_{m=1}^M$有$M$种“文档-主题”分布，即对每篇文档都会产生一个$\vec\vartheta_m$分布。每篇文档$m$中有$n$个词，每个词$w_{m,n}$都有一个主题$z_{m,n}$，该词实际是由$\vec\varphi_{z_{m,n}}$产生。具体过程下面再说。

贝叶斯估计

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

$$\rm posterior = {likelihood \cdot prior \over evidence}$$
即
$$\rm p(X|\theta) = {p(\theta|X) \cdot p(\theta) \over p(X)}$$
在MLE和MAP中，由于是要求函数最大值时的参数，所以都不会考虑evidence。但在贝叶斯估计中，不再直接取极值，所以还会考虑evidence，下面的这个积分也是通常贝叶斯估计中最难处理的部分：
$$\rm p(X) = \int p(\theta|X) \cdot p(\theta)d\theta$$

共轭先验

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

由于pLSA中参数“主题-词项”$p(w|z)$和“文档-主题”$p(z|d)$都服从多项分布，因此选择Dirichlet分布作为他们的共轭先验。在pLSA中，假设这俩多项分布是确定的，我们在已经假设确定的分布下，选择具体的主题和词项；但在LDA中，这俩多项分布的参数是被当作随机变量，服从另一层先验分布，所以，Dirichlet分布也称作分布的分布，其定义如下：

$$\begin{align} \mathrm{Dir}(\vec{p}|\vec{\alpha}) &\triangleq {\Gamma \left(\sum_{k=1}^K \alpha_k \right) \over \prod_{k=1}^K \Gamma(\alpha_k)} \prod_{k=1}^K p_k^{\alpha_k -1} \tag 1 \\ &\triangleq \frac{1}{\Delta(\vec{\alpha})} \prod_{k=1}^K p_k^{\alpha_k -1} \end{align}$$
其中，$\vec{p}$是要猜测的随机向量，$\vec{\alpha}$是超参数，$\Delta(\vec{\alpha})$称作Delta函数，是Dirichlet分布的归一化系数。
多项分布定义：
$$\mathrm{Mult}(\vec{n}|\vec{p}, N) \triangleq {N \choose \vec{n}} \prod_{k=1}^K p_k^{n_k} \tag 2$$
其中，$\vec{p}$和$\vec{n}$ 服从约束$\sum_k p_k=1$和$\sum_k n_k = N$。而$\Gamma(x+1) = x!$就是阶乘在实数集上的扩展。因此，公式(1)和(2)有相同的形式，所以，这俩分布也称作Dirichlet-Multinomail共轭。

接着我们再来看看，对于非负实数$\alpha$和$\beta$，我们有如下关系：

$$Beta(p|\alpha,\beta)+BioCount(m_1,m_2)=Beta(\alpha+m_1,\beta+m_2)$$
其中$BioCount(m_1,m_2)$对应的是二项分布$B(m_1+m_2,p)$的计数。从以上我们可以看出，Beta分布中的超参数$\alpha$和$\beta$可以理解为物理计数，也被称为伪计数(pseudo-count)。同理推广到多项式分布有：
$$Dir(\vec p|\vec\alpha)+MultCount(\vec m)=Dir(\vec p|\vec\alpha+\vec m)$$
其中$MultCount(\vec m)$对应的是多项分布$Mult(\vec m,\vec p)$的计数。同理，Dir分布中的超参数$\vec\alpha$也可以理解为伪计数(pseudo-count)

现在考虑Dirichlet分布的特殊情况，超参数$\vec\alpha$中的值相等，即$\vec\alpha=\alpha\vec I$，变量$\alpha$也称为聚集参数。下图就是$\alpha$取不同值时候对称Dirichlet分布的图形：

很明显可以看到，当$\alpha=1$时，Dirichlet分布退化为均匀分布。当$\alpha>1$时，先验概率会更倾向于$p_1=p_2=...=p_k$的情况。当$\alpha<1$时，先验概率会更倾向于$p_i=1,p_{other}=0$的情况。换句话说，当聚合参数$\alpha=1$时，先验退化为均匀分布。当大于1时，LDA就聚合在更多的主题，相反小于1，则聚合在少数几个主题。

生成过程

在LDA中，“文档-主题”向量$\vec{\vartheta}_m$由超参数为$\vec\alpha$的Dirichlet分布生成，“主题-词项”向量$\vec{\varphi}_k$由超参数为$\vec\beta$ 的Dirichlet分布生成。

符号解释：

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

根据概率图，整个样本集合的生成过程如下：

1. 对所有的主题$k\in[1,K]$(即生成$\Phi_{K\times N}$，$K$个主题，对应词表$V$共有$N$个词)：
   
   • 生成$K$个主题的“主题-词项”分布$\vec{\varphi}_k \sim \mathrm{Dir}(\vec{\beta})$($N$维矩阵，对应词表$V$中的每个词项的概率)
2. 对所有的文档$m \in [1,M]$：
   
   • 生成当前文档$m$相应的“文档-主题”分布$\vec{\vartheta}_m \sim \mathrm{Dir}(\vec{\alpha})$ ($K$维向量，即第$m$篇文档对应的每个主题的概率)
   • 生成当前文档$m$的长度$N_m \sim \mathrm{Poiss(\xi)}$
   • 对当前文档$m$中的所有词$n \in [1, N_m]$：
     
     • 生成当前位置的词的所属主题$z_{m,n} \sim \mathrm{Mult}(\vec\vartheta_m)$
     • 根据之前生成的主题分布$\Phi_{K\times N}$，生成当前位置的词的相应词项$w_{m,n} \sim \mathrm{Mult}(\vec\varphi_{z_{m,n}})$

上述过程可以分解为

• $\vec\alpha\rightarrow \vec\theta_m \rightarrow z_{m,n}$, 这个过程表示在生成第$m$篇文档的时候，先从“文档-主题”分布$\vec\vartheta_m$中选出第$n$个词的topic编号$z_{m,n}$；
• $\vec\beta \rightarrow \vec\varphi_k \rightarrow w_{m,n} | k=z_{m,n}$, 这个过程表示用如下动作生成语料中第$m$篇文档的第$n$个词：在上帝手头的$K$个“主题-词项”分布$\Phi_{K\times N}$中，挑选主题为$k=z_{m,n}$的那个$\vec\varphi_{z_{m,n}}$分布生成 $w_{m,n}$；

理解 LDA最重要的就是理解这两个物理过程。 LDA 模型在基于$K$个 topic 生成语料中的$M$篇文档的过程中， 由于是 bag-of-words模型，有一些物理过程是相互独立可交换的。由此， LDA 生成模型中，$M$篇文档会对应于$M$个独立的 Dirichlet-Multinomial 共轭结构；$K$个 topic 会对应于 K 个独立的Dirichlet-Multinomial共轭结构。所以理解LDA所需要的所有数学就是理解Dirichlet-Multiomail共轭，其它都就是理解物理过程。现在我们进入细节， 来看看 LDA 模型是如何被分解为 $M$+$K$ 个Dirichlet-Multinomial 共轭结构的。

由第一个物理过程，我们知道 $\vec\alpha\rightarrow \vec\theta_m \rightarrow z_{m,n}$ 表示生成第 $m$ 篇文档中的所有词对应的topics，显然 $\vec\alpha\rightarrow \vec\theta_m$ 对应于 Dirichlet 分布， $\vec\theta_m \rightarrow z_{m,n}$ 对应于 Multinomial 分布， 所以整体是一个 Dirichlet-Multinomial 共轭结构；

$$\begin{align} p(\vec{z}_m |\vec{\alpha}) & = \int p(\vec{z}_m|\vec\theta_m) p(\vec\theta_m|\vec{\alpha})d\vec\theta_m \notag \\ & = \int \prod_{k=1}^V p_k^{n_k} Dir(\vec\theta_m|\vec{\alpha}) d\vec\theta_m \notag \\ & = \int \prod_{k=1}^V p_k^{n_k} \frac{1}{\Delta(\vec{\alpha})} \prod_{k=1}^V p_k^{\alpha_k -1} d\vec{p} \notag \\ & = \frac{1}{\Delta(\vec{\alpha})} \int \prod_{k=1}^V p_k^{n_k + \alpha_k -1} d\vec{p} \notag \\ & = \frac{\Delta(\vec{n}_m+\vec{\alpha})}{\Delta(\vec{\alpha})} \end{align}$$
其中 $\vec n_m = (n_{m}^{(1)}, \cdots, n_{m}^{(K)})$ 表示第 $m$ 篇文档中第 $k$ 个topic 产生的词的个数。由于语料中 $M$ 篇文档的 topics 生成过程相互独立，所以我们得到 $M$ 个相互独立的 Dirichlet-Multinomial 共轭结构，从而我们可以得到整个语料中 topics 生成概率
$$\begin{align} p(\vec{\mathbf{z}} |\vec{\alpha}) & = \prod_{m=1}^M p(\vec{z}_m |\vec{\alpha}) \notag \\ &= \prod_{m=1}^M \frac{\Delta(\vec{n}_m+\vec{\alpha})}{\Delta(\vec{\alpha})} \end{align}$$

目前为止，我们由 $M$ 篇文档得到了 $M$ 个 Dirichlet-Multinomial 共轭结构，还有额外 $K$ 个 Dirichlet-Multinomial共轭结构在哪儿呢？为了形象话说明LDA生成文档的过程，我们认为假设文档是通过抛骰子生成的。在上帝按照之前的规则玩 LDA 游戏的时候，上帝是先完全处理完成一篇文档，再处理下一篇文档。文档中每个词的生成都要抛两次骰子，第一次抛一个doc-topic骰子得到 topic, 第二次抛一个topic-word骰子得到 word，每次生成每篇文档中的一个词的时候这两次抛骰子的动作是紧邻轮换进行的。如果语料中一共有 N 个词，则上帝一共要抛 2N次骰子，轮换的抛doc-topic骰子和 topic-word骰子。但实际上有一些抛骰子的顺序是可以交换的，我们可以等价的调整2N次抛骰子的次序：前N次只抛doc-topic骰子得到语料中所有词的 topics，然后基于得到的每个词的 topic 编号，后N次只抛topic-word骰子生成 N 个word。于是上帝在玩 LDA 游戏的时候，可以等价的按照如下过程进行：

1. 上帝有两个大坛子的骰子，第一个坛子装的是doc-topic骰子，第二个坛子装的是topic-word骰子。
2. 上帝随机从第二坛子中独立选取 $K$ 个topic-word骰子，编号1到K。
3. 每次生成一篇新文档前，上帝先从第一个坛子中随机抽取一个doc-topic骰子， 然后重复投掷这个doc-topic骰子，为每个词生成一个topic编号$z$。重复如上过程处理每篇文档，为每篇文档每个词生成topic编号，但是词还未生成。
4. 从头到尾，对语聊中每篇文档的每个topic编号$z$，选择 $K$ 个topic-word骰子中编号为$z$的那个，投掷这个骰子，于是生成对应的词语，直到生成完所有词语。

以上游戏是先生成了语料中所有词的 topic, 然后对每个词在给定 topic 的条件下生成 word。在语料中所有词的 topic 已经生成的条件下，任何两个 word 的生成动作都是可交换的。于是我们把语料中的词进行交换，把具有相同 topic 的词放在一起

$$\begin{align*} \vec{\mathbf{w}}’ &= (\vec{w}_{(1)}, \cdots, \vec{w}_{(K)}) \\ \vec{\mathbf{z}}’ &= (\vec{z}_{(1)}, \cdots, \vec{z}_{(K)}) \end{align*}$$
其中，$\vec w_{(k)}$ 表示这些词都是由第 k 个 topic 生成的， $\vec z_{(k)}$ 对应于这些词的 topic 编号，所以$\vec z_{(k)}$中的分量都是k。

对应于概率图中的第二个物理过程 $\vec\beta \rightarrow \vec\varphi_k \rightarrow w_{m,n} | k=z_{m,n}$，在 $k=z_{m,n}$ 的限制下，语料中任何两个由 topic k 生成的词都是可交换的，即便他们不再同一个文档中，所以我们此处不再考虑文档的概念，转而考虑由同一个 topic 生成的词。考虑如下过程 $\vec\beta \rightarrow \vec\varphi_k \rightarrow w_{(k)}$ ，容易看出， 此时 $\vec\beta \rightarrow \vec\varphi_k$ 对应于 Dirichlet 分布， $\vec\varphi_k \rightarrow w_{(k)}$ 对应于 Multinomial 分布， 所以整体也还是一个 Dirichlet-Multinomial 共轭结构；

同样的，我们可以得到

$$p(\vec{w}_{(k)} |\vec{\beta}) = \frac{\Delta(\vec{n}_k+\vec{\beta})}{\Delta(\vec{\beta})}$$
其中 $\vec n_k = (n_{k}^{(1)}, \cdots, n_{k}^{(V)})$ 表示中第 $k$ 个topic 产生的词中 word $t$的个数。而语料中 $K$ 个 topics 生成words 的过程相互独立，所以我们得到 $K$ 个相互独立的 Dirichlet-Multinomial 共轭结构，从而我们可以得到整个语料中词生成概率
$$\begin{align} p(\vec{\mathbf{w}} |\vec{\mathbf{z}},\vec{\beta}) &= p(\vec{\mathbf{w}}’ |\vec{\mathbf{z}}’,\vec{\beta}) \notag \\ &= \prod_{k=1}^K p(\vec{w}_{(k)} | \vec{z}_{(k)}, \vec{\beta}) \notag \\ &= \prod_{k=1}^K \frac{\Delta(\vec{n}_k+\vec{\beta})}{\Delta(\vec{\beta})} \end{align}$$

因此我们可以得到联合概率分布$p(\vec{w},\vec{z})$：

$$\begin{align} p(\vec{\mathbf{w}},\vec{\mathbf{z}} |\vec{\alpha}, \vec{\beta}) &= p(\vec{\mathbf{w}} |\vec{\mathbf{z}}, \vec{\beta}) p(\vec{\mathbf{z}} |\vec{\alpha}) \notag \\ &= \prod_{k=1}^K \frac{\Delta(\vec{n}_k+\vec{\beta})}{\Delta(\vec{\beta})} \prod_{m=1}^M \frac{\Delta(\vec{n}_m+\vec{\alpha})}{\Delta(\vec{\alpha})} \end{align}$$

因此，由该生成过程可知，第 $m$ 篇文档中第 $n$ 个词 $t$ 的生成概率：

$$p(w_{m,n}=t|\vec{\vartheta}_m, \Phi) = \sum_{k=1}^K p(w_{m,n}=t|\vec{\varphi}_k) p(z_{m,n}=k|\vec{\vartheta}_m)$$

LDA Gibbs Sampler

有了联合分布$p(\vec{\mathbf{w}},\vec{\mathbf{z}})$, 万能的 MCMC 算法就可以发挥作用了！于是我们可以考虑使用 Gibbs Sampling 算法对这个分布进行采样。通常均匀分布Uniform(0,1)的样本，即我们熟悉的类rand()函数，可以由线性同余发生器生成；而其他的随机分布都可以在均匀分布的基础上，通过某种函数变换得到，比如，正态分布可以通过Box-Muller变换得到。然而，这种变换依赖于计算目标分布的积分的反函数，当目标分布的形式很复杂，或者是高维分布时，很难简单变换得到，而随机模拟(MCMC)就是求解近似解的一种强有力的方法。随机模拟的核心就是对一个分布进行抽样(Sampling)，更多关于随机抽样可以参考这篇博客。

当然由于 $\vec{\mathbf{w}}$ 是观测到的已知数据，只有 $\vec{\mathbf{z}}$ 是隐含的变量，所以我们真正需要采样的是分布 $p(\vec{\mathbf{z}}|\vec{\mathbf{w}})$。语料库 $\vec{\mathbf{w}}$ 中的第$i$个词我们记为$w_i$, 其中$i=(m,n)$是一个二维下标，对应于第$m$篇文档的第 $n$个词，我们用 $\neg{i}$ 表示去除下标为i的词。那么按照 Gibbs Sampling 算法的要求，我们要求得任一个坐标轴 i 对应的条件分布 $p(z_i=k|\vec{z}_{\neg{i}},\vec{\mathbf{w}})$ 。假设已经观测到的词 $w_i=t$， 则由贝叶斯法则，我们容易得到

$$\begin{align*} p(z_i = k|\vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}}) \propto p(z_i = k, w_i = t |\vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}}_{\neg i}) \\ \end{align*}$$
由于$z_i=k,w_i=t$ 只涉及到第 m 篇文档和第 k 个 topic，所以上式的条件概率计算中, 实际上也只会涉及到如下两个Dirichlet-Multinomial 共轭结构
1. $\vec\alpha\rightarrow \vec\theta_m \rightarrow z_{m,n}$,
2. $\vec\beta \rightarrow \vec\varphi_k \rightarrow w_{m,n} | k=z_{m,n}$,
其它的 M+K-2 个 Dirichlet-Multinomial 共轭结构和$z_i=k,w_i=t$是独立的。因此语料去掉第 i 个词对应的 $(z_i,w_i)$，并不改变我们之前讨论的 M+K 个 Dirichlet-Multinomial 共轭结构，只是某些地方的计数会减少。所以 $\vec{\theta}_m, \vec{\varphi}_k$ 的后验分布都是 Dirichlet:
$$\begin{align*} p(\vec{\theta}_m|\vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}}_{\neg i}) &= Dir(\vec{\theta}_m| \vec{n}_{m,\neg i} + \vec{\alpha}) \\ p(\vec{\varphi}_k|\vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}}_{\neg i}) &= Dir( \vec{\varphi}_k| \vec{n}_{k，\neg i} + \vec{\beta}) \end{align*}$$

使用上面两个式子，把以上想法综合一下，我们就得到了如下的 Gibbs Sampling 公式的推导

$$\begin{align*} p(z_i = k|\vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}}) & \propto p(z_i = k, w_i = t |\vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}}_{\neg i}) \\ &= \int p(z_i = k, w_i = t, \vec{\theta}_m,\vec{\varphi}_k | \vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}}_{\neg i}) d \vec{\theta}_m d \vec{\varphi}_k \\ &= \int p(z_i = k, \vec{\theta}_m|\vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}}_{\neg i}) \cdot p(w_i = t, \vec{\varphi}_k | \vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}}_{\neg i}) d \vec{\theta}_m d \vec{\varphi}_k \\ &= \int p(z_i = k |\vec{\theta}_m) p(\vec{\theta}_m|\vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}}_{\neg i}) \cdot p(w_i = t |\vec{\varphi}_k) p(\vec{\varphi}_k|\vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}}_{\neg i}) d \vec{\theta}_m d \vec{\varphi}_k \\ &= \int p(z_i = k |\vec{\theta}_m) Dir(\vec{\theta}_m| \vec{n}_{m,\neg i} + \vec{\alpha}) d \vec{\theta}_m \\ & \hspace{0.2cm} \cdot \int p(w_i = t |\vec{\varphi}_k) Dir( \vec{\varphi}_k| \vec{n}_{k,\neg i} + \vec{\beta}) d \vec{\varphi}_k \\ &= \int \theta_{mk} Dir(\vec{\theta}_m| \vec{n}_{m,\neg i} + \vec{\alpha}) d \vec{\theta}_m \cdot \int \varphi_{kt} Dir( \vec{\varphi}_k| \vec{n}_{k,\neg i} + \vec{\beta}) d \vec{\varphi}_k \\ &= E(\theta_{mk}) \cdot E(\varphi_{kt}) \\ &= \hat{\theta}_{mk} \cdot \hat{\varphi}_{kt} \\ \end{align*}$$

$z_i=k,w_i=t$的概率只和两个 Dirichlet-Multinomial 共轭结构关联。而最终得到的 $\hat\theta_{mk},\ \hat\varphi_{kt}$ 就是对应的两个 Dirichlet 后验分布在贝叶斯框架下的参数估计。借助于前面介绍的Dirichlet 参数估计的公式 ，我们有

$$\begin{align*} \hat{\theta}_{mk} &= \frac{n_{m,\neg i}^{(k)} + \alpha_k}{\sum_{k=1}^K (n_{m,\neg i}^{(k)} + \alpha_k)} \tag 1\\ \hat{\varphi}_{kt} &= \frac{n_{k,\neg i}^{(t)} + \beta_t}{\sum_{t=1}^V (n_{k,\neg i}^{(t)} + \beta_t)} \tag 2 \end{align*}$$
于是，我们最终得到了 LDA 模型的 Gibbs Sampling 公式
$$\begin{equation} p(z_i = k|\vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}}) \propto \frac{n_{m,\neg i}^{(k)} + \alpha_k}{\sum_{k=1}^K (n_{m,\neg i}^{(k)} + \alpha_k)} \cdot \frac{n_{k,\neg i}^{(t)} + \beta_t}{\sum_{t=1}^V (n_{k,\neg i}^{(t)} + \beta_t)} \tag 3 \end{equation}$$

Training

有了 LDA 模型，当然我们的目标有两个

• 估计模型中的参数 $\vec{\varphi}_1, \cdots, \vec{\varphi}_K$ 和 $\vec{\theta}_1, \cdots, \vec{\theta}_K$
• 对于新来的一篇文档 $doc_{new}$，我们能够计算这篇文档的 topic 分布 $\vec{\theta}_{new}$。

有了 Gibbs Sampling 公式， 我们就可以基于语料训练 LDA 模型，并应用训练得到的模型对新的文档进行 topic 语义分析。训练的过程就是获取语料中的 $(z,w)$ 的样本，而模型中的所有的参数都可以基于最终采样得到的样本进行估计。统计量$n_m^{(z)}$和$n_z^{(t)}$分别是M×K和K×V矩阵，它们每行的和分别是M维向量$n_m$(文档长度)和K维向量$n_k$。Gibbs sampling算法有三个阶段：初始化、burn-in和sampling。具体算法如下：

• 算法：$\mathrm{LdaGibbs}({\vec{w},\alpha,\beta,K})$
• 输入：单词向量$\vec{w}$ ，超参数$\alpha$和$\beta$，主题数K
• 全局变量：统计量${n_m^{(k)}}$、${n_k^{(t)}}$，以及它们的总数${n_m}、{n_k}$，全部条件概率数组$p(z_i|\cdot)$
• 输出：主题向量${\vec{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$
    • 增加“主题-词项”总数：$nk+=1$

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

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

由这个topic-word 频率矩阵我们可以计算每一个p(word|topic)概率，从而算出模型参数$\vec{\varphi}_1, \cdots, \vec{\varphi}_K$， 这就是上帝用的 K 个 topic-word 骰子。当然，语料中的文档对应的骰子参数 $\vec{\theta}_1, \cdots, \vec{\theta}_M$， 在以上训练过程中也是可以计算出来的，只要在 Gibbs Sampling 收敛之后，统计每篇文档中的 topic 的频率分布，我们就可以计算每一个 p(topic|doc) 概率，于是就可以计算出每一个$\vec{\theta}_m$。由于参数$\vec{\theta}_m$ 是和训练语料中的每篇文档相关的，对于我们理解新的文档并无用处，所以工程上最终存储 LDA 模型时候一般没有必要保留。通常，在 LDA 模型训练的过程中，我们是取 Gibbs Sampling 收敛之后的 n 个迭代的结果进行平均来做参数估计，这样模型质量更高。

Inference

有了 LDA 的模型，对于新来的文档 $doc_{new}$, 我们如何做该文档的 topic 语义分布的计算呢？基本上 inference 的过程和 training 的过程完全类似。对于新的文档， 我们只要认为 Gibbs Sampling 公式中的 $\hat\varphi_{kt}$ 部分是稳定不变的，是由训练语料得到的模型提供的，所以采样过程中我们只要估计该文档的 topic 分布 $\vec{\theta}_{new}$ 就好了

[随机初始化]

• 设置全局变量${n_m^{(k)}}$、${n_m}$为零
• 对所有文档 $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$

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

• 对所有文档 $m \in [1, M]$：
  
  • 对文档 m 中的所有单词 $m \in [1, M]$：
    
    • 删除该文档的主题计数：$n_m^{(k)}-=1; n_m-=1;$
    • 通过Training得到的$\hat\varphi_{kt}$和迭代得到的$\hat\varphi_{kt}$，带入公式(3)采样出该文档的新主题：$\tilde{k} \sim p(z_i|\vec{z}_{\neg{i}}, \vec{w})$
    • 增加该文档的新主题计数：$n_m^{(k)}+=1; n_m+=1;$
• 如果Markov链收敛：
  
  • 根据公式(2)生成文档-主题分布 $\Theta$

参考资料

1. http://www.engr.uvic.ca/~seng474/svd.pdf
2. http://blog.jqian.net/post/plsa.html
3. http://blog.tomtung.com/2011/10/plsa/
4. http://blog.jqian.net/post/lda.html#toc_0