介绍

Mikolov在Distributed Representations of Words and Phrases and their Compositionality中,提到高频词的subsampling问题,以下我对相关节选进行的简单中文翻译:

在非常大的语料中,最高频的词汇,可能出现上百万次(比如:in, the, a这些词)。这样的词汇通常比其它罕见词提供了更少的信息量。例如,当skip-gram模型,通过观察”France”和”Paris”的共现率(co-occurrence),Skip-gram模型可以从中获益;而”France”和”the”的共现率则对模型贡献很少;而几乎每个词都常在句子中与”the”共同出现过。该思路也常用在相反的方向上;高频词的向量表示,在上百万样本训练完后不会出现大变化。

为了度量这种罕见词与高频词间存在不平衡现象,我们使用一个简单的subsampling方法:训练集中的每个词\(w_i\),以下面公式计算得到的概率进行抛弃:

\[P(w_i)=1-\sqrt{\frac{t}{f(w_i)}}\]

\(f(w_i)\)是\(w_i\)的词频,t为选中的一个阀值,通常为\(10^{-5}\)周围(1e-5=0.00001)。我们之所以选择该subsampling公式,是因为:它可以很大胆的对那些词频大于t的词进行subsampling,并保留词频的排序(ranking of the frequencies)。尽管subsampling公式的选择是拍脑袋出来的(启发式的heuristically),我们发现它在实践中很有效。它加速了学习,并极大改善了罕见词的学习向量的准确率(accuracy)。

具体实现

在当中的描述是:

\[P(w_i)=1-(\sqrt{\frac{sample}{freq(w_i)}}+\frac{sample}{freq(w_i)})\]

实际代码如下:

// 进行subsampling,随机丢弃常见词,保持相同的频率排序ranking.
// The subsampling randomly discards frequent words while keeping the ranking same
if (sample > 0) {

  // 计算相应的抛弃概率ran.
  real ran = (sqrt(vocab[word].cn / (sample * train_words)) + 1) * (sample * train_words) / vocab[word].cn;

  // 生成一个随机数next_random.
  next_random = next_random * (unsigned long long)25214903917 + 11;

  // 如果random/65535 - ran > 0, 则抛弃该词,继续
  if (ran < (next_random & 0xFFFF) / (real)65536) 
      continue;
}

为此,我简单地写了一个python程序,做了个测试。程序托管在github上,点击下载

词频越大,ran值就越小。subsampling进行抽样时被抽到的概率就越低。

1.png

参考:

Distributed Representations of Words and Phrases and their Compositionality

介绍

如果你在大学期间学过信息论或数据压缩相关课程,一定会了解Huffman Tree。建议首先在wikipedia上重新温习下Huffman编码与Huffman Tree的基本概念: Huffman编码wiki

简单的说(对于文本中的字母或词),Huffman编码会将出现频率较高(也可认为是:权重较大)的字(或词)编码成短码,而将罕见的字(或词)编码成长码。对比长度一致的编码,能大幅提升压缩比例。

而Huffman树指的就是这种带权路径长度最短的二叉树。权指的就是权重W(比如上面的词频);路径指的是:从根节点到叶子节点的路径长度L;带权路径指的是:树中所有的叶结点的权值乘上其到根结点的路径长度。WPL=∑W*L,它是最小的。

word2vec的Huffman-Tree实现

为便于word2vec的Huffman-Tree实现,我已经将它单独剥离出来,具体代码托管到github上: huffman_tree代码下载。示例采用的是wikipedia上的字母:

即:F:2, O:3, R:4, G:4, E:5, T:7

这里有两个注意事项:

  • 1.对于单个节点分枝的编码,wikipedia上的1/0分配是定死的:即左为0,右为1(但其实分左右是相对的,左可以调到右,右也可以调到左)。而word2vec采用的方法是,两侧中哪一侧的值较大则为1,值较小的为0。当然你可以认为(1为右,0为左)。
  • 2.word2vec会对词汇表中的词汇预先从大到小排好序,然后再去创建Huffman树。在CreateBinaryTree()调用后,会生成最优的带权路径最优的Huffman-Tree。

最终生成的图如下:

此图中,我故意将右边的T和RG父结节调了下左右,这样可以跳出上面的误区(即左为0,右为1。其实是按大小来判断0/1)

相应的数据结构为:

/**
 * word与Huffman树编码
 */
struct vocab_word {
  long long cn;     // 词在训练集中的词频率
  int *point;       // 编码的节点路径
  char *word,       // 词
       *code,       // Huffman编码,0、1串
       codelen;     // Huffman编码长度
};

最后,上面链接代码得到的结果:

1
2
3
4
5
6
word=T	cn=7	codelen=2	code=10	point=4-3-
word=E	cn=5	codelen=2	code=01	point=4-2-
word=G	cn=4	codelen=3	code=111	point=4-3-1-
word=R	cn=4	codelen=3	code=110	point=4-3-1-
word=O	cn=3	codelen=3	code=001	point=4-2-0-
word=F	cn=2	codelen=3	code=000	point=4-2-0-

整个计算过程设计的比较精巧。使用了三个数组:count[]/binary[]/parent_node[],这三个数组构成相应的Huffman二叉树。有vocab_size个叶子节点。最坏情况下,每个节点下都有一个叶子节点,二叉树的总节点数为vocab_size * 2 - 1就够用了。代码使用的是 vocab_size * 2 + 1。

当然,如果你有兴趣关注下整棵树的构建过程的话,也可以留意下这部分输出:

count[]:	7 5 4 4 3 2 1000000000000000 1000000000000000 1000000000000000 1000000000000000 1000000000000000 1000000000000000
binary[]:	0 0 0 0 0 0 0 0 0 0 0 0
parent[]:	0 0 0 0 0 0 0 0 0 0 0 0
	
--step 1:
count[]:	7 5 4 4 3 2 5 1000000000000000 1000000000000000 1000000000000000 1000000000000000 1000000000000000
binary[]:	0 0 0 0 1 0 0 0 0 0 0 0
parent[]:	0 0 0 0 6 6 0 0 0 0 0 0
	
--step 2:
count[]:	7 5 4 4 3 2 5 8 1000000000000000 1000000000000000 1000000000000000 1000000000000000
binary[]:	0 0 1 0 1 0 0 0 0 0 0 0
parent[]:	0 0 7 7 6 6 0 0 0 0 0 0
	
--step 3:
count[]:	7 5 4 4 3 2 5 8 10 1000000000000000 1000000000000000 1000000000000000
binary[]:	0 1 1 0 1 0 0 0 0 0 0 0
parent[]:	0 8 7 7 6 6 8 0 0 0 0 0
	
--step 4:
count[]:	7 5 4 4 3 2 5 8 10 15 1000000000000000 1000000000000000
binary[]:	0 1 1 0 1 0 0 1 0 0 0 0
parent[]:	9 8 7 7 6 6 8 9 0 0 0 0
	
--step 5:
count[]:	7 5 4 4 3 2 5 8 10 15 25 1000000000000000
binary[]:	0 1 1 0 1 0 0 1 0 1 0 0
parent[]:	9 8 7 7 6 6 8 9 10 10 0 0

参考

1.Huffman编码

Y Bengio在2003年发布了paper《A Neural Probabilistic Language Model》,即NNLM。我们来看一下这篇10多年前的paper中的主要内容:

1.介绍

语言建模和其它学习问题很难,因为存在一个基础问题:维度灾难(curse of dimensionality)。当你想建模许多离散随机变量间(比如:一个句子中的词,或者在一个数据挖掘任务中的离散变量)的联合分布时,这个现象特别明显。例如,如果你想建模自然语言中在一个size=10w的词汇表V的10个连续词的联合分布,潜在会有\(100000^{10}-1=10^{50}-1\)个自由参数(free parameters)。当建模连续变量时,我们可以更轻易地获得泛化(例如:使用多层神经网络或高斯混合模型的平滑分类函数),因为学到的函数可以具有一些局部平滑属性。对于离散空间,泛化结果并不明显:这些离散变量的任何变化,都可能对要估计的函数值具有一个剧烈的影响;当离散变量的值可取的数目很大时,大多数观察到的目标相互之间在hamming距离上差距很大

受非参数化密度估计的启发,对不同学习算法如何泛化进行可视化的一种有效方法是,认为将初始集中在training points(例如:训练句子)上的概率质量(probability mass)在一个更大的空间(volume)中分散是的(distributed)。在高维中,将概率质量(probability mass)进行分散(distribute)很重要,而非围绕在每个training point周围的所有方向上均匀。我们在本paper会展示提出的这种方法与state-of-art统计语言建模方法的不同。

统计语言模型可以通过下面的公式进行表示,给定所有之前出现的词,对下一个词的条件概率:

\[\hat{P}(W_1^T) = \prod\limits_{1}^{T} \hat{P}(w_t | w_t^{t-1})\]

其中,\(w_t\)是第t个词,子序列可以写成:\(w_i^j = (w_i, w_{i+1}, \cdots, w_{j-1}, w_j)\)。这样的统计语言模型在许多应用中很有用,比如:语音识别,语言翻译,信息检索。统计语言模型对这样的应用有极大影响。

当构建自然语言的统计模型时,可以通过使用词序(word order)来减小建模问题的难度,实际上在词序列中时序更接近的词在统计学上相互依赖更强。因而,对于一个大数目的上下文(比如:最近n-1个词的组合)中的每个词,可以近似使用n-gram模型来构建对于下一个词的条件概率表

\[\hat{P}(w_t | w_1^{t-1}) \approx \hat{P}(w_t | w_{t-n+1}^{t-1})\]

我们只考虑这些在训练语料中实际出现、或者发生足够频繁的连续词的组合。当n个词的一个新组合出现在训练语料中,会发生什么?我们不想分配零概率到这样cases中,因为这些新组合很可能会出现,而且他们可能对于更大的上下文size会出现的更频繁。一个简单的答案是,使用一个更小的上下文size来看下预测概率,比如:back-off trigram模型、或者smoothed trigram模型。在这样的模型中,从在训练语料中的词序列获取的的模型如何泛化到新的词序列上呢?一种方式是,给这样的插值或back-off ngram模型对应一个生成模型。本质上,一个新的词序列可以通过“粘合(gluning)”在训练语料中频繁出现的非常短或长度为1,2…,n的重叠块来生成。在特别的back-off或插值n-gram算法中,获得下一块(piece)的概率的规则是隐式的。通常,研究者们使用n=3(比如:trigrams),并获取state-of-art结果。。。。

1.1 使用分布式表示来解决维度灾难

简单的,提出的方法可以如下进行总结:

  • 1.将词表中的每个词与一个distributed word feature vector(在\(R^m中\)一个real-valued vector)进行关联。
  • 2.采用在序列中的词的feature vectors来表示词序列的联合概率函数.
  • 3.同时学习word vectors和概率函数的参数

feature vector可以表示词的不同方面:每个词与一个向量空间中的一个点(point)相关联。features的数目(比如:实验中采用m=30, 60或100)比词汇表的size要小很多。概率函数被表示成:给定之前的词,一个关于下一个词的条件概率乘积(例如:使用一个multi-layer NN来预测)。该函数具有这样的参数,它以最大化训练数据的log似然来进行迭代调参。feature vectors是通过学习得到的,但它们可以使用语义特征的先验知识进行初始化。

为什么会有效?在之前的示例中,如果我们知道,dog和cat扮演着相似的角色(语义上和结构上),对于(the,a), (bedroom, room), (is, was), (running, walking)是相似的,我们可以很自然地对下面的句子进行泛化:

The cat is walking in the bedroom

泛化成:

A dog was running in a room

或者:

The cat is running in a room

A dog is walking in a bedroom

The dog was walking in the room

以及许多其它组合。在我们提出的模型中,它是可以泛化的,因为“相似”的词被认为是具有一个相似的feature vector,因为概率函数是一个关于这些feature values的平滑函数,在特征上的微小变化在概率上只会引入很小的变化。因些,在训练数据中上述句子之一的出现将会增加概率,不仅仅是那些句子,而且包括在句子空间中它们的“邻居”的组合数目。

1.2 之前工作

略。

2.一个神经模型

训练集是关于\(w_t \in V\)的一个序列\(w_1, \cdots, w_T\),其中词汇表V是一个有限的大集合。学习目标是学习一个好的模型,使得:

\[f(w_t, \cdots, w_{t-n+1}) = \hat{P}(w_t \mid w_1^{t-1})\]

也就是说给出很高的out-of-sample likelihood。下面,我们会上报关于\(1/\hat{P}(w_t \mid w_1^{t-1})\)的几何平均,它被称为“困惑度(perplexity)”,它也是平均负log似然的指数。模型的唯一限制是,对于\(w_1^{t-1}\)的任意选择,\(\sum\limits_{i=1}^{\mid V\mid} f(i, w_{t-1}, \cdots, w_{t-n+1}) = 1\),其中\(f>0\)。通过将这些条件概率相乘,可以获得一个关于这些词序列的联合概率的模型。

我们将函数\(f(w_t, \cdots, w_{t-n+1})= \hat{P}(w_1^{t-1})\)解耦成两部分:

  • 1.一个映射函数C,它将V的任意元素i映射到一个真实向量\(C(i) \in R^m\)。它表示与词典中每个词相关的分布式特征向量(distributed feature vectors)。实际上,C被表示成一个关于自由参数的\(\mid V \mid \times m\)的矩阵。
  • 2.词上的概率函数,由C进行表示:对于在上下文中的词,\((C(w_{t-n+1}), \cdots, C(w_{t-1}))\),会使用一个函数g将一个关于feature vectors的输入序列映射到一个关于下一个词\(w_t \in V\)的条件概率分布上。g的输出是一个vector,它的第i个元素可以估计概率\(\hat{P}(w_t = i \mid w_1^{t-1})\),如图1所示。
\[f(i, w_{t-1}, \cdots, w_{t-n+1}) = g(i, C(w_{t-1}), \cdots, C(w_{t-n+1}))\]

1.png

图1: 神经网络结构:\(f(i, w_{t-1}, \cdots, w_{t-n+1}) = g(i, C(w_{t-1}), \cdots, C(W_{t-n+1}))\),其中g是神经网络,C(i)是第i个word feature vector。

函数f是这两个mappings(C和g)的一个组合,其中C对于在上下文中所有词是共享的(shared)。这两部分每一个都与一些参数有关。mapping C的参数就简单的是feature vectors自身,通过一个\(\mid V\mid \times m\)的矩阵C进行表示,其中第i行表示第i个词的feature vector C(i)。函数g通过一个feed-forward或RNN或其它参数化函数(parametrized function)来实现,它使用参数\(\mathcal{W}\)。

训练通过搜索\(\theta\)来完成,它会在训练语料上最大化penalized log-likelihood:

\[L = \frac{1}{T} \sum\limits_{t} log f(w_t, w_{t-1}, \cdots, w_{t-n+1}; \theta) + R(\theta)\]

其中\(R(\theta)\)是一个正则项。例如,在我们的实验中,R是一个权重衰减项(weight cecay penalty),它只应用于神经网络的weights和矩阵C上,而非应用在biases上。

在上述模型中,自由参数(free parameters)的数目与词汇表中词数目V成线性比例。它也只与阶数n成线性比例:如果介入更多共享结构,例如,使用一个time-decay神经网络或一个RNN(或者一个关于两种网络的组合),比例因子可以被减小到sub-linear。

在以下大多数实验上,除word features mapping外,神经网络具有一个hidden layer,可选的,还有从word features到output上的直接连接(direct connections)。因此,实际有两个hidden layers:

  • 1.shared word features layer C:它是线性的
  • 2.普通的双曲正切( hyperbolic tangent) hidden layer

更精准的,神经网络会计算如下的函数,它使用一个softmax output layer,它会保证正例概率(positive probabilities)求和为1:

\[\hat{P}(w_t | w_{t-1}, \cdots, w_{t-n+1}) = \frac{e^{y_{w_t}}}{\sum\limits_{i} e^{y_i}}\]

其中,\(y_i\)是对于每个output word i的未归一化log概率,它会使用参数b, W, U, d和H,具体计算如下:

\[y = b + W x + U tanh(d+Hx)\]

…(1)

其中,双曲正切tanh以element-by-element的方式进行应用,W可选为0(即没有直接连接),x是word features layer activation vector,它是从矩阵C的input word features的拼接(concatenation):

\[x = (C(w_{t-1}), C(W_{t-2}), \cdots, C(W_{t-n+1}))\]

假设h是hidden units的数目,m是与每个词相关的features的数目。当从word features到outputs上没有直接连接(direct connections)时,矩阵W被设置为0。该模型的自由参数是output biases b(具有\(\mid V \mid\)个元素),hidden layer biases d(具有h个元素),hidden-to-output weights U(一个\(\mid V \mid \times h\)的矩阵),word features到output weights W(一个\(\mid V\mid \times (n-1) m\)的矩阵),hidden layer weights H(一个 \(h \times (n-1) m\)的矩阵),word features C(一个\(\mid V\mid \times m\)的矩阵):

\[\theta = (b,d,W,U,H,C)\]

自由参数的数目是 \(\mid V\mid (1 + nm + h) + h(1 + (n-1) m)\)。主导因子是\(\mid V\mid (nm+h)\)。注意,在理论上,如果在weights W和H上存在一个weight decay(C上没有),那么W和H可以收敛到0,而C将会增长。实际上,当使用SGA(随机梯度上升)训练时,我们不会观察到这样的行为。

在神经网络上进行SGA包含,包含在执行以下的在训练语料的第t个词之后的迭代更新中:

\[\theta \leftarrow \theta + \epsilon \frac{\partial log \hat{P}(w_t | w_{t-1}, \cdots, w_{t-n+1})}{\partial \theta}\]

其中,\(\epsilon\)是learning rate。注意,大部分参数不需要更新、或者在每个样本之后训练:\(C(j)\)的word features,不会出现在input window中。

混合模型。在我们的实验中,我们已经发现,通过将神经网络的概率预测与一个interpolated trigram model相组合,可以提升效果,使用一个固定的权重0.5,一个可学习权重(在验证集上达到最大似然),或者一个weights集合,是基于在context的频率为条件的(使用相似的过程:在interpolated trigram上组合trigram, bigram, and unigram)。

其它

略.

参考:

http://www.jmlr.org/papers/volume3/bengio03a/bengio03a.pdf

mllib中lda源码分析(一)中,我们描述了LDA的原理、以及变分推断batch算法和online算法的推导。这一节会再描述下spark MLlib中的ml实现。

spark MLlib的实现,是基于变分推断算法实现的,后续的版本会增加Gibbs sampling。本文主要关注online版本的lda方法。(代码主要以1.6.1为主)

一、Batch算法

二、Online算法

ml中的lda相当于一层wrapper,更好地适配ml中的架构(tranformer/pipeline等等),调用的实际上是mllib中的lda实现。

1.LDA.run(主函数)

  def run(documents: RDD[(Long, Vector)]): LDAModel = {
    // 初始化(em=>em, online=>online).
    val state = ldaOptimizer.initialize(documents, this)

    // 迭代开始. 
    var iter = 0
    val iterationTimes = Array.fill[Double](maxIterations)(0)
    while (iter < maxIterations) {
      val start = System.nanoTime()
      
      // optimzier对应的迭代函数.
      state.next()
      
      val elapsedSeconds = (System.nanoTime() - start) / 1e9
      iterationTimes(iter) = elapsedSeconds
      iter += 1
    }

    // 获取最终完整的的LDAModel. 
    state.getLDAModel(iterationTimes)
  }

OnlineOptimizer.next

  override private[clustering] def next(): OnlineLDAOptimizer = {
    // 先进行sample
    val batch = docs.sample(withReplacement = sampleWithReplacement, miniBatchFraction,
      randomGenerator.nextLong())
    
    // 如果为空,返回
    if (batch.isEmpty()) return this

    // 提交batch,进行lda迭代.
    submitMiniBatch(batch)
  }

submitMinibatch

/**
   * 提供语料的minibatch给online LDA模型.为出现在minibatch中的terms它会自适应更新topic distribution。
   * 
   * Submit a subset (like 1%, decide by the miniBatchFraction) of the corpus to the Online LDA
   * model, and it will update the topic distribution adaptively for the terms appearing in the
   * subset.
   */
  private[clustering] def submitMiniBatch(batch: RDD[(Long, Vector)]): OnlineLDAOptimizer = {
    iteration += 1
    val k = this.k
    val vocabSize = this.vocabSize
    val expElogbeta = exp(LDAUtils.dirichletExpectation(lambda)).t      // β ~ Dirichlet(λ)
    val expElogbetaBc = batch.sparkContext.broadcast(expElogbeta)       // 
    val alpha = this.alpha.toBreeze
    val gammaShape = this.gammaShape

    // step 1: 对每个partition做map,单独计算E-Step. stats(stat, gammaPart)
    val stats: RDD[(BDM[Double], List[BDV[Double]])] = batch.mapPartitions { docs =>

      // 1.过滤掉空的.
      val nonEmptyDocs = docs.filter(_._2.numNonzeros > 0)

      // 2.stat 是一个DenseMatrix:  k x vocabSize
      val stat = BDM.zeros[Double](k, vocabSize)

      // 
      var gammaPart = List[BDV[Double]]()

      // 3.遍历partition上的所有document:执行EM算法.
      // E-step可以并行计算.
      // M-step需要reduce,作合并,然后再计算.
      nonEmptyDocs.foreach { case (_, termCounts: Vector) =>
        // 3.1  取出document所对应的词的id。(支持dense/sparse).
        val ids: List[Int] = termCounts match {
          case v: DenseVector => (0 until v.size).toList
          case v: SparseVector => v.indices.toList
        }

        // 3.2 更新状态 E-step => gammad ().
        val (gammad, sstats) = OnlineLDAOptimizer.variationalTopicInference(
          termCounts, expElogbetaBc.value, alpha, gammaShape, k)

        // 3.3 根据所对应id取出对应列. 更新sstats(对应主题状态) 
        stat(::, ids) := stat(::, ids).toDenseMatrix + sstats

        // 3.4 更新该partition上每个文档的gammad的gammaPart列表中.
        gammaPart = gammad :: gammaPart
      }

      // 4.mapPartition返回iterator,每个partition返回一个迭代器(stat, gammaPart)
      // stat: k x V matrix. 针对该partition上的文档,所更新出的每个词在各主题上的分布.
      // gammaPart: list[vector[K]] 该分区上每个文档的gammad列表.
      Iterator((stat, gammaPart))
    }

    // step 2: 对mini-batch的所有partition上的stats(stat,gammaPart)中的stat进行求总和.
    val statsSum: BDM[Double] = stats.map(_._1).reduce(_ += _)
    expElogbetaBc.unpersist()

    // step 3: 对mini-batch的所有partition上的stats(stat,gammaPart)中的gammaPart进行更新.
    val gammat: BDM[Double] = breeze.linalg.DenseMatrix.vertcat(
      stats.map(_._2).reduce(_ ++ _).map(_.toDenseMatrix): _*)

    // step 4: 更新batchResult ( K x V), 每个元素上乘以 E[log(β)]
    val batchResult = statsSum :* expElogbeta.t

    // M-step:
    // 更新λ.
    // Note that this is an optimization to avoid batch.count
    updateLambda(batchResult, (miniBatchFraction * corpusSize).ceil.toInt)

    // 如果需要优化DocConentration,是否更新alpha
    if (optimizeDocConcentration) updateAlpha(gammat)
    this
  }

variationalTopicInference

里面比较核心的一个方法是variationalTopicInference:

  private[clustering] def variationalTopicInference(
      termCounts: Vector,
      expElogbeta: BDM[Double],
      alpha: breeze.linalg.Vector[Double],
      gammaShape: Double,
      k: Int): (BDV[Double], BDM[Double]) = {

    // step 1: 词的id和cnt: (ids, cnts) 
    val (ids: List[Int], cts: Array[Double]) = termCounts match {
      case v: DenseVector => ((0 until v.size).toList, v.values)
      case v: SparseVector => (v.indices.toList, v.values)
    }

    // step 2: 初始化: gammad ~ Γ(100, 0.01), 100维
    // gammd: mini-batch的变分分布: q(θ | γ) 
    // expElogthetad: paper中的exp(E[logθ]), 其中θ ~ Dirichlet(γ)
    // expElogbetad:  paper中的exp(E[logβ]), 其中β : V * K
    // phiNorm:  ϕ ∝ exp{E[logθ] + E[logβ]},ϕ ~ θ*β. 非零.
    // Initialize the variational distribution q(theta|gamma) for the mini-batch
    val gammad: BDV[Double] =
      new Gamma(gammaShape, 1.0 / gammaShape).samplesVector(k)                   // K
    val expElogthetad: BDV[Double] = exp(LDAUtils.dirichletExpectation(gammad))  // K
    val expElogbetad = expElogbeta(ids, ::).toDenseMatrix                        // ids * K

    // 加上一个很小的数,让它非零.
    val phiNorm: BDV[Double] = expElogbetad * expElogthetad :+ 1e-100            // ids
    var meanGammaChange = 1D
    val ctsVector = new BDV[Double](cts)                                         // ids

    // 单个mini-batch里的loop.
    // 迭代,直到 β 和 ϕ 收敛. paper中是0.000001,此处用的是1e-3.
    // Iterate between gamma and phi until convergence
    while (meanGammaChange > 1e-3) {
      val lastgamma = gammad.copy

      // breeze操作:https://github.com/scalanlp/breeze/wiki/Universal-Functions
      // ":"为elementwise操作符;其中的:=,对象相同,内容赋值
      // 计算:γ_dk, 先对ctsVector进行归一化,再乘以转置E(log(β)]^T,再pointwise乘E(log(θ)).
      // 各矩阵或向量的维度: K为topic, ids:为词汇表V
      // gammad: Vector(K),单个文档上每个主题各有一γ值.
      // expElogthetad: Vector(K),单个文档上每个主题各有一θ值.
      // expElogbetad: Matrix(ids*K),每个词,每个主题都有一β值.
      // ctsVector: Vector(ids),每个词上有一个词频.
      // phiNorm: Vector(ids),每个词上有一个ϕ分布值。
      // 
      // 
      //        K                  K * ids               ids
      gammad := (expElogthetad :* (expElogbetad.t * (ctsVector :/ phiNorm))) :+ alpha
      
      // 更新 exp(E[logθ]), expElogbetad不需要更新
      expElogthetad := exp(LDAUtils.dirichletExpectation(gammad))
      
      // 更新 phiNorm, 
      // TODO: Keep more values in log space, and only exponentiate when needed.
      phiNorm := expElogbetad * expElogthetad :+ 1e-100

      // 平均变化量.
      meanGammaChange = sum(abs(gammad - lastgamma)) / k
    }

    // sstats: mini-batch下每个主题下的词分布.
    // n 
    // sstatsd: k x 1 * 1 x ids => k x ids
    val sstatsd = expElogthetad.asDenseMatrix.t * (ctsVector :/ phiNorm).asDenseMatrix
    (gammad, sstatsd)
  }

ok。关于alpha的更新等,此处不再解释。详见源码。

在对MLlib中的lda进行源码分析之前,我们先熟悉下blei的LDA paper以及相关的概念。

介绍

首先,先做下约定。

  • 一个词(word)是离散数据的基本单位,被定义成词汇表中的一个item,由{1,…,V}进行索引。我们将词表示成one-hot格式。使用上标来表示components,在词汇表中第v个词被表示成一个V维的向量w,其中$ w^v=1, w^u=0(u\neq{v}) $。
  • 一个文档(document)是一个N个词组成的序列,表示成W=(w1,w2,…,wN),其中wn是序列上的第n个词。
  • 一个语料(corpus)是一个关于M个文档的集合,被表示成D={W1,W2,…,Wm}.

(blei论文中用黑体的w表示文档,为便于区分,此处用大写W)。

我们希望找到一个关于该语料的概率模型,它不仅可以将高概率分配给该语料的组成文档,还可以将高概率分配给其它“相似”文档。

LDA

LDA是一个关于语料的生成概率模型。基本思想是,文档可以被表示成在隐主题(latent topics)上的随机混合,其中每个主题(topic)都以一个在词(words)上的分布进行表示。

对于在语料D中的每个文档W,LDA假设如下的生成过程:

  • 1.选择 N ~ Poisson(ξ)
  • 2.选择 θ ~ Dir(α)
  • 3.对于每个文档(N个词汇wn):
    • (a) 选择一个主题$z_n$ ~ Multinomial(θ)
    • (b) 以主题$z_n$为条件,使用多项分布条件概率(multinomial probability):$ P(w_n \mid z_n, \beta) $选中一个词$w_n$

在该基础模型上,做了一些简化假设。

  • 首先,Dirichlet分布的维度k(主题变量z的维度),假设是已知并且确定的。
  • 第二,词概率由一个k x V的矩阵β进行参数化,其中$ \beta_{ij}=p(w^j=1 \mid z^i=1) $, 目前当成是一个待估计的确定量。
  • 最后,Poisson猜想是不严格的,可以使用更多的实际文档长度分布。注意N对于所有其它数据生成变量(θ和z)是独立的。它是个辅助变量,我们通常忽略它的随机性。

一个k维的Dirichlet随机变量θ,它的取值为(k-1)-simplex,(一个k维向量θ,它在泛化三角形(k-1)-simplex之内,其中:$ \theta_i \geq 0, \sum_{i=1}^k\theta_i=1 $) ,在该simplex上具有以下的概率密度:

\(p(\theta|\alpha)=\frac{\Gamma(\sum_{i=1}^{k}\alpha_{i})}{\prod_{i=1}^{k}\Gamma(\alpha_{i})}{\theta_{1}}^{\alpha_{1}-1}...{\theta_{k}}^{\alpha_{k}-1}\) ……(1)

其中,参数α是一个k维的vector,相应的components上: αi > 0, 其中 Γ(x)为Gamma函数。Dirichlet是在simplex上的一个合适分布——它在指数族内,具有有限维的充分统计量(finite dimensional sufficient statistics),与multinomial分布共轭。在paper第5部分,这些属性将有助于LDA的inference和parameter estimation算法。

给定参数 α 和 β,我们可以给出关于一个主题混合θ,一个关于N个主题的z集合(每个词都有一个主题),一个N个词W的的联合分布

\(p(\theta,z,W|\alpha,\beta)=p(\theta|\alpha) \prod_{n=1}^{N}p(z_n|\theta)p(w_n|z_n,\beta)\) ……(2)

其中$p(z_n \mid \theta) $可以简单认为:θi对于唯一的i,有$ {z_{n}}^{i}=1 $。在θ上积分(θ上连续),并在z上进行求和(z上离散),我们可以得到一个关于文档的边缘分布(marginal distribution)

\(p(W|\alpha,\beta)=\int p(\theta|\alpha)(\prod_{n=1}^{N}\sum_{z_n}p(z_n|\theta)p(w_n|z_n,\beta))d\theta\) ……(3)

最后,将所有单个文档(document)的边缘分布进行连乘运算,我得到整个语料(corpus)的概率

\[p(D|\alpha,\beta)=\prod_{d=1}^{M}\int p(\theta_{d}|\alpha)(\prod_{n=1}^{N}\sum_{z_{dn}}p(z_n|\theta)p(w_n|z_n,\beta))d\theta_{d}\]

LDA可以表示成图1所示的概率图模型。有三个级别的LDA表述。

  • 语料级参数:α 和 β是语料级参数(corpus-level parameters),假设在生成一个语料的过程中只抽样一次。
  • 文档级变量:θd是文档级变量(document-level variable),每个文档抽样一次。
  • 词级别变量:$z_{dn}$和$w_{dn}$是词级别变量(word-level variables),对于每个文档中的每个词抽样一次。

将LDA与一个简单的Dirichlet-multinomial clustering模型相区分很重要。一个经典的clustering模型将涉及到一个两层模型(two-level),对于一个语料只抽样一次Dirichlet;对于语料中的每个文档,只选一次multinomial clustering变量;对于在基于该cluster变量条件的文档上只选择一个词集合。有了多个clustering models后,这样的模型会限制一个文档只与单个主题相关联。而LDA涉及三层(three-level),尤其是主题节点(topic node)在单个文档中被重复抽样。在该模型下,文档可以与多个主题相关联。

图一:LDA的图形化模型表示。”plates”表示可重复。outer plate表示文档,inner plate表示在单个文档中关于主题和词汇的可重复的选择。

通常在贝叶斯统计建模(Bayesian statistical modeling)中学到与图1相似的结构。它们被称为层级模型(hierarchical models),或更精确地称为:条件独立层级模型(con-ditionally independent hierarchical models)(Kass and Steffey, 1989).这样的模型经常作为参数期望贝叶斯模型(parametric empirical Bayes models)被提到,它也被用于参数估计中。在第5部分,我们将采用经验贝叶斯方法(empirical Bayes approach),来估计在LDA的简单实现中的α 和 β, 我们也会考虑更完整的Bayesian方法。

2.1 LDA与可交换性

随机变量 {z1,…,zN}的一个有限集合,如果联合分布在置换(permutation)后不变,我们称之为是可交换的(exchangeable)。如果π是从整数1到N的一个置换(permutation):

\[p(z_1,...,z_N)=p(z_{\pi(1)},...,z_{\pi(N)})\]

一个无限的随机变量序列,如果每个有限的子序列是可交换的,那么就是无限可交换的(infinitely exchangeable)。

De Finetti’s representation理论声称:一个无限可交换的随机变量序列的联合概率,如果一个随机参数从这些分布中抽取,那么问题中的随机变量,在其于该参数条件下是独立同分布的(iid: independent and identically distributed)。

在LDA中,我们假设,词由主题(确定的条件分布)生成,这些主题在一个文档中是无限可交换的(infinitely exchangeable)。由de Finetti的理论,一个词和主题的序列的概率,必须具有以下形式:

\[p(W,z)=\int p(\theta) (\prod_{n=1}^{N} p(z_n|\theta) p(w_n|z_n)) d\theta\]

其中,θ是在这些主题上的一个服从多项分布的随机变量。我们可以在等式(3)中获取LDA分布,通过对主题变量边缘化,并让θ服从一个Dirichlet分布。

2.2 其它模型

unigram模型

在unigram模型中,每个文档中的词都从一个单独的multinomial分布独立抽取得到:

\[p(W)=\prod_{n=1}^{N} p(w_n)\]

Mixture of unigrams

如果unigram模型和一个离散随机主题变量z混合,我们就得到了一个mixture of unigrams model。在该混合模型下,每个文档的生成:首先选中一个主题z,然后以条件多项概率$p(w \mid z)$独立生成N个词。一个文档的概率为:

\[p(W)=\sum_{z} p(z) \prod_{n=1}^{N} p(w_n|z)\]

当从一个语料估计时,词分布可以看成是:假设每个文档只有一个主题,在该假设下的主题表示。该假设对于大语料集来说是一个很受限的模型。

相反的,LDA允许文档展现多个主题。这会额外带来另一个参数形销:有k-1参数与mixture of unigrams中的p(z)有关,而LDA模型中有k个参数与$ p(\theta \mid \alpha)$有关。

pLSI

pLSI是另一个广泛使用的文档模型(Hofmann 1999)。pLSI模型认为,一个文档d,词wn,对于一个未观察到的主题z来说是条件独立的:

\[p(d,w_n)=p(d)\sum_z p(w_n|z) p(z|d)\]

pLSI模型尝试放宽maxture of unigrams模型中作出的简化猜想:每个文档只由一个主题生成。在某种意义上,它会捕获概率:一个文档可以包含多个主题,因为p(z | d)可看成是对于一个特定文档d的主题权重的混合。然后,注意,d是一个dummy index,指向在训练集中文档列表。这样,d是一个多项分布随机变量,值尽可能与训练文档一样多,模型学到的主题混合 p(z | d)只对应于被训练的文档。出于该原因,pLSI并不是一种定义良好的(well-defined)文档生成模型;天然不支持使用它来分配概率给一个之前未见过的文档。

图3: 不同模型的图表示

2.3 几何学解释

LDA与其它隐主题模型不同,可以通过隐空间(latent space)的几何表示来解释。来看下在每种模型下,一个文档在几何上是如何表示的。

上面所描述的所有4种模型(unigram, mixture of unigrams, pLSI, LDA),都是在词的空间分布上操作。每个这样的分布,可以看成是在(V-1)-simplex上的一个点,我们称之为word simplex。simplex是泛化三角形,关于simplex,参见simplex的wiki

unigram模型是在word simplex上找到单个点,假定语料中的所有词都来自该分布。隐变量模型(letent variable models)会考虑在word simplex上的k个点,并基于这些点形成一个子simplex(sub-simplex),我们称之为topic simplex。注意,在topic simplex上的任意点,在word simplex也是一个点。使用topic simplex的不同隐变量模型,以不同的方式生成一个文档。

  • 1-gram混合模型(mixture of unigram model):假设对于每个文档,在word simplex上的k个点(也就是说:topic simplex的其中一个角落)被随机选中,该文档的所有词汇都从这些点的分布上抽取得到。
  • pLSI模型:假定一个训练文档的每个词都来自一个随机选中的topic。这些主题(topics)本身从一个在这些主题上的指定文档分布上抽取到,例如,在topic simplex上的一个点。对于每个文档都有一个这样的分布;训练文档的集合定义了在topic simplex上的一个经验分布。
  • LDA: 假定,已观察到和未观察到的文档上的每个词,都由一个随机选中的topic生成,这个topic从一个随机选中的参数的分布上抽取。该参数从来自topic simplex的一个平滑分布上的每个文档中抽样得到。

图4: 嵌在包含三个词的word simplex上的三个主题的topic simplex。word simplex的角(corners)对应于三个分布,每个词各自都具有一个概率分布。topic simplex的三个顶点对应于在词上的三个不同分布。The mixture of unigrams模型会将每个文档放到topic simplex上其中一个角(corners)上。pLSI模型会引入根据x的topic simplex的一个经验分布。LDA则在由等高线表示的topic simplex上放置一个平滑分布。

3.推断与参与估计

3.1 inference

为了使用LDA,我们需要解决的核心推断问题是,对于一个给定文档,计算这些隐变量的后验分布(主题分布)

\[p(\theta,z|W,\alpha,\beta)=\frac{p(\theta,z,W|\alpha,\beta)}{p(W|\alpha,\beta)}\]

不幸的是,该分布很难计算。为了归一化该分布,我们对隐变量边缘化(marginalize),并将等式(3)表示成模型参数的形式:

\[p(W|\alpha,\beta)=\frac{\Gamma(\sum_i\alpha_i)}{\prod_{i}\Gamma(\alpha_i)} \int (\prod_{i=1}^{k}\theta_{i}^{\alpha_i-1}) (\prod_{n=1}^{N}\sum_{i=1}^{k}\prod_{j=1}^{V}(\theta_i\beta_{ij})^{w_n^j}) d\theta\]

该函数很难解,因为θ 和 β在隐主题的求和上相耦合(Dickey, 1983)。Dickey展示了该函数是一个在Dirichlet分布(可以表示成一个超几何函数)的特定扩展下的期望。它被用在一个Beyesian上下文中…

尽管后验分布很难精准推断(inference),有许多近似推断算法可以用于LDA,包括Laplace近似,变分近似(variational approximation),马尔可夫链蒙特卡罗法MCMC(jordan,1999)。本节我们描述了一种在LDA中简单的凸变分推断算法(convexity-based variational algorithm for inference)。其它方法在paper 第8部分讨论。

3.2 变分推断

凸变分推断算法的基本思想是,充分利用Jensen不等式来获取一个在log likelihood上的可调下界(jordan et al.,1999)。本质上,可以考虑一个关于下界的家族(a family of lower bounds),由一个变分参数(variational parameters)集合进行索引。变分参数由一个优化过程进行选择,它会尝试找到最紧可能的下界(tightest possible lower bound)。

获取一个可控下界家族的一个简单方法是,考虑将原始的图示模型进行简单修改,移去一些边(edge)和节点(node)。将图1所示的LDA模型进行修改, θ 和 β之间的耦合,由于边θ, z, 和 W之间存在边(edge)。通过抛弃这些边以及W节点,生成一个更简单的图示模型,它有自由的变分参数,我们可以在隐变量上获取一个分布族。该分布族由以下的变分分布组成:

\(q(\theta,z|\gamma,\phi)=q(\theta|\gamma)\prod_{n=1}^{N}q(z_n|\phi_n)\) ……(4)

其中,Dirichlet分布参数γmultinomial分布参数(φ1 , . . . , φN),是自由变分参数。

图4: (左)LDA的图示模型 (右)采用变分分布来近似LDA后验的图示模型

在指定了一个简化版本的概率分布族后,下一步是建立一个优化问题,来确定变分参数γ 和 φ的值。正如在附录A中展示的,找到一个在log likelihood上的紧下限,可以直接翻译成下面的优化问题:

\((\gamma^{*},\phi^{*})=arg min_{\gamma,\phi}^{} D(q(\theta,z|\gamma,\phi)||p(\theta,z|W,\alpha,\beta))\)……(5)

变分参数的最优值,可以通过对变分分布$ q(\theta,z \mid \gamma,\phi) $和真实后验概率$ p(\theta,z \mid W,\alpha,\beta) $间的KL散度进行最小化得到。该最小化可以通过一个迭代型的定点方法(fixed-point method)完成。特别的,我们在附录A.3中展示了:通过计算KL散度的导数,将它们等于0, 可以得到以下更新等式的pair:

\(\phi_{ni} \propto \beta_{iw_n}exp\{E_q[log(\theta_i)|\gamma]\}\)……(6)

\(\gamma_{i}=\alpha_{i}+\sum_{n=1}^N\phi_{ni}\) ……(7)

在附录A.1中,多项分布的期望更新可以以如下方式计算:

\(E_q[log(\theta_i)|\gamma]=\Phi(\gamma_i)-\Phi(\sum_{j=1}^{k}\gamma_{j})\) ……(8)

其中,Ψ 是logΓ函数通过Taylor展开式的首个导数项。

等式(6)和(7)具有一个吸引人的直觉解释。Dirichlet分布更新(Dirichlet update)是一个在给定变分分布($ E[z_n \mid \phi_n] $>)下给定的期望观察的后验Dirichlet。多项分布更新(multinomial update)类似于使用贝叶斯理论,$ p(z_n \mid w_n) \propto p(w_n \mid z_n) $,其中p(zn)是在变分分布下的期望值的log的指数近似。

注意,变分分布实际上是一个条件分布,由一个关于W的函数区分。因为在等式(5)中的优化问题,由确定的W所管理,因而,最优解($ \gamma^{\ast},\phi^{\ast} $)是一个关于W的函数。我们可以将产生的变分分布写成:$ q(\theta,z \mid \gamma^{\ast}(W),\phi^{\ast}(W)) $,其中,我们已经在W上显式地作出了独立性。这样,变分分布可以看成是一个对后验分布$ p(\theta,z \mid W,\alpha,\beta)$的近似。

在文本语言中,最优参数($ \gamma^{\ast}(W), \phi^{\ast}(W) $)是由文档决定的(document-specific)。特别的,我们可以将Dirichlet参数$ \gamma^{\ast}(W) $看成是在topic simplex中提供一个文档表示。

图6: LDA的variational inference算法

图6展示了变分推断过程的总结,对于 γ 和 φn具有合适的起始点。从该伪代码我们可以知道:LDA的变分推断算法,每次迭代需要O((N+1)k)次操作。经验上,我们可以找到对于单个文档的迭代数,与文档中词的数目相近似。因而总的操作数与$ N^2k $相近。

3.3 参数估计

本节描述了一种在LDA模型中用于参数估计的经验贝叶斯方法(empirical Bayes method)。特别的,给定一个文档集的语料库:D={W1,W2,…WM}。我们希望找到α 和 β,来最大化整个语料数据的(marginal)log likelihood:

\[l(\alpha,\beta)=\sum_{d=1}^{M}logp(W_d|\alpha,\beta)\]

如上所述,p(W | α,β)的计算很困难。然而,变分推断提供给我们一个在log likelihood上的下界,我们可以分别根据α和β进行最大化。我们找到了近似经验贝叶斯估计,通过一种交替变分EM过程(alternating variational EM procedure),来分别对变分参数γ 和 φ,最大化下界。接着,变分参数确定下来后,再分别根据模型参数α 和 β,最大化下界。

我们在附录A.4中提供了一个用于LDA的变分EM算法的导数的详细推导。推导以下面的迭代形式描述:

  • 1.(E-step): 对于每个文档,找到变分参数{$ \gamma_d^{\ast},\phi_d^{\ast}: d \in D $}。这个在上一节已描述。
  • 2.(M-step): 分别根据模型参数α 和 β,最大化在log likelihood上产生的下界。这相当于,在E-step中计算出的合适的后验概率下,为每个文档使用期望的足够统计量,找到极大似然估计。

这两个step会一直循环,直到log likelihood的下界收敛为止。

在附录A.4,我们展示了对于条件多项分布参数β进行M-step更新(M-step update),可以写成:

\(\beta_{ij} \propto \sum_{d=1}^{M}\sum_{n=1}^{N_d}\phi_{dni}^{*}w_{dn}^j\) ……(9)

我们会进一步展示Dirichlet参数α的M-step更新,它可以使用一种有效在线性时间上增长的Newton-Raphson方法,其中的Hessian是可逆的。

3.4 Smoothing

对于许多文档的大语料,具有很大的词汇表size,通常会有严重的sparsity问题。对于一个新文档,它包含的词汇很可能没有出现在训练语料中。对于这些词来说,多项分布参数(multinomial parameters)的最大似然估计会分配概率为0,因而,这些新文档的概率为0。应付该问题的方法是,对多项分布参数进行”平滑(smooth)”,对于所有词汇表item都分配一个正的概率,不管它们是否出现在训练集当中。最常用的方法是Laplace smoothing;这本质上会产生在多项分布参数上的均匀Dirichlet先验分布下,对后验分布求平均。

不幸的是,在混合模型设置中,简单的Laplace smoothing不再适合最大化后验方法。事实上,在多项分布参数上放置一个Dirichlet先验,我们可以在混合模型设置上得到一个难解的后验,这与在基本的LDA模型上得到一个难解的后验的原因基本类似。我们提出解决该问题的方法是,简单地应用变分推断方法到扩展模型上,它包括了在多项分布参数上的Dirichlet smoothing

图7:带平滑的LDA图示模型

在LDA设置中,我们可以获得如图7所示的扩展模型。我们将β看成是一个k x V的随机矩阵(每行表示一个mixture component),其中我们假设每行都是从一个可交换的Dirichlet分布上独立抽取的。我们现在扩展我们的inference过程,将βi看成是随机变量,它天然具有一个基于训练语料条件的后验分布。此时我们来看下这个更完整的LDA Bayesian方法。

将一个变分过程看成是一个Bayesian inference,它在随机变量β, θ 和 z上放置一个独立分布(Attias,2000):

\[q(\beta_{1:k},z_{1:M},\theta_{1:M}|\lambda,\phi,\gamma)=\prod_{i=1}^{k}Dir(\beta_i|\lambda_i) \prod_{d=1}^M q_d(\theta_d,z_d|\phi_d,\gamma_d)\]

其中 $ q_d(\theta,z \mid \phi,\gamma) $是在等式(4)中定义的变分分布。很容易验证,产生的变分推断过程会产生等式(6)和(7)来分别作为对变分参数φ 和 γ更新等式,同时也会对一个新的变分参数λ做一个额外更新:

\[\lambda_{ij}=\eta + \sum_{d=1}^{M} \sum_{n=1}^{N_d} \phi_{dni}^{*} w_{dn}^j\]

迭代这些等式直到收敛,会产生一个在β, θ 和 z上的一个合适后验分布。

现在我们将超参数η用于可交换的Dirichlet,如同前面提到的超参数α。设置这些超参数的方法是:我们使用变分EM来寻找基于边缘似然上的这些参数的最大似然估计。这些过程在附录A.4描述。

4.online算法

对于主题模型,后验概率是很难计算的,可以通过近似后验推断来计算。当前的近似后验推导算法分为两大类:抽样方法(sampling approaches)和优化解法(optimization approaches)。抽样方法有MCMC,优化解法有VB(Variational Bayes)。VB在经验上比MCMC更快,与MCMC准确率相当,对于大数据集来说,VB是更吸引人。

但是,使用VB在大数据上计算很困难。标准的“batch”VB算法,会迭代分析每个观察(observation),并迭代更新数据集范围的变分参数。batch算法的每次迭代开销,对于非常大的数据集来说不切实际。在主题建模(topic modeling)应用中,这一点特别重要。主题建模会对不能手工标注的大型文档集合进行归纳隐式结构。

为此,blei等提出了一种online VB算法,这种算法基于online随机优化(online stochastic optimization),在大数据集上,它比batch算法更快地产生好的参数估计。Online LDA可以方便地分析大量文档集合,不需要本地存储,或者收集文档,每个文档可以以流式(streaming)到达,看完后抛弃。

以下的等式重新标号。

假设有K个主题。每个主题定义了一个在词汇表上的多项分布,假设这些主题从一个Dirichlet上抽取得到: βk ∼ Dirichlet(η)。给定主题后,LDA为每个文档d假设了以下的生成过程。

  • 首先,从主题$ \theta_d \sim Dirichlet(\alpha) $上抽取一个分布;
  • 接着,对于在文档中的每个词i,从主题权重$ z_{di} \sim \theta_d $ 上抽取一个主题索引$ z_{di} \in {1, . . . , K }$,从选中的主题上抽取观察词$w_{di}$,$w_{di} \sim \beta_{z_{di}}$。出于简洁性,我们假设在θ 和 β上具有对称先验,但该假设很容易放宽。

注意,如果我们在主词分配z上进行求和,接着,我们得到:$p(w_{di} \mid \theta_d,\beta)=\sum_k\theta_{dk}\beta_{kw}$。这会导致LDA的”multinomial PCA”解释。我们可以将LDA看成是将一个词频n的矩阵(其中$n_{dw}$是词w出现在文档d中的次数),概率分解(probabilistic factorization)成一个关于主题权重θ的矩阵,以及一个主题字典β。我们的工作可以看成是在线矩阵分解技术的一个扩展,它会最优化squared error来得到更泛化的概率公式。

  • 主题(the topics): β
  • 主题比例(topic proportions): θ
  • 主题分配(topic assignments): z

4.1 Batch算法

在VB推断中,真实的后验由一个更简单的分布q(z,θ,β)来近似,它由一个自由参数集进行索引。这些参数用于优化:最大化置信下界(ELBO:Evidence Lower BOund)

\(logp(W|\alpha,\eta) \geq \mathcal{L} (W,\phi,\gamma,\lambda) \triangleq E_q[logp(W,z,\theta,\beta|\alpha,\eta)] - E_q[logq(z,\theta,\beta)]\) …… (1)

最大化ELBO,相当于最小化q(z,θ,β)与后验$p(z,\theta,\beta \mid W,\alpha,\eta)$间的KL散度。我们选择一个完整的因子分布q(fully factorized distribution):

\[q(z_{di}=k)=\phi_{dw_{di}k}\] \[q(\theta_d)=Dirichlet(\theta_d;\gamma_d)\] \[q(\beta_k)=Dirichlet(\beta_k;\lambda_k)\]

……(2)

每个词的主题分配z的后验,可以由φ参数化,在每个文档上的主题权重θ可以由γ进行参数化,在主题上的后验β可以由λ进行参数化。为便于记忆,我们将λ作为”主题(the topics)”。等式1分解为:

\[\mathcal{L} (W,\phi,\gamma,\lambda) = \sum_{d} \{ E_q[log p(w_d | \theta_d,z_d,\beta)] \\ + E_q [log p(z_d | \theta_d)] - E [logq(z_d)] + E [logp(\theta_d|\alpha)] - E [log q(\theta_d)] \\ + ( E [logp(\beta |\eta)] - E [log q(\beta)])/D \}\]

……(3)

注意,我们在对文档求和中引入了每个语料的项,可以通过文档D的数目进行划分。该步可以帮助我们对online VB算法求导。

接着,我们将上述期望展开成变分参数的函数。变分目标函数只依赖于$ n_dw $, 词w在文档中的出现次数。当使用VB时,文档可以由词数(word counts)来进行归纳:

\[\mathcal{L} = \sum_d \sum_w n_{dw} \sum_k \phi_{dwk}( E_q[log\theta_{dk}] + E[log\beta_{kw}] - log\phi_{dwk}) \\ - log\Gamma(\sum_k\gamma_{dk}) + \sum_k(\alpha-\gamma_{dk})E[log\theta_{dk}] + log\Gamma(\gamma_{dk}) \\+ (\sum_k -log\Gamma(\sum_w\lambda_{kw}) + \sum_w(\eta-\lambda_{kw})E[log\beta_{kw}] + log\Gamma(\lambda_{kw}))/D \\ + log\Gamma(K\alpha) - K log\Gamma(\alpha) + (log\Gamma(\mathcal{V}\eta)-\mathcal{V}log\Gamma(\eta))/D \\ \triangleq \sum_d \ell(n_d,\phi_d,\gamma_d,\lambda)\]

……(4)

其中V是词汇表的size,D是是文档数目。$ \ell(n_d,\phi_d,\gamma_d,\lambda) $ 表示文档d对ELBO的贡献度。

L可以使用在变分参数φ, γ, λ 上的坐标上升法进行优化:

\[\phi_{dwk} \propto exp \{ E_q[log\theta_{dk}] + E_q[log\beta_{kw}]\}\] \[\gamma_{dk}=\alpha + \sum_{w}n_{dw}\phi_{dwk}\] \[\lambda_{kw}=\eta+\sum_d n_{dw} \phi_{dwk}\]

……(5)

在q分布下logθ和logβ的期望为:

\[E[log\theta_{dk}]=\Phi(\gamma_{dk})- \Phi(\sum_{i=1}^K \gamma_{di})\] \[E[log\beta_{kw}]=\Phi(\gamma_{kw})- \Phi(\sum_{i=1}^V \gamma_{ki})\]

……(6)

其中Ψ表示digamma函数。

等式5的更新可以保证ELBO收敛到一个稳定点。通过EM算法(Expectation-Maximization)来完成,我们可以将这些更新划分到”E-step”:它会保持λ固定,迭代更新γ 和 φ 直到收敛,接着会进行”M-step”:由给定φ更新λ。实例上,如果在每个E-step上重新初始化γ 和 φ,该算法会收敛到一个更好的解。算法1就是batch VB for LDA:

Online算法

算法1具有常数量的内存开销,经验上,会比batch collpased Gibbs sampling收敛更快。然而,它仍需要在每次迭代时,将整个语料库传进去。因而,如果对于非常大的数据集会很慢,天然不适合不断到来的新数据。我们提出了一种在线变分推断算法来拟合λ,该参数是在主题分布β上的变分后验。我们的算法几科与batch VB算法一样简单,但对于大数据集收敛更快。

要想让主题参数λ的设置更好,就要让在算法1中的E-step中拟合得到每个文档的变分参数γ和φ之后,衡量ELBO的L要尽可能地高。将γ(nd, λ)和φ(nd, λ)设置为由E-step产生的γd 和 φd。我们的学习目标是设置λ,使下面目标最大化:

\(\mathcal{L}(n,\lambda) \triangleq \sum_d \ell (n_d,\gamma(n_d,\lambda), \phi(n_d,\lambda), \lambda)\) ……(7)

$ \ell(n_d,\gamma_d,\phi_d,\lambda) $是第d个文档对等式(4)的变分上界的贡献度。这类似于最小平方矩阵分解的目标函数,尽管LDA的ELBO比简单的squared loss function更不方便。

Online VB for LDA在算法2中描述。

第t个词频向量$n_t$被观察看,保持λ固定,我们执行一个E-step来找到局部最优解$\gamma_{t}$和$\phi_{t}$。接着,我们计算$\hat{\lambda}$,假如整个语料由单个文档$n_t$重复D次组成,那么λ的设置就是最优的(给定φt)。D是提供给算法的唯一文档数目,比如:一个语料的size。(在真实online的例子中:D->∞,对应于β的经验Bayes估计)。我们接着更新λ,使用一个λ的先前值和$\hat{\lambda}$进行加权平均得到。$\hat{\lambda}$的权重为:$\rho_t \triangleq (\tau_0 + t)^{-\kappa}$,其中$\kappa \in (0.5,1]$控制着$\hat{\lambda}$的旧值被遗忘的比率,τ0 ≥ 0会减缓算法的早期迭代。条件κ ∈ (0.5, 1]需要保证收敛。我们会在下一节展示,online LDA对应于在变分目标函数L上的一个随机自然梯度算法(stochastic natural gradient algorithm)。

该算法很像paper[16]中提出在online VB,在模型上有隐数据——最重要的区别是,我们使用一个近似的E-step来优化γt 和 φt, 因为我们不能精确计算条件分布$ p(z_t,\theta_t \mid \beta,n_t,\alpha) $。

Mini-batches: 在随机学习中的常用技术是,在每次更新时考虑多个观察(observations)以降噪。在online LDA中,这意味着计算$\hat{\lambda}$使用S>1的观察:

\(\hat{\lambda}_{kw}= \eta + \frac{D}{S} \sum_{S} n_{tsk} \phi_{tskw}\) ……(8)

其中$n_{ts}$是在mini-batch t中的第s个文档。该文档的变分参数$\phi_{ts}$和$\gamma_{ts}$是使用普通E-step进行拟合。注意,当S=D 和 κ = 0时,online VB就成了batch VB。

超参数估计:在batch 变分LDA中,超参数α 和 η的点估计,可以拟合给出的γ 和 λ,它使用线性时间的Newton-Raphson法。我们同样将α 和 η的更新并入到online LDA:

\[\alpha \leftarrow \alpha - \rho_{t} \hat{\alpha}(\gamma_t)\] \[\eta \leftarrow \rho_{t} \hat{\eta}(\lambda)\]

……(9)

其中,$\hat{\alpha}(\gamma_t)$是Hessian逆矩阵乘以梯度$\triangledown_{\alpha}l(n_t,\gamma_t,\phi_t,\lambda$,$ \hat{\eta}(\lambda))$是Hessian的逆矩阵乘以梯度$\triangledown_{\eta}\mathcal{L}$,$\rho_{t} \triangleq (\tau_0 + t)^{-\kappa}$。

5.3 收敛分析

此处不详说,见paper。