Spark MLlib — EMLDA

  LDA(Latent Dirichlet allocation)是一种主题模型，它可以将文档集中每篇文档的主题按照概率分布的形式给出，也即根据给定的一篇文档，推测其主题分布。同时它是一种无监督学习算法，在训练时不需要手工标注的训练集，需要的仅仅是文档集以及指定主题的数量k即可。此外LDA的另一个优点则是，对于每一个主题均可找出一些词语来描述它。本文主要介绍LDA涉及的数学知识以及Spark MLlib中基于Graphx的实现方式。

1.理论基础

  这部分内容主要参考《LDA数学八卦》和《通俗理解LDA主题模型》，涉及到数学内容真是挺多的，网上大多数博文的阐述过程大多一上来就摆出一大堆理论基础，自己学习的过程中也会觉得一开始就很吃力，而且不能把众多理论和最终的模型联系起来，所以一直想找到一种易于理解的方式来阐述该模型。

1.1 LDA引入— Dirichlet先验分布

  在NLP领域，文本的表现形式通常是有序的词序列,即d=(ω1,ω2,...,ωn)。文本建模的实质就是为文本对应的词序列建模。《LDA数学八卦》中以上帝掷骰子的游戏形象的阐述了各种模型生成文本的过程，这里简单概括下各种模型的本质，并假设m篇文本的总词数为N，词典中单词总数为V：

1.1.1 主题无关模型

• Unigram Model：投掷一个V面的骰子，重复N次生成词序列，则这N个随机变量（即生成的N个单词）是独立同分布的，该分布的概率函数即为骰子每一面出现的概率为p⃗ ={p1,p2,...pV}。则词典中每个单词出现次数（记为n⃗ =n1,n2,...,nV）的联合分布满足Multinomial分布(多项分布)，其概率密度函数为：
  p(n⃗ )=Mult(n⃗ |p⃗ ,N)=(Nn⃗ )∏k=1Vpnkk=N!n1!n2!...nV!pn11pn22...pnVV(1.1)
• 贝叶斯Unigram Model：在Unigram Model的基础上为p⃗ 加入了先验分布，即需要先以一定的概率选出骰子。这里的先验分布选择了Multinomial式分布的共轭先验分布Dirichlet分布，Dirichlet分布的一般表现形式如下：
  Dir(p⃗ |α⃗ )=Γ(∑Kk=1αk)∏Kk=1Γ(αk)∏k=1Kpnkk(1.2)
  
  其中Γ(x)=∫∞0tx−1e−tdt=(n−1)!为gamma函数。
  共轭先验分布的本质为：Dirichlet先验+Multinomial分布—>后验分布也为Dirichlet分布 即：
  Dir(p⃗ |α⃗ )+Mult(m⃗ )=Dir(p⃗ |α⃗ +m⃗ )(1.3)
  
  关于gamma函数的更多性质以及Dirichlet先验和Multinomial分布的共轭关系的证明待有空专门整理。

1.1.2 主题相关模型

  文本的生成过程过程依然是生成单词序列的过程，只不过引入主题后，需要先确定主题，再根据主题生成单词的过程，这也更符合人类平常的语言习惯。

• PLSA模型：先投掷一个K面的doc-topic骰子，确定主题编号z，再从K个V面的topic-word骰子中选择编号为z的骰子投掷，得到一个单词，重复该过程生成文档。假设文档dm={ω1,ω2,...,ωn}都有其主题序列概率向量θ⃗ m，编号为k的topic-word骰子的单词序列概率向量为φ⃗ k，则文档dm中每个单词w的生成概率为：
  p(ω|dm)=∑z=1Kp(ω|z)p(z|dm)=∑z=1Kφzwθmz(1.4)
  
  则整篇文档dm的生成概率为：
  p(ω⃗ |dm)=∏i=1n∑z=1Kp(ω|z)p(z|dm)=∏i=1n∑z=1Kφzwθmz(1.5)
• LDA模型：相当于是为 PLSA模型中的主题序列向量θ⃗ m和单词序列向量为φ⃗ k都加入了先验分布。从Unigram Model的投掷过程可知θ⃗ m和φ⃗ k都是满足Multinomial分布的，因此先验分布的选择都是Dirichlet分布，分别为α⃗ 和β⃗ 。

1.1.3 LDA物理图模型

Figure 1. LDA物理图模型

根据该概率图可以得到生成文档的过程可表述为：
(1) 从Dirichlet分布α⃗ 中采样生成文档m的主题分布θ⃗ m
(2) 从主题的Multinomial分布θ⃗ m中采样生成文档m中第 n个词的主题zm,n
(3) 从Dirichlet分布β⃗ 中采样生成主题zm,n(其编号为k) 对应的词语分布φ⃗ k
(4) 从词语的Multinomial分布φ⃗ k中采样生成最终的词语ωm,n
其中α⃗ →θ⃗ m→z⃗ m表示生成第m篇文档中的所有词对应的topics，显然α⃗ →θ⃗ m对应于Dirichlet分布，θ⃗ m→z⃗ m对应于Multinomial分布,所以整体是一个Dirichlet-Multinomial共轭结构：

同理β⃗ →φ⃗ k→w⃗ (k)也是一个Dirichlet-Multinomial共轭结构：

1.2 隐藏变量求解

PLSA模型和LDA模型模型都引入了隐藏变量，即主题z。隐藏变量的参数估计问题常采用的方式为EM算法和Gibbs Sampling。

1.2.1 EM算法

  EM算法的本质也是求解最大似然估计，只不过这里的似然函数是可观察样本x={x(1),x(2),...,x(m)}和隐藏变量z的联合分布概率p(x,z)的似然函数，即：

ℓ(θ)=∑i=1mlog(p;θ)=∑i=1mlog∑zp(x,z;θ)(1.6)

而最大化的似然函数的求解问题通过迭代两步来完成：
(1) E步：在参数确定的情况下(第一次迭代时会初始化参数)通过Jensen不等式【对凹函数而言，期望的函数值大于等于函数值的期望】找到ℓ(θ)下界
(2) M步：通过参数调整使得下界不断上升以此来逼近ℓ(θ)
详细的算法求解可参考《JerryLead—EM算法》。PLSA就是采用EM算法去求解“文档-主题”矩阵Θ=(θ⃗ 1,θ⃗ 2,...,θ⃗ M)T和“主题-词项”矩阵Φ=(φ⃗ 1,φ⃗ 2,...,φ⃗ K)T

1.2.2 Gibbs Sampling

  LDA模型中隐藏变量的参数求解可以采用变分EM算法或Gibbs采样，Spark MLlib中采用的变分EM算法和在线学习的方法，其中EM LDA基于gibbs采样原理来估计参数的。这里只是简单概述下Gibbs采样的本质
  Gibbs 采样解决的主要问题就是如何采样(采样主题及该主题下的单词)使得的采样的样本符合其联合概率的p(x,z)，从而根据这些样本去计算后验概率实现隐藏变量参数估计。 这又要追溯到了马尔科夫稳态的相关知识。
  在《LDA数学八卦》中给出了一个关于后代经济状态迭代的例子，即在给定初代收入阶级(下层，中层，上层)分布π0以及从父代到子代，收入阶层的变化的转移概率P的情况下，计算第n代后代的收入阶级分布πn。这里的状态转移过程是个一阶马尔科夫链，也即当前状态只与其上一个状态有关，而与更早之前的状态无关，即满足如下条件：

P(Xt+1=x|Xt,Xt−1,...)=P(Xt+1=x|Xt)(1.7)

  因此，显然有πn=πn−1P=...=π0Pn。由计算结果可以看出从某一代起，分布结果就收敛了，而且收敛的结果是一个和初代收入阶级分布无关的结果。我们称之为马尔科夫链能收敛到平稳分布。从而引出以下思考：
  对于给定的概率分布p(x)，我们希望能有一种便捷的生成其样本。如果我们能构造一个转移矩阵为P的马氏链，使得该马氏链的平稳分布恰好是p(x)，那么我们从任何一个初始状态x0出发沿着马尔科夫链转移，得到一个转移序列x0,x1,...,xn,xn+1,...,如果马尔科夫链在第n步已经收敛了，于是我们就得到了概率分布满足p(x)的样本xn,xn+1...这里给出Gibbs Sampling计算LDA的迭代公式如下：

p(zi=k|z¬i→,w⃗ )∝n(k)m,¬i+αk∑Kk=1(n(k)m,¬i+αk)⋅n(t)k,¬i+βt∑Vt=1(n(t)k,¬i+βt)∝(n(k)m,¬i+αk)⋅(n(t)k,¬i+βt)∑Vt=1(n(t)k,¬i+βt)(1.8)(1.9)

  其中：¬i表示在语料中去除(zi,wi)用词i之外的单词来重新计算该单词对应的主题概率；
  n(k)m,¬i表示文档m中除了单词i之外归属于主题k的词的数目；
  n(t)k,¬i表示除了单词i之外主题k拥有的词t的总数目

GraphX基础知识

  Spark GraphX的原理细节方面的知识可参考shijinkui/spark_graphx_analyze，这里只介绍GraphX本算法设计到的基本操作的概念

2.1 图的构造

  从上面的分析可知，我们需要维护包含隐藏变量(主题)属性的M∗K规模的矩阵n(k)m以及K∗V规模的矩阵n(t)k，目标是要求出文档—主题分布矩阵θ⃗ m以及主题—词分布矩阵φ⃗ k。为了能够统一地分布式存储n(k)m和n(t)k矩阵，Graphx用了一种巧妙的方法将两个矩阵分布在不同类型的顶点数据上。它把n(t)k矩阵转置为n(k)t矩阵，使得n(k)m和n(k)t矩阵的行统一为K维向量。n(k)m中每篇文档对应一个文档节点，该节点数据为一个K维向量，向量中每个元素表示该主题再该文档中出现的次数（也即该文档中归属于该主题的词频）。n(k)t中每个单词对应一个单词节点，该节点数据也为一个K维向量，向量中每个元素表示该单词在该主题中出现的次数。如下图所示为文档总数为2，语料单词总数为4的情况下LDA模型生成图为：

Figure 2. LDA模型生成图
其中m1,m2为文档节点，t1,t2,t3,t4为词顶点，顶点数据n⃗ 是一个K维向量。n(k)m矩阵以nm1→和nm2→两个向量的形式分别存储在m1,m2两个文档节点上。n(t)k矩阵以nt1→,nt2→,nt3→,nt4→四个向量的形式分别存储t1,t2,t3,t4四个词节点上。与此同时，GraphX以文档到词作为边，以词频作为边数据，把语料库构造成图。上图中文档节点m1和t1之间对应一条有向边，边数据为单词t1在文档m1中的词频。

2.2 图的分布式存储

  Grpahx中构建完图时默认会使用“partitionBy(PartitionStrategy.EdgePartition1D)”方法采用vertexcut(点分割)方式存储图(见Figure 3.)。这种存储方式特点 是任何一条边只会出现在一台机器上，每个点有可能分布到不同的机器上。当点被分割到不同机器上时，是相同的镜像，但是有一个点作为主点(master), 其他的点作为虚点(ghost)，当点B的数据发生变化时,先更新点B的master的数据，然后将所有更新好的数据发送到B的ghost所在的所有机 器，更新B的ghost。这样做的好处是在边的存储上是没有冗余的，而且对于某个点与它的邻居的交互操作，只要满足交换律和结合律，比如求邻居权重的和， 求点的所有边的条数这样的操作，可以在不同的机器上并行进行，只要把每个机器上的结果进行汇总就可以了，网络开销也比较小。代价是每个点可能要存储多份， 更新点要有数据同步开销。

Figure 3. 点分割存储

2.3 图的聚合操作

  LDA算法的迭代过程中需要对每篇文档的每个单词重新选取主题，即重新计算p(zi=k|z¬i→,w⃗ )。这里又利用了GraphX的边操作的特性，它将语料库中每篇文档中的每个词的操作转化为在图中每条边上的操作；而每个单词又要归纳包含该单词的所有文章(即和该单词节点相邻的所有文档节点)对其主题概率分布(即该单词节点上存储的数据—K维向量)的影响。这整个过程就是图的聚合操作，图的聚合操作是一个Map-Reduce操作，如Figure 4所示。

Figure 4. 图的聚合操作

Spark MLlib实践

  这里先给出官方使用LDA模型的例子JavaLatentDirichletAllocationExample.java中的部分关键代码，以此作为源码分析的切入点：

其中line33~46是加载原始文档—词频数据并将每行数据解析为Vector的结构。其中原始数据如下，每一行表示一个文档，每一列表示一个单词，每一个元素Dm,w表示第m篇文档中单词w的词频

  line60中首先new LDA()中指定了默认使用变分EM算法(EMLDA)来学习模型，setK(3)指定了主题个数K为3(不指定的话默认为10)，也对应上文中提到的最终构成的LDA图中文档顶点和图节点的数据内容是3维向量

org.apache.spark.mllib.clustering.LDA.scaladef this() = this(k = 10, maxIterations = 20, docConcentration = Vectors.dense(-1),    topicConcentration = -1, seed = Utils.random.nextLong(), checkpointInterval = 10,    ldaOptimizer = new EMLDAOptimizer)

下面重点分析run方法，这里面完成了LDA模型图的构建以及参数的迭代。

org.apache.spark.mllib.clustering.LDA.scaladef run(documents: RDD[(Long, Vector)]): LDAModel = {    val state = ldaOptimizer.initialize(documents, this)    var iter = 0    val iterationTimes = Array.fill[Double](maxIterations)(0)    while (iter < maxIterations) {      val start = System.nanoTime()      state.next()      val elapsedSeconds = (System.nanoTime() - start) / 1e9      iterationTimes(iter) = elapsedSeconds      iter += 1    }    state.getLDAModel(iterationTimes)  }

(1) 其中initialize()方法完成图的构建（文档节点，词顶点以及文档节点指向词顶点的边）,这里给出部分核心代码

org.apache.spark.mllib.clustering.EMLDAOptimizer.scalaoverride private[clustering] def initialize(      docs: RDD[(Long, Vector)],      lda: LDA): EMLDAOptimizer = {    val docConcentration = lda.getDocConcentration //也即先验参数 alpha    val topicConcentration = lda.getTopicConcentration //也即先验参数 beta    val k = lda.getK //也即主题个数    ......    ......    this.docConcentration = if (docConcentration == -1) (50.0 / k) + 1.0 else docConcentration    this.topicConcentration = if (topicConcentration == -1) 1.1 else topicConcentration    val randomSeed = lda.getSeed    // 为每篇文档中的每个词创建一条(文档->单词)的边，该边以三元组(文档Id，单词Id，词频)的形式存储    val edges: RDD[Edge[TokenCount]] = docs.flatMap { case (docID: Long, termCounts: Vector) =>      // filter()方法过滤掉词频为0的单词      termCounts.asBreeze.activeIterator.filter(_._2 != 0.0).map { case (term, cnt) =>        /*这里的term2index为单词编号，在词频矩阵作为输入条件的情况下，每一列代表一个单词的词频信息，因此为列编号即为单词编号，方法实现里是按照(-1,-2,...)依次编号矩阵的列*/        Edge(docID, term2index(term), cnt)      }    }    /* 根据边生成顶点，edge.srcId边的起点即文档顶点，edge.dstId边的中点即词顶点(顶点数据均为为K维数据)，edge.attr边数据即词频信息。*/    val docTermVertices: RDD[(VertexId, TopicCounts)] = {      val verticesTMP: RDD[(VertexId, TopicCounts)] =        edges.mapPartitionsWithIndex { case (partIndex, partEdges) =>          val random = new Random(partIndex + randomSeed)          partEdges.flatMap { edge =>             //先随机生成K维的主题分布参数            val gamma = normalize(BDV.fill[Double](k)(random.nextDouble()), 1.0)            //为K维随机参数*词频信息作为这条边为其两个顶点分配的K维顶点数据信息            val sum = gamma * edge.attr            Seq((edge.srcId, sum), (edge.dstId, sum))          }        }      verticesTMP.reduceByKey(_ + _)//reduce操作汇总所有边分配给顶点的数据信息    }    //根据顶点和边信息生成图，并采用点分割模式进行分布式存储    this.graph = Graph(docTermVertices, edges).partitionBy(PartitionStrategy.EdgePartition1D)    this.k = k    this.vocabSize = docs.take(1).head._2.size    this.checkpointInterval = lda.getCheckpointInterval    this.graphCheckpointer = new PeriodicGraphCheckpointer[TopicCounts, TokenCount](      checkpointInterval, graph.vertices.sparkContext)    this.graphCheckpointer.update(this.graph)    this.globalTopicTotals = computeGlobalTopicTotals() //计算所有词的主题分布概率和    this  }private def computeGlobalTopicTotals(): TopicCounts = {    val numTopics = k    /*filter方法从所有的图顶点中过滤出词顶点，在前面生成边过程中对单词编号(term2index)时,是按照(-1,-2,...)编号的，所以这里的过滤其实就是查看顶点编号是否小于0*/    graph.vertices.filter(isTermVertex).values.fold(BDV.zeros[Double](numTopics))(_ += _)  }

(2) 接着在循环体里面调用next()方法，迭代计算模型参数

org.apache.spark.mllib.clustering.EMLDAOptimizer.scalaoverride private[clustering] def next(): EMLDAOptimizer = {    val eta = topicConcentration  // 也即先验参数 beta    val W = vocabSize            // 单词总数，也即输入矩阵中列数    val alpha = docConcentration // 也即先验参数 alpha    val N_k = globalTopicTotals   //所有词的主题分布概率和    val sendMsg: EdgeContext[TopicCounts, TokenCount, (Boolean, TopicCounts)] => Unit =      (edgeContext) => {        // N_{wj} 词汇w在文档中的频次        val N_wj = edgeContext.attr        // E-STEP: 计算gamma_{wjk}:词汇w在文档j中分配给主题k的概率,参考公式(2.9)        val scaledTopicDistribution: TopicCounts =          computePTopic(edgeContext.srcAttr, edgeContext.dstAttr, N_k, W, eta, alpha) *= N_wj         //将计算出来的gamma_{wjk}发送消息给边的源顶点和目的顶点        edgeContext.sendToDst((false, scaledTopicDistribution))        edgeContext.sendToSrc((false, scaledTopicDistribution))      }    // 顶点合并消息，用于Map阶段，每个分区中每个点收到的消息合并，以及reduce阶段，合并不同分区的消息    val mergeMsg: ((Boolean, TopicCounts), (Boolean, TopicCounts)) => (Boolean, TopicCounts) =      (m0, m1) => {        val sum =          if (m0._1) {            m0._2 += m1._2          } else if (m1._1) {            m1._2 += m0._2          } else {            m0._2 + m1._2          }        (true, sum)      }    // M-STEP: 每个节点通过收集邻居数据来更新主题权重数据     val docTopicDistributions: VertexRDD[TopicCounts] =      graph.aggregateMessages[(Boolean, TopicCounts)](sendMsg, mergeMsg)        .mapValues(_._2)    // 根据最新顶点数据更新图     val newGraph = Graph(docTopicDistributions, graph.edges)    graph = newGraph    graphCheckpointer.update(newGraph)    globalTopicTotals = computeGlobalTopicTotals()    this  }

最后给出这个Example的输出，主题—词矩阵，如下：