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∈W 在每篇文档 d∈D 中出现的次数 n(w,d)。pLSA 还假设对于每对出现的 (d,w) 都对应着一个表示“主题”的隐藏变量 z∈Z。pLSA 是一个生成模型，它假设 d、w 和 z 之间的关系用贝叶斯网络表示是这样的（从 (Blei et al., 2003) 偷的图）：

实心的节点 d 和 w 表示我们能观察到的文档和单词，空心的 z 表示我们观察不到的隐藏变量，用来表示隐含的主题。图中用了所谓的“盘子记法”，即用方框表示随机变量的重复。这里方框右下角的字母 M 和 N 分别表示有 M 篇文档，第 j 篇文档有 Nj 个单词。每条有向边表示随机变量间的依赖关系。也就是说，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)：利用前者我们可以知道每篇文章中各主题所占的比重，利用后者我们则能知道各单词在各主题中出现的概率，从而进一步找出各主题的“关键词”。记 θ=(P(z|d),P(w|z))，表示我们希望估计的模型参数。当然 θ 不仅仅代表两个数，而是对于每对 (w(j),z(k)) 和 (d(i),z(k)) ，我们都要希望知道 P(z(k)|d(i)) 和 P(w(j)|z(k)) 的值。也就是说，模型中共有 |Z|⋅|D|+|W|⋅|Z| 个参数。我们还知道：

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

接下来 M-step 就是要求 argmaxθQt(θ) 了。利用基本的微积分工具 ^[2]，可以分别对每对 (w(j),z(k))和 (d(i),z(k)) 求出：

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 似然度确实是单调上升的.

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

参考文献

(Hofmann, 1999) 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.
(Gildea & Hofmann, 1999) Gildea, D. and Hofmann, T. 1999. Topic-based language models using EM. Proceedings of the 6th European Conference on Speech (1999), 2167-2170.
(Brants, 2005) Brants, T. 2005. Test Data Likelihood for PLSA Models. Information Retrieval. (2005), 181-196.
(Blei et al., 2003) Blei, D.M. et al. 2003. Latent Dirichlet Allocation. Journal of Machine Learning Research. 3, 4-5 (2003), 993-1022.

1. ^[1] 这里 Hofmann 自己在 (Hofmann, 1999) 和 (Gildea & Hofmann, 1999) 中使用了不同的形式。本文和 (Gildea & Hofmann, 1999)、(Brants, 2005) 一样选择不去理会 P(d)。因为正如 (Brants, 2005) 中指出、(Blei et al., 2003) 及很多其它文献吐槽的那样，(Hofmann, 1999) 中的模型算出的 P(d) 实在坑爹，当 d 不在训练集中时 P(d) 就一律为0，没什么意义，还不如别估计它呢。另外，(Hofmann, 1999) 中额外引入了一个参数 β来“解决”过度拟合问题，但 (Brants, 2005) 中指出这一问题实际并不存在，因此本文也对此忽略不提。 ↩
2. ^[2] 具体而言，这里要求的是 Qt(θ) 在 ∑wP(w|z)=1 和 ∑zP(z|d)=1 约束条件下的极值。根据拉格朗日乘数法，解：

    

    

   ∇θ(Q(θ)+∑zαz(∑wP(w|z)−1)+∑dβd(∑zP(z|d)−1))=0

    

   ↩

3. ^[3] 完整的程序和数据在这里。 ↩
4. ^[4] 吐槽：用 Matlab 做简单字符串处理怎么都那么恶心！长度不同的字符串竟然算是不同类型的！Cell array 怎么那么难用！ ↩