1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30 package org.openimaj.ml.gmm;
31
32 import java.util.Arrays;
33 import java.util.EnumSet;
34
35 import org.apache.commons.math.util.MathUtils;
36 import org.openimaj.math.matrix.MatrixUtils;
37 import org.openimaj.math.statistics.MeanAndCovariance;
38 import org.openimaj.math.statistics.distribution.AbstractMultivariateGaussian;
39 import org.openimaj.math.statistics.distribution.DiagonalMultivariateGaussian;
40 import org.openimaj.math.statistics.distribution.FullMultivariateGaussian;
41 import org.openimaj.math.statistics.distribution.MixtureOfGaussians;
42 import org.openimaj.math.statistics.distribution.MultivariateGaussian;
43 import org.openimaj.math.statistics.distribution.SphericalMultivariateGaussian;
44 import org.openimaj.ml.clustering.DoubleCentroidsResult;
45 import org.openimaj.ml.clustering.kmeans.DoubleKMeans;
46 import org.openimaj.util.array.ArrayUtils;
47 import org.openimaj.util.pair.IndependentPair;
48
49 import Jama.Matrix;
50 import gnu.trove.list.array.TDoubleArrayList;
51
52
53
54
55
56
57
58
59
60
61
62 public class GaussianMixtureModelEM {
63
64
65
66
67
68
69 public static enum CovarianceType {
70
71
72
73
74 Spherical {
75 @Override
76 protected void setCovariances(MultivariateGaussian[] gaussians, Matrix cv) {
77 double mean = 0;
78
79 for (int i = 0; i < cv.getRowDimension(); i++)
80 for (int j = 0; j < cv.getColumnDimension(); j++)
81 mean += cv.get(i, j);
82 mean /= (cv.getColumnDimension() * cv.getRowDimension());
83
84 for (final MultivariateGaussian mg : gaussians) {
85 ((SphericalMultivariateGaussian) mg).variance = mean;
86 }
87 }
88
89 @Override
90 protected MultivariateGaussian[] createGaussians(int ngauss, int ndims) {
91 final MultivariateGaussian[] arr = new MultivariateGaussian[ngauss];
92 for (int i = 0; i < ngauss; i++) {
93 arr[i] = new SphericalMultivariateGaussian(ndims);
94 }
95
96 return arr;
97 }
98
99 @Override
100 protected void mstep(EMGMM gmm, GaussianMixtureModelEM learner, Matrix X, Matrix responsibilities,
101 Matrix weightedXsum,
102 double[] norm)
103 {
104 final Matrix avgX2uw = responsibilities.transpose().times(X.arrayTimes(X));
105
106 for (int i = 0; i < gmm.gaussians.length; i++) {
107 final Matrix weightedXsumi = new Matrix(new double[][] { weightedXsum.getArray()[i] });
108 final Matrix avgX2uwi = new Matrix(new double[][] { avgX2uw.getArray()[i] });
109
110 final Matrix avgX2 = avgX2uwi.times(norm[i]);
111 final Matrix mu = ((AbstractMultivariateGaussian) gmm.gaussians[i]).mean;
112 final Matrix avgMeans2 = MatrixUtils.pow(mu, 2);
113 final Matrix avgXmeans = mu.arrayTimes(weightedXsumi).times(norm[i]);
114 final Matrix covar = MatrixUtils.plus(avgX2.minus(avgXmeans.times(2)).plus(avgMeans2),
115 learner.minCovar);
116
117 ((SphericalMultivariateGaussian) gmm.gaussians[i]).variance = MatrixUtils.sum(covar)
118 / X.getColumnDimension();
119 }
120 }
121 },
122
123
124
125 Diagonal {
126 @Override
127 protected void setCovariances(MultivariateGaussian[] gaussians, Matrix cv) {
128 for (final MultivariateGaussian mg : gaussians) {
129 ((DiagonalMultivariateGaussian) mg).variance = MatrixUtils.diagVector(cv);
130 }
131 }
132
133 @Override
134 protected MultivariateGaussian[] createGaussians(int ngauss, int ndims) {
135 final MultivariateGaussian[] arr = new MultivariateGaussian[ngauss];
136 for (int i = 0; i < ngauss; i++) {
137 arr[i] = new DiagonalMultivariateGaussian(ndims);
138 }
139
140 return arr;
141 }
142
143 @Override
144 protected void mstep(EMGMM gmm, GaussianMixtureModelEM learner, Matrix X, Matrix responsibilities,
145 Matrix weightedXsum,
146 double[] norm)
147 {
148 final Matrix avgX2uw = responsibilities.transpose().times(X.arrayTimes(X));
149
150 for (int i = 0; i < gmm.gaussians.length; i++) {
151 final Matrix weightedXsumi = new Matrix(new double[][] { weightedXsum.getArray()[i] });
152 final Matrix avgX2uwi = new Matrix(new double[][] { avgX2uw.getArray()[i] });
153
154 final Matrix avgX2 = avgX2uwi.times(norm[i]);
155 final Matrix mu = ((AbstractMultivariateGaussian) gmm.gaussians[i]).mean;
156 final Matrix avgMeans2 = MatrixUtils.pow(mu, 2);
157 final Matrix avgXmeans = mu.arrayTimes(weightedXsumi).times(norm[i]);
158
159 final Matrix covar = MatrixUtils.plus(avgX2.minus(avgXmeans.times(2)).plus(avgMeans2),
160 learner.minCovar);
161
162 ((DiagonalMultivariateGaussian) gmm.gaussians[i]).variance = covar.getArray()[0];
163 }
164 }
165 },
166
167
168
169 Full {
170 @Override
171 protected MultivariateGaussian[] createGaussians(int ngauss, int ndims) {
172 final MultivariateGaussian[] arr = new MultivariateGaussian[ngauss];
173 for (int i = 0; i < ngauss; i++) {
174 arr[i] = new FullMultivariateGaussian(ndims);
175 }
176
177 return arr;
178 }
179
180 @Override
181 protected void setCovariances(MultivariateGaussian[] gaussians, Matrix cv) {
182 for (final MultivariateGaussian mg : gaussians) {
183 ((FullMultivariateGaussian) mg).covar = cv.copy();
184 }
185 }
186
187 @Override
188 protected void mstep(EMGMM gmm, GaussianMixtureModelEM learner, Matrix X, Matrix responsibilities,
189 Matrix weightedXsum,
190 double[] norm)
191 {
192
193
194 final int nfeatures = X.getColumnDimension();
195 for (int c = 0; c < learner.nComponents; c++) {
196 final Matrix post = responsibilities.getMatrix(0, X.getRowDimension() - 1, c, c).transpose();
197
198 final double factor = 1.0 / (ArrayUtils.sumValues(post.getArray()) + 10 * MathUtils.EPSILON);
199
200 final Matrix pXt = X.transpose();
201 for (int i = 0; i < pXt.getRowDimension(); i++)
202 for (int j = 0; j < pXt.getColumnDimension(); j++)
203 pXt.set(i, j, pXt.get(i, j) * post.get(0, j));
204
205 final Matrix argcv = pXt.times(X).times(factor);
206 final Matrix mu = ((FullMultivariateGaussian) gmm.gaussians[c]).mean;
207
208 ((FullMultivariateGaussian) gmm.gaussians[c]).covar = argcv.minusEquals(mu.transpose().times(mu))
209 .plusEquals(Matrix.identity(nfeatures, nfeatures).times(learner.minCovar));
210 }
211 }
212 },
213
214
215
216
217 Tied {
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269 @Override
270 protected void setCovariances(MultivariateGaussian[] gaussians,
271 Matrix cv)
272 {
273 for (final MultivariateGaussian mg : gaussians) {
274 ((FullMultivariateGaussian) mg).covar = cv;
275 }
276 }
277
278 @Override
279 protected MultivariateGaussian[] createGaussians(int ngauss, int ndims) {
280 final MultivariateGaussian[] arr = new MultivariateGaussian[ngauss];
281 final Matrix covar = new Matrix(ndims, ndims);
282
283 for (int i = 0; i < ngauss; i++) {
284 arr[i] = new FullMultivariateGaussian(new Matrix(1, ndims), covar);
285 }
286
287 return arr;
288 }
289
290 @Override
291 protected void mstep(EMGMM gmm, GaussianMixtureModelEM learner, Matrix X, Matrix responsibilities,
292 Matrix weightedXsum, double[] norm)
293 {
294
295 final int nfeatures = X.getColumnDimension();
296
297 final Matrix avgX2 = X.transpose().times(X);
298 final double[][] mudata = new double[gmm.gaussians.length][];
299 for (int i = 0; i < mudata.length; i++)
300 mudata[i] = ((FullMultivariateGaussian) gmm.gaussians[i]).mean.getArray()[0];
301 final Matrix mu = new Matrix(mudata);
302
303 final Matrix avgMeans2 = mu.transpose().times(weightedXsum);
304 final Matrix covar = avgX2.minus(avgMeans2)
305 .plus(Matrix.identity(nfeatures, nfeatures).times(learner.minCovar))
306 .times(1.0 / X.getRowDimension());
307
308 for (int i = 0; i < learner.nComponents; i++)
309 ((FullMultivariateGaussian) gmm.gaussians[i]).covar = covar;
310 }
311 };
312
313 protected abstract MultivariateGaussian[] createGaussians(int ngauss, int ndims);
314
315 protected abstract void setCovariances(MultivariateGaussian[] gaussians, Matrix cv);
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334 protected abstract void mstep(EMGMM gmm, GaussianMixtureModelEM learner, Matrix X,
335 Matrix responsibilities, Matrix weightedXsum, double[] inverseWeights);
336 }
337
338
339
340
341
342
343
344 public static enum UpdateOptions {
345
346
347
348 Means,
349
350
351
352 Weights,
353
354
355
356 Covariances
357 }
358
359 protected static class EMGMM extends MixtureOfGaussians {
360 EMGMM(int nComponents) {
361 super(null, null);
362
363 this.weights = new double[nComponents];
364 Arrays.fill(this.weights, 1.0 / nComponents);
365 }
366 }
367
368 private static final double DEFAULT_THRESH = 1e-2;
369 private static final double DEFAULT_MIN_COVAR = 1e-3;
370 private static final int DEFAULT_NITERS = 100;
371 private static final int DEFAULT_NINIT = 1;
372
373 CovarianceType ctype;
374 int nComponents;
375 private double thresh;
376 private double minCovar;
377 private int nIters;
378 private int nInit;
379
380 private boolean converged = false;
381 private EnumSet<UpdateOptions> initOpts;
382 private EnumSet<UpdateOptions> iterOpts;
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409 public GaussianMixtureModelEM(int nComponents, CovarianceType ctype, double thresh, double minCovar,
410 int nIters, int nInit, EnumSet<UpdateOptions> iterOpts, EnumSet<UpdateOptions> initOpts)
411 {
412 this.ctype = ctype;
413 this.nComponents = nComponents;
414 this.thresh = thresh;
415 this.minCovar = minCovar;
416 this.nIters = nIters;
417 this.nInit = nInit;
418 this.iterOpts = iterOpts;
419 this.initOpts = initOpts;
420
421 if (nInit < 1) {
422 throw new IllegalArgumentException("GMM estimation requires at least one run");
423 }
424 this.converged = false;
425 }
426
427
428
429
430
431
432
433
434
435 public GaussianMixtureModelEM(int nComponents, CovarianceType ctype) {
436 this(nComponents, ctype, DEFAULT_THRESH, DEFAULT_MIN_COVAR, DEFAULT_NITERS, DEFAULT_NINIT, EnumSet
437 .allOf(UpdateOptions.class), EnumSet.allOf(UpdateOptions.class));
438 }
439
440
441
442
443
444
445
446
447
448
449 public boolean hasConverged() {
450 return converged;
451 }
452
453
454
455
456
457
458
459
460
461
462 public MixtureOfGaussians estimate(Matrix X) {
463 return estimate(X.getArray());
464 }
465
466
467
468
469
470
471
472
473
474
475 public MixtureOfGaussians estimate(double[][] X) {
476 final EMGMM gmm = new EMGMM(nComponents);
477
478 if (X.length < nComponents)
479 throw new IllegalArgumentException(String.format(
480 "GMM estimation with %d components, but got only %d samples", nComponents, X.length));
481
482 double max_log_prob = Double.NEGATIVE_INFINITY;
483
484 for (int j = 0; j < nInit; j++) {
485 gmm.gaussians = ctype.createGaussians(nComponents, X[0].length);
486
487 if (initOpts.contains(UpdateOptions.Means)) {
488
489 final DoubleKMeans km = DoubleKMeans.createExact(nComponents);
490 final DoubleCentroidsResult means = km.cluster(X);
491
492 for (int i = 0; i < nComponents; i++) {
493 ((AbstractMultivariateGaussian) gmm.gaussians[i]).mean.getArray()[0] = means.centroids[i];
494 }
495 }
496
497 if (initOpts.contains(UpdateOptions.Weights)) {
498 gmm.weights = new double[nComponents];
499 Arrays.fill(gmm.weights, 1.0 / nComponents);
500 }
501
502 if (initOpts.contains(UpdateOptions.Covariances)) {
503
504 final Matrix cv = MeanAndCovariance.computeCovariance(X);
505
506 ctype.setCovariances(gmm.gaussians, cv);
507 }
508
509
510 final TDoubleArrayList log_likelihood = new TDoubleArrayList();
511
512
513 converged = false;
514 double[] bestWeights = null;
515 MultivariateGaussian[] bestMixture = null;
516 for (int i = 0; i < nIters; i++) {
517
518 final IndependentPair<double[], double[][]> score = gmm.scoreSamples(X);
519 final double[] curr_log_likelihood = score.firstObject();
520 final double[][] responsibilities = score.secondObject();
521 log_likelihood.add(ArrayUtils.sumValues(curr_log_likelihood));
522
523
524 if (i > 0 && Math.abs(log_likelihood.get(i) - log_likelihood.get(i - 1)) < thresh) {
525 converged = true;
526 break;
527 }
528
529
530 mstep(gmm, X, responsibilities);
531
532
533 if (nIters > 0) {
534 if (log_likelihood.getQuick(i) > max_log_prob) {
535 max_log_prob = log_likelihood.getQuick(i);
536 bestWeights = gmm.weights;
537 bestMixture = gmm.gaussians;
538 }
539 }
540
541
542
543 if (Double.isInfinite(max_log_prob) && nIters > 0) {
544 throw new RuntimeException(
545 "EM algorithm was never able to compute a valid likelihood given initial " +
546 "parameters. Try different init parameters (or increasing n_init) or " +
547 "check for degenerate data.");
548 }
549
550 if (nIters > 0) {
551 gmm.gaussians = bestMixture;
552 gmm.weights = bestWeights;
553 }
554 }
555 }
556
557 return gmm;
558 }
559
560 protected void mstep(EMGMM gmm, double[][] X, double[][] responsibilities) {
561 final double[] weights = ArrayUtils.colSum(responsibilities);
562 final Matrix resMat = new Matrix(responsibilities);
563 final Matrix Xmat = new Matrix(X);
564
565 final Matrix weighted_X_sum = resMat.transpose().times(Xmat);
566 final double[] inverse_weights = new double[weights.length];
567 for (int i = 0; i < inverse_weights.length; i++)
568 inverse_weights[i] = 1.0 / (weights[i] + 10 * MathUtils.EPSILON);
569
570 if (iterOpts.contains(UpdateOptions.Weights)) {
571 final double sum = ArrayUtils.sumValues(weights);
572 for (int i = 0; i < weights.length; i++) {
573 gmm.weights[i] = (weights[i] / (sum + 10 * MathUtils.EPSILON) + MathUtils.EPSILON);
574 }
575 }
576
577 if (iterOpts.contains(UpdateOptions.Means)) {
578
579 final double[][] wx = weighted_X_sum.getArray();
580
581 for (int i = 0; i < nComponents; i++) {
582 final double[][] m = ((AbstractMultivariateGaussian) gmm.gaussians[i]).mean.getArray();
583
584 for (int j = 0; j < m[0].length; j++) {
585 m[0][j] = wx[i][j] * inverse_weights[i];
586 }
587 }
588 }
589
590 if (iterOpts.contains(UpdateOptions.Covariances)) {
591 ctype.mstep(gmm, this, Xmat, resMat, weighted_X_sum, inverse_weights);
592 }
593 }
594
595 @Override
596 public GaussianMixtureModelEM clone() {
597 try {
598 return (GaussianMixtureModelEM) super.clone();
599 } catch (final CloneNotSupportedException e) {
600 throw new RuntimeException(e);
601 }
602 }
603 }