本文主要译自paper 1.

Abstract

神经网络概率语言模型(NPLM)与过去广泛使用的n-gram语言模型相互竞争,甚至偶尔优于后者。NPLM最大的缺点是训练和测试时间超长。Morin和Bengio提出了一种层次化语言模型,它会构建一棵关于词的二叉树,它比非层次化的使用专家知识的模型快2个数量级。我们引入了一种快速的层次化语言模型,它使用一种简单的基于特征的算法,可以从数据中自动构建word树。我们接着展示,产生的模型胜过非层次化神经网络模型、n-gram模型。

介绍

统计语言模型的关注点是,构建基于词序列的概率模型。这样的模型可以被用于区分可能的序列与不可能的序列,广泛用于语音识别,信息检索,机器翻译。大多数统计语言模型都基于Markov猜想:一个词的分布只依赖于 在它之前出现的固定数目的词。该猜想明显是错误的,但它很简洁方便,因为它将固定长度的词序列的概率分布建模问题,缩减成:给定前面固定数目的词(称为上下文:context),建模下一个词的分布。此处,我们使用:P(wn|w1:n-1)来表示分布,wn是下一个词, w1:n-1表示上下文(w1,..,wn-1)。

目前n-gram模型是最流行的统计语言模型,因为简单,并且性能很好。这些模型的条件概率表,通过统计训练数据中的n元组,并进行归一化进行估计。因为n元组的数据是n的指数级,对raw counts进行平滑会达到很好的性能。n-gram模型有许多平滑方法,详见paper 2. 尽管n-gram模型有许多复杂的平滑方法,n-gram模型很难利用大量上下文,因为数据的稀疏性问题很严重。这种现像的主要原因是,经典的n-gram模型是本质上有界的条件概率表,里面的条目都是相互独立的。这些模型不会利用这样的事实:即相似的词常出现在相似的上下文中,因为它们没有相似性的概率。基于分类的n-gram模型(paper 3)主要解决该问题,通过将词和上下文,基于使用模式聚类到分类中,使用这些分类信息来提升泛化。它会提升n-gram的性能,这种方法引入了严格死板的相似性,因为每个词很典型,都属于确定的某个类。

另一种可选的、更灵活的抵消数据稀疏性问题的方法是,将每个词使用一个real-valued的特征向量,它会捕获该特性,以便相似的上下文中的词,具有相似的特征向量。接着,下一个词的条件概率可以被建模成一个关于上下文词和下一个词的平滑函数。这种方法提供了自动平滑,因为对于一个给定的上下文,相似的词可以保证分配到相似的概率。同样的,相似的上下文现在也可能具有相似的表示,并会生成对下一个词相似的预测。大多数基于该方法的模型,使用一个前馈神经网络(feed-forwark NN),将上下文词汇的特征向量映射到下一个词的分布上。这种方法的可能最好模型是神经网络语言模型NPLM(paper 4),在100w词级别的语料上,它胜过n-gram模型。

层次化神经网络语言模型

NPLM的主要缺点是,这样的相类似模型,训练和测试时间很慢。因为下一个词的概率计算,需要显式对词汇表中所有词进行归一化,对于给定下一个词的概率计算开销,以及对于在下一个词之上的所有分布的计算开销,事实上两者几乎一样:它们会花费关于词汇表size的线性时间开销。由于在这样的模型上,计算精确的梯度,需要重复计算给定上下文的下一个词的概率,通过增加概率来更新模型参数,训练时间与词汇表size成线性关系。通常,自然语言数据集包含着上万的词汇,这意味着即使以最简单方式训练这种类NPLM模型,在实际中的计算开销仍过大。一种加速该过程的方法是,使用专门的重要性sampling过程,来逼近要学习所需的梯度(paper 5)。然而,该方法可以本质上加速训练时间,测试时间的开销依然很快。

我们引入了层次化NPLM(paper 6),对比普通的NPLM,它在训练和测试的时间复杂度均做到了指数级的衰减。通过二叉树替代NPLM中非结构化的词汇表,可以表示成一个层次化的词汇表词簇。每个词对应于树上的叶子节点,可以由顶点到叶子节点的路径唯一确定。如果N是词汇表的词数,并且树是一棵平衡树,任何词都可以由一串O(log N)二分决策来指定,它会指定两个子节点,当前节点可以访问到下一节点。该过程将N-way选择替换成一串O(log N)的二分选择。在概率术语中,N-way归一化,可以替换成一串O(log N)的局部(二分)归一化。结果上,词汇表中的词分布,可以通过提供访问左子节点的概率来指定。在层次化NPLM中,这样NPLM方式的局部概率的计算:采用一个特征向量作为上下文词汇,同时给当前节点的一个特征向量作为输入。下一个词的概率,由该词所对应路径的二分决策的概率指定。

当数据集包含上百万的词时,该模型的表现优于基于分类的3-gram,但比paper 6中的NPLM表现差。这种层次化NPLM模型,比普通的NPLM快2个数量级。这种方法的主要限制,主要是用于构建word树的过程。该树可以从WordNet IS-A分类体系开始,通过结合人工和数据驱动处理,并将它转换成一个二叉树。我们的目标是,将该过程替换成从训练数据中自动构建树,不需要任何专家知识。我们也探索了使用树(里面的词汇至少出现一次)的性能优点。

log-bilinear model

我们使用log-bilinear语言模型(LBL:Paper 7)作为我们层次化模型的基础,因为它的性能十分出色,并且很简洁。类似于所有神经网络语言模型,LBL模型将每个词表示成一个real-valued的特征向量。我们将词w的特征向量表示成:,所有词的特征向量组成矩阵R. 模型的目标:给定上下文w1:n-1,预测下一个词wn。我们将要预测下一个词的的特征向量表示成,它可以通过对上下文词汇的特征向量进行线性组合得到:

公式一:

其中,Ci是权重矩阵,它与位置i的上下文有关。接着,要预测的特征向量,和词汇表中每个词的特征向量,两者间的相似度通过内积计算。该相似度接着进行指数化,归一化,从而获取下一个词的分布:

公式二:

这里的bw是词w的偏置bias,它被用于捕获上下文独立的词频。

注意,LBL模型可以被解释成一种特殊的前馈神经网络,它只有一个线性隐层,以及一个softmax输出层。该网络的输入是上下文词汇的特征向量,而从隐层到输出层的权重矩阵,则可以简单通过特征向量矩阵R指定。该向量隐单元的激活,对应于下一个词的预测特征向量。不同于NPLM,LBL模型需要计算隐层激活,每次预测只计算一次,隐层不存在非线性。虽然LBL模型很简单,但表现很好,在相当大的数据集上,优于NPLM和n-gram模型。

4.层次化log-bilinear模型

我们的层次化语言模型基于该模型. 该模型的特点是,使用log-bilinear语言模型来计算每个节点的概率,并处理树中每个词多次出现率。注意,使用多个词出现率的思想,由papar 6提出,但它并没有实现。

HLBL的第一个组件是一棵关于词作为叶子的二叉树。接着,我们假设词汇表中每个词对应一个叶子结点。接着,每个词可以通过从树顶点到叶子节点的路径唯一指定。该路径本身可以编码成一串二进制字符串d(每个节点会做二分决策),对应于:访问当前节点的左子节点。例如,字符串”10”表示的是:从顶点,到左子节点,然后到右子节点。这样,每个词可以被表示成一串二进制码。

HLBL模型的第二个组件是,每个节点会做出二分决策的概率模型,在我们的case中,使用了一个LBL模型的改进版本。在HLBL模型中,和非层次化的LBL类似,上下文词汇表示成real-valued特征向量。树上的每个非叶子节点,也具有一个特征向量与它相关联,它用于区分词是在该节点的左子树还是在右子树。不同于上下文词汇,被预测的词,被表示成使用它们的二进制码。然而,这种表示相当灵活,因为编码中的每种二进制数字,相当于该节点做出的决策,它依赖于节点的特征向量。

在HLBL模型中,给定上下文,下一个词是w的概率,就是由该词编码指定的一串二分决策的概率。因为在一个节点上做出一个决策的概率,只取决于要预测的特征向量(由上下文决定),以及该节点的特征向量,我们可以将下一个词的概率表示成二分决策的概率乘积:

公式三:

di是词汇w的第i位数字编码,qi是编码所对应路径上第i个节点的特征向量。每个决策的概率为:

公式四:

其中,是logistic函数,而是使用公式一计算得到的要预测的特征向量,在该公式中的bi是该节点的偏置bias,它可以捕获当离开该节点到它的左子节点上的上下文独立的趋势。

D(w)表示词汇w所对应的一串编码集合。每个词允许多个编码,可以更好地预测词,具有多种含义或者多种使用模式。每个词使用多种编码,也可以使它更容易将多种不同的词层次结构,组合成单一的结构。这也反映了一个事实:没有单一的层次结构可以表示词之间的所有关系。

使用LBL模型(而非NPLM)来计算局部概率,允许我们避免计算隐层的非线性,这可以让我们的层次模型,比层次化NPLM模型更快作出预测。更重要的是,对于O(log N)个决策中的每一个,层次化NPLM都需要计算一次隐层激活,而HLBL模型只在每次预测时计算一次要预测的特征向量。然而,在LBL模型中,单个二分决策的概率计算的时间复杂性,仍然是特征向量维度D的二次方,这会使高维特征向量的计算开销很高昂。我们可以使时间复杂度与D是线性关系,通过将权重矩阵Ci限制为对角阵。注意,对于上下文size=1的情况,该限制条件不会减小模型表示的幂(power),我们相信,这种损失(loss),多于可以通过使用更复杂的树更快的训练时间的补偿。

HLBL模型可以通过最大化(带罚项的)log似然进行训练。因为下一个词的概率只依赖于上下文权重、上下文词汇的特征向量、以及从顶点到词汇所在叶子节点路径上的节点的特征向量,在每个训练case上,只有一小部分(log级别)的参数需要更新。

参考

介绍

word2vec中的logistic function计算使用的是快速算法,初始先算好值,然后在神经网络训练迭代中直接查表加快计算和训练速度。设计也还算精巧。

代码

在初始化时,有代码片段一,会对expTable做专门的初始化:

#define EXP_TABLE_SIZE 1000
#define MAX_EXP 6



// step 4: 分配logistic查表.
expTable = (real *)malloc((EXP_TABLE_SIZE + 1) * sizeof(real));
  
// 初始化: 预先计算好指数运算表. 
for (i = 0; i < EXP_TABLE_SIZE; i++) {
    expTable[i] = exp((i / (real)EXP_TABLE_SIZE * 2 - 1) * MAX_EXP); // Precompute the exp() table
    expTable[i] = expTable[i] / (expTable[i] + 1);                   // Precompute f(x) = x / (x + 1)
  }

expTable数据的大小为1000,存了1000个值用于查表。

而在模型训练阶段,代码片段二:

// 前向传播: f=∑ neu1*syn1
// Propagate hidden -> output
for (c = 0; c < layer1_size; c++) {
  f += neu1[c] * syn1[c + l2];
}

// 范围判断,查表
if (f <= -MAX_EXP) 
  continue;
else if (f >= MAX_EXP) 
  continue;
else 
  f = expTable[(int)((f + MAX_EXP) * (EXP_TABLE_SIZE / MAX_EXP / 2))];

将代码片段二代入到代码片段一,我们可以发现很神奇的现象:

将上式简单做个替换,把f替成z,整体输入看成i:

而再做一次运算,即知该表的值刚好与 logit函数 y = 1/(1+e^-z) 相同。

为什么要这么设计呢?

通过查表法,当然是为了加快指数运算。如果在每次迭代时,都去做指数运算,cpu开销蛮大。这样的设计,可以优化训练时间。

我来来看前前一段代码中,expTable[0 … EXP_TABLE_SIZE]中具体对应的值是哪些?这里,我提供了它的python实现:代码下载 ,绘制出实际的取值曲线(logistic regression的某一段取值),详见上面链接提供的代码:

另外,还有一点,就是这样的分割:EXP_TABLE_SIZE=1000, MAX_EXP=6? 把数字代入,即知:

为什么要做这样的trick进行分解呢?上图展示的是从输入z的角度去理解。要真正理解为什么这样取,还要结合实际代码中Hierarchical softmax代码里的情况去理解。

里面的f输入范围在(-MAX_EXP, MAX_EXP), 即(-6, 6)。

第一阶段:(f + MAX_EXP)/MAX_EXP/2 * EXP_TABLE_SIZE 所做的操作:

  • f + MAX_EXP: 平移,(0, 12)
  • (f + MAX_EXP)/MAX_EXP: 压缩,(0,2)
  • (f + MAX_EXP)/MAX_EXP/2: 归一化,(0,1)
  • (f + MAX_EXP)/MAX_EXP/2 * EXP_TABLE_SIZE: 即将归一化后的f映射到[0,1000)个expTable的格子里,即是第一式的输入i

第二阶段:(i / (real)EXP_TABLE_SIZE * 2 - 1) * MAX_EXP 所做的操作:

  • i / EXP_TABLE_SIZE: 将上面第一阶段的i,重新压缩到(0,1)
  • i / (real)EXP_TABLE_SIZE * 2: 拉伸成(0,2)
  • i / (real)EXP_TABLE_SIZE * 2 - 1: 平移为(-1,1)
  • (i / (real)EXP_TABLE_SIZE * 2 - 1) * MAX_EXP: 放大为(-6,6)

ok,这时候,我们就会明白,为啥取(-6,6)了。因为:

  • logit(6)=0.997527376843
  • logit(-6)=0.00247262315663

这个范围内的值足够了。

但如果你细发现,两步合在一起,我们发现仿佛又回到了原点,于是可能开始会怀疑,为什么不干脆直接用f,省去中间这么多的变换,直接原封不动地输入expTable呢?

比如类似这样的实现:

def logit(x):
    e = math.exp(x)
    return e/(e+1)

我们可以推导下,如果让你设计一个需要1000个值的数组,存储logit函数值。假设数组名为expTable,即:

  • 输入x为(0,1000)
  • expTable[x]的输出值对应概率值(0,1)

由于logit(0)=0.5,常用的一个实现是:

从0开始刚好为0,输入(0,+inf),输出(0,1)。如果想让x刚好对应索引,刚在原基础上除以1000即可。即:

再在原基础上做一下范围控制,就是我们上面的:

如果希望得到更精准的值,将EXP_TABLE_SIZE和MAX_EXP调大些即可以。

介绍

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方法:训练集中的每个词wi,以下面公式计算得到的概率进行抛弃:

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

具体实现

有道之前的中,提到的该subsampling描述也不准确。在当中的描述是:

而实际中,采用的是:

部分代码:

// 进行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上,点击下载

下面提供了三种方法,最终生成的图可以看下。对于前两种方法,基本能做到与原先词频正相关,最后使用时,需要设置一个阀值,砍掉高频词。而最后一种方法,效果也不错(虽然偶有会存留高频词,或者低频词也同样被砍掉)。而word2vec采用的正是第3种方法(大于0的采样点会被抛弃掉)。

参考:

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编码

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的更新等,此处不再解释。详见源码。