Gradient-based training of Gaussian Mixture Models in High-Dimensional Spaces
Authors: Alexander Gepperth
Benedikt Pfülb
Venue:  (2020)
Keywords: GMM
SGD

Abstract: We present an approach for efficiently training Gaussian Mixture Models (GMMs) with Stochastic Gradient Descent (SGD) on large amounts of high-dimensional data (e.g., images). In such a scenario, SGD is strongly superior in terms of execution time and memory usage, although it is conceptually more complex than the traditional Expectation-Maximization (EM) algorithm.
For enabling SGD training, we propose three novel ideas:
First, we show that minimizing an upper bound to the GMM log likelihood instead of the full one is feasible and numerically much more stable way in high-dimensional spaces.
Secondly, we propose a new regularizer that prevents SGD from converging to pathological local minima.
And lastly, we present a simple method for enforcing the constraints inherent to GMM training when using SGD.
We also propose an SGD-compatible simplification to the full GMM model based on local principal directions, which avoids excessive memory use in high-dimensional spaces due to quadratic growth of covariance matrices.
Experiments on several standard image datasets show the validity of our approach, and we provide a publicly available TensorFlow implementation.

Under review as a conference paper at ICLR 2020
GRADIENT-BASED TRAINING OF GAUSSIAN MIXTURE MODELS IN HIGH-DIMENSIONAL SPACES
Anonymous authors Paper under double-blind review
ABSTRACT
We present an approach for efficiently training Gaussian Mixture Models (GMMs) with Stochastic Gradient Descent (SGD) on large amounts of high-dimensional data (e.g., images). In such a scenario, SGD is strongly superior in terms of exe- cution time and memory usage, although it is conceptually more complex than the traditional Expectation Expectation-Maxiiizzaiii Maximization (EM) algorithm. For enabling SGD train- ing, we propose three novel ideas: First, we show that minimizing an upper bound to the GMM log likelihood instead of the full one feasible and numerically much more stable way in high-dimensional spaces. Secondly, we propose new an- nealing procedure that prevents SGD from converging to pathological local min- ima. We also propose an SGD-compatible simplification to the full GMM model based on local principal directions, which avoids excessive memory in high- dimensional spaces due to quadratic growth of covariance matrices. Experiments several standard image datasets show the validity of our approach, and we pro- vide a publicly available TensorFlow implementation.
INTRODUCTION
This contribution is in the context of Gaussian Mixture Models (GMMs), which is probabilistic unsupervised method for clustering and data modeling. GMMs have been used in a wide range of scenarios, e.g., Melnykov al. (2010). Traditionally, free parameters of GMM are estimated (EM) algorithm, which has the appealing property of
using the so-called poppectation-Maxiiiann requiring learning rates and which automatically enforces all constraints that GMMs impose.
1.1 MOTIVATION
GMMs have several properties that make their application to image data appealing, like the ability to generate samples from the modeled distribution, or their intrinsic outlier detection ability. However, for such high-dimensional and data-intensive scenarios which require high degree of paralleliza- tion to efficiently solvable, the traditional way of training GMMs by EM is reaching its limits, both in terms of memory consumption and execution speed. Memory requirements skyrocket mainly because GMM intrinsically batch-type algorithm that needs to store and process all samples in memory in order to be parallelizable. In addition, computation time becomes excessive because EM involves covariance matrix inversions, which become prohibitive for high data dimensionality. Here, an online mini-batch type of optimization such Stochastic Gradient Descent (SGD) has strong advantages w.r.t. the memory consumption in particular, but also w.r.t. computation time since it re- quires no matrix inversions. Therefore, we seek develop a Stochastic Gradient Descent algorithm for GMMs, making it possible to efficiently train GMMs on large collections of images. Addition- ally, since the number of free GMM parameters depends quadratically on the data dimensionality D, we seek to develop simplified GMM models with a less severe impact of D.
1.2 RELATED WORK
There are several proposals for what is termed "incremental" or "online" EM for training GMMs
Vlassis & Likas ) Engel & Heinen (2010); Pinto & Engel (2015); Cederborg al. (2010)
Song & Wang (2005) Kristan et al. (2008). allowing GMMs to be trained by providing one sample
time. However, none of these models reproduce the original GMM algorithm faithfully, which is

---

Under review as a conference paper at ICLR 2020
why most of them additionally resort to component pruning merging heuristies for components, leading to learning dynamics that are difficult to understand. This is also the case for the Locally Weighted Projection Regression (LWPR) algorithm Vijayakumar et al. (2005) which can be inter- preted as GMM model adapted for online learning. To our knowledge, there is only a single work proposing train GMMs by SGD Hosseini & Sra (2013). Here, constraint enforcement ensured by using manifold optimization techniques and carefully chosen set of regularizers, leading to rather complex model with several additional and hard-to-choose hypee-parameters. The idea of in- cluding regularization//mmaaiin into GMMs (in the context of EM) was proposed in Verbeek et al. (2005); Ormoneit & Tresp (1998), although their regularizes significantly differs from ours. Ap- proximating the full GMM log-likelihood is used in Verbeek et al. (2003); Pinheiro & Bates (1995); Dognin et al. (2099,, although together with EM training. There has been significant work on simpli- fying GMM models (e.g.. Jakvoljvvvv
leading to what is called parsimonious GMM models, which is reviewed in Bouveyron & Saumee- (2014).
1.3 CONTRIBUTIONS
The main novel contributions of this article are:
generic and novel method for training GMMs using standard SGD that enforces all constraints and is numerically stable for high dimensional data (e.g., D > 500)
annealing procedure that ensures convergence from wide range of initial conditions in streaming setting (excluding, e.g., centroid initialization by k-means)
simplification of the basic GMM model that is particularly suited for applying GMMs to high-dimensional data
an analysis of the link between SGD-trained GMM models and self-organizing maps
Apart from these novel contributions, we provide publicly available TensorFlow implementation
2 DATASETS
We use two datasets for all experiments. MNIST (LeCun et al. 1998) consists of gray scale images of handwritten digits (0-9) and is common benchmark for computer vision systems and classifica- tion problems. SVHN (Netzer et al. 2011) is a 10-class benchmark based on color images of house numbers (0-9). We use the cropped and centered digit format. For both benchmarks, we use either the full images (32 x 32 x 3 for SVHN, 28 x 28 x for MNIST), or crop of the central 6 x 6 patch of each image. The class information not used since GMM is purely unsupervised method.
METHODS: STOCHASTIC GRADIENT DESCENT FOR GMMS
Gaussian Mixture Models (GMMs) normally formulated in terms of K component probabilities,
modeled by multi-variate Gaussians. It is assumed that each data sample æ, represented as a single row of the data matrix X, has been sampled from single Gaussian component k, selected with a priori probability тą. Sampling is performed from the Gaussian conditional probability k())k), whose parameters are the centroids f4 and covariance matrices Σ₁: p(z)k) N(x;; μ,, E). A GMM aims to explain the data by minimizing the log-likelihood function
L
log TTP(w,,
Free parameters to be adapted are the component probabilities Tee the component centroids μ4 and the component covariance matrices Σg.
3.1 MAX-COMPONENT APPROXIMATION
As a first step to implement Stochastic Gradient Descent (SGD) for optimizing the log-likelihood (see Eq. [}. we observe that the component weights T. and the conditional probabilities p(æ)k) are
https: http::///gtthbbb.ccmm/ooooooooe icl r200GMM......
2

---

Under review as a conference paper at ICLR 2020
positive by definition. It is therefore evident that any single component of the inner sum the components & is a lower bound for the whole inner sum. The largest of these K lower bounds is given by the maximum over the components, thus giving
Σ log Σ Txp(anl() S LMC
Σ N
1
log max; (TTP(2n)1))
N
1
Σ max; log (Tap(æn)):)).
N
This is what we term the max-component approximation of Eq. Since LMC 2 L. we can decrease by minimizing LMC- The advantage of LMC that not affected by numerical instabilities as L is. An easy way to see this is by considering the case of a single Gaussian component with weight Tą = centroid μt = and whose covariance matrix Dg has only diagonal entries which are all equal to 0.1. In this case, an 1000ddimensional input vector æ; 1 Vi would give conditional probability of p(ar|k [k) = (2x0.1) (0001)-5000 e 5000 which will, depending on floating point in- accuracy, produce NaN infinity for the first factor due to overflow, and zero for the second factor due to underflow, resulting in an undefined behavior.
3.2 CONSTRAINT ENFORCEMENT
GMMs require the weights to be normalized: Σg TA= and elements of the diagonal matrices to be non-negative: Eig 2 O Vi, j. In traditional GMM training via Expectation-Maxiizzatinn these constraints are taken care of by the method of Lagrangian multipliers. In SGD, they must be we adopt the approach proposed in Hosseini & Sra (2015), which replaces them by other free parameters & from which the тTş are computed such that normalization is ensured. One such computational scheme is
enforced in different manner. For the weights
exp(Ex)
TA Σj exp(E,)
(3)
For ensuring non- negativeeves of the covariances, identify all covariances inferior to certain minimal value ymm after each gradient descent step, and subsequently clip them to min
3.3 ANNEALING PROCEDURE
A major issue for SGD optimization of GMMs are local minima that obviously correspond to unde- sirable solutions. Prominently among these is what we term the single-component solution: here, a single component k" has weight close 1, with its centroid and covariance matrix being given by the mean and covariance of the data:
(4) Another undesirable local minimum is the degenerate solution where all components have the same weight, centroid and covariance matrix, the latter being identical to the Lngle-componen solution. See Fig. (left) for a visualization of the latter.
T 1, P1 = E[]] Σx- = (ov(X).
Such solutions usually evolve from an incorrect initialization of the GMM model, and it is easy to show that the gradient w.r.t. all model parameters is zero in both of these scenarios. This is easy to show; As the gradients read (with Tin denoting the responsibility of component / for sample n)
ØC
1
@pm
=
N
Σ
Ti (Irm
p1m)
(5)
Oc
Σ
2 Enn
(Inm Pmm )2
(6)
N
it is easy to see that when either all or single centroid/covariance matrix are equal to the total mean and variance of the data (J4m y From and Σim (arm )2m,))) obtain zero gradients for
degenerate solutions (T,, = 7/4 = 1/ K) and single-cone compoeet solutions (Elm,, : 7/- = 1), respectively.
Our approach for avoiding these undesirable solutions in principle is punish their characteristic response patterns by an appropriate modification the loss function that is minimized, i.c., LMc- The following guidelines should be adhered to:
3

---

Under review as a conference paper at ICLR 2020
single-comppnent and degenerate solution should be punished by the annealing procedure
mainly wish to ensure good starting point for optimization, annealing should be active only at the beginning of the training process
as training advances, annealed loss should smoothly transition into the original loss
the number of free parameters introduced by the annealing procedure should be low
These requirements are met by what we term the smoothed max-component log-likelihood LMC:
Σ max: ( Σ ghi log (m)p(cnlli)) N
1
LMC
Here, we assign normalized coefficient vector gt. to each Gaussian mixture component k. The entries of g are computed in the following fashion:
Assume that the K Gaussian components are arranged on a ID grid of dimensions (1 or 2D grid of dimensions (vK), > Each linear component index k has thus an unique associated ID or 2D coordinate c(k).
Assume that the vector gl of length K is actually representing a ID structure dimension (1, K) or 2D structure dimension (v K" VKK). Each linear vector index g has thus unique associated ID or 2D coordinate c(j).
The entries of the vector 9k computed
9kj = exp (-(cL) - e(k:))
(8)
20
and subsequently normalize to have unit sum. Thus, Eq. essentially represents convolution of the probabilities k)))k), arranged in an ID or 2D grid, with Gaussian convolution filter, thus amounting smoothing operation. The ID or 2D variance in Eq. iss parameter that must be set as a function of the grid size such that Gaussians are neither homogeneous nor delta peaks. Thus, the loss function in Eq. is maximized if the log probabilities follow an unimodal Gaussian profile of variance a, whereas single-component and degenerate solutions are punished.
is trivial to see that the annealed loss function in Eq. reduces the non-anneeled form Eq. 2 in the limit where a - O. This is because the vectors It approach Kronecker deltas in this case, with only single entry of value 1, which essentially removes the inner sum Eq. 7 By making time-dependent, starting at intermediate value of a (to)) =0 oo and then let it approach small final value a (too) = doc. we can smoothly transition the annealed loss function Eq. Eq. to the original max-component log- lig-likeldoood Eq. Time dependency of (t)) thus be chosen be:
do
t < to >> loa
(tt)
Too To exp (TT)
t0 < 1 < Loo
where the time constant in the exponential is chosen as = log ao-ax
too to ensure a smooth transition.
3.4 FULL GMM MODEL TRAINED BY SGD AND ITS HYPER- PARAMETERS
Putting everything together, training Gaussian Mixture Models with SGD performed by minimiz- ing the smoothed max-component log-likelibe lodliikel of Eq. 7 enforcing the constraints on the com- ponent weights and covariances as detailed in Sec. 3.2 and transitioning from the smoothed to the "bare" loss function detailed in Sec. Sec. In order to avoid convergence to other undesirable local minima by SGD, we introduce three weighting constants X z and Ay, all the [0, range, for controlling the relative adaptation speeds for the three groups of free parameters in the GMM models
4

---

Under review as a conference paper at ICLR 2020
(setting, e.g., X 0 would disable the adaptation of weights). We have thus the update rules:
DCmc
(x + 1) (t( - c. AM
MLMC and
Na (1 mb(()) - A
(10)
0pp
dLMC
Ex + Ex(t) x 0
Centroids are initialized random values in range I-p', weights are chosen equiprobable, that is to say 1/K, for all components, and covariance matrix entries are uniformly initialized to very small positive values E (a good ad-hoc choice is 0.15). Last but not least, a learning rate is required for performing SGD. We find that our algorithm depends very weakly <= 0.01 is usually a good choice. We summarize good practices for choosing nypee-parameters of the proposed SGD approach for GMMs in App. [ component centroids by k-means as is usually recommended when training GMMs by EM.
and that fixed value of
Please not that it is not required to initialize the
3.5 SIMPLIFICATION AND MEMORY/PERFFOMMNT ISSUES
For the full GMM model, regardless of whether it is trained by EM SGD, the memory quirements a GMM with K components depend strongly on the data dimensionality D. A good measure for memory requirements is the number of free model parameters M computed as M K KD2 + D + 1). To reduce this number, we propose several strategies that are compatible with training by SGD (see also Fig. D
Figure I: Two possibilities for simplifying GMMs, illustrated for dataset with two clusters and two Gaussian components. Green dots represent data samples and arrows represent normalized local principal directions along which variances (represented by ellipses) covariance matrix, where it is assumed that clusters vary only in the coordinate directions (not the case here). Right: non-diagonal covariance matrix where only limited number (here one) of principal directions is considered for variance computation.
computed. Left: diagonal
Diagonal Covariance Matrix A well-known simplification of GMMs consists of using diago- nal covariance matrix instead of a full one. This reduces the number of trainable parameters to K(2 1), where D is the dimensionality of the data vectors and K the number of components.
Local Principal Directions A potentially very powerful simplification consist of using a diagonal covariance matrix of S D diagonal entries, and letting SGD adapt the local principal directions along which these variances are measured for each component (in addition to the other free param- eters, course). The expression of the conditional component probabilities now includes, for each component k, the S normalized principal direction vectors dks; = 1, S and reads:
x dkx)
kk
exp
(11)
V (2 STI
Ens
32
This leads value of M = ((SS +)D+1), and thus achieves satisfactory model simplification only if S can be chosen very small and is high. order for the model to fully determined, another constraint needs to be enforeed after each SGD step. namely the orthogonality of the local
principal directions for each component A: P Pr = 6im VI, m < S. For now, the matrix D, whose the S vectors dys is subjected QR decomposition D Q R2 and setting D Qi-
5

---

Under review as a conference paper at ICLR 2020
4 EXPERIMENTS
Unless otherwise stated, the experiments in this section will always be conducted with the follow- ing parameter values (see App. C for justification of these choices): total iterations T- 1 1000
v VK. to 0.3T, &o 0.8T, mini- size B = μ 0.01, dg 1.2, x 0.011, min 0.15 and K = 5 5. The adaptation strengths Ag. A and X are all set to 1.0. Except for Sec. 3.5 we will use diagonal covariance matrix for each of the components. Training/test data are taken from MNIST or SVHN, either using full images or the central 6 x 6 patches.
0.011,
7
2
6 8
9
3 5
3
Exemplary results for learned centroids, using normal experimental conditions, trained on variants of the two datasets (from left to right): MNIST, MNIST patches, SVHN, SVHN patches.
Figure 2:
4.1 VALIDITY THE MAX-COMPONENT APPROXIMATION AND COMPARISON EM
Since we minimize only an approximation to the GMM log-likelihood, first important question is about the quality this approximation. this end, we plot both energy functions over time when training on full MNIST and SVHN images, see Fig. 3 This figure shows that the log- log-likelihood obtained by EM-based training (using the implementation of skleam) with the same configuration (using 10 iterations and k-means-based initialization) achieves very similar bog-likelihood, which confirms that SGD indeed comparable in performance.
apprctuatiion Boood
M-C apoadeuciio
lng-Ibdlhedd
KS1
1000
30
Jiees
201
180
10
120
SGD iteration
SGD iteration
SGD iteration
Figure 3: Comparing the max nax-component approximation to the negative log-likelihooe during SGD training, shown for MNIST (left), MNIST patches (center) and SVHN patches (right).
4.2 EFFECTS AND BENEFITS OF ANNEALING
In order to demonstrate the beneficial effects the annealing procedure, we conduct the basic experiment on full MNIST with annealing effectively turned off, which is achieved by setting 0==0 + 0.001 0.011. As a result, in Fig. 4 we observe typical single-component solutions when visualizing the centroids and weights, along with markedly inferior (higher) loss value of s 155, in contrast to == 138 when annealing is turned on. As side effect, annealing enforces a topo- logical ordering of centroids in the manner of SOM (see also App. for mathematical analysis), which is only of aesthetic value for the moment, but might be exploited in the future.
4.3 BASIC FEASIBILITY AND PARAMETER SPACE
Here, we put everything together and train an GMM using SGD on various datasets very high and intermediate dimensionality (results: Fig. 2 and App. D) In particular, analyze how the choice of the relative adaptation strengths \y, Ay and number of components affects the final log-likelihood. The results presented in Tab. T suggest that increasing K is always beneficial, and

---

Under review as a conference paper at ICLR 2020
9
2
/
9
3
5
8
7
O
9
6
2
8
4
Figure 4: Effects of annealing SGD convergence to centroids MNIST. Left: single-componen
solution (no annealing). Right: regular solution (with annealing).
that the covariances and weights should be adapted more slowly than the centroids. The batch size B had an effect convergence speed but significant impact on the final loss value (not shown here).
Table I: Results for various hyper-parameter parameter values, always given as the final max- nax-component log- likelihood Lvc- The full log- ligllikll values were always within 0.5 of the given value, and smaller by construction. Each cell contains value for the MNIST and one for the SVHN dataset.
Parameter
K 25
K 36
49
K 64
Ax Ag = 1.0,
B ==1
138.1, 735.2
136.3, 704.7
133.3, 693.5
134.1, 682.9
Ax Ag 0.1,
==1
138.3. 715.6
134.6.
706.5
691.3 133.5,
130.1, 691.4
Ax, xg = 0.0,
==1
138.1, 722.5
135.5, 714.3
133.2, 689.1
130.8, 688.4
4.4 ROBUSTNESS TO INITIAL CONDITIONS
MBC appooduutiia lopilizzhood
36C apeodaatto loplisthoo
24C1 appostration
lopela
EM
EM
0
300
140
100
2000
20,0
0000
2000 SCD itcration
0000
SGD iteration
SGD iteration
Figure 5: Comparing different initialization ranges for centroids: [-0.5, 0.5] (left), [-0.7, 0.7] (center) and [-1, 1] (right). The orange line represents the same value in all diagrams, the seeming difference is due to scaling as initial loss values some experiments are much higher (i.e., worse).
To test whether SGD converges as function of the initial conditions, we train three GMMs on MNIST using three different initializations: j² E {1.0, 0.7, All of these choices led non- convergence of the EM-based GMM implementation of sklearn. Differently from the other exper- iments, we choose a learning rate learning. The alternative is to double the training time which works well. Fig. 5 shows the de- velopment of the loss function in each of the three cases, and we observe that all three cases indeed converge by comparing the final loss values to Fig. F (middle) and to the EM baseline. This behavior persists across all datasets, but only MNIST is shown due to space limitations.
0.02 since the original value does not lead to convergent
4.5 SGD TRAINING OF SIMPLIFIED GMMS WITH PRINCIPAL LOCAL DIRECTIONS
We train a GMM with E (26, 50, 100} local principal directions on MNIST and SVHN patches and log the final loss function values in order to assess model quality. We find that learned centroids and final loss values approach those of the diagonal model for > 20 for both datasets. Please see App. for visualization of the final centroids as function of S. For the full datasets, find that

---

Under review as a conference paper at ICLR 2020
MNIST and SVHN require values of S 2 100 and S 2 300, respectively for achieving similar loss values as the diagonal model.
5 DISCUSSION
We showed that GMMs can be trained on high-dimensionn image data by SGD, at very small batch sizes, in a purely online fashion. Memory and execution time requirements are very modest, and can be reduced even more by using intelligent simplifications of the GMM model. Here, we would like to discuss few noteworthy points and suggest some avenues for improving the model:
Robust Convergence Unlike EM approaches which require a careful selection of starting point, the presented SGD scheme is very robust to initial conditions due the proposed annealing scheme. In all experiments with high-dimensional image data, we found no initializations that did not, given reasonable time, converge to a regular solution.
Numerical Stability While we are not optimizing the full log-likelihood here, the max-component approximation is actually good one, and has the advantage of absolute numerical stability. Clearly, clever ways might be found to avoid instabilities in the log-likelihood function computation, but when it comes to automatically computed gradients this is no longer possible, and indeed we found that it was mainly the (inaccessible) gradient computations where NaN values happened mostly.
Complex annealing procedure The complex way of regularizing the model, which is similar in spirit to simulated annealing approaches, may seem daunting, although the time constants can be determined by simple heuristics, see App. [ A next step will be to replace this procedure by an automated scheme based on evaluation of the loss function: when it is stationary, can be decreased slightly. Repeating this strategy until a small minimal value is reached should remove the need to fix annealing parameters except for time constants which can be chosen function of dataset size.
Free Parameters Apart from the annealing mechanism, the model contains the same parameters any SGD approach would, the learning rate and the weighting constants for the model parameters: Xx, An and Ax. We found no experimental evidence suggesting that the latter require values differ- ent from 1, although from general principles it might be sensible to adapt the variances and weights slower than the centroids. Further, is required understand how these constants can be used to obtain better solutions. Interesting is that standard DNN optimizers like Adam Kingma & Ba (2014) did not seem to produce satisfactory solutions which also requires looking into.
Approximate Gradient Descent The procedure described in Sec. 3.4 performs gradient descent on an approximation to the full log-likelihood Eq. [T] This might seem unsatisfactory, so we would like to point out that the original EM procedure does nothing else: it optimizes lower bound of the log-likelihoode In particular, the M-step is not guaranteed to improve the log-likelihood all: all we know is that it will never make it worse. In our SGD approach, there is simple way to fix this shortcoming, which we will tackle as an next step: we can replace the max-operation in Eq. b a softmax function with initially large steepness parameter which will make it indistinguishable from a discrete max operation. After convergence, w can be relaxed, which will result in a smooth transition to optimizing the full log- likeeiiooood but from an already well-converged state.
CONCLUSION AND OUTLOOK
On more abstract level, we developed the presented method because of our interest in continual or incremental learning, which essentially about eliminating the well-known catastrophic forgetting effect. We believe GMMs are essential building block for continual learning models since their parameter updates are purely local in the sense that only components that are close to the currently best-matching component are updated. Next steps will consist of constructing deep network classi- fiers from convolutional GMM layers, with a readout layer on top (all layers being trained by SGD). Further, we will investigate how to sample efficiently from such hierarchical convolutional GMMs, allow generating large batches of samples in a replay-based architecture for continual learning.
8

---

Under review as a conference paper at ICLR 2020
REFERENCES
Charles Bouveyron and Camille
Model-based clustering of high high--dimensional data: A review. Computational Statisties & Data Analysis, 71:52-78, 78, 2014.
Thomas Cederborg, Ming Li, Adrien Baranes, and Pierre Yves Oudeyer. Incremental local online gaussian mixture regression for imitation learning of multiple tasks. In 2010 JEEE/RSJ Interna- tional Conference on Intelligent Robots and Systems, pp. 267-274, IEEE, 2010.
Pierre I. Dognin, Vaibhava Goel, John R Hershey, and Peder A Olsen. A fast, accurate approxi- mation to log likelihood of gaussian mixture models. In 2009 IEEE International Conference on Acoustics, Speech and Signal Processing, pp. 3817-3820. IEEE, 2009.
Paulo Martins Engel and Milton Roberto Heinen. Incremental learning of multivariate gaussian mixture models. In Antônio Carlos da Rocha Costa, Rosa Maria Vicari, and Flavio Tonidandel (eds.), Advances in Artificial Intelligence SBIA 2010, pp. 82-91, Berlin, Heidelberg, 2010. Springer Berlin Heidelberg. ISBN 978-3-642-16138-4
Tom Heskes. Energy functions for self-organizing maps. Kohonen maps, pp. 303-315. Elsevier, 1999.
Reshad Hosseini and Suvrit Sra. Matrix manifold optimization for gaussian mixtures. In Advances in Neural Information Processing Systems, pp. 910-918, 2015.
Nikša M Jakovljević. Gaussian mixture model with precision matrices approximated by sparsely represented eigenvectors. In 2014 22nd Telecommunications Forum Telfor (TELFOR), pp. 435- 440. IEEE, 2014.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization, 2014. URL http::/arxiv. wwtgabs..///oooooooot.../o..eest /0bs/1410 6980 cite arxiv: 6980ommmttt Published as con- ference paper at the 3rd International Conference for Learning Representations, San Diego, 2015.
Teuvo Kohonen. The self-organizing map. Proceedings of the IEEE, 78(9):1464 1990.
Matej Kristan, Danijel Skocaj, and Aleš Leonardis. Incremental learning with gaussian mixture models. In Computer Vision Winter Workshop, pp. 25-32, 2008.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11 (6(11):2277-2322 1998.
Volodymyr Melnykov. Ranjan Maitra, et al. Finite mixture models and model-based clustering.
Statistics Surveys, 416, 2010.
Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, and Andrew Ng. Reading digits in natural images with unsupervised feature learning. In NIPS workshop on deep learning and unsupervised feature learning, volume 2011, pp. 2011.
Dirk Ormoneit and Volker Tresp. Averaging, maximum penalized likelihood and bayesian estima- tion for improving gaussian mixture probability density estimates. IEEE Transactions on Neural Networks, 9(4):639-650, 1998.
Pfülb and A Gepperth. A Comprehensive, Application-orienan Study of Catastrophic Forgetting in DNNs. In International Conference on Learning Representations (ICIR), 2019.
José C Pinheiro and Douglas M Bates. Approximations to the likeliiid function in the nonlin- car mixed-effects model. Journal of computational and Graphical Statistics, 12-3: 1995.
Rafael Coimbra Pinto and Paulo Martins Engel. A fast incremental gaussian mixture model. PloS one, 10(10):e01399)) 139901,, 2015.
Mingzhou Song and Hongbin Wang. Highly efficient incremental estimation of gaussian mixture models for online data stream clustering. In Intelligent Computing: Theory and Applications III, volume 5803, pp. 183. International Society for Opties and Photonies, 2005.

---

Under review as a conference paper at ICLR 2020
Jakob J Verbeek, Nikos Vlassis, and Ben Kröse. Efficient greedy learning of gaussian mixture models. Neural computation, 15(2):469 485, 2003.
Jakob J Verbeek, Nikos Vlassis, and Ben JA Kröse. Self-organizing mixture models. Neurocomput ing, 63:99-123, 2005,
Sethu Vijayakumar, Aaron Dosuuza, and Stefan Schaal. Incremental online learning high dimen- sions. Neural computation, 17(12):2602-2663 2005.
Nikos Vlassis and Aristidis Likas. A greedy em algorithm for gaussian mixture learning. Neural Processing Letters 15 Feb 2002. ISSN 1573-773X, doi: 10.1023/A1013384411111 URL /ttt::///ooll.ooo///000000011111
10

---

Under review as a conference paper at ICLR 2020
A LOCAL PRINCIPAL DIRECTIONS
Here, we give visualization of the resulting centroids for various choices of when training on MNIST patches.
Figure 6: Comparing different values of S when training GMMs with the local principal directions
simplification. Shown are centroids for S = 25 (left), S = 50 (middle) and S = 100 (right).
B RIGOROUS LINK TO SELF-ORGANIZING MAPS
is known that the oelf-organizing map (SOM, Kohonen 1990)) rule has no energy function is minimizing. However, some modifications (see Heskes 1999)) have been proposed that ensure the existence 8 energy function. These energy-baset SOM models reproduce all features of the original model and very closely. In the notation of this article, SOMs model the data through K prototypes μı and K neighborhood functions gs defined a periodic 2D grid, and their energy function is written
learning rule that the original SOM algorithm is actually approximating
maxk Σ Dky Eaglla |& 1 Σ N
LsOM
(12)
Discarding constant terms, we find that this is actually identical to the max-component log- likelihood approximation given Eq. [ with constant equiprobable weights тı and constant diagonal Σ with equal diagonal entries. To our knowledge, this is the first time that rigorous link between SOMs and GMMs has been established based on comparison energy functions, showing that SOMs actually implementing annealing-type -ype approximation to the full GMM model, with component weights and variances chosen in a special way.
C RULES-OF-THUND FOR SGD-TRAINING OF GMMS ON IMAGES
When training DNNs by SGD, several parameters need be according to heuristics. Here, we present some rules-of-thumb for performing this selection with GMMs. Generally, we can always our experiments, it was always feasible However, tuning these parameters
set the batch size 1, Mmin 0.15 and 0.01. to set to to half epoch, and too to to half epoch, can speed up training significantly. The number of prototypes critical choice, but follows simple "the more the better" rule. From the mathematical foundations of GMMs, it is evident that more components must always be able to reach lower loss (except for possible local minima). The relative adaptation strengths A, Ag and A can always be set (at least we did not find any examples where that did not work).
1. Choose the number of iterations an increasing function of K: the more components, the more free parameters, the more time to train is needed
2. Flatten images, and put them as rows with random indices (as usual for SGD) into a tensor X
3. Choose quadratic number of components n
4. Initialize to n²/3.
5. Preliminary train the GMM with constant dg- Select to as the first time step where loss does no longer decrease. It never hurts if lo is chosen large, just takes longer.

---

Under review as a conference paper at ICLR 2020
6. Select foe by trial and error the train set: start with very large values and use the smallest value that gives the same final loss value. Larger values for to never hurt but take longer.
7. Train with the chosen parameters!
D MORE RESULTS
Visualization of prototypes training on different datasets (extracted random classes, see Pfülb & Gepperth (2019) for details about the datasets). Each result (in the form of centroid visualizations) recorded after training for epoch on given dataset with K (see) Figure Jand Figure 8
Figure 7: Visualization of prototypes training different datasets (from left to right): FashionM- NIST (15000 iterations), NotMNIST (132275 iterations) and MADBase (15000 iterations).
Figure 8: Visualization of prototypes training different datasets (from left to right): Devanagari (4500 iterations), EMNIST (86250 iterations, rotated and mirrored MNIST classes) and Fruits (1475 iterations).
12


Reviews:
[{"review_id": "rkgY7q_4cS", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper tackles the problem of online learning for GMMs, in the context of high-dimensional data.  Specifically, the authors limit the scope to SGD-like approaches and EM-like optimization.  They also offer a TF implementation.\n\nI feel that this work is largely incremental, but more importantly indicating the authors' lack of understanding of the very long history of (online) EM.  While the authors do acknowledge some of the online EM work, they go on to develop what is a rather ad-hoc approach to online EM.\n\nThe max-component approximation in Sec. 3.1 is claimed to address the issue of numerical stability.  The authors do not appear to resort to the log-sum-exp \"trick\", which tackles such problems.  (In fact, their max approx is of this type.)\n\nSec. 3.2 uses a very standard representation of multinouli in terms of its natural parameters, which the authors again do not refer to.\n\nThe \"smoothing\" in Sec. 3.3 is hard to justify and difficult to understand, esp. the gridding approach.  Why not use hierarchical priors instead?\n\nIn Sec. 3.4, additional smoothing is accomplished using a subspace approach, which requires QR decomposition.  How will this affect computational efficiency, if the subspace needs to be recomputed?\n\nFinally, I have strong concerns about the experimental evaluation.  The authors choose datasets where the sample ambient space is at most 28x28, which is not exactly (very) high-dimensional.  \n\nI am mostly concerned about the evaluation.  ", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "HkeSb3jgqS", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper proposes a new method based on Stochastic Gradient Descent to train Gaussian Mixture Models, and studies this method especially in the context of high dimension. \nThe method is based on optimizing max-lower bound to the log-likelihood, together with a regularization term, using stochastic gradient descent. \nThe method is applied to two image datasets, and seem to produce sensible results. \n\nHowever, in my opinion, the results and presentation do not seem at the level suitable for publication.\n\nThere are no guarantees for convergence/other performance measures for the new method, in contrast to recent methods based on moment matching (e.g Ge, Huand and Kakade, 2015). Therefore, the new method and paper should provide excellent empirical results to be worthy of publication. \n\nThe authors write in the 'related work' section that GMM with regularization was proposed by [Verbeek et al. 2005], but it is an older idea - for example [Ormoneit&Tresp 1998] \n\nIn Section 3.1, the motivation for the maximization in eq. (2) is unclear. \nWhy is it easier to compute the gradient this way rather than keep the original likelihood?\nMoreover, the max operation is non-smooth and can cause problems with the definition of the gradient at some points. \nThe authors point to a problem of underflow/overflow when evaluating the gradient of the \nfull log-likelihood because the densities p(x | k) can be very small - but it is standard practice in EM to keep all probabilities multiplied by say their maximum p_max and keep log(p_max) separately, to avoid underflow problems. \n\nSection 3.2: I don't understand the requirement \\Sigma_{i,j} >=0. Is it a requirement for each entry of the covariance matrix? (which covariance matrix? there are K such matrices). The requirement should be that each matrix is positive-definite, not the entries.\n\nSection 3.3: The first sentence is wrong - EM can also suffer from the problem of local minima. \nAlso, the single-component solution doesn't seem like a local minimum - but rather the log-likelihood is unbounded here\nsince you can put another component on one of the data points with infinitely small variance. \nThe degenerate solution where all Gaussians have equal weights, mean and variance does not seem like a local minimum\nof the log-likelihood. Say the data really comes from 2 distinct Gaussians - then separating the two Gaussians means a bit would increase the log-likelihood. is not a local minimum. I'm not even sure if the gradient is zero at this point - the authors should show this. Maybe the authors mean their modified loss L_MC - this should be stated clearly. \n\nThe change of regularization during the training process seems like a heuristic that worked well for the authors, but it is thus unclear what optimization problem is the optimization solving. The regularization is thus included for optimization reasons, and not in the usual sense of regularization. \n\n\nThe expression for \\tau at the end of Section 3.3 seems wrong to me. I don't see how plugging it into eq. (7) gives a continuous \\sigma(t) function. \n\nSection 3.4: What is \\mu^i? I didn't see a definition \n\nI don't understand the local principal directions covariance structure. The authors write 'a diagonal covariance matrix of S < D entries). But what about the other D-S coordinates? are they all zero? or can have any values? \nThe parameters count lists (S+1)*D+1 for each Gaussian so I'm assuming S*D parameters are used for the covariance matrix, but it is unclear how. Eq. (9) has the parameters vectors d_{ks} for the principal directions, together with the \\Sigma_ss scalar values - it would be good to relate them to the mean and variance of the Gaussians. \n\n\nSection 4: The paragraph that describes the experimental details at the beginning is repeated twice. \n\nThe experimental results are not very convincing. The images in Figure 2,4 were picked by an unknown (manual?) criteria. \n\nIn the comparison to EM in Figure 3 there are missing details - which log-likelihood is used? L, L_MC? or different ones for different methods? is this test set log-likelihood? what fraction of the data was used? \nThere is also no comparison of running time between the two methods. \n\n\n", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "ByxznYJAFB", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper describes in detail a proper implementation of SGD for learning GMMs. GMMs are admittedly one of the basic models in unsupervised learning, so the topic is relevant to ICLR, even though they are not particularly a hot topic.\n\nThe paper is overall clear and well-written. The main contributions are an effective learning of GMMs from random initialization that is competitive (in terms of final loss) to EM training with K-means initialization. The authors discuss a max-component loss instead of standard likelihood for numerical stability, using a softmax reparametrization for component weights to ensure they are in the simplex, and an annealed smoothing of the max-likelihood based on arbitrarily embedding components indexes on a 2D regular grid. Experiments are shown on MNIST and SVHN comparing different hyper parameter settings to a baseline EM implementation from scikit-learn.\n\n- max-component: the use of the log(max p_k) instead of log(sum_k p_k) is sufficiently motivated to avoid well-known numerical problems that arise from directly computing the p_k rather than in log space. However, it is also a standard trick (e.g., in logsoftmax computations) to compute log(sum_k p_k) from log-probabilities using (for k* in argmax_k p_k):\n\nlog(sum_k p_k) = log(p_k*) + log(sum_k exp(log p_k - log p_k*))\n\nwhich is numerically more stable. Maybe I'm missing something, but I do not understand the advantage of the max-component compared to this, so it seems to me that the max-component trick is a fairly limited contribution in itself.\n\n- the main novelty seems to be the regularizer. The authors present experiments to show the effectiveness of it to avoid single-component solutions that may arise from random initialization, which is an interesting point of the paper. The motivation for smoothing on the 2D grid is still somewhat mysterious to me, even though the relationship with SOMs of Appendix B is interesting.\n\n- the paper falls a bit short in showing any improvement compared to other baselines. The experiments describe some ablations, but does not really answer some questions: is regularization important if one uses K-means initialization? are there any practical advantages (e.g. learning speed, hyperparamter tuning) compared to standard EM? The authors say that the selection of the starting point in EM is important. This is a fair point, but it also seems solved by K-means. While the authors describe a \"tutorial\" for choosing the hyper parameters of their method, it still seems fairly manual. So the practical advantage of the method (which probably exists) would benefit from more comparisons. \n\n- one of the difficulties in training GMMs comes from learning covariance matrices. While the authors discuss some way to train low-rank models, the only successful results seem to be with diagonal covariance matrices, which seems much easier. For instance, at least a toy example in which the low-rank version is useful would be interesting.\n\nOverall, it seems to me that the work is serious, and describes possibly interesting how-tos for training GMMs. The main contribution is to describe a simple method to learn GMMs from random initialization. The main technical novelty seems the regularizer, which seems to work. The method has the advantage of simplicity, but successful results have only been shown with diagonal covariance matrices and it is unclear exactly what is gained over EM+K-means initialization. \n\n\nother comments:\n- negative log-likelihood is used in the results section. It would be good to clarify it somewhere since the paper only mentions \"log-likelihood\" but report a loss that should be minimized\n\n- Section 4.4 \"Differently from the other experiments, we choose a learning rate of \u000f = 0.02 since the original value does not lead to convergent learning. The alternative is to double the training time which works as well.\" -> in the first sentence, I suppose it is more a matter of \"slow learning\" than \"non-convergent learning\"\n\n", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper presents an SGD-based learning of a Gaussian mixture model, designed to match a data streaming setting.

The reviews state that the paper contains some quite good points, such as
* the simplicity and scalability of the method, and its robustness w.r.t. the initialization of the approach;
* the SOM-like approach used to avoid degenerated solutions;

Among the weaknesses are
* an insufficient discussion wrt the state of the art, e.g. for online EM;
* the description of the approach seems yet not mature (e.g., the constraint enforcement boils down to considering that the $\pi_k$ are obtained using softmax; the discussion about the diagonal covariance matrix vs the use of local principal directions is not crystal clear);
* the fact that experiments need be strengthened.

I thus encourage the authors to rewrite and polish the paper, simplifying the description of the approach and better positioning it w.r.t. the state of the art (in particular, mentioning the data streaming motivation from the start). Also, more evidence, and a more thorough analysis thereof, must be provided to back up the approach and understand its limitations.
Decision: Reject