Spark2.0机器学习系列之8：聚类分析(K-Means,Bisecting K-Means,LDA，高斯混合模型）

在写这篇文章之前，先说一些题外话。

许多机器学习算法（如后面将要提到的LDA）涉及的数学知识太多，前前后后一大堆，理解起来不是那么容易。

面对复杂的机器学习模型，尤其是涉及大量数学知识的模型，我们往往要花费大量的时间和精力去推导数学算法（公式），如果过分沉湎于此会忽略了很多背后也许更重要的东西，正所谓只见树木，不见森林，而这是缺乏远见，是迷茫的。

我们需要深入理解模型背后的逻辑和所蕴含的或简或繁的思想。某些思想甚至可能是很美的思想，很伟大的思想。这些理解，使得面对复杂的问题时候，面对陌生问题时，会左右我们的选择。我们会首先想起谁，谁是我们工具箱里最得心应手的工具，就像自己的孩子一样亲切。而有些方法我们注定不喜欢，也用不太懂。毕竟我们不是每个人都会去写复杂的机器学习算法code，尤其是复杂的分布式code，但也不意味着我们不需要关注模型的数学原理，因为没这些数学原理的支撑，也许你永远理解不了这个模型。但是仅仅掌握公式的推导甚至从头到尾能写出完整代码，也还是远远不够的。比如说我们都知道为了控制over-fitting，需要正则化，为什么在看似很精确的公式里面硬生生的多加入一个正则项，就能够达到控制over-fitting的目的呢？又如何判断所选择的正则化参数是合理的？不大也不小，刚好合适。公式我们都能完整的推导吧，但也许有些事我们想的还是不那么明白，用的时候自然也就会有一些疑虑。

所以呢，要做到深入浅出真的很难，这也是我们学习ML需要努力的方向，理解这些ML模型背后的逻辑，所蕴含的美，所表达的的思想，甚至是哲理，也许能助你我成长到一个新的高度。

好了，言归正传。在Spark2.0版本中（不是基于RDD API的MLlib），共有四种聚类方法：

（1）K-means

（2）Latent Dirichlet allocation (LDA)

（3）Bisecting k-means（二分k均值算法）

（4）Gaussian Mixture Model (GMM)。

基于RDD API的MLLib中，共有六种聚类方法：

（1）K-means

（2）Gaussian mixture

（3）Power iteration clustering (PIC)

（4）Latent Dirichlet allocation (LDA)**

（5）Bisecting k-means

（6）Streaming k-means

多了Power iteration clustering (PIC)和Streaming k-means两种。

本文将介绍LDA主题建模，其它方法在我的整个系列中都会有介绍。

什么是LDA主题建模？

隐含狄利克雷分配（LDA，Latent Dirichlet Allocation）是一种主题模型(Topic Model，即从所收集的文档中推测主题)。甚至可以说LDA模型现在已经成为了主题建模中的一个标准，是实践中最成功的主题模型之一。那么何谓“主题”呢？，就是诸如一篇文章、一段话、一个句子所表达的中心思想。不过从统计模型的角度来说，我们是用一个特定的词频分布来刻画主题的，并认为一篇文章、一段话、一个句子是从一个概率模型中生成的。也就是说在主题模型中，主题表现为一系列相关的单词，是这些单词的条件概率。形象来说，主题就是一个桶，里面装了出现概率较高的单词(参见下面的图)，这些单词与这个主题有很强的相关性。

LDA可以用来识别大规模文档集（document collection）或语料库（corpus）中潜藏的主题信息。它采用了词袋（bag of words）的方法，这种方法将每一篇文档视为一个词频向量，从而将文本信息转化为了易于建模的数字信息。但是词袋方法没有考虑词与词之间的顺序，这简化了问题的复杂性，同时也为模型的改进提供了契机。每一篇文档代表了一些主题所构成的一个概率分布，而每一个主题又代表了很多单词所构成的一个概率分布。

LDA可以被认为是如下的一个聚类过程：

（1）各个主题（Topics）对应于各类的“质心”，每一篇文档被视为数据集中的一个样本。

（2）主题和文档都被认为存在一个向量空间中，这个向量空间中的每个特征向量都是词频（词袋模型）

（3）与采用传统聚类方法中采用距离公式来衡量不同的是，LDA使用一个基于统计模型的方程，而这个统计模型揭示出这些文档都是怎么产生的。

Latent Dirichlet allocation (LDA) is a topic model which infers topics from a collection of text documents. LDA can be thought of as a clustering algorithm as follows:

（1）Topics correspond to cluster centers, and documents correspond to examples (rows) in a dataset.

（2）Topics and documents both exist in a feature space, where feature vectors are vectors of word counts (bag of words).

（3）Rather than estimating a clustering using a traditional distance, LDA uses a function based on a statistical model of how text documents are generated.

下面的几段文字来源于：http://www.tuicool.com/articles/reaIra6

它基于一个常识性假设：文档集合中的所有文本均共享一定数量的隐含主题。基于该假设，它将整个文档集特征化为隐含主题的集合，而每篇文本被表示为这些隐含主题的特定比例的混合。

LDA的这三位作者在原始论文中给了一个简单的例子。比如给定这几个主题：Arts、Budgets、Children、Education，在这几个主题下，可以构造生成跟主题相关的词语，如下图所示：

然后可以根据这些词语生成如下图所示的一篇文章（其中不同颜色的词语分别对应上图中不同主题下的词）

表面上理解LDA比较简单，无非就是：当看到一篇文章后，我们往往喜欢推测这篇文章是如何生成的，我们可能会认为某个作者先确定这篇文章的几个主题，然后围绕这几个主题遣词造句，表达成文。

前面说了这么多，在推导模型前，总结几条核心思想：

（1）隐含主题，形象的说就是一个桶，里面装了出现概率较高的单词，从聚类的角度来说，各个主题（Topics）对应于各类的“质心”，主题和文档都被认为存在于同一个词频向量空间中。（2）在文档集合中的所有文本均共享一定数量的隐含主题的假设下，我们将寻找一个基于统计模型的方程。

LDA的核心公式如下：

p(w|d)=∑k=1Kp(w|zk)?p(zk|d)

上面式子的左边，就是文档的词频，是很容易统计得到的。如果一篇文章中有很多个词语，那么就有很多个等式了。再如果我们收集了很多的文档，那么就有更多的等式了。这时候就是一个矩阵了，等式左边的矩阵是已知的，右边其实就是我们要求解的目标-与隐含主题相关，图示如下：

LDA建模的过程及算法优化

看一段很概述性的文字：

LDA是一个层次贝叶斯模型，把模型的参数也看作随机变量，从而可以引入控制参数的参数，实现彻底的“概率化”。LDA模型的Dirichlet的先验分布，文档d上主题的多项式分布。目前，参数估计是LDA最重要的任务，主要有两种方法：Gibbs抽样法（计算量大，但相对简单和精确）和变分贝叶斯推断法（计算量小，精度度弱）。

看来关于LDA的数学求解，还有很多知识要学习，本文篇幅太长了，我将再写一篇文章详细介绍贝塔分布，多项式、Dirichlet的先验分布、Gibbs抽样法，参数估计等相关的知识。更好的建议是阅读rickjin科普文章《LDA数学八卦》，仔细看完这个，那么就没有什么不明白的了。

应该说掌握了上述概率相关的知识以后，不难理解下面：LDA模型中一篇文档生成的方式：

(1)从狄利克雷分布α中取样生成文档i的主题分布θi

(2)从主题的多项式分布θi中取样生成文档i第j个词的主题zi,j

(3)从狄利克雷分布β中取样生成主题zi,j的词语分布?zi,j

(4)从词语的多项式分布?zi,j中采样最终生成词语wi,j

至此为之，我们要去考虑，怎么去计算这两个矩阵，怎么去优化的问题了。Spark采用的两种优化算法：

（1）EMLDAOptimizer 通过在likelihood函数上计算最大期望EM，提供较全面的结果。

（2）OnlineLDAOptimizer 通过在小批量数据上迭代采样实现online变分推断，比较节省内存。在线变分预测是一种训练LDA模型的技术，它以小批次增量式地处理数据。由于每次处理一小批数据，我们可以轻易地将其扩展应用到大数据集上。MLlib按照 Hoffman论文里最初提出的算法实现了一种在线变分学习算法。

LDA supports different inference algorithms via setOptimizer function. EMLDAOptimizer learns clustering using expectation-maximization on the likelihood function and yields comprehensive results, while OnlineLDAOptimizer uses iterative mini-batch sampling for online variational inference and is generally memory friendly.

一些数学铺垫

参数估计（parameter estimation）

用系统理论的语言来看：在已知系统模型结构时，用系统的输入和输出数据计算系统模型参数的过程。18世纪末德国数学家C.F.高斯首先提出参数估计的方法，他用最小二乘法计算天体运行的轨道。20世纪60年代，随着电子计算机的普及，参数估计有了飞速的发展。参数估计有多种方法：矩估计、极大似然法、一致最小方差无偏估计、最小风险估计、同变估计、最小二乘法、贝叶斯估计、极大验后法、最小风险法和极小化极大熵法等。最基本的方法是最小二乘法和极大似然法（来自百度百科，百度百科倒是经常贴“高大上”的描述方法，但是有些东西怎么看却都是前言不接后语，连个公式渲染都没有，好歹简单学学wikipedia，让中国人上网查个数学公式贼费劲，还乱七八糟的都不标准不统一，吐槽一下，呵呵）。

从系统理论角度出发比较“高大上”，说的也比较抽象，但是考虑到任何机器学习模型也可以理解为一个系统，上面那段话好好理解一下也挺好，我们可以提取出下面一些理解：

（1）参数估计并不是数理统计里面一个特有的术语，一上来就看数理统计里面的参数估计，就太狭隘了，也就没有那么深入的理解。但是也不可否认，参数估计很多方法都是数理统计的方法，但也不全是，如著名的最小二乘法（在那里都能见到这哥们的影子）。

（2）既然都开始估计参数了，那么模型的结构肯定是已知的，也就是说：机器学习选取什么算法，决策树？还是SVM？还是Beyes概率的，深度学习的？，你自己其实已经有主意了的，手头上有一些样本数据（观测值），任务是怎么去求解（优化）系统（模型）的参数。

正是由于最小二乘法和极大似然法也是ML里面最基本的使用广泛的参数估计方法，我们需要重点掌握。

另外从数理统计学的语言讲：参数估计是根据从总体中抽取的样本估计总体分布中包含的未知参数的方法。人们常常需要根据手中的数据，分析或推断数据反映的本质规律。即根据样本数据如何选择统计量去推断总体的分布或数字特征等。统计推断是数理统计研究的核心问题。

可以看出：

(1) 正是由于我们往往只能获得一少部分样本，却被要求拿来估计复杂系统的真实参数，所以数学家才会根据不同的模型，应用于不同的数据情况下，设计出如此多的巧妙办法，尽量使估计“无偏，一致，有效”。

（2）数理统计学里面的参数估计目的更加明确，即推断总体的分布，如总体概率密度函数的类型及参数，总体样本的统计参数（均值、方差等等）。显然机器学习的目标更加远大，我们需要估计各类复杂模型的参数，当然也包括一部分以概率统计理论为基础的ML模型（比如本文后面要介绍的GMM和LDA）。

联合概率分布 Joint probability distribution

首先推荐一篇博文：zouxy09的CSDN 在里面举例的帮助下，比较容易看懂很多概念：

从最大似然到EM算法浅解 http://blog.csdn.net/zouxy09/article/details/8537620/

概率方法的很多公式都很长，看起来挺吓人，很多时候其实都是由于采用的是联合概率分布，对他的形式把握透彻了，这些公式其实没有那么难理解。

很多时候（比如说朴素Beyes分类器），会简化问题，首先做独立性假设，虽然可能并不完全满足独立性假设，也都这样做，而且发现有很好的鲁棒性。

独立性假设就是：一个事件发生概率不受其他事件影响的事件.两(多)个独立事件同时发生的概率等于两(多)者分别发生的概率之积.假如A、B互相独立,那么A、B同时发生的概率是P(AB)=P(A)×P(B)

如果我们抽取了N个样本，那么可以认为是发生了N次独立事件，所以有下面比较严谨的表述。

设x1,x2,...,xN是从概率密度函数p(x;θ)中随机抽取的N个样本，从而得到联合概率密度函数p(X;θ),其中X={x1,x2,...,xN}是样本集。假设不同样本之间具有统计的独立性，则有：

p(X;θ)=p(x1,x2,...,xN;θ)=∏k=1k=Np(xk,θ)，θ∈Θ

p(xk,θ)代表概率密度函数的参数为θ时，xk发生的概率。这里的参数θ是一个向量（参数向量）。当然θ是可以变化（不断变化θ代表不同的概率密度分布函数，只是具有相同的类型），因此θ是参数向量集合的Θ的一个元素。

最大似然估计

未完待续…

Spark 代码分析、参数设置及结果评价

SPARK中可选参数

（1）K：主题数量（或者说聚簇中心数量）

（2）maxIterations：EM算法的最大迭代次数，设置足够大的迭代次数非常重要，在早期的迭代过程中，EM通常包含无用的主题，但是随着迭代次数的增多这种情况会得到显著的改善，至少需要设置20次的迭代，50-100次是更合理的设置，取决于你的数据集。

（3）docConcentration（Dirichlet分布的参数α)：文档在主题上分布的先验参数（超参数α)。当前必须大于1，值越大，推断出的分布越平滑。默认为-1，自动设置。

（4）topicConcentration（Dirichlet分布的参数β)：主题在单词上的先验分布参数。当前必须大于1，值越大，推断出的分布越平滑。默认为-1，自动设置。

（5）checkpointInterval：检查点间隔。maxIterations很大的时候，检查点可以帮助减少shuffle文件大小并且可以帮助故障恢复。

SPARK中模型的评估

（1）perplexity是一种信息理论的测量方法，b的perplexity值定义为基于b的熵的能量（b可以是一个概率分布，或者概率模型），通常用于概率模型的比较

LDA takes in a collection of documents as vectors of word counts and the following parameters (set using the builder pattern):

（1）k: Number of topics (i.e., cluster centers)

（2）：optimizer: Optimizer to use for learning the LDA model, either EMLDAOptimizer or OnlineLDAOptimizer

（3）docConcentration: Dirichlet parameter for prior over documents’ distributions over topics. Larger values encourage smoother inferred distributions.

（4）topicConcentration: Dirichlet parameter for prior over topics’ distributions over terms (words). Larger values encourage smoother inferred distributions.

（5）：maxIterations: Limit on the number of iterations. It is important to do enough iterations. In early iterations, EM often has useless topics, but those topics improve dramatically after more iterations. Using at least 20 and possibly 50-100 iterations is often reasonable, depending on your dataset.

（6）checkpointInterval: If using checkpointing (set in the Spark configuration), this parameter specifies the frequency with which checkpoints will be created. If maxIterations is large, using checkpointing can help reduce shuffle file sizes on disk and help with failure recovery.

All of spark.mllib’s LDA models support:

（1）describeTopics: Returns topics as arrays of most important terms and term weights

（2）topicsMatrix: Returns a vocabSize by k matrix where each column is a topic

详细代码注释

import org.apache.spark.sql.SparkSession
import org.apache.log4j.{Level, Logger}
import org.apache.spark.ml.clustering.LDA

object myClusters {
  def main(args:Array[String]){

    //屏蔽日志
    Logger.getLogger("org.apache.spark").setLevel(Level.ERROR)
    Logger.getLogger("org.eclipse.jetty.server").setLevel(Level.OFF)

    val warehouseLocation = "/Java/Spark/spark-warehouse"
    val spark=SparkSession
            .builder()
            .appName("myClusters")
            .master("local[4]")
            .config("spark.sql.warehouse.dir",warehouseLocation)
            .getOrCreate();   

    val dataset_lpa=spark.read.format("libsvm")
                          .load("/spark-2.0.0-bin-hadoop2.6/data/mllib/sample_lda_libsvm_data.txt")
    //------------------------------------1 模型训练-----------------------------------------
    /**
     * k: 主题数，或者聚类中心数
     * DocConcentration：文章分布的超参数(Dirichlet分布的参数)，必需>1.0，值越大，推断出的分布越平滑
     * TopicConcentration：主题分布的超参数(Dirichlet分布的参数)，必需>1.0，值越大，推断出的分布越平滑
     * MaxIterations：迭代次数，需充分迭代，至少20次以上
     * setSeed：随机种子
     * CheckpointInterval：迭代计算时检查点的间隔
     * Optimizer：优化计算方法，目前支持"em", "online" ，em方法更占内存，迭代次数多内存可能不够会抛出stack异常
     */
    val lda=new LDA()
                .setK(3)
                .setTopicConcentration(3)
                .setDocConcentration(3)
                .setOptimizer("online")
                .setCheckpointInterval(10)
                .setMaxIter(100)

    val model=lda.fit(dataset_lpa)   

    /**生成的model不仅存储了推断的主题，还包括模型的评价方法。*/
    //---------------------------------2 模型评价-------------------------------------

    //模型的评价指标：ogLikelihood，logPerplexity
    //（1）根据训练集的模型分布计算的log likelihood，越大越好。
    val ll = model.logLikelihood(dataset_lpa)

    //（2）Perplexity评估，越小越好
    val lp = model.logPerplexity(dataset_lpa)

    println(s"The lower bound on the log likelihood of the entire corpus: $ll")
    println(s"The upper bound bound on perplexity: $lp")

    //---------------------------------3 模型及描述------------------------------
    //模型通过describeTopics、topicsMatrix来描述 

    //（1）描述各个主题最终的前maxTermsPerTopic个词语（最重要的词向量）及其权重
    val topics=model.describeTopics(maxTermsPerTopic=2)
    println("The topics described by their top-weighted terms:")
    topics.show(false)

    /**主题    主题包含最重要的词语序号                     各词语的权重
        +-----+-------------+------------------------------------------+
        |topic|termIndices  |termWeights                               |
        +-----+-------------+------------------------------------------+
        |0    |[5, 4, 0, 1] |[0.21169509638828377, 0.19142090510443274]|
        |1    |[5, 6, 1, 2] |[0.12521929515791688, 0.10175547561034966]|
        |2    |[3, 10, 6, 9]|[0.19885345685860667, 0.18794498802657686]|
        +-----+-------------+------------------------------------------+
     */    

    //（2） topicsMatrix: 主题-词分布，相当于phi。
    val topicsMat=model.topicsMatrix
    println("topicsMatrix")
    println(topicsMat.toString())
     /**topicsMatrix
        12.992380082908886  0.5654447550856024  16.438154549631257
        10.552480038361052  0.6367807085306598  19.81281695100224
        2.204054885551135   0.597153999004713   6.979803589429554
    *
    */

    //-----------------------------------4 对语料的主题进行聚类---------------------
    val topicsProb=model.transform(dataset_lpa)
    topicsProb.select("label", "topicDistribution")show(false)

    /** 文档序号 各主题的权重
        +-----+--------------------------------------------------------------+
        |label|topicDistribution                                             |
        +-----+--------------------------------------------------------------+
        |0.0  |[0.523730754859981,0.006564444943344147,0.46970480019667477]  |
        |1.0  |[0.7825074858166653,0.011001204994496623,0.206491309188838]   |
        |2.0  |[0.2085069748527087,0.005698459472719417,0.785794565674572]   |
    */
  }
}

对参数进行调试

//对迭代次数进行循环
for(i<-Array(5,10,20,40,60,120,200,500)){
    val lda=new LDA()
                .setK(3)
                .setTopicConcentration(3)
                .setDocConcentration(3)
                .setOptimizer("online")
                .setCheckpointInterval(10)
                .setMaxIter(i)
    val model=lda.fit(dataset_lpa) 

    val ll = model.logLikelihood(dataset_lpa)
    val lp = model.logPerplexity(dataset_lpa)

    println(s"$i $ll")
    println(s"$i $lp")
 }

可以得到如下的结果：logPerplexity在减小，LogLikelihood在增加，最大迭代次数需要设置50次以上，才能收敛：