PLSA

PLSA模型

PLSA和LDA很像，都属于主题模型，即它们都认为上帝在写文章时先以一定概率选择了一个主题，然后在这主题下以一定概率选择了一个词，重复这个过程就完成了一篇文章，即$p(d_i,w_j)=p(z_k|d_i)p(w_j|z_k)$，其中$d$表示文章，$w$表示词，$z$表示主题。

模型求解

模型求解即求出所有的$p(z_k|d_i)$和$p(w_j|z_k)$，这样就可以生成任意篇文章了。

这里有有必要补充个基础概念--条件概率和后验概率。所谓条件概率就是“由因得果”，在$p(z_k|d_i)$中$d_i$是因，$z_k$是果，所以$p(z_k|d_i)$就是个条件概率，同样$p(w_j|z_k)$也是条件概率。所谓后验概率就是“执果寻因”，即观察到了系统的输出和输出，去探寻系统内部的运作机理。对应到PLSA模型中就是观察到$d_i$中出现了一个词$w_j$（文档和词都是观察变量），去探寻连接$d_i$和$w_j$的是哪个主题（主题是隐含变量），如下图所示，其实就是求$p(z_k|d_i,w_j)$。

图1. PLSA的概率图模型

下面用EM算法求解模型参数$p(z_k|d_i)$和$p(w_j|z_k)$。

E-Step

E是Expection(期望)的意思，即根据上一轮得到的模型参数求隐含变量的期望，对应到PLSA模型中就是根据上轮得到的模型参数$p(z_k|d_i)$和$p(w_j|z_k)$计算每篇文档中每个词背后对应的主题的概率$p(z_k|d_i,w_j)$。回头看看图1，从$d_i$到$w_j$一共有$K$条路径，途经$z_k$的概率为

\begin{equation}p(z_k|d_i,w_j)=\frac{p(z_k|d_i)p(w_j|z_k)}{\sum_k{p(z_k|d_i)p(w_j|z_k)}}\label{post1}\end{equation}

 这里的条件概率$p(z_k|d_i)$和$p(w_j|z_k)$是由上一轮的M-Step得到的，初始时$p(z_k|d_i)$和$p(w_j|z_k)$由随机赋值得到。

如果完全由贝叶斯公式推导是这样的

\begin{equation}p(z_k|d_i,w_j)=\frac{p(z_k)p(d_i,w_j|z_k)}{p(d_i,w_j)}=\frac{p(z_k)p(z_k|d_i)p(w_j|z_k)}{\sum_k{p(z_k|d_i)p(w_j|z_k)}}\label{post2}\end{equation}

跟公式\ref{post1}相比，公式\ref{post2}分子中多了个$p(z_k)$。于是计算$p(z_k|d_i,w_j)$就出现了两个不同的版本，两种版本的代码我都见过，但是PLSA原创作者使用的是公式\ref{post1}。

M-Step

M是极大似然估计（Maximum Likelihood Estimate，MLE）的意思，在已知后验概率的情况下通过MLE的方法求条件概率。

当我们已知所有的 $p(z_k|d_i,w_j)$时，统计一下在所有文章中由$z_k$到$w_j$的次数，再统计一下在所有文章中由$z_k$到任意$w$的次数，两个次数相除就得到了$p(w_j|z_k)$

\begin{equation}p(w_j|z_k)=\frac{\sum_i{p(z_k|d_i,w_j)}}{\sum_i{\sum_j{p(z_k|d_i,w_j)}}}\label{cond1}\end{equation}

同样，统计一下在文章$d_i$当中主题$z_k$出现的次数，再统计一下文章$d_i$中所有主题$z$的出现次数，两者相除就得到了$p(z_k|d_i)$

 \begin{equation}p(z_k|d_i)=\frac{\sum_j{p(z_k|d_i,w_j)}}{\sum_j{\sum_k{p(z_k|d_i,w_j)}}}\label{cond2}\end{equation}

且慢，不是说M-Step是用MLE的方法求条件概率吗？这种简单地统计频数，让两个频数相除跟MLE有什么关系呢？其实频数相除就是由MLE推导出来的，我们举一个简单的例子来证明MLE和频数相除是等价的（直接拿PLSA的例子来证明会比较复杂，中间还牵涉到拉格朗日数乘法）。投了10次硬币，6次正面向上，4次反面向上，问这枚硬币正面向上的概率是多少。用频数相除的方法可以很容易地得到正面向上的概率是$\frac{6}{10}$。如果是用MLE求解，先设正面向上的概率为$p$，则似然函数为$p^6{(1-p)^4}$，对数似然函数为$ln{p^6}+ln{(1-p)^4}=6ln{p}+4ln{(1-p)}$，为求对数似然函数的极大值点我们令其导数为0，$\frac{6}{p}-\frac{4}{1-p}=0$，得$p=\frac{6}{10}$。所以两种方法等价。

公式\ref{cond1}和\ref{cond2}没有考虑到一个词出现在文章的不同位置其权重实际上是不一样的，比如一个词出现在正文里我们算作1次出现，如果出现在标题里就应该算作1.5次出现。于是改进后的条件概率计算公式为

\begin{equation}p(w_j|z_k)=\frac{\sum_i{weight_{ij}\cdot p(z_k|d_i,w_j)}}{\sum_i{\sum_j{weight_{ij}\cdot p(z_k|d_i,w_j)}}}\label{cond3}\end{equation}

\begin{equation}p(z_k|d_i)=\frac{\sum_j{weight_{ij}\cdot p(z_k|d_i,w_j)}}{\sum_j{\sum_k{weight_{ij}\cdot p(z_k|d_i,w_j)}}}\label{cond4}\end{equation}

$weight_{ij}$是$w_j$在$d_i$中的权重。

PLSA用于推荐

PLSA是一个词袋模型（BOW, Bag Of Word），它不考虑词在文档中出现的顺序，但可以把词在文档中的权重考虑进来。我们把这些概念平行推广到推荐系统中来，一个用户的购买记录看作是一个文档，购买的每一件商品看作是一个词，用户对商品的评分看作是词在文档中的权重。套用PLSA算法就可以得到用户在各个隐含主题上的向量表示$p(z_k|d_i)$，基于这个向量再去计算相似用户，接着套用协同过滤算法给用户推荐商品。

Java实现

PLSA.java

package plsa;

import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.File;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;

/**
 * 最初的代码来自于https://code.google.com/archive/p/mltool4j/，源代码在计算p(z|d,w)时使用了p(z)，但是在传统的PLSA算法中p(z)根本就没有出现过，所以我对源代码做了改动。
 * 
 * @author orisun
 * @date 2016年7月13日
 */
public class PLSA {
    private Dataset dataset = null;
    private Posting[][] invertedIndex = null;
    private int M = -1; //文档数
    private int V = -1; //词汇数
    private int K = -1; //主题数

    public boolean doPLSA(String datafilePath, int ntopics, int iters) {
        try {
            this.dataset = new Dataset(datafilePath);
        } catch (IOException e) {
            e.printStackTrace();
            return false;
        }
        this.M = this.dataset.size();
        this.V = this.dataset.getFeatureNum();
        this.K = ntopics;

        //建立term-->doc的倒排索引，在计算p(w|z)时可以提高速度
        this.buildInvertedIndex(this.dataset);
        this.runEM(iters);
        return true;
    }

    /**
     * 建立term-->doc的倒排索引，在计算p(w|z)时可以提高速度
     * @param ds
     * @return
     */
    @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;
                list[w].add(new Posting(m, position));
            }
        }
        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 runEM(int iters) {
        // p(z|d), size: M x K
        double[][] Pz_d = new double[this.M][this.K];

        // p(w|z), size: K x V
        double[][] Pw_z = new double[this.K][this.V];

        // p(z|d,w), size: M x K x doc.size()
        double[][][] Pz_dw = new double[this.M][this.K][];

        // L: log-likelihood value
        double L = -1;

        //初始时，随机初始化参数
        this.init(Pz_d, Pw_z, Pz_dw);
        for (int it = 0; it < iters; it++) {
            System.out.println("iteration " + it);
            // E-step
            if (!this.Estep(Pz_d, Pw_z, Pz_dw)) {
                System.out.println("EM,  in E-step");
            }

            // M-step
            if (!this.Mstep(Pz_dw, Pw_z, Pz_d)) {
                System.out.println("EM, in M-step");
            }

            File modelPath = new File("model");
            if (modelPath.exists()) {
                if (modelPath.isFile()) {
                    modelPath.delete();
                    modelPath.mkdirs();
                }
            } else {
                modelPath.mkdirs();
            }
            //进入最后几轮迭代时，保存参数
            if (it > iters - 10) {
                L = calcLoglikelihood(Pz_d, Pw_z);
                System.out.println("[" + it + "]" + "\tlikelihood: " + L);
                outputPzd(Pz_d, "model/doc_topic." + it);//即文档向量
                outputPwz(Pw_z, "model/topic_word." + it);
            }
        }

        return false;
    }

    /**
     * 拿计算好的文档向量，去计算所有文档跟第1篇文档的相似度。以此来验证PLSA得到的文档向量是合理的。
     */
    public void test(String docVecFile) {
        BufferedReader br = null;
        try {
            br = new BufferedReader(new FileReader(docVecFile));
            String line = br.readLine();
            if (line == null) {
                return;
            }
            String[] arr = line.split("\\s+");
            if (arr.length < 1 + this.K) {
                System.err.println("1st doc vector's length is less than " + this.K);
                return;
            }
            double[] vec1 = new double[this.K];
            double norm1 = 0.0;//向量模长
            for (int i = 1; i < 1 + this.K; i++) {
                vec1[i - 1] = Double.parseDouble(arr[i]);
                norm1 += vec1[i - 1] * vec1[i - 1];
            }
            norm1 = Math.sqrt(norm1);
            Map<String, Double> simMap = new HashMap<String, Double>();
            while ((line = br.readLine()) != null) {
                arr = line.split("\\s+");
                if (arr.length == 1 + this.K) {
                    String docName = arr[0];
                    double[] vec2 = new double[this.K];
                    double norm2 = 0.0;//向量模长
                    double prod = 0.0;//向量内积
                    for (int i = 1; i < 1 + this.K; i++) {
                        vec2[i - 1] = Double.parseDouble(arr[i]);
                        norm2 += vec2[i - 1] * vec2[i - 1];
                        prod += vec1[i - 1] * vec2[i - 1];
                    }
                    norm2 = Math.sqrt(norm2);
                    double sim = prod / (norm1 * norm2);
                    simMap.put(docName, sim);
                }
            }

            //按相似度从大到小排序
            List<Entry<String, Double>> simList = new ArrayList<Entry<String, Double>>(
                simMap.entrySet());
            Collections.sort(simList, new Comparator<Entry<String, Double>>() {
                @Override
                public int compare(Entry<String, Double> o1, Entry<String, Double> o2) {
                    if (o1.getValue() > o2.getValue()) {
                        return -1;
                    } else if (o1.getValue() < o2.getValue()) {
                        return 1;
                    } else {
                        return 0;
                    }
                }
            });
            //输出前100个与文档1最相似的文档
            for (int i = 0; i < 100 && i < simList.size(); i++) {
                System.out.println(simList.get(i).getKey() + "\t" + simList.get(i).getValue());
            }
        } catch (IOException e) {
            e.printStackTrace();
        } finally {
            try {
                br.close();
            } catch (IOException e) {
            }
        }

    }

    private boolean init(double[][] Pz_d, double[][] Pw_z, double[][][] Pz_dw) {
        // p(z|d), size: M x K
        for (int m = 0; m < this.M; m++) {
            double norm = 0.0;
            for (int z = 0; z < this.K; z++) {
                Pz_d[m][z] = Math.random();
                norm += Pz_d[m][z];
            }

            for (int z = 0; z < this.K; z++) {
                Pz_d[m][z] /= 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: M x K x doc.size()
        for (int m = 0; m < this.M; m++) {
            for (int z = 0; z < this.K; z++) {
                Pz_dw[m][z] = new double[this.dataset.getDataAt(m).size()];
            }
        }
        return false;
    }

    private boolean Estep(double[][] Pz_d, 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_d[m][z] * Pw_z[z][w];
                    Pz_dw[m][z][position] = val;
                    norm += val;
                }
                // 当前文档中的当前词，在各个主题上的概率分布进行归一化
                for (int z = 0; z < this.K; z++) {
                    Pz_dw[m][z][position] /= norm;
                }
            }
        }
        return true;
    }

    private boolean Mstep(double[][][] Pz_dw, double[][] Pw_z, double[][] Pz_d) {
        // p(z|d)
        for (int m = 0; m < this.M; m++) {
            double norm = 0.0;
            for (int z = 0; z < this.K; z++) {
                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[m][z][position];
                }
                Pz_d[m][z] = sum;
                norm += sum;
            }

            // normalization
            for (int z = 0; z < this.K; z++) {
                Pz_d[m][z] /= norm;
            }
        }

        // 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[m][z][position];
                }
                Pw_z[z][w] = sum;
                norm += sum;
            }
            // normalization
            for (int w = 0; w < this.V; w++) {
                Pw_z[z][w] /= norm;
            }
        }

        return true;
    }

    private double calcLoglikelihood(double[][] Pz_d, 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_d[m][z] * Pw_z[z][w];
                }
                L += n * Math.log10(sum);
            }
        }
        return L;
    }

    /**
     * 输出每篇文档在各个主题上的概率分布
     * 
     * @param outFile
     */
    private void outputPzd(double[][] Pz_d, String outFile) {
        BufferedWriter bw = null;
        try {
            bw = new BufferedWriter(new FileWriter(outFile));
            for (int i = 0; i < this.M; i++) {
                String docName = this.dataset.getDataAt(i).docName;
                bw.write(docName);
                for (int j = 0; j < this.K; j++) {
                    bw.write("\t");
                    bw.write(String.valueOf(Pz_d[i][j]));
                }
                bw.newLine();
            }
        } catch (IOException e) {
            e.printStackTrace();
        } finally {
            if (bw != null) {
                try {
                    bw.close();
                } catch (IOException e) {
                }
            }
        }
    }

    /**
     * 输出每个主题下的top100的词
     * 
     * @param outFile
     */
    private void outputPwz(double[][] Pw_z, String outFile) {
        BufferedWriter bw = null;
        try {
            bw = new BufferedWriter(new FileWriter(outFile));
            for (int i = 0; i < Pw_z.length; i++) {
                Map<String, Double> wordWeight = new HashMap<String, Double>();//词在该主题下的权重
                for (int j = 0; j < Pw_z[i].length; j++) {
                    String word = this.dataset.features.get(j);
                    wordWeight.put(word, Pw_z[i][j]);
                }
                List<Entry<String, Double>> wordWeightList = new ArrayList<Entry<String, Double>>(
                    wordWeight.entrySet());
                Collections.sort(wordWeightList, new Comparator<Entry<String, Double>>() {
                    @Override
                    public int compare(Entry<String, Double> o1, Entry<String, Double> o2) {
                        if (o1.getValue() > o2.getValue()) {
                            return -1;
                        } else if (o1.getValue() < o2.getValue()) {
                            return 1;
                        } else {
                            return 0;
                        }
                    }
                });
                for (int j = 0; j < wordWeightList.size() && j < 100; j++) {
                    bw.write(wordWeightList.get(j).getKey() + ":" + wordWeightList.get(j).getValue()
                             + "\t");
                }
                bw.newLine();
            }
        } catch (IOException e) {
            e.printStackTrace();
        } finally {
            if (bw != null) {
                try {
                    bw.close();
                } catch (IOException e) {
                }
            }
        }
    }

    public static void main(String[] args) {
        int nTopic = 50;
        int nIter = 100;
        PLSA plsa = new PLSA();
        if (args.length < 1) {
            System.err.println("train data in docs/user2vec");
            plsa.doPLSA("docs/user2vec", nTopic, nIter);
        } else {
            System.out.println("train data in " + args[0]);
            if (args.length >= 2) {
                nTopic = Integer.parseInt(args[1]);
            }
            if (args.length >= 3) {
                nIter = Integer.parseInt(args[2]);
            }
            plsa.doPLSA(args[0], nTopic, nIter);
        }
        System.out.println("end PLSA");

        String docVecFile = "model/doc_topic." + (nIter - 1);
        plsa.test(docVecFile);
    }
    //nohup java -cp .:plsa.jar plsa.PLSA /data/orisun/cf/data/user_graph.txt 50 100 &
}

Dataset.java

package plsa;

import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.io.IOException;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

/**
 * 文档集合
 * 
 * @author orisun
 * @date 2016年7月10日
 */
public class Dataset {

    /** 文档集合 **/
    List<Data> datas = new ArrayList<Data>();
    /** 记录每个词的编号 **/
    Map<String, Integer> featureIndex = new HashMap<String, Integer>();
    List<String> features = new ArrayList<String>();

    int size() {
        return datas.size();
    }

    int getFeatureNum() {
        return featureIndex.size();
    }

    Data getDataAt(int i) {
        return datas.get(i);
    }

    /**
     * 
     * @param dataDir
     *            如果dataDir是文档集所在的目录。文档格式：每行存储一个词及词在文件中的权重，空格分隔。每篇文档中词可以有重复。<br>
     *            如果所有文档都放在dataDir这一个文件里面，则文件每行的格式为:文件名\t词:权重\t词:权重……
     * @throws IOException
     */
    Dataset(String dataDir) throws IOException {
        File path = new File(dataDir);
        if (path.exists()) {
            int featureNum = 0;
            if (path.isDirectory()) {
                File[] files = path.listFiles();
                for (File file : files) {
                    Data data = new Data();
                    data.docName = file.getName();
                    BufferedReader br = new BufferedReader(new FileReader(file));
                    String line = null;
                    while ((line = br.readLine()) != null) {
                        String[] arr = line.trim().split("\\s+");
                        if (arr.length == 2) {
                            String word = arr[0];
                            double weight = Double.parseDouble(arr[1]);
                            Integer index = featureIndex.get(word);
                            if (index == null) {
                                featureIndex.put(word, featureNum);
                                features.add(word);
                                index = featureNum;
                                featureNum++;
                            }
                            Feature feature = new Feature(index, weight);
                            data.features.add(feature);
                        }
                    }
                    br.close();
                    datas.add(data);
                }
            } else if (path.isFile()) {
                BufferedReader br = new BufferedReader(new FileReader(path));
                String line = null;
                while ((line = br.readLine()) != null) {
                    String[] arr = line.trim().split("\\s+");
                    if (arr.length >= 2) {
                        Data data = new Data();
                        data.docName = arr[0];
                        for (int i = 1; i < arr.length; i++) {
                            String[] brr = arr[i].split(":");
                            if (brr.length == 2) {
                                String word = brr[0];
                                double weight = Double.parseDouble(brr[1]);
                                Integer index = featureIndex.get(word);
                                if (index == null) {
                                    featureIndex.put(word, featureNum);
                                    features.add(word);
                                    index = featureNum;
                                    featureNum++;
                                }
                                Feature feature = new Feature(index, weight);
                                data.features.add(feature);
                            }
                        }
                        datas.add(data);
                    }
                }
                br.close();
            }
        }
    }

}

Data.java

package plsa;

import java.util.ArrayList;
import java.util.List;

/**
 * 文档
 * 
 * @author orisun
 * @date 2016年7月10日
 */
public class Data {

    /** 文档中的所有词 **/
    List<Feature> features = new ArrayList<Feature>();
    /** 文档名称 **/
    String docName;

    int size() {
        return features.size();
    }

    Feature getFeatureAt(int i) {
        return features.get(i);
    }
}

Feature.java

package plsa;

/**
 * 词
 * 
 * @author orisun
 * @date 2016年7月10日
 */
public class Feature {

    /** 该词在所有词中的编号 **/
    int dim;
    /** 该词在指定文档中的权重 **/
    double weight;

    Feature(int index, double weight) {
        this.dim = index;
        this.weight = weight;
    }
}

Posting.java

package plsa;

/**
 * 倒排索引
 * 
 * @author orisun
 * @date 2016年7月10日
 */
public class Posting {

    /** 文档编号 **/
    int docID;
    /** 词在文档中的位置 **/
    int pos;

    Posting(int docID, int pos) {
        this.docID = docID;
        this.pos = pos;
    }
}