LSA，pLSA原理及其代碼實現

一. LSA

1. LSA原理

LSA(latent semantic analysis)潛在語義分析，也被稱為 LSI(latent semantic index)，是 Scott Deerwester, Susan T. Dumais 等人在 1990 年提出來的一種新的 索引和檢索方法。該方法和傳統向量空間模型(vector space model)一樣使用向量來表示詞(terms)和文檔(documents)，並通過向量間的關系(如夾角)來判斷詞及文檔間的關系；不同的是，LSA 將詞和文檔映射到 潛在語義空間，從而去除了原始向量空間中的一些“噪音”，提高了信息檢索的精確度。

引用吳軍老師在 “矩陣計算與文本處理中的分類問題” 中的總結：

三個矩陣有非常清楚的物理含義。第一個矩陣 U 中的每一行表示意思相關的一類詞，其中的每個非零元素表示這類詞中每個詞的重要性（或者說相關性），數值越大越相關。最后一個矩陣 V 中的每一列表示同一主題一類文章，其中每個元素表示這類文章中每篇文章的相關性。中間的矩陣 D 則表示類詞和文章類之間的相關性。因此，我們只要對關聯矩陣 X 進行一次奇異值分解，我們就可以同時完成了近義詞分類和文章的分類。（同時得到每類文章和每類詞的相關性）。

傳統向量空間模型使用精確的詞匹配，即精確匹配用戶輸入的詞與向量空間中存在的詞，無法解決 一詞多義(polysemy)和 一義多詞(synonymy)的問題。實際上在搜索中，我們實際想要去比較的不是詞，而是 隱藏在詞之后的意義和概念。

LSA 的核心思想是將詞和文檔映射到 潛在語義空間，再比較其相似性。

舉個簡單的栗子，對一個 Term-Document 矩陣做SVD分解，並將左奇異向量和右奇異向量都取后2維（之前是3維的矩陣），投影到一個平面上（潛在語義空間），可以得到：

在圖上，每一個紅色的點，都表示一個詞，每一個藍色的點，都表示一篇文檔，這樣我們可以對這些詞和文檔進行聚類，比如說 stock 和 market 可以放在一類，因為他們老是出現在一起，real 和 estate 可以放在一類，dads，guide 這種詞就看起來有點孤立了，我們就不對他們進行合並了。按這樣聚類出現的效果，可以提取文檔集合中的近義詞，這樣當用戶檢索文檔的時候，是用語義級別（近義詞集合）去檢索了，而不是之前的詞的級別。這樣一減少我們的檢索、存儲量，因為這樣壓縮的文檔集合和PCA是異曲同工的，二可以提高我們的用戶體驗，用戶輸入一個詞，我們可以在這個詞的近義詞的集合中去找，這是傳統的索引無法做到的。

2. LSA的優點

1）低維空間表示可以刻畫同義詞，同義詞會對應着相同或相似的主題。

2）降維可去除部分噪聲，是特征更魯棒。

3）充分利用冗余數據。

4）無監督/完全自動化。

5）與語言無關。

3. LSA的缺點

1）LSA可以處理向量空間模型無法解決的一義多詞(synonymy)問題，但不能解決 一詞多義(polysemy)問題。因為LSA將每一個詞映射為潛在語義空間中的一個點，也就是說一個詞的多個意思在空間中對於的是同一個點，並沒有被區分。

2）SVD的優化目標基於L-2 norm 或者 Frobenius Norm  的，這相當於隱含了對數據的高斯分布假設。而 term 出現的次數是非負的，這明顯不符合 Gaussian 假設，而更接近 Multi-nomial 分布。

3）特征向量的方向沒有對應的物理解釋。

4）SVD的計算復雜度很高，而且當有新的文檔來到時，若要更新模型需重新訓練。

5）沒有刻畫term出現次數的概率模型。

6）對於count vectors 而言，歐式距離表達是不合適的（重建時會產生負數）。

7）維數的選擇是ad-hoc的。

 

二. pLSA

首先，我們可以看看日常生活中人是如何構思文章的。如果我們要寫一篇文章，往往是先確定要寫哪幾個主題。譬如構思一篇自然語言處理相關的文章，可能40%會談論語言學，30%談論概率統計，20%談論計算機，還有10%談論其它主題。

對於語言學，容易想到的詞包括：語法，句子，主語等；對於概率統計，容易想到的詞包括：概率，模型，均值等；對於計算機，容易想到的詞包括：內存，硬盤，編程等。我們之所以能想到這些詞，是因為這些詞在對應的主題下出現的概率很高。我們可以很自然的看到，一篇文章通常是由多個主題構成的，而每一個主題大概可以用與該主題相關的頻率最高的一些詞來描述。以上這種想法由Hofmann於1999年給出的pLSA模型中首先進行了明確的數學化。Hofmann認為一篇文章（Doc）可以由多個主題（Topic）混合而成，而每個Topic都是詞匯上的概率分布，文章中的每個詞都是由一個固定的Topic生成的。下圖是英語中幾個Topic的例子。

 

pLSA的建模思路分為兩種。

1. 第一種思路

以的概率從文檔集合中選擇一個文檔

以的概率從主題集合中選擇一個主題

以的概率從詞集中選擇一個詞

有幾點說明：

• 以上變量有兩種狀態：observed ( & ) 和 latent ()
• 來自文檔，但同時是集合（元素不重復），相當於一個詞匯表

直接的，針對observed variables做建立likelihood function：

其中，為pair出現的次數。為加以區分，之后使用與標識對應文檔與詞匯數量。兩邊取，得：

其中，倒數第二步旨在將暴露出來。由於likelihood function中與存在latent variable，難以直接使用MLE求解，很自然想到用E-M算法求解。E-M算法主要分為Expectation與Maximization兩步。

  Step 1: Expectation

假設已知與，求latent variable的后驗概率

Step 2: Maximization

求關於參數和的Complete data對數似然函數期望的極大值，得到最優解。帶入E步迭代循環。

由式可得：

此式后部分為常量。故令：

建立以下目標函數與約束條件：

只有等式約束，使用Lagrange乘子法解決：

對與求駐點，得：

令，得：

，故有：

同理，有：

將與回代Expectation：

，循環迭代。

pLSA的建模思想較為簡單，對於observed variables建立likelihood function，將latent variable暴露出來，並使用E-M算法求解。其中M步的標准做法是引入Lagrange乘子求解后回代到E步。

 

總結一下使用EM算法求解pLSA的基本實現方法：

(1)E步驟：求隱含變量Given當前估計的參數條件下的后驗概率。

(2)M步驟：最大化 Complete data對數似然函數的期望，此時我們使用E步驟里計算的隱含變量的后驗概率，得到新的參數值。

兩步迭代進行直到收斂。

2. 第二種思路

這個思路和上面思路的區別就在於對P(d,w)的展開公式使用的不同，思路二使用的是3個概率來展開的，如下：

這樣子我們后面的EM算法的大致思路都是相同的，就是表達形式變化了，最后得到的EM步驟的更新公式也變化了。當然，思路二得到的是3個參數的更新公式。如下：

 

 你會發現，其實多了一個參數是P(z)，參數P（d|z）變化了（之前是P(z|d)），然后P(w|z)是不變的，計算公式也相同。

給定一個文檔d，我們可以將其分類到一些主題詞類別下。

PLSA算法可以通過訓練樣本的學習得到三個概率，而對於一個測試樣本，其中P(w|z)概率是不變的，但是P(z)和P(d|z)都是需要重新更新的，我們也可以使用上面的EM算法，假如測試樣本d的數據，我們可以得到新學習的P(z)和P(d|z)參數。這樣我們就可以計算：

為什么要計算P(z|d)呢？因為給定了一個測試樣本d，要判斷它是屬於那些主題的，我們就需要計算P(z|d)，就是給定d，其在主題z下成立的概率是多少，不就是要計算嗎。這樣我們就可以計算文檔d在所有主題下的概率了。

這樣既可以把一個測試文檔划歸到一個主題下，也可以給訓練文檔打上主題的標記，因為我們也是可以計算訓練文檔它們的的。如果從這個應用思路來說，思路一說似乎更加直接，因為其直接計算出來了。

3. pLSA的優勢

1）定義了概率模型，而且每個變量以及相應的概率分布和條件概率分布都有明確的物理解釋。

2）相比於LSA隱含了高斯分布假設，pLSA隱含的Multi-nomial分布假設更符合文本特性。

3）pLSA的優化目標是是KL-divergence最小，而不是依賴於最小均方誤差等准則。

4）可以利用各種model selection和complexity control准則來確定topic的維數。

4. pLSA的不足

 

1）概率模型不夠完備：在document層面上沒有提供合適的概率模型，使得pLSA並不是完備的生成式模型，而必須在確定document i的情況下才能對模型進行隨機抽樣。

2）隨着document和term 個數的增加，pLSA模型也線性增加，變得越來越龐大。

3）EM算法需要反復的迭代，需要很大計算量。

 

針對pLSA的不足，研究者們又提出了各種各樣的topic based model, 其中包括大名鼎鼎的Latent Dirichlet Allocation (LDA)。

 

三. pLSA的Python代碼實現

1. preprocess.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
class Preprocess:
    def __init__(self, fname, fsw):
        self.fname = fname
        # doc info
        self.docs = []
        self.doc_size = 0
        # stop word info
        self.sws = []
        # word info
        self.w2id = {}
        self.id2w = {}
        self.w_size = 0
        # stop word list init
        with open(fsw, 'r') as f:
            for line in f:
                self.sws.append(line.strip())
    def __work(self):
        with open(self.fname, 'r') as f:
            for line in f:
                line_strip = line.strip()
                self.doc_size += 1
                self.docs.append(line_strip)
                items = line_strip.split()
                for it in items:
                    if it not in self.sws:
                        if it not in self.w2id:
                            self.w2id[it] = self.w_size
                            self.id2w[self.w_size] = it
                            self.w_size += 1
        self.w_d = np.zeros([self.w_size, self.doc_size], dtype=np.int)
        for did, doc in enumerate(self.docs):
            ws = doc.split()
            for w in ws:
                if w in self.w2id:
                    self.w_d[self.w2id[w]][did] += 1
    def get_w_d(self):
        self.__work()
        return self.w_d
    def get_word(self, wid):
        return self.id2w[wid]
if __name__ == '__main__':
    fname = './data.txt'
    fsw = './stopwords.txt'
    pp = Preprocess(fname, fsw)

2. plsa.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
import time
import logging
def normalize(vec):
    s = sum(vec)
    for i in range(len(vec)):
        vec[i] = vec[i] * 1.0 / s
def llhood(w_d, p_z, p_w_z, p_d_z):
    V, D = w_d.shape
    ret = 0.0
    for w, d in zip(*w_d.nonzero()):
        p_d_w = np.sum(p_z * p_w_z[w,:] * p_d_z[d,:])
        if p_d_w > 0:
            ret += w_d[w][d] * np.log(p_d_w)
    return ret
class PLSA:
    def __init__(self):
        pass
    def train(self, w_d, Z, eps):
        V, D = w_d.shape
        # create prob array, p(d|z), p(w|z), p(z)
        p_d_z = np.zeros([D, Z], dtype=np.float)
        p_w_z = np.zeros([D, Z], dtype=np.float)
        p_z = np.zeros([Z], dtype=np.float)
        # initialize
        p_d_z = np.random.random([D, Z])
        for d_idx in range(D):
            normalize(p_d_z[d_idx])
        p_w_z = np.random.random([V, Z])
        for w_idx in range(V):
            normalize(p_w_z[w_idx])
        p_z = np.random.random([Z])
        normalize(p_z)
        # iteration until converge
        step = 1
        pp_d_z = p_d_z.copy()
        pp_w_z = p_w_z.copy()
        pp_z = p_z.copy()
        while True:
            logging.info('[ iteration ] step %d' % step)
            step += 1
            p_d_z *= 0.0
            p_w_z *= 0.0
            p_z *= 0.0
            # run EM algorithm
            for w_idx, d_idx in zip(*w_d.nonzero()):
                #print '[ EM ] >>>>>> E step : '
                p_z_d_w = pp_z * pp_d_z[d_idx,:] * pp_w_z[w_idx,:]
                normalize(p_z_d_w)
                #print '[ EM ] >>>>>> M step : '
                tt = w_d[w_idx, d_idx] * p_z_d_w
                p_w_z[w_idx,:] += tt
                p_d_z[d_idx,:] += tt
                p_z += tt
            normalize(p_w_z)
            normalize(p_d_z)
            p_z = p_z / w_d.sum()
            # check converge
            l1 = llhood(w_d, pp_z, pp_w_z, pp_d_z)
            l2 = llhood(w_d, p_z, p_w_z, p_d_z)
            diff = l2 - l1
            logging.info('[ iteration ] l2-l1  %.3f - %.3f = %.3f ' % (l2, l1, diff))
            if abs(diff) < eps:
                logging.info('[ iteration ] End EM ')
                return (l2, p_d_z, p_w_z, p_z)
            pp_d_z = p_d_z.copy()
            pp_w_z = p_w_z.copy()
            pp_z = p_z.copy()

3. main.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from preprocess import Preprocess as PP
from plsa import PLSA
import numpy as np
import logging
import time
def main():
    # setup logging --------------------------
    logging.basicConfig(filename='plsa.log',
                        level=logging.INFO,
                        format='%(asctime)s %(filename)s[line:%(lineno)d] %(levelname)s %(message)s',
                        datefmt='%a, %d %b %Y %H:%M:%S')
    #console = logging.StreamHandler()
    #console.setLevel(logging.INFO)
    #logging.getLogger('').addHandler(console)
    # some basic configuration ---------------
    fname = './data.txt'
    fsw = './stopwords.txt'
    eps = 20.0
    key_word_size = 10
    # preprocess -----------------------------
    pp = PP(fname, fsw)
    w_d = pp.get_w_d()
    V, D = w_d.shape
    logging.info('V = %d, D = %d' % (V, D))
    # train model and get result -------------
    pmodel = PLSA()
    for z in range(3, (D+1), 10):
        t1 = time.clock()
        (l, p_d_z, p_w_z, p_z) = pmodel.train(w_d, z, eps)
        t2 = time.clock()
        logging.info('z = %d, eps = %f, time = %f' % (z, l, t2-t1))
        for itz in range(z):
            logging.info('Topic %d' % itz)
            data = [(p_w_z[i][itz], i) for i in range(len(p_w_z[:,itz]))]
            data.sort(key=lambda tup:tup[0], reverse=True)
            for i in range(key_word_size):
                logging.info('%s : %.6f ' % (pp.get_word(data[i][1]), data[i][0]))
if __name__ == '__main__':
    main()

 

版權聲明：

   本文由笨兔勿應所有，發布於http://www.cnblogs.com/bentuwuying。如果轉載，請注明出處，在未經作者同意下將本文用於商業用途，將追究其法律責任。