機率語言模型及其變形系列(1)-PLSA及EM演算法_em演算法

本系列博文介紹常見機率語言模型及其變形模型，主要總結PLSA、LDA及LDA的變形模型及參數Inference方法。初步計劃內容如下

第一篇：PLSA及EM演算法

第二篇：LDA及Gibbs Samping

第三篇：LDA變形模型-Twitter LDA，TimeUserLDA，ATM，Labeled-LDA，MaxEnt-LDA等

第四篇：基於變形LDA的paper分類總結

第五篇：LDA Gibbs Sampling的JAVA實現

第一篇 PLSA及EM演算法

[本文PDF版本下載地址 PLSA及EM演算法-yangliuy]

前言：本文主要介紹PLSA及EM演算法，首先給出LSA（隱性語義分析）的早期方法SVD，然後引入基於機率的PLSA模型，其參數學習採用EM演算法。接著我們分析如何運用EM演算法估計一個簡單的mixture unigram 語言模型和混合高斯模型GMM的參數，最後總結EM演算法的一般形式及運用關鍵點。對於改進PLSA，引入hyperparameter的LDA模型及其Gibbs Sampling參數估計方法放在本系列後面的文章LDA及Gibbs Samping介紹。

1 LSA and SVD

LSA(隱性語義分析)的目的是要從文本中發現隱含的語義維度-即“Topic”或者“Concept”。我們知道，在文檔的空間向量模型（VSM）中，文檔被表示成由特徵詞出現機率組成的多維向量，這種方法的好處是可以將query和文檔轉化成同一空間下的向量計算相似性，可以對不同詞項賦予不同的權重，在文本檢索、分類、聚類問題中都得到了廣泛應用，在基於貝葉斯演算法及KNN演算法的newsgroup18828文本分類器的JAVA實現和基於Kmeans演算法、MBSAS演算法及DBSCAN演算法的newsgroup18828文本聚類器的JAVA實現系列文章中的分類聚類演算法大多都是採用向量空間模型。然而，向量空間模型沒有能力處理一詞多義和一義多詞問題，例如同義字也分別被表示成獨立的一維，計算向量的餘弦相似性時會低估使用者期望的相似性；而某個詞項有多個詞義時，始終對應同一維度，因此計算的結果會高估使用者期望的相似性。

LSA方法的引入就可以減輕類似的問題。基於SVD分解，我們可以構造一個原始向量矩陣的一個低秩逼近矩陣，具體的做法是將詞項文檔矩陣做SVD分解

  其中是以詞項(terms)為行, 文檔(documents)為列做一個大矩陣. 設一共有t行d列,  矩陣的元素為詞項的tf-idf值。然後把的r個對角元素的前k個保留（最大的k個保留）, 後面最小的r-k個奇異值置0, 得到；最後計算一個近似的分解矩陣

則在最小二乘意義下是的最佳逼近。由於最多包含k個非零元素，所以的秩不超過k。通過在SVD分解近似，我們將原始的向量轉化成一個低維隱含語義空間中，起到了特徵降維的作用。每個奇異值對應的是每個“語義”維度權重，將不太重要的權重設為0，只保留最重要的維度資訊，去掉一些資訊“nosie”,因而可以得到文檔的一種更優表示形式。將SVD分解降維應用到文檔聚類的JAVA實現可參見此文。

IIR中給出的一個SVD降維的執行個體如下：

左邊是原始矩陣的SVD分解，右邊是只保留權重最大2維，將原始矩陣降到2維後的情況。

2 PLSA

儘管基於SVD的LSA取得了一定的成功，但是其缺乏嚴謹的數理統計基礎，而且SVD分解非常耗時。Hofmann在SIGIR'99上提出了基於機率統計的PLSA模型，並且用EM演算法學習模型參數。PLSA的機率圖模型如下

其中D代表文檔，Z代表隱含類別或者主題，W為觀察到的單詞，表示單詞出現在文檔的機率，表示文檔中出現主題下的單詞的機率，給定主題出現單詞的機率。並且每個主題在所有詞項上服從Multinomial 分布，每個文檔在所有主題上服從Multinomial 分布。整個文檔的產生過程是這樣的：

(1) 以的機率選中文檔；

(2) 以的機率選中主題；

(3) 以的機率產生一個單詞。

我們可以觀察到的資料就是對，而是隱含變數。的聯合分布為

而和分布對應了兩組Multinomial 分布，我們需要估計這兩組分布的參數。下面給出用EM演算法估計PLSA參數的詳細推導過程。

3 Estimate parameters in PLSA  by EM

（註：本部分主要參考Tomas Hoffman, Unsupervised Learning by Probabilistic Latent Semantic Analysis.）

如文本語言模型的參數估計-最大似然估計、MAP及貝葉斯估計一文所述，常用的參數估計方法有MLE、MAP、貝葉斯估計等等。但是在PLSA中，如果我們試圖直接用MLE來估計參數，就會得到似然函數

其中是term 出現在文檔中的次數。注意這是一個關於和的函數，一共有N*K + M*K個自變數(注意這裡M表示term的總數，一般文獻習慣用V表示)，如果直接對這些自變數求偏導數，我們會發現由於自變數包含在對數和中，這個方程的求解很困難。因此對於這樣的包含“隱含變數”或者“缺失資料”的機率模型參數估計問題，我們採用EM演算法。

EM演算法的步驟是：

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

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

兩步迭代進行直到收斂。

先解釋一下什麼是Incomplete data和complete data。Zhai老師在一篇經典的EM演算法Notes中講到，當未經處理資料的似然函數很複雜時，我們通過增加一些隱含變數來增強我們的資料，得到“complete data”,而“complete data”的似然函數更加簡單，方便求極大值。於是，原始的資料就成了“incomplete data”。我們將會看到，我們可以通過最大化“complete data”似然函數的期望來最大化"incomplete data"的似然函數，以便得到求似然函數最大值更為簡單的計算途徑。

針對我們PLSA參數估計問題，在E步驟中，直接使用貝葉斯公式計算隱含變數在當前參數取值條件下的後驗機率，有

在這個步驟中，我們假定所有的和都是已知的，因為初始時隨機賦值，後面迭代的過程中取前一輪M步驟中得到的參數值。

在M步驟中，我們最大化Complete data對數似然函數的期望。在PLSA中，Incomplete data 是觀察到的，隱含變數是主題，那麼complete data就是三元組，其期望是

注意這裡是已知的，取的是前面E步驟裡面的估計值。下面我們來最大化期望，這又是一個多元函數求極值的問題，可以用拉格朗日乘數法。拉格朗日乘數法可以把條件極值問題轉化為無條件極值問題，在PLSA中目標函數就是,約束條件是

由此我們可以寫出拉格朗日函數

這是一個關於和的函數，分別對其求偏導數，我們可以得到

注意這裡進行過方程兩邊同時乘以和的變形，聯立上面4組方程，我們就可以解出M步驟中通過最大化期望估計出的新的參數值

解方程組的關鍵在於先求出,其實只需要做一個加和運算就可以把的係數都化成1，後面就好計算了。

然後使用更新後的參數值，我們又進入E步驟，計算隱含變數 Given當前估計的參數條件下的後驗機率。如此不斷迭代，直到滿足終止條件。

注意到我們在M步驟中還是使用對Complete Data的MLE，那麼如果我們想加入一些先驗知識進入我們的模型，我們可以在M步驟中使用MAP估計。正如文本語言模型的參數估計-最大似然估計、MAP及貝葉斯估計中投硬幣的二項分布實驗中我們加入“硬幣一般是兩面均勻的”這個先驗一樣。而由此計算出的參數的估計值會在分子分母中多出關於先驗參數的preduo counts,其他步驟都是一樣的。具體可以參考Mei Qiaozhu 的Notes。

PLSA的實現也不難，網上有很多實現code。

例如這個PLSA的EM演算法實現 http://ezcodesample.com/plsaidiots/PLSAjava.txt

主要的類如下(作者Andrew Polar)

//The code is taken from://http://code.google.com/p/mltool4j/source/browse/trunk/src/edu/thu/mltool4j/topicmodel/plsa//I noticed some difference with original Hofmann concept in computation of P(z). It is //always even and actually not involved, that makes this algorithm non-negative matrix //factoring and not PLSA.//Found and tested by Andrew Polar. //My version can be found on semanticsearchart.com or ezcodesample.com

class ProbabilisticLSA{private Dataset dataset = null;    private Posting[][] invertedIndex = null;    private int M = -1; // number of data    private int V = -1; // number of words    private int K = -1; // number of topics    public ProbabilisticLSA()    {    }    public boolean doPLSA(String datafileName, int ntopics, int iters)    {    File datafile = new File(datafileName);        if (datafile.exists())        {        if ((this.dataset = new Dataset(datafile)) == null)            {        System.out.println("doPLSA, dataset == null");                return false;            }            this.M = this.dataset.size();            this.V = this.dataset.getFeatureNum();            this.K = ntopics;                         //build inverted index            this.buildInvertedIndex(this.dataset);            //run EM algorithm            this.EM(iters);            return true;                    }        else        {        System.out.println("ProbabilisticLSA(String datafileName), datafile: " + datafileName + " doesn't exist");            return false;        }    }    //Build the inverted index for M-step fast calculation. Format:    //invertedIndex[w][]: a unsorted list of document and position which word w    // occurs.     //@param ds the dataset which to be analysis    @SuppressWarnings("unchecked")    private boolean buildInvertedIndex(Dataset ds)    {    ArrayList<Posting>[] list = new ArrayList[this.V];        for (int k=0; k<this.V; ++k) {        list[k] = new ArrayList<Posting>();        }                for (int m = 0; m < this.M; m++)        {        Data d = ds.getDataAt(m);         for (int position = 0; position < d.size(); position++)         {        int w = d.getFeatureAt(position).dim;         // add posting         list[w].add(new Posting(m, position));        }        }        // convert to array        this.invertedIndex = new Posting[this.V][];        for (int w = 0; w < this.V; w++)        {        this.invertedIndex[w] = list[w].toArray(new Posting[0]);        }        return true;    }    private boolean EM(int iters)    {    // p(z), size: K        double[] Pz = new double[this.K];        // p(d|z), size: K x M        double[][] Pd_z = new double[this.K][this.M];        // p(w|z), size: K x V        double[][] Pw_z = new double[this.K][this.V];        // p(z|d,w), size: K x M x doc.size()        double[][][] Pz_dw = new double[this.K][this.M][];         // L: log-likelihood value         double L = -1;         // run EM algorithm         this.init(Pz, Pd_z, Pw_z, Pz_dw);         for (int it = 0; it < iters; it++)         {          // E-step         if (!this.Estep(Pz, Pd_z, Pw_z, Pz_dw))         {         System.out.println("EM,  in E-step");             }             // M-step             if (!this.Mstep(Pz_dw, Pw_z, Pd_z, Pz))             {             System.out.println("EM, in M-step");             }             L = calcLoglikelihood(Pz, Pd_z, Pw_z);             System.out.println("[" + it + "]" + "\tlikelihood: " + L);         }                         //print result         for (int m = 0; m < this.M; m++)         {         double norm = 0.0;         for (int z = 0; z < this.K; z++) {         norm += Pd_z[z][m];             }         if (norm <= 0.0) norm = 1.0;         for (int z = 0; z < this.K; z++) {         System.out.format("%10.4f", Pd_z[z][m]/norm);             }             System.out.println();        }         return false;    }       private boolean init(double[] Pz, double[][] Pd_z, double[][] Pw_z, double[][][] Pz_dw)    {    // p(z), size: K        double zvalue = (double) 1 / (double) this.K;        for (int z = 0; z < this.K; z++)        {        Pz[z] = zvalue;        }        // p(d|z), size: K x M        for (int z = 0; z < this.K; z++)        {        double norm = 0.0;            for (int m = 0; m < this.M; m++)            {            Pd_z[z][m] = Math.random();                norm += Pd_z[z][m];            }            for (int m = 0; m < this.M; m++)            {            Pd_z[z][m] /= norm;            }        }        // p(w|z), size: K x V        for (int z = 0; z < this.K; z++)        {        double norm = 0.0;            for (int w = 0; w < this.V; w++)            {            Pw_z[z][w] = Math.random();                norm += Pw_z[z][w];            }            for (int w = 0; w < this.V; w++)            {            Pw_z[z][w] /= norm;            }        }        // p(z|d,w), size: K x M x doc.size()        for (int z = 0; z < this.K; z++)        {        for (int m = 0; m < this.M; m++)            {         Pz_dw[z][m] = new double[this.dataset.getDataAt(m).size()];            }        }        return false;    }    private boolean Estep(double[] Pz, double[][] Pd_z, double[][] Pw_z, double[][][] Pz_dw)    {    for (int m = 0; m < this.M; m++)        {    Data data = this.dataset.getDataAt(m);            for (int position = 0; position < data.size(); position++)            {            // get word(dimension) at current position of document m                int w = data.getFeatureAt(position).dim;                double norm = 0.0;                for (int z = 0; z < this.K; z++)                {                double val = Pz[z] * Pd_z[z][m] * Pw_z[z][w];                    Pz_dw[z][m][position] = val;                    norm += val;                }                // normalization                for (int z = 0; z < this.K; z++)                {                Pz_dw[z][m][position] /= norm;                }            }        }        return true;    }    private boolean Mstep(double[][][] Pz_dw, double[][] Pw_z, double[][] Pd_z, double[] Pz)    {    // p(w|z)        for (int z = 0; z < this.K; z++)        {        double norm = 0.0;            for (int w = 0; w < this.V; w++)            {            double sum = 0.0;                Posting[] postings = this.invertedIndex[w];                for (Posting posting : postings)                {                int m = posting.docID;                    int position = posting.pos;                    double n = this.dataset.getDataAt(m).getFeatureAt(position).weight;                    sum += n * Pz_dw[z][m][position];                }                Pw_z[z][w] = sum;                norm += sum;            }            // normalization            for (int w = 0; w < this.V; w++)            {            Pw_z[z][w] /= norm;            }        }        // p(d|z)        for (int z = 0; z < this.K; z++)        {        double norm = 0.0;            for (int m = 0; m < this.M; m++)            {            double sum = 0.0;                Data d = this.dataset.getDataAt(m);                for (int position = 0; position < d.size(); position++)                {                double n = d.getFeatureAt(position).weight;                    sum += n * Pz_dw[z][m][position];                }                Pd_z[z][m] = sum;                norm += sum;            }            // normalization            for (int m = 0; m < this.M; m++)            {            Pd_z[z][m] /= norm;            }        }        //This is definitely a bug        //p(z) values are even, but they should not be even        double norm = 0.0;        for (int z = 0; z < this.K; z++)        {        double sum = 0.0;            for (int m = 0; m < this.M; m++)            {            sum += Pd_z[z][m];            }            Pz[z] = sum;            norm += sum;       }        // normalization        for (int z = 0; z < this.K; z++)        {        Pz[z] /= norm;        //System.out.format("%10.4f", Pz[z]);  //here you can print to see        }        //System.out.println();        return true;    }    private double calcLoglikelihood(double[] Pz, double[][] Pd_z, double[][] Pw_z)    {    double L = 0.0;        for (int m = 0; m < this.M; m++)        {        Data d = this.dataset.getDataAt(m);            for (int position = 0; position < d.size(); position++)            {            Feature f = d.getFeatureAt(position);                int w = f.dim;                double n = f.weight;                double sum = 0.0;                for (int z = 0; z < this.K; z++)                {                sum += Pz[z] * Pd_z[z][m] * Pw_z[z][w];                }                L += n * Math.log10(sum);            }        }        return L;    }}public class PLSA {public static void main(String[] args) {ProbabilisticLSA plsa = new ProbabilisticLSA();//the file is not used, the hard coded data is used instead, but file name should be valid,//just replace the name by something valid.plsa.doPLSA("C:\\Users\\APolar\\workspace\\PLSA\\src\\data.txt", 2, 60);        System.out.println("end PLSA");    }}

4 Estimate parameters in a simple mixture unigram language model by EM

在PLSA的參數估計中，我們使用了EM演算法。EM演算法經常用來估計包含“缺失資料”或者“隱含變數”模型的參數估計問題。這兩個概念是互相聯絡的，當我們的模型中有“隱含變數”時，我們會認為未經處理資料是“不完全的資料”，因為隱含變數的值無法觀察到；反過來，當我們的資料incomplete時，我們可以通過增加隱含變數來對“缺失資料”建模。

為了加深對EM演算法的理解，下面我們來看如何用EM演算法來估計一個簡單混合unigram語言模型的參數。本部