Probabilistic latent semantic analysis (概率潜在语义分析,pLSA) 是一种 Topic model,在99年被 Thomas Hofmann 提出。它和随后提出的 LDA 使得 Topic Model 成为了研究热点,其后的模型大都是建立在二者的基础上的。

我们有时会希望在数量庞大的文档库中自动地发现某些结构。比如我们希望在文档库发现若干个“主题”,并将每个主题用关键词的形式表现出来。 我们还希望知道每篇文章中各个主题占得比重如何,并据此判断两篇文章的相关程度。而 pLSA 就能完成这样的任务。

我之前取了 Wikinews 中的 1000 篇新闻,试着用 pLSA 在其中发现 K=15 个主题。比如一篇关于 Wikileaks 的阿萨奇被保释消息的新闻,算法以 100% 的概率把它分给了主题 9,其关键词为:

media phone hacking wikileaks assange australian stated information investigation murdoch

可以看到这个主题的发现还是非常靠谱的。又比如这条中国人民的老朋友威胁要大打打核战争 的新闻。 算法把它以 97.7% 的概率分给了主题 3,2.3% 的概率分给了主题 7。主题 3 的关键词是:

south north court china military death tornado service million storm

主题 7 的关键词是:

nuclear plant power japan million carbon radiation china water minister

可以看到这条新闻和主题 3 中的“南北”、“军事”、“中国”、“死亡”这些信息联系在一起,和主题 7 中的“核”、“中国”联系在一起。 应该是因为我的数据集中与北朝鲜核问题相关的新闻只有不到 10 条,而 10 个主题的划分并不够细致,所以关于“朝核问题”或者“核武器”的这样的主题并没能被分离出来。但可以看到即使是这样结果也是很 make sense 的。

那我们就来看看 pLSA 模型是怎么回事吧。和很多模型一样,pLSA 遵从 bag-of-words 假设, 即只考虑一篇文档中单词出现的次数,忽略单词的先后次序关系,且每个单词的出现都是彼此独立的。 这样一来,我们观察到的其实就是每个单词 \(w \in W\) 在每篇文档 \(d \in D\) 中出现的次数 \(n(w,d)\)。 pLSA 还假设对于每对出现的 \((d,w)\) 都对应着一个表示“主题”的隐藏变量 \(z \in Z\)。 pLSA 是一个生成模型,它假设 \(d\)\(w\)\(z\) 之间的关系用贝叶斯网络表示是这样的(从 [Blei03] 偷的图):

image0

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

  1. \(P(d)\) 的先验概率选择一篇文档 \(d\)
  2. 选定 \(d\) 后,以 \(P(z|d)\) 的概率选中主题 \(z\)
  3. 选中主题 \(z\) 后,以 \(P(w|z)\) 的概率选中单词 \(w\)

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

\begin{align*} P(d,w) &= P(d)P(w|d) \\ P(w|d) &= \sum_z P(w|z)P(z|d) \end{align*}

根据最大log似然估计法,我们要求的就是

\begin{align*} \arg\max_\theta L(\theta) &= \arg\max_\theta \sum_{d,w} n(d,w)\log P(d,w;\theta) \\ &= \arg\max_\theta \sum_{d,w} n(d,w)\log P(w|d;\theta)P(d) \\ &= \arg\max_\theta \left\{ \sum_{d,w} n(d,w)\log P(w|d;\theta) + \sum_{d,w} n(d,w)\log P(d) \right\} \end{align*}

这里由于 \(\sum_{d,w} n(d,w)\log P(d)\) 这一项与 \(\theta\) 无关,在 \(\arg\max_\theta\) 中可以被直接扔掉。 [1]

因此

\begin{align*} \arg\max_\theta L(\theta) &= \arg\max_\theta \sum_{d,w} n(d,w)\log P(w|d;\theta) \\ &= \arg\max_\theta \sum_{d,w} n(d,w)\log \sum_z P(w|z)P(z|d) \end{align*}

这里出现了 \(\log\)\(\sum\) 的形式,导致很难直接拿它做最大似然。但假如能观察到 \(z\),问题就很简单了。于是我们想到根据 EM 算法 (参见我的上篇笔记),可以用下式迭代逼近 \(\arg\max_\theta L(\theta)\)

\begin{equation*} \arg\max_\theta Q_t(\theta) = \arg\max_\theta \sum_{d,w} n(d,w) E_{z|d,w;\theta_t}[\log P(w, z|d;\theta)] \end{equation*}

其中

\begin{align*} E_{z|d,w;\theta_t}[\log P(w, z|d;\theta)] &= \sum_z P(z|d,w;\theta_t) \log P(w, z|d;\theta) \\ &= \sum_z P(z|d,w;\theta_t) [\log P(w|z) + \log P(z|d)] \\ \end{align*}

在 E-step 中,我们需要求出 \(Q_t(\theta)\) 中除 \(\theta\) 外的其它未知量,也就是说对于每组 \((d^{(i)}, w^{(j)}, z^{(k)})\) 我们都需要求出 \(P(z^{(k)}|d^{(i)},w^{(j)};\theta_t)\)。 根据贝叶斯定理贝叶斯定理,我们知道:

\begin{equation*} 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)} \end{equation*}

\(P_t(z|d)\)\(P_t(w|z)\) 就是上轮迭代求出的 \(\theta_t\)。这样就完成了 E-step。

接下来 M-step 就是要求 \(\arg\max_\theta Q_t(\theta)\) 了。利用基本的微积分工具 [2],可以分别对每对 \((w^{(j)}, z^{(k)})\)\((d^{(i)}, z^{(k)})\) 求出:

\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*}

以上就是 pLSA 算法了。最后贴个我用 MATLAB 写的实现 [3]

function [p_w_z, p_z_d, Lt] = pLSA(n_dw, n_z, iter_num)
% PLSA  Fit a pLSA model on given data
%       in which n_dw(d,w) is the number of occurrence of word w
%       in document d, d, n_z is the number of topics to be discovered
%

% pre-allocate space
[n_d, n_w] = size(n_dw); % max indices of d and w
p_z_d = rand(n_z, n_d); % p(z|d)
p_w_z = rand(n_w, n_z); % p(w|z)
n_p_z_dw = cell(n_z, 1); % n(d,w) * p(z|d,w)
for z = 1:n_z
    n_p_z_dw{z} = sprand(n_dw);
end

p_dw = sprand(n_dw); % p(d,w)
Lt = []; % log-likelihood
for i = 1:iter_num
    %disp('E-step');
    for d = 1:n_d
        for w = find(n_dw(d,:))
            for z = 1:n_z
                n_p_z_dw{z}(d,w) = p_z_d(z,d) * p_w_z(w,z) * ...
                    n_dw(d,w) / p_dw(d, w);
            end
        end
    end

    %disp('M-step');
    %disp('update p(z|d)')
    concat = cat(2, n_p_z_dw{:}); % make n_p_z_dw{:}(d,:)) possible
    for d = 1:n_d
        for z = 1:n_z
            p_z_d(z,d) = sum(n_p_z_dw{z}(d,:));
        end
        p_z_d(:,d) = p_z_d(:,d) / sum(concat(d,:));
    end

    %disp('update p(w|z)')
    for z = 1:n_z
        for w = 1:n_w
            p_w_z(w,z) = sum(n_p_z_dw{z}(:,w));
        end
        p_w_z(:,z) = p_w_z(:,z) / sum(n_p_z_dw{z}(:));
    end

    % update p(d,w) and calculate likelihood
    L = 0;
    for d = 1:n_d
        for w = find(n_dw(d,:))
            p_dw(d,w) = 0;
            for z = 1:n_z
                p_dw(d,w) = p_dw(d,w) + p_w_z(w,z) * p_z_d(z,d);
            end
            L = L + n_dw(d,w) * log(p_dw(d, w));
        end
    end

    Lt = [Lt; L];
    %plot(Lt); ylim([2*median(Lt)-L-0.1 L+(L-median(Lt))/2+0.1]);
    %drawnow; pause(0.1)
end

end

第一次拿 Mablab 写程序,比较丑…… [4]

下图是 Log 似然度随迭代收敛的情况。可以看到收敛速度还是相对较快的。 而且由于是 EM 算法的缘故,Log 似然度确实是单调上升的.

image1

最后,pLSA 的问题是在文档的层面上没有一个概率模型,每篇文档的 P(d|z) 都是需要拟合的模型参数。 这就导致参数的数目会随文档数目线性增长、不能处理训练集外的文档这样的问题。所以02 David Blei、Andrew Ng(就是正在 ml-class.org 里上公开课的那位) 和 Michael Jordan 又提出了一个更为简洁的模型:LDA。有时间的话下次再写了。

[1]这里 Hofmann 自己在 [Hofmann99][Gildea99] 中使用了不同的形式。本文和 Gildea99[Brants05] 一样选择不去理会 \(P(d)\)。因为正如 Brants05 中指出、Blei03 及很多其它文献吐槽的那样,Hofmann99 中的模型算出的 \(P(d)\) 实在坑爹,当 \(d\) 不在训练集中时 \(P(d)\) 就一律为0,没什么意义,还不如别估计它呢。另外 (Hofmann, 1999) 中额外引入了一个参数 \(\beta\) 来“解决”过度拟合问题,但 Brants05 中指出这一问题实际并不存在,因此本文也对此忽略不提。
[2]

具体而言,这里要求的是 \(Q_t(\theta)\)\(\sum_w P(w|z) = 1\)\(\sum_z P(z|d) = 1\) 约束条件下的极值。根据拉格朗日乘数法,解:

\begin{equation*} \nabla_\theta \left( Q(\theta) + \sum_z \alpha_z (\sum_w P(w|z) -1) + \sum_d \beta_d (\sum_z P(z|d) -1) \right) = \mathbf{0} \end{equation*}
[3]完整的程序和数据在这里
[4]吐槽:用 Matlab 做简单字符串处理怎么都那么恶心!长度不同的字符串竟然算是不同类型的!Cell array 怎么那么难用!
[Blei03]Blei, D.M. et al. 2003. Latent Dirichlet Allocation. Journal of Machine Learning Research. 3, 4-5 (2003), 993-1022.
[Hofmann99]Hofmann, T. 1999. Probabilistic latent semantic indexing. Proceedings of the 22nd annual international ACM SIGIR conference on Research and development in information retrieval SIGIR 99. pages, (1999), 50-57.
[Gildea99]Gildea, D. and Hofmann, T. 1999. Topic-based language models using EM. Proceedings of the 6th European Conference on Speech (1999), 2167-2170.
[Brants05]Brants, T. 2005. Test Data Likelihood for PLSA Models. Information Retrieval. (2005), 181-196.