Repository: spark
Updated Branches:
  refs/heads/master f377431a5 -> 986977340


SPARK-5400 [MLlib] Changed name of GaussianMixtureEM to GaussianMixture

Decoupling the model and the algorithm

Author: Travis Galoppo <[email protected]>

Closes #4290 from tgaloppo/spark-5400 and squashes the following commits:

9c1534c [Travis Galoppo] Fixed invokation instructions in comments
d848076 [Travis Galoppo] SPARK-5400 Changed name of GaussianMixtureEM to 
GaussianMixture to separate model from algorithm


Project: http://git-wip-us.apache.org/repos/asf/spark/repo
Commit: http://git-wip-us.apache.org/repos/asf/spark/commit/98697734
Tree: http://git-wip-us.apache.org/repos/asf/spark/tree/98697734
Diff: http://git-wip-us.apache.org/repos/asf/spark/diff/98697734

Branch: refs/heads/master
Commit: 986977340d0d02dbd0346bd233dbd93b8c8e74c9
Parents: f377431
Author: Travis Galoppo <[email protected]>
Authored: Fri Jan 30 15:32:25 2015 -0800
Committer: Xiangrui Meng <[email protected]>
Committed: Fri Jan 30 15:32:25 2015 -0800

----------------------------------------------------------------------
 .../examples/mllib/DenseGaussianMixture.scala   |  67 +++++
 .../spark/examples/mllib/DenseGmmEM.scala       |  67 -----
 .../mllib/clustering/GaussianMixture.scala      | 251 +++++++++++++++++++
 .../mllib/clustering/GaussianMixtureEM.scala    | 251 -------------------
 .../GMMExpectationMaximizationSuite.scala       |  83 ------
 .../mllib/clustering/GaussianMixtureSuite.scala |  83 ++++++
 6 files changed, 401 insertions(+), 401 deletions(-)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/spark/blob/98697734/examples/src/main/scala/org/apache/spark/examples/mllib/DenseGaussianMixture.scala
----------------------------------------------------------------------
diff --git 
a/examples/src/main/scala/org/apache/spark/examples/mllib/DenseGaussianMixture.scala
 
b/examples/src/main/scala/org/apache/spark/examples/mllib/DenseGaussianMixture.scala
new file mode 100644
index 0000000..df76b45
--- /dev/null
+++ 
b/examples/src/main/scala/org/apache/spark/examples/mllib/DenseGaussianMixture.scala
@@ -0,0 +1,67 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *    http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.spark.examples.mllib
+
+import org.apache.spark.{SparkConf, SparkContext}
+import org.apache.spark.mllib.clustering.GaussianMixture
+import org.apache.spark.mllib.linalg.Vectors
+
+/**
+ * An example Gaussian Mixture Model EM app. Run with
+ * {{{
+ * ./bin/run-example mllib.DenseGaussianMixture <input> <k> <convergenceTol>
+ * }}}
+ * If you use it as a template to create your own app, please use 
`spark-submit` to submit your app.
+ */
+object DenseGaussianMixture {
+  def main(args: Array[String]): Unit = {
+    if (args.length < 3) {
+      println("usage: DenseGmmEM <input file> <k> <convergenceTol> 
[maxIterations]")
+    } else {
+      val maxIterations = if (args.length > 3) args(3).toInt else 100
+      run(args(0), args(1).toInt, args(2).toDouble, maxIterations)
+    }
+  }
+
+  private def run(inputFile: String, k: Int, convergenceTol: Double, 
maxIterations: Int) {
+    val conf = new SparkConf().setAppName("Gaussian Mixture Model EM example")
+    val ctx  = new SparkContext(conf)
+    
+    val data = ctx.textFile(inputFile).map { line =>
+      Vectors.dense(line.trim.split(' ').map(_.toDouble))
+    }.cache()
+      
+    val clusters = new GaussianMixture()
+      .setK(k)
+      .setConvergenceTol(convergenceTol)
+      .setMaxIterations(maxIterations)
+      .run(data)
+    
+    for (i <- 0 until clusters.k) {
+      println("weight=%f\nmu=%s\nsigma=\n%s\n" format 
+        (clusters.weights(i), clusters.gaussians(i).mu, 
clusters.gaussians(i).sigma))
+    }
+    
+    println("Cluster labels (first <= 100):")
+    val clusterLabels = clusters.predict(data)
+    clusterLabels.take(100).foreach { x =>
+      print(" " + x)
+    }
+    println()
+  }
+}

http://git-wip-us.apache.org/repos/asf/spark/blob/98697734/examples/src/main/scala/org/apache/spark/examples/mllib/DenseGmmEM.scala
----------------------------------------------------------------------
diff --git 
a/examples/src/main/scala/org/apache/spark/examples/mllib/DenseGmmEM.scala 
b/examples/src/main/scala/org/apache/spark/examples/mllib/DenseGmmEM.scala
deleted file mode 100644
index de58be3..0000000
--- a/examples/src/main/scala/org/apache/spark/examples/mllib/DenseGmmEM.scala
+++ /dev/null
@@ -1,67 +0,0 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements.  See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License.  You may obtain a copy of the License at
- *
- *    http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-package org.apache.spark.examples.mllib
-
-import org.apache.spark.{SparkConf, SparkContext}
-import org.apache.spark.mllib.clustering.GaussianMixtureEM
-import org.apache.spark.mllib.linalg.Vectors
-
-/**
- * An example Gaussian Mixture Model EM app. Run with
- * {{{
- * ./bin/run-example org.apache.spark.examples.mllib.DenseGmmEM <input> <k> 
<covergenceTol>
- * }}}
- * If you use it as a template to create your own app, please use 
`spark-submit` to submit your app.
- */
-object DenseGmmEM {
-  def main(args: Array[String]): Unit = {
-    if (args.length < 3) {
-      println("usage: DenseGmmEM <input file> <k> <convergenceTol> 
[maxIterations]")
-    } else {
-      val maxIterations = if (args.length > 3) args(3).toInt else 100
-      run(args(0), args(1).toInt, args(2).toDouble, maxIterations)
-    }
-  }
-
-  private def run(inputFile: String, k: Int, convergenceTol: Double, 
maxIterations: Int) {
-    val conf = new SparkConf().setAppName("Gaussian Mixture Model EM example")
-    val ctx  = new SparkContext(conf)
-    
-    val data = ctx.textFile(inputFile).map { line =>
-      Vectors.dense(line.trim.split(' ').map(_.toDouble))
-    }.cache()
-      
-    val clusters = new GaussianMixtureEM()
-      .setK(k)
-      .setConvergenceTol(convergenceTol)
-      .setMaxIterations(maxIterations)
-      .run(data)
-    
-    for (i <- 0 until clusters.k) {
-      println("weight=%f\nmu=%s\nsigma=\n%s\n" format 
-        (clusters.weights(i), clusters.gaussians(i).mu, 
clusters.gaussians(i).sigma))
-    }
-    
-    println("Cluster labels (first <= 100):")
-    val clusterLabels = clusters.predict(data)
-    clusterLabels.take(100).foreach { x =>
-      print(" " + x)
-    }
-    println()
-  }
-}

http://git-wip-us.apache.org/repos/asf/spark/blob/98697734/mllib/src/main/scala/org/apache/spark/mllib/clustering/GaussianMixture.scala
----------------------------------------------------------------------
diff --git 
a/mllib/src/main/scala/org/apache/spark/mllib/clustering/GaussianMixture.scala 
b/mllib/src/main/scala/org/apache/spark/mllib/clustering/GaussianMixture.scala
new file mode 100644
index 0000000..5c626fd
--- /dev/null
+++ 
b/mllib/src/main/scala/org/apache/spark/mllib/clustering/GaussianMixture.scala
@@ -0,0 +1,251 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *    http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.spark.mllib.clustering
+
+import scala.collection.mutable.IndexedSeq
+
+import breeze.linalg.{DenseVector => BreezeVector, DenseMatrix => 
BreezeMatrix, diag, Transpose}
+
+import org.apache.spark.mllib.linalg.{Matrices, Vector, Vectors, DenseVector, 
DenseMatrix, BLAS}
+import org.apache.spark.mllib.stat.distribution.MultivariateGaussian
+import org.apache.spark.mllib.util.MLUtils
+import org.apache.spark.rdd.RDD
+import org.apache.spark.util.Utils
+
+/**
+ * This class performs expectation maximization for multivariate Gaussian
+ * Mixture Models (GMMs).  A GMM represents a composite distribution of
+ * independent Gaussian distributions with associated "mixing" weights
+ * specifying each's contribution to the composite.
+ *
+ * Given a set of sample points, this class will maximize the log-likelihood 
+ * for a mixture of k Gaussians, iterating until the log-likelihood changes by 
+ * less than convergenceTol, or until it has reached the max number of 
iterations.
+ * While this process is generally guaranteed to converge, it is not guaranteed
+ * to find a global optimum.  
+ * 
+ * @param k The number of independent Gaussians in the mixture model
+ * @param convergenceTol The maximum change in log-likelihood at which 
convergence
+ * is considered to have occurred.
+ * @param maxIterations The maximum number of iterations to perform
+ */
+class GaussianMixture private (
+    private var k: Int, 
+    private var convergenceTol: Double, 
+    private var maxIterations: Int,
+    private var seed: Long) extends Serializable {
+  
+  /** A default instance, 2 Gaussians, 100 iterations, 0.01 log-likelihood 
threshold */
+  def this() = this(2, 0.01, 100, Utils.random.nextLong())
+  
+  // number of samples per cluster to use when initializing Gaussians
+  private val nSamples = 5
+  
+  // an initializing GMM can be provided rather than using the 
+  // default random starting point
+  private var initialModel: Option[GaussianMixtureModel] = None
+  
+  /** Set the initial GMM starting point, bypassing the random initialization.
+   *  You must call setK() prior to calling this method, and the condition
+   *  (model.k == this.k) must be met; failure will result in an 
IllegalArgumentException
+   */
+  def setInitialModel(model: GaussianMixtureModel): this.type = {
+    if (model.k == k) {
+      initialModel = Some(model)
+    } else {
+      throw new IllegalArgumentException("mismatched cluster count (model.k != 
k)")
+    }
+    this
+  }
+  
+  /** Return the user supplied initial GMM, if supplied */
+  def getInitialModel: Option[GaussianMixtureModel] = initialModel
+  
+  /** Set the number of Gaussians in the mixture model.  Default: 2 */
+  def setK(k: Int): this.type = {
+    this.k = k
+    this
+  }
+  
+  /** Return the number of Gaussians in the mixture model */
+  def getK: Int = k
+  
+  /** Set the maximum number of iterations to run. Default: 100 */
+  def setMaxIterations(maxIterations: Int): this.type = {
+    this.maxIterations = maxIterations
+    this
+  }
+  
+  /** Return the maximum number of iterations to run */
+  def getMaxIterations: Int = maxIterations
+  
+  /**
+   * Set the largest change in log-likelihood at which convergence is 
+   * considered to have occurred.
+   */
+  def setConvergenceTol(convergenceTol: Double): this.type = {
+    this.convergenceTol = convergenceTol
+    this
+  }
+  
+  /**
+   * Return the largest change in log-likelihood at which convergence is
+   * considered to have occurred.
+   */
+  def getConvergenceTol: Double = convergenceTol
+
+  /** Set the random seed */
+  def setSeed(seed: Long): this.type = {
+    this.seed = seed
+    this
+  }
+
+  /** Return the random seed */
+  def getSeed: Long = seed
+
+  /** Perform expectation maximization */
+  def run(data: RDD[Vector]): GaussianMixtureModel = {
+    val sc = data.sparkContext
+    
+    // we will operate on the data as breeze data
+    val breezeData = data.map(u => u.toBreeze.toDenseVector).cache()
+    
+    // Get length of the input vectors
+    val d = breezeData.first().length
+    
+    // Determine initial weights and corresponding Gaussians.
+    // If the user supplied an initial GMM, we use those values, otherwise
+    // we start with uniform weights, a random mean from the data, and
+    // diagonal covariance matrices using component variances
+    // derived from the samples    
+    val (weights, gaussians) = initialModel match {
+      case Some(gmm) => (gmm.weights, gmm.gaussians)
+      
+      case None => {
+        val samples = breezeData.takeSample(withReplacement = true, k * 
nSamples, seed)
+        (Array.fill(k)(1.0 / k), Array.tabulate(k) { i => 
+          val slice = samples.view(i * nSamples, (i + 1) * nSamples)
+          new MultivariateGaussian(vectorMean(slice), initCovariance(slice)) 
+        })  
+      }
+    }
+    
+    var llh = Double.MinValue // current log-likelihood 
+    var llhp = 0.0            // previous log-likelihood
+    
+    var iter = 0
+    while(iter < maxIterations && Math.abs(llh-llhp) > convergenceTol) {
+      // create and broadcast curried cluster contribution function
+      val compute = sc.broadcast(ExpectationSum.add(weights, gaussians)_)
+      
+      // aggregate the cluster contribution for all sample points
+      val sums = breezeData.aggregate(ExpectationSum.zero(k, 
d))(compute.value, _ += _)
+      
+      // Create new distributions based on the partial assignments
+      // (often referred to as the "M" step in literature)
+      val sumWeights = sums.weights.sum
+      var i = 0
+      while (i < k) {
+        val mu = sums.means(i) / sums.weights(i)
+        BLAS.syr(-sums.weights(i), 
Vectors.fromBreeze(mu).asInstanceOf[DenseVector],
+          Matrices.fromBreeze(sums.sigmas(i)).asInstanceOf[DenseMatrix])
+        weights(i) = sums.weights(i) / sumWeights
+        gaussians(i) = new MultivariateGaussian(mu, sums.sigmas(i) / 
sums.weights(i))
+        i = i + 1
+      }
+   
+      llhp = llh // current becomes previous
+      llh = sums.logLikelihood // this is the freshly computed log-likelihood
+      iter += 1
+    } 
+    
+    new GaussianMixtureModel(weights, gaussians)
+  }
+    
+  /** Average of dense breeze vectors */
+  private def vectorMean(x: IndexedSeq[BreezeVector[Double]]): 
BreezeVector[Double] = {
+    val v = BreezeVector.zeros[Double](x(0).length)
+    x.foreach(xi => v += xi)
+    v / x.length.toDouble 
+  }
+  
+  /**
+   * Construct matrix where diagonal entries are element-wise
+   * variance of input vectors (computes biased variance)
+   */
+  private def initCovariance(x: IndexedSeq[BreezeVector[Double]]): 
BreezeMatrix[Double] = {
+    val mu = vectorMean(x)
+    val ss = BreezeVector.zeros[Double](x(0).length)
+    x.map(xi => (xi - mu) :^ 2.0).foreach(u => ss += u)
+    diag(ss / x.length.toDouble)
+  }
+}
+
+// companion class to provide zero constructor for ExpectationSum
+private object ExpectationSum {
+  def zero(k: Int, d: Int): ExpectationSum = {
+    new ExpectationSum(0.0, Array.fill(k)(0.0), 
+      Array.fill(k)(BreezeVector.zeros(d)), 
Array.fill(k)(BreezeMatrix.zeros(d,d)))
+  }
+  
+  // compute cluster contributions for each input point
+  // (U, T) => U for aggregation
+  def add(
+      weights: Array[Double], 
+      dists: Array[MultivariateGaussian])
+      (sums: ExpectationSum, x: BreezeVector[Double]): ExpectationSum = {
+    val p = weights.zip(dists).map {
+      case (weight, dist) => MLUtils.EPSILON + weight * dist.pdf(x)
+    }
+    val pSum = p.sum
+    sums.logLikelihood += math.log(pSum)
+    val xxt = x * new Transpose(x)
+    var i = 0
+    while (i < sums.k) {
+      p(i) /= pSum
+      sums.weights(i) += p(i)
+      sums.means(i) += x * p(i)
+      BLAS.syr(p(i), Vectors.fromBreeze(x).asInstanceOf[DenseVector],
+        Matrices.fromBreeze(sums.sigmas(i)).asInstanceOf[DenseMatrix])
+      i = i + 1
+    }
+    sums
+  }  
+}
+
+// Aggregation class for partial expectation results
+private class ExpectationSum(
+    var logLikelihood: Double,
+    val weights: Array[Double],
+    val means: Array[BreezeVector[Double]],
+    val sigmas: Array[BreezeMatrix[Double]]) extends Serializable {
+  
+  val k = weights.length
+  
+  def +=(x: ExpectationSum): ExpectationSum = {
+    var i = 0
+    while (i < k) {
+      weights(i) += x.weights(i)
+      means(i) += x.means(i)
+      sigmas(i) += x.sigmas(i)
+      i = i + 1
+    }
+    logLikelihood += x.logLikelihood
+    this
+  }  
+}

http://git-wip-us.apache.org/repos/asf/spark/blob/98697734/mllib/src/main/scala/org/apache/spark/mllib/clustering/GaussianMixtureEM.scala
----------------------------------------------------------------------
diff --git 
a/mllib/src/main/scala/org/apache/spark/mllib/clustering/GaussianMixtureEM.scala
 
b/mllib/src/main/scala/org/apache/spark/mllib/clustering/GaussianMixtureEM.scala
deleted file mode 100644
index 899fe5e..0000000
--- 
a/mllib/src/main/scala/org/apache/spark/mllib/clustering/GaussianMixtureEM.scala
+++ /dev/null
@@ -1,251 +0,0 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements.  See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License.  You may obtain a copy of the License at
- *
- *    http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-package org.apache.spark.mllib.clustering
-
-import scala.collection.mutable.IndexedSeq
-
-import breeze.linalg.{DenseVector => BreezeVector, DenseMatrix => 
BreezeMatrix, diag, Transpose}
-
-import org.apache.spark.mllib.linalg.{Matrices, Vector, Vectors, DenseVector, 
DenseMatrix, BLAS}
-import org.apache.spark.mllib.stat.distribution.MultivariateGaussian
-import org.apache.spark.mllib.util.MLUtils
-import org.apache.spark.rdd.RDD
-import org.apache.spark.util.Utils
-
-/**
- * This class performs expectation maximization for multivariate Gaussian
- * Mixture Models (GMMs).  A GMM represents a composite distribution of
- * independent Gaussian distributions with associated "mixing" weights
- * specifying each's contribution to the composite.
- *
- * Given a set of sample points, this class will maximize the log-likelihood 
- * for a mixture of k Gaussians, iterating until the log-likelihood changes by 
- * less than convergenceTol, or until it has reached the max number of 
iterations.
- * While this process is generally guaranteed to converge, it is not guaranteed
- * to find a global optimum.  
- * 
- * @param k The number of independent Gaussians in the mixture model
- * @param convergenceTol The maximum change in log-likelihood at which 
convergence
- * is considered to have occurred.
- * @param maxIterations The maximum number of iterations to perform
- */
-class GaussianMixtureEM private (
-    private var k: Int, 
-    private var convergenceTol: Double, 
-    private var maxIterations: Int,
-    private var seed: Long) extends Serializable {
-  
-  /** A default instance, 2 Gaussians, 100 iterations, 0.01 log-likelihood 
threshold */
-  def this() = this(2, 0.01, 100, Utils.random.nextLong())
-  
-  // number of samples per cluster to use when initializing Gaussians
-  private val nSamples = 5
-  
-  // an initializing GMM can be provided rather than using the 
-  // default random starting point
-  private var initialModel: Option[GaussianMixtureModel] = None
-  
-  /** Set the initial GMM starting point, bypassing the random initialization.
-   *  You must call setK() prior to calling this method, and the condition
-   *  (model.k == this.k) must be met; failure will result in an 
IllegalArgumentException
-   */
-  def setInitialModel(model: GaussianMixtureModel): this.type = {
-    if (model.k == k) {
-      initialModel = Some(model)
-    } else {
-      throw new IllegalArgumentException("mismatched cluster count (model.k != 
k)")
-    }
-    this
-  }
-  
-  /** Return the user supplied initial GMM, if supplied */
-  def getInitialModel: Option[GaussianMixtureModel] = initialModel
-  
-  /** Set the number of Gaussians in the mixture model.  Default: 2 */
-  def setK(k: Int): this.type = {
-    this.k = k
-    this
-  }
-  
-  /** Return the number of Gaussians in the mixture model */
-  def getK: Int = k
-  
-  /** Set the maximum number of iterations to run. Default: 100 */
-  def setMaxIterations(maxIterations: Int): this.type = {
-    this.maxIterations = maxIterations
-    this
-  }
-  
-  /** Return the maximum number of iterations to run */
-  def getMaxIterations: Int = maxIterations
-  
-  /**
-   * Set the largest change in log-likelihood at which convergence is 
-   * considered to have occurred.
-   */
-  def setConvergenceTol(convergenceTol: Double): this.type = {
-    this.convergenceTol = convergenceTol
-    this
-  }
-  
-  /**
-   * Return the largest change in log-likelihood at which convergence is
-   * considered to have occurred.
-   */
-  def getConvergenceTol: Double = convergenceTol
-
-  /** Set the random seed */
-  def setSeed(seed: Long): this.type = {
-    this.seed = seed
-    this
-  }
-
-  /** Return the random seed */
-  def getSeed: Long = seed
-
-  /** Perform expectation maximization */
-  def run(data: RDD[Vector]): GaussianMixtureModel = {
-    val sc = data.sparkContext
-    
-    // we will operate on the data as breeze data
-    val breezeData = data.map(u => u.toBreeze.toDenseVector).cache()
-    
-    // Get length of the input vectors
-    val d = breezeData.first().length
-    
-    // Determine initial weights and corresponding Gaussians.
-    // If the user supplied an initial GMM, we use those values, otherwise
-    // we start with uniform weights, a random mean from the data, and
-    // diagonal covariance matrices using component variances
-    // derived from the samples    
-    val (weights, gaussians) = initialModel match {
-      case Some(gmm) => (gmm.weights, gmm.gaussians)
-      
-      case None => {
-        val samples = breezeData.takeSample(withReplacement = true, k * 
nSamples, seed)
-        (Array.fill(k)(1.0 / k), Array.tabulate(k) { i => 
-          val slice = samples.view(i * nSamples, (i + 1) * nSamples)
-          new MultivariateGaussian(vectorMean(slice), initCovariance(slice)) 
-        })  
-      }
-    }
-    
-    var llh = Double.MinValue // current log-likelihood 
-    var llhp = 0.0            // previous log-likelihood
-    
-    var iter = 0
-    while(iter < maxIterations && Math.abs(llh-llhp) > convergenceTol) {
-      // create and broadcast curried cluster contribution function
-      val compute = sc.broadcast(ExpectationSum.add(weights, gaussians)_)
-      
-      // aggregate the cluster contribution for all sample points
-      val sums = breezeData.aggregate(ExpectationSum.zero(k, 
d))(compute.value, _ += _)
-      
-      // Create new distributions based on the partial assignments
-      // (often referred to as the "M" step in literature)
-      val sumWeights = sums.weights.sum
-      var i = 0
-      while (i < k) {
-        val mu = sums.means(i) / sums.weights(i)
-        BLAS.syr(-sums.weights(i), 
Vectors.fromBreeze(mu).asInstanceOf[DenseVector],
-          Matrices.fromBreeze(sums.sigmas(i)).asInstanceOf[DenseMatrix])
-        weights(i) = sums.weights(i) / sumWeights
-        gaussians(i) = new MultivariateGaussian(mu, sums.sigmas(i) / 
sums.weights(i))
-        i = i + 1
-      }
-   
-      llhp = llh // current becomes previous
-      llh = sums.logLikelihood // this is the freshly computed log-likelihood
-      iter += 1
-    } 
-    
-    new GaussianMixtureModel(weights, gaussians)
-  }
-    
-  /** Average of dense breeze vectors */
-  private def vectorMean(x: IndexedSeq[BreezeVector[Double]]): 
BreezeVector[Double] = {
-    val v = BreezeVector.zeros[Double](x(0).length)
-    x.foreach(xi => v += xi)
-    v / x.length.toDouble 
-  }
-  
-  /**
-   * Construct matrix where diagonal entries are element-wise
-   * variance of input vectors (computes biased variance)
-   */
-  private def initCovariance(x: IndexedSeq[BreezeVector[Double]]): 
BreezeMatrix[Double] = {
-    val mu = vectorMean(x)
-    val ss = BreezeVector.zeros[Double](x(0).length)
-    x.map(xi => (xi - mu) :^ 2.0).foreach(u => ss += u)
-    diag(ss / x.length.toDouble)
-  }
-}
-
-// companion class to provide zero constructor for ExpectationSum
-private object ExpectationSum {
-  def zero(k: Int, d: Int): ExpectationSum = {
-    new ExpectationSum(0.0, Array.fill(k)(0.0), 
-      Array.fill(k)(BreezeVector.zeros(d)), 
Array.fill(k)(BreezeMatrix.zeros(d,d)))
-  }
-  
-  // compute cluster contributions for each input point
-  // (U, T) => U for aggregation
-  def add(
-      weights: Array[Double], 
-      dists: Array[MultivariateGaussian])
-      (sums: ExpectationSum, x: BreezeVector[Double]): ExpectationSum = {
-    val p = weights.zip(dists).map {
-      case (weight, dist) => MLUtils.EPSILON + weight * dist.pdf(x)
-    }
-    val pSum = p.sum
-    sums.logLikelihood += math.log(pSum)
-    val xxt = x * new Transpose(x)
-    var i = 0
-    while (i < sums.k) {
-      p(i) /= pSum
-      sums.weights(i) += p(i)
-      sums.means(i) += x * p(i)
-      BLAS.syr(p(i), Vectors.fromBreeze(x).asInstanceOf[DenseVector],
-        Matrices.fromBreeze(sums.sigmas(i)).asInstanceOf[DenseMatrix])
-      i = i + 1
-    }
-    sums
-  }  
-}
-
-// Aggregation class for partial expectation results
-private class ExpectationSum(
-    var logLikelihood: Double,
-    val weights: Array[Double],
-    val means: Array[BreezeVector[Double]],
-    val sigmas: Array[BreezeMatrix[Double]]) extends Serializable {
-  
-  val k = weights.length
-  
-  def +=(x: ExpectationSum): ExpectationSum = {
-    var i = 0
-    while (i < k) {
-      weights(i) += x.weights(i)
-      means(i) += x.means(i)
-      sigmas(i) += x.sigmas(i)
-      i = i + 1
-    }
-    logLikelihood += x.logLikelihood
-    this
-  }  
-}

http://git-wip-us.apache.org/repos/asf/spark/blob/98697734/mllib/src/test/scala/org/apache/spark/mllib/clustering/GMMExpectationMaximizationSuite.scala
----------------------------------------------------------------------
diff --git 
a/mllib/src/test/scala/org/apache/spark/mllib/clustering/GMMExpectationMaximizationSuite.scala
 
b/mllib/src/test/scala/org/apache/spark/mllib/clustering/GMMExpectationMaximizationSuite.scala
deleted file mode 100644
index 198997b..0000000
--- 
a/mllib/src/test/scala/org/apache/spark/mllib/clustering/GMMExpectationMaximizationSuite.scala
+++ /dev/null
@@ -1,83 +0,0 @@
-/*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements.  See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License.  You may obtain a copy of the License at
- *
- *    http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-
-package org.apache.spark.mllib.clustering
-
-import org.scalatest.FunSuite
-
-import org.apache.spark.mllib.linalg.{Vectors, Matrices}
-import org.apache.spark.mllib.stat.distribution.MultivariateGaussian
-import org.apache.spark.mllib.util.MLlibTestSparkContext
-import org.apache.spark.mllib.util.TestingUtils._
-
-class GMMExpectationMaximizationSuite extends FunSuite with 
MLlibTestSparkContext {
-  test("single cluster") {
-    val data = sc.parallelize(Array(
-      Vectors.dense(6.0, 9.0),
-      Vectors.dense(5.0, 10.0),
-      Vectors.dense(4.0, 11.0)
-    ))
-    
-    // expectations
-    val Ew = 1.0
-    val Emu = Vectors.dense(5.0, 10.0)
-    val Esigma = Matrices.dense(2, 2, Array(2.0 / 3.0, -2.0 / 3.0, -2.0 / 3.0, 
2.0 / 3.0))
-
-    val seeds = Array(314589, 29032897, 50181, 494821, 4660)
-    seeds.foreach { seed =>
-      val gmm = new GaussianMixtureEM().setK(1).setSeed(seed).run(data)
-      assert(gmm.weights(0) ~== Ew absTol 1E-5)
-      assert(gmm.gaussians(0).mu ~== Emu absTol 1E-5)
-      assert(gmm.gaussians(0).sigma ~== Esigma absTol 1E-5)
-    }
-  }
-  
-  test("two clusters") {
-    val data = sc.parallelize(Array(
-      Vectors.dense(-5.1971), Vectors.dense(-2.5359), Vectors.dense(-3.8220),
-      Vectors.dense(-5.2211), Vectors.dense(-5.0602), Vectors.dense( 4.7118),
-      Vectors.dense( 6.8989), Vectors.dense( 3.4592), Vectors.dense( 4.6322),
-      Vectors.dense( 5.7048), Vectors.dense( 4.6567), Vectors.dense( 5.5026),
-      Vectors.dense( 4.5605), Vectors.dense( 5.2043), Vectors.dense( 6.2734)
-    ))
-  
-    // we set an initial gaussian to induce expected results
-    val initialGmm = new GaussianMixtureModel(
-      Array(0.5, 0.5),
-      Array(
-        new MultivariateGaussian(Vectors.dense(-1.0), Matrices.dense(1, 1, 
Array(1.0))),
-        new MultivariateGaussian(Vectors.dense(1.0), Matrices.dense(1, 1, 
Array(1.0)))
-      )
-    )
-    
-    val Ew = Array(1.0 / 3.0, 2.0 / 3.0)
-    val Emu = Array(Vectors.dense(-4.3673), Vectors.dense(5.1604))
-    val Esigma = Array(Matrices.dense(1, 1, Array(1.1098)), Matrices.dense(1, 
1, Array(0.86644)))
-    
-    val gmm = new GaussianMixtureEM()
-      .setK(2)
-      .setInitialModel(initialGmm)
-      .run(data)
-      
-    assert(gmm.weights(0) ~== Ew(0) absTol 1E-3)
-    assert(gmm.weights(1) ~== Ew(1) absTol 1E-3)
-    assert(gmm.gaussians(0).mu ~== Emu(0) absTol 1E-3)
-    assert(gmm.gaussians(1).mu ~== Emu(1) absTol 1E-3)
-    assert(gmm.gaussians(0).sigma ~== Esigma(0) absTol 1E-3)
-    assert(gmm.gaussians(1).sigma ~== Esigma(1) absTol 1E-3)
-  }
-}

http://git-wip-us.apache.org/repos/asf/spark/blob/98697734/mllib/src/test/scala/org/apache/spark/mllib/clustering/GaussianMixtureSuite.scala
----------------------------------------------------------------------
diff --git 
a/mllib/src/test/scala/org/apache/spark/mllib/clustering/GaussianMixtureSuite.scala
 
b/mllib/src/test/scala/org/apache/spark/mllib/clustering/GaussianMixtureSuite.scala
new file mode 100644
index 0000000..c2cd56e
--- /dev/null
+++ 
b/mllib/src/test/scala/org/apache/spark/mllib/clustering/GaussianMixtureSuite.scala
@@ -0,0 +1,83 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *    http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.spark.mllib.clustering
+
+import org.scalatest.FunSuite
+
+import org.apache.spark.mllib.linalg.{Vectors, Matrices}
+import org.apache.spark.mllib.stat.distribution.MultivariateGaussian
+import org.apache.spark.mllib.util.MLlibTestSparkContext
+import org.apache.spark.mllib.util.TestingUtils._
+
+class GaussianMixtureSuite extends FunSuite with MLlibTestSparkContext {
+  test("single cluster") {
+    val data = sc.parallelize(Array(
+      Vectors.dense(6.0, 9.0),
+      Vectors.dense(5.0, 10.0),
+      Vectors.dense(4.0, 11.0)
+    ))
+    
+    // expectations
+    val Ew = 1.0
+    val Emu = Vectors.dense(5.0, 10.0)
+    val Esigma = Matrices.dense(2, 2, Array(2.0 / 3.0, -2.0 / 3.0, -2.0 / 3.0, 
2.0 / 3.0))
+
+    val seeds = Array(314589, 29032897, 50181, 494821, 4660)
+    seeds.foreach { seed =>
+      val gmm = new GaussianMixture().setK(1).setSeed(seed).run(data)
+      assert(gmm.weights(0) ~== Ew absTol 1E-5)
+      assert(gmm.gaussians(0).mu ~== Emu absTol 1E-5)
+      assert(gmm.gaussians(0).sigma ~== Esigma absTol 1E-5)
+    }
+  }
+  
+  test("two clusters") {
+    val data = sc.parallelize(Array(
+      Vectors.dense(-5.1971), Vectors.dense(-2.5359), Vectors.dense(-3.8220),
+      Vectors.dense(-5.2211), Vectors.dense(-5.0602), Vectors.dense( 4.7118),
+      Vectors.dense( 6.8989), Vectors.dense( 3.4592), Vectors.dense( 4.6322),
+      Vectors.dense( 5.7048), Vectors.dense( 4.6567), Vectors.dense( 5.5026),
+      Vectors.dense( 4.5605), Vectors.dense( 5.2043), Vectors.dense( 6.2734)
+    ))
+  
+    // we set an initial gaussian to induce expected results
+    val initialGmm = new GaussianMixtureModel(
+      Array(0.5, 0.5),
+      Array(
+        new MultivariateGaussian(Vectors.dense(-1.0), Matrices.dense(1, 1, 
Array(1.0))),
+        new MultivariateGaussian(Vectors.dense(1.0), Matrices.dense(1, 1, 
Array(1.0)))
+      )
+    )
+    
+    val Ew = Array(1.0 / 3.0, 2.0 / 3.0)
+    val Emu = Array(Vectors.dense(-4.3673), Vectors.dense(5.1604))
+    val Esigma = Array(Matrices.dense(1, 1, Array(1.1098)), Matrices.dense(1, 
1, Array(0.86644)))
+    
+    val gmm = new GaussianMixture()
+      .setK(2)
+      .setInitialModel(initialGmm)
+      .run(data)
+      
+    assert(gmm.weights(0) ~== Ew(0) absTol 1E-3)
+    assert(gmm.weights(1) ~== Ew(1) absTol 1E-3)
+    assert(gmm.gaussians(0).mu ~== Emu(0) absTol 1E-3)
+    assert(gmm.gaussians(1).mu ~== Emu(1) absTol 1E-3)
+    assert(gmm.gaussians(0).sigma ~== Esigma(0) absTol 1E-3)
+    assert(gmm.gaussians(1).sigma ~== Esigma(1) absTol 1E-3)
+  }
+}


---------------------------------------------------------------------
To unsubscribe, e-mail: [email protected]
For additional commands, e-mail: [email protected]

Reply via email to