|
| 1 | +/* |
| 2 | + * Licensed to the Apache Software Foundation (ASF) under one or more |
| 3 | + * contributor license agreements. See the NOTICE file distributed with |
| 4 | + * this work for additional information regarding copyright ownership. |
| 5 | + * The ASF licenses this file to You under the Apache License, Version 2.0 |
| 6 | + * (the "License"); you may not use this file except in compliance with |
| 7 | + * the License. You may obtain a copy of the License at |
| 8 | + * |
| 9 | + * http://www.apache.org/licenses/LICENSE-2.0 |
| 10 | + * |
| 11 | + * Unless required by applicable law or agreed to in writing, software |
| 12 | + * distributed under the License is distributed on an "AS IS" BASIS, |
| 13 | + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 14 | + * See the License for the specific language governing permissions and |
| 15 | + * limitations under the License. |
| 16 | + */ |
| 17 | + |
| 18 | +package org.apache.spark.mllib.clustering |
| 19 | + |
| 20 | +import scala.collection.mutable.IndexedSeq |
| 21 | + |
| 22 | +import breeze.linalg.{DenseVector => BreezeVector, DenseMatrix => BreezeMatrix, diag, Transpose} |
| 23 | +import org.apache.spark.rdd.RDD |
| 24 | +import org.apache.spark.mllib.linalg.{Matrices, Vector, Vectors} |
| 25 | +import org.apache.spark.mllib.stat.impl.MultivariateGaussian |
| 26 | +import org.apache.spark.mllib.util.MLUtils |
| 27 | + |
| 28 | +/** |
| 29 | + * This class performs expectation maximization for multivariate Gaussian |
| 30 | + * Mixture Models (GMMs). A GMM represents a composite distribution of |
| 31 | + * independent Gaussian distributions with associated "mixing" weights |
| 32 | + * specifying each's contribution to the composite. |
| 33 | + * |
| 34 | + * Given a set of sample points, this class will maximize the log-likelihood |
| 35 | + * for a mixture of k Gaussians, iterating until the log-likelihood changes by |
| 36 | + * less than convergenceTol, or until it has reached the max number of iterations. |
| 37 | + * While this process is generally guaranteed to converge, it is not guaranteed |
| 38 | + * to find a global optimum. |
| 39 | + * |
| 40 | + * @param k The number of independent Gaussians in the mixture model |
| 41 | + * @param convergenceTol The maximum change in log-likelihood at which convergence |
| 42 | + * is considered to have occurred. |
| 43 | + * @param maxIterations The maximum number of iterations to perform |
| 44 | + */ |
| 45 | +class GaussianMixtureEM private ( |
| 46 | + private var k: Int, |
| 47 | + private var convergenceTol: Double, |
| 48 | + private var maxIterations: Int) extends Serializable { |
| 49 | + |
| 50 | + /** A default instance, 2 Gaussians, 100 iterations, 0.01 log-likelihood threshold */ |
| 51 | + def this() = this(2, 0.01, 100) |
| 52 | + |
| 53 | + // number of samples per cluster to use when initializing Gaussians |
| 54 | + private val nSamples = 5 |
| 55 | + |
| 56 | + // an initializing GMM can be provided rather than using the |
| 57 | + // default random starting point |
| 58 | + private var initialModel: Option[GaussianMixtureModel] = None |
| 59 | + |
| 60 | + /** Set the initial GMM starting point, bypassing the random initialization. |
| 61 | + * You must call setK() prior to calling this method, and the condition |
| 62 | + * (model.k == this.k) must be met; failure will result in an IllegalArgumentException |
| 63 | + */ |
| 64 | + def setInitialModel(model: GaussianMixtureModel): this.type = { |
| 65 | + if (model.k == k) { |
| 66 | + initialModel = Some(model) |
| 67 | + } else { |
| 68 | + throw new IllegalArgumentException("mismatched cluster count (model.k != k)") |
| 69 | + } |
| 70 | + this |
| 71 | + } |
| 72 | + |
| 73 | + /** Return the user supplied initial GMM, if supplied */ |
| 74 | + def getInitialModel: Option[GaussianMixtureModel] = initialModel |
| 75 | + |
| 76 | + /** Set the number of Gaussians in the mixture model. Default: 2 */ |
| 77 | + def setK(k: Int): this.type = { |
| 78 | + this.k = k |
| 79 | + this |
| 80 | + } |
| 81 | + |
| 82 | + /** Return the number of Gaussians in the mixture model */ |
| 83 | + def getK: Int = k |
| 84 | + |
| 85 | + /** Set the maximum number of iterations to run. Default: 100 */ |
| 86 | + def setMaxIterations(maxIterations: Int): this.type = { |
| 87 | + this.maxIterations = maxIterations |
| 88 | + this |
| 89 | + } |
| 90 | + |
| 91 | + /** Return the maximum number of iterations to run */ |
| 92 | + def getMaxIterations: Int = maxIterations |
| 93 | + |
| 94 | + /** |
| 95 | + * Set the largest change in log-likelihood at which convergence is |
| 96 | + * considered to have occurred. |
| 97 | + */ |
| 98 | + def setConvergenceTol(convergenceTol: Double): this.type = { |
| 99 | + this.convergenceTol = convergenceTol |
| 100 | + this |
| 101 | + } |
| 102 | + |
| 103 | + /** Return the largest change in log-likelihood at which convergence is |
| 104 | + * considered to have occurred. |
| 105 | + */ |
| 106 | + def getConvergenceTol: Double = convergenceTol |
| 107 | + |
| 108 | + /** Perform expectation maximization */ |
| 109 | + def run(data: RDD[Vector]): GaussianMixtureModel = { |
| 110 | + val sc = data.sparkContext |
| 111 | + |
| 112 | + // we will operate on the data as breeze data |
| 113 | + val breezeData = data.map(u => u.toBreeze.toDenseVector).cache() |
| 114 | + |
| 115 | + // Get length of the input vectors |
| 116 | + val d = breezeData.first.length |
| 117 | + |
| 118 | + // Determine initial weights and corresponding Gaussians. |
| 119 | + // If the user supplied an initial GMM, we use those values, otherwise |
| 120 | + // we start with uniform weights, a random mean from the data, and |
| 121 | + // diagonal covariance matrices using component variances |
| 122 | + // derived from the samples |
| 123 | + val (weights, gaussians) = initialModel match { |
| 124 | + case Some(gmm) => (gmm.weight, gmm.mu.zip(gmm.sigma).map { case(mu, sigma) => |
| 125 | + new MultivariateGaussian(mu.toBreeze.toDenseVector, sigma.toBreeze.toDenseMatrix) |
| 126 | + }) |
| 127 | + |
| 128 | + case None => { |
| 129 | + val samples = breezeData.takeSample(true, k * nSamples, scala.util.Random.nextInt) |
| 130 | + (Array.fill(k)(1.0 / k), Array.tabulate(k) { i => |
| 131 | + val slice = samples.view(i * nSamples, (i + 1) * nSamples) |
| 132 | + new MultivariateGaussian(vectorMean(slice), initCovariance(slice)) |
| 133 | + }) |
| 134 | + } |
| 135 | + } |
| 136 | + |
| 137 | + var llh = Double.MinValue // current log-likelihood |
| 138 | + var llhp = 0.0 // previous log-likelihood |
| 139 | + |
| 140 | + var iter = 0 |
| 141 | + while(iter < maxIterations && Math.abs(llh-llhp) > convergenceTol) { |
| 142 | + // create and broadcast curried cluster contribution function |
| 143 | + val compute = sc.broadcast(ExpectationSum.add(weights, gaussians)_) |
| 144 | + |
| 145 | + // aggregate the cluster contribution for all sample points |
| 146 | + val sums = breezeData.aggregate(ExpectationSum.zero(k, d))(compute.value, _ += _) |
| 147 | + |
| 148 | + // Create new distributions based on the partial assignments |
| 149 | + // (often referred to as the "M" step in literature) |
| 150 | + val sumWeights = sums.weights.sum |
| 151 | + var i = 0 |
| 152 | + while (i < k) { |
| 153 | + val mu = sums.means(i) / sums.weights(i) |
| 154 | + val sigma = sums.sigmas(i) / sums.weights(i) - mu * new Transpose(mu) // TODO: Use BLAS.dsyr |
| 155 | + weights(i) = sums.weights(i) / sumWeights |
| 156 | + gaussians(i) = new MultivariateGaussian(mu, sigma) |
| 157 | + i = i + 1 |
| 158 | + } |
| 159 | + |
| 160 | + llhp = llh // current becomes previous |
| 161 | + llh = sums.logLikelihood // this is the freshly computed log-likelihood |
| 162 | + iter += 1 |
| 163 | + } |
| 164 | + |
| 165 | + // Need to convert the breeze matrices to MLlib matrices |
| 166 | + val means = Array.tabulate(k) { i => Vectors.fromBreeze(gaussians(i).mu) } |
| 167 | + val sigmas = Array.tabulate(k) { i => Matrices.fromBreeze(gaussians(i).sigma) } |
| 168 | + new GaussianMixtureModel(weights, means, sigmas) |
| 169 | + } |
| 170 | + |
| 171 | + /** Average of dense breeze vectors */ |
| 172 | + private def vectorMean(x: IndexedSeq[BreezeVector[Double]]): BreezeVector[Double] = { |
| 173 | + val v = BreezeVector.zeros[Double](x(0).length) |
| 174 | + x.foreach(xi => v += xi) |
| 175 | + v / x.length.toDouble |
| 176 | + } |
| 177 | + |
| 178 | + /** |
| 179 | + * Construct matrix where diagonal entries are element-wise |
| 180 | + * variance of input vectors (computes biased variance) |
| 181 | + */ |
| 182 | + private def initCovariance(x: IndexedSeq[BreezeVector[Double]]): BreezeMatrix[Double] = { |
| 183 | + val mu = vectorMean(x) |
| 184 | + val ss = BreezeVector.zeros[Double](x(0).length) |
| 185 | + x.map(xi => (xi - mu) :^ 2.0).foreach(u => ss += u) |
| 186 | + diag(ss / x.length.toDouble) |
| 187 | + } |
| 188 | +} |
| 189 | + |
| 190 | +// companion class to provide zero constructor for ExpectationSum |
| 191 | +private object ExpectationSum { |
| 192 | + def zero(k: Int, d: Int): ExpectationSum = { |
| 193 | + new ExpectationSum(0.0, Array.fill(k)(0.0), |
| 194 | + Array.fill(k)(BreezeVector.zeros(d)), Array.fill(k)(BreezeMatrix.zeros(d,d))) |
| 195 | + } |
| 196 | + |
| 197 | + // compute cluster contributions for each input point |
| 198 | + // (U, T) => U for aggregation |
| 199 | + def add( |
| 200 | + weights: Array[Double], |
| 201 | + dists: Array[MultivariateGaussian]) |
| 202 | + (sums: ExpectationSum, x: BreezeVector[Double]): ExpectationSum = { |
| 203 | + val p = weights.zip(dists).map { |
| 204 | + case (weight, dist) => MLUtils.EPSILON + weight * dist.pdf(x) |
| 205 | + } |
| 206 | + val pSum = p.sum |
| 207 | + sums.logLikelihood += math.log(pSum) |
| 208 | + val xxt = x * new Transpose(x) |
| 209 | + var i = 0 |
| 210 | + while (i < sums.k) { |
| 211 | + p(i) /= pSum |
| 212 | + sums.weights(i) += p(i) |
| 213 | + sums.means(i) += x * p(i) |
| 214 | + sums.sigmas(i) += xxt * p(i) // TODO: use BLAS.dsyr |
| 215 | + i = i + 1 |
| 216 | + } |
| 217 | + sums |
| 218 | + } |
| 219 | +} |
| 220 | + |
| 221 | +// Aggregation class for partial expectation results |
| 222 | +private class ExpectationSum( |
| 223 | + var logLikelihood: Double, |
| 224 | + val weights: Array[Double], |
| 225 | + val means: Array[BreezeVector[Double]], |
| 226 | + val sigmas: Array[BreezeMatrix[Double]]) extends Serializable { |
| 227 | + |
| 228 | + val k = weights.length |
| 229 | + |
| 230 | + def +=(x: ExpectationSum): ExpectationSum = { |
| 231 | + var i = 0 |
| 232 | + while (i < k) { |
| 233 | + weights(i) += x.weights(i) |
| 234 | + means(i) += x.means(i) |
| 235 | + sigmas(i) += x.sigmas(i) |
| 236 | + i = i + 1 |
| 237 | + } |
| 238 | + logLikelihood += x.logLikelihood |
| 239 | + this |
| 240 | + } |
| 241 | +} |
0 commit comments