Fitting large mixture models using stochastic component selection
Authors: Milan Papez
Tomáš Pevný
Vaclav Smidl
Venue: NeurIPS 2021 Submitted (2021)
Keywords: Mixture models
expectation-maximization
sum-product networks
Metropolis-Hastings
unsupervised learning

Abstract: Traditional methods for unsupervised learning of finite mixture models require evaluating the likelihood of all components of the mixture. This becomes computationally prohibitive when the number of components is large, as it is in the sum-product (transform) networks. Therefore, we propose an approach combining the expectation-maximization and the Metropolis-Hastings algorithm to evaluate only a lower number of, stochastically sampled, components, thus substantially reducing the computational cost. We put emphasis on the generality of our method, equipping it with the ability to train both shallow and deep mixture models which involve complex, and possibly nonlinear, transformations. The performance of our method is illustrated in a variety of synthetic and real-data contexts, considering deep models, such as mixtures of normalizing flows and sum-product (transform) networks.

Fitting large mixture models using stochastic component selection
Anonymous Author(s)
Affiliation Address email
Abstract
Traditional methods for unsupervised learning of finite mixture models require to evaluate the likelihood of all components of the mixture. This becomes computa- tionally prohibitive when the number of components is large, as it is, for example, in the sum-product (transform) networks. As remedy, we propose an approach combining the expectation maximization and the evaluate only small number of, stochastically sampled, components, thus substantially reducing the computational cost. We put emphasis on generality of our method, equipping with the ability train both shallow and deep mixture models which involve complex, and possibly nonlinear, transformations. The performance of our method is illustrated in variety of synthetic and real-data contexts, considering deep models, such as mixtures of normalizing flows and sum-product (transform) networks.
algorithm
10
11
12
Introduction
14 Finite mixture models [40] constitute a fundamental class of density estimation models. They have been successfully applied in diverse fields, including bioinformatics [49], econometrics [1]] engineering [33], etc. A mixture model relies on a weighted sum of probability distributions- here referred to components-to cluster N unlabelled datapoints into K categories. The traditional maximum likelihood techniques train the model by optimizing either (i) the marginal likelihood via
16
18
gradient-descent (50) or (ii) the evidence lower bound via variational methods [4), including the
20
expectation expectation-maxinee maximization (EM) [S]. The dependence structure among approximate, variational,
21
distributions then ranges from the fully independent (mean-field) [25 to fully dependent [30] The
22
sampling-basee techniques target the posterior distribution using sequential Monte Carlo or
Markov chain Monte Carlo [5]. e.g. via the Gibbs 34 Metropolis-Faatingg sampling 38 The computational cost of these methods typically scales with OTK operations, where N and K
24
25
are defined above, T' is the number of iterations and D the dimension of data.
26 Various methods to decrease the computational cost via any factor OTTK have been proposed. T can be lowered by proper initialization, e.g. the optimal seeding [5); an efficient step-size schedule, 28 e.g. the line-search [58]; or increased estimation precision, e.g. the variance reduction [8. is often reduced using the coreset methods, which approximate the original dataset by weighted dataset
30 such that the exact and approximate marginal likelihoods are close. The weighted variants of the
variational 17 59 6 and sampling-based [39 methods then process the coresets. Reducing D relies
32 on the compression of data into smaller representations via random projections (53, 2], which is
achieved in ways: (i) each data item is projected into individual representation [1] (ii) all
34
data items are projected into an overall representation, commonly referred to as sketch 28 22
Nevertheless, all the aforementioned techniques-including those with reduced computational cost-
36 evaluate all K components. This is very demanding for large models, and the problem is even more
Submitted to Conference on Neural Information Processing Systems (NeurIPS 2021). Do not distribute.

---

Table 1: The computational features of various EM algorithms. We compare whether the methods (i) perform the computations with reduced number of data (minibatching), (ii) update a lower number of statistics, (iii) make less evaluations of the conditional likelihood, and (iv) are suitable for training of deep models. Here, EM, SA, S, MC and MH stand for expectation expecattin-mmxxiiizztion stochastic approximation, sparse, truncated, Monte Carlo and Metropolis-Hasstnns,, respectively.
Feature / Algorithm
EM
SAEM
SSAEM
TSAEM
MCSAEM
MHSAEM
13
44
24
18
11
(ours)
N datapoints
X
M K statistics
X
M < K likelihoods
X
X
X
X
deep models
X
X
X
X
X
37
severe for mixtures involving intricate models, such as neural networks [1 42 Gaussian processes
38 [57]. normalizing flows [48]; deep mixtures, including sum-product (transform) networks 45 47), deep Gaussian mixture models [55], etc. In spite of this, a little attention has been paid to the design 40 of algorithms which does not evaluate all K components. The notable exceptions are the sparse EM algorithm (24) and the truncated variational EM algorithm (18), see Table and Section for details. Moreover, the methods are mostly tailored for specific class of mixture models, e.g. the Gaussian
39
41
42
43
mixture models.
In this paper, we make the following contributions:
We propose EM-based algorithm which relies the MH sampler stochastically evaluate less components in mixture models, substantially reducing the computational cost.
48
47
We design our method to enable optimization of fairly generic EM objective functions, making suitable for training of both shallow and deep mixture models.
48
We apply approach to Gaussian mixture mdoels (GMMs) and their generalizations: product-transform networks (SPTNs) and mixtures of real-valued non-volume preserving (real NVP) flows [15], reaching approximately 100x speed-up compared state-of-the-art methods.
50
51
Problem formulation
52
A finite mixture model characterizes the relation between an observed (known) variable, & E X C R,,, 54 and latent (unknown) variable, 2 E = {l...., .....,KK K}. via the marginal (incomplete-datt likelihood
58
55 in the following form:
po Σ Pmo (x/2 k)) k),
56 where 8 := (Ti, ................................ TK nk) E θ are unknown parameters. Here, 7) are the parameters of the 57 conditional likelihood, Pu. (z)2). and Ty is the weight which parameterizes the prior, Px, (2) π2.
58 and satisfies F for each & E Z and ΣK mi 1.
Given a set of independent and identically distributed data, = (z.) goal is to learn the unknown parameters the marginal log-likelihood,
60
((0)= = log
log
Σ Pme (1g/22 = k)) (2i k). 1=
(2)
61 The marginalization in (2) is tractable for almost all forms of Pr. (z|2). Indeed, we consider Pa. (2|2) 82 belong to an arbitrary family of r). -i.ddifferentiable probability distributions. However, we assume that K is high, making the marginalization in 2 computationally costly, thus rendering the optimization
64 objective presumably intractable. Therefore, we want to design computationally efficient algorithm, requiring only M < K evaluations of p., (z|2) ()) at each iteration.

---

66
3
The EM algorithm
67
The maximum likelihood estimation seeks the parameters maximizing the marginal log- likelihoooood 68 GML = arg maxage L(0). The traditional EM algorithm [3] addresses this task indirectly, i.e. by
69 optimizing the evidence lower bound (ELBO),
(())
2(0) H(0) = ELBO(
(3)
where H((0) == (x) (z/x)] is the differential entropy an estimate, 0 O, and
70
(0) ppg(x)x) (oggoo
po(z:
k. xzi)
(4)
71 the EM objective function. Here, (z( x) is the joint (complete-data) likelihood, and pa (2|x)) is the posterior distribution over the latent variables == (z)) Given initial value, lo the EM
algorithm produces sequence of estimates, (6, (8,)z=1. by alternating between the expectation (E) and maximization (M) steps,
74
E-step: 2:-1(0);
(5)
M-step: 0, == arg max 16)
(6)
BE
This sequence is guaranteed to monotonically tighten the ELBO, arriving local optimum of under mild regularity assumptions [56].
76
77
The EM algorithm is computationally expensive, since (4) evaluates po(z;, xi)) for each E Z and
78 i E (I...., .......... .N.. This has be performed for all E (I....; T) in (5). Albeit the marginal factor,
79 Pn. (2), is just the cheap categorical distribution, the conditional factor, Pn: (21[2). typically involves
80 high-ding nigh-dimensional operations (e.g.. the inversion the full D D-dimensional covariance matrices
81 in the GMMs). Moreover, the M-step (6) is also expensive for large K. This holds despite that (6)
82 can be reduced to closed-form updates of expected sufficient statistics for Pq. (x|=) belonging to the
88 exponential family [4]] (again, due to high D). All in all, the computational complexity of the EM
84
algorithm scales with (TDNK).
85
If cannot be computed under a closed-form solution, one can resort to direct gradient-descent optimization of Q(0). where arg max replaced by one (or more) step(s) of a gradient descent technique. The EM algorithm is then referred to as the generalized EM algorithm [56).
86
87
4 The generalized MHSAEM algorithm
89
We design a version of the generalized EM algorithm suitable for scenarios where 4 can represent deep, discrete, latent variable models, thus being parameterized by possibly complex nonlinear transformations. We particularly focus on decreasing the the number of operations in the generalized EM algorithm from O(TDNK) OTDDBM, where B N and M < K.
81
88
4.1
E-step
34
We reduce the cost of evaluating the EM objective function (4) by combining the minibatching (as used many times before) and the Monte Carlo sampling. Namely, the specifie application of the latter to generic mixture models is the key contribution of this paper.
85
96
Minibatching. At each iteration, t, compute the conditional expectation in 4 only for subset- here referred to as minibatch-oo the original full dataset, i.c. (xr))ez. Here, / is a set of B < N indices, i, sampled uniformly without replacement from (1, .I..... N). This substantially decreases the necessary computations compared to the full sweep over all datapoints [23]
96
100
101
Monte Carlo sampling. For each i I, we want to draw M < K random samples from pe (z)) in order to obtain a Monte Carlo estimate of (4). The straightforward way to do this would to draw the samples directly from Pa 2; (z).. However, direct sampling from po (r.;) does not lead to any substantial decrease in the number of operations. This is caused by the fact that even for a single
102
104
sample of 24, we have to first compute the normalizing factor, pe (xx),, to obtaining the posterior,
3

---

106
po;z)la;). This requires K expensive evaluations of pa (2(2,.20 24), which is precisely what we want to avoid. Our approach is resort to the Markov chain Monte Carlo (MCMC), which allows us to sample from po(2i)o;), with the computational complexity decreasing to only single evaluation of po(2,, ..) per single sample of 2;-
107
108
109
110
MCMC methods obviate the computation of the normalizing factor in p=(2,(24) by simulating
Markov chain, (=s,))==1, from a transition kernel, ži, -(.), which leaves po(z; (a;) as its unique stationary (invariant) distribution, starting from an initial value 24.0- The specific form of P determines the structure of an MCMC method. We chose the Metropolis-Hastinng (MH) sampler.
111
112
113
114
which represents p(2,,e-1, 25,e) as follows: given E = 2,, draw sample from the proposal
115
distribution 2, ~ q(- (t), compute the acceptance ratio,
a 2.) min min Pa=,.t-1 (2. 19
Pos,,--1 (2s. -q
and, if u (22 23) where is drawn from uniform distribution, Uniform 10,1)-accepp accept the 117 sample and set 2i,t = z; otherwise, set 2i,t = 2, For each i E I and t E (I,...,7). we repeat Therefore, at every current iteration, we continue extend the chain from the point where we left at the previous iteration, 1, by 120 taking & 2M - Under mild regularity assumptions (52), the chain passes the transition period
this process M times, construing set Z,,t = - 2
118
(the burn-in phase), and the samples can then be used to approximate the conditional expectation in
122
(4 as follows:
1 Σ Σ
2
log Pr. (2:(2)
(8)
M
iET 2EZ...
123
Note that, to ensure this approach is truly efficient, we have to draw only M < K samples at each iteration, t; otherwise, for M = K, we may rather compute the exact marginalizatton in (4), since tractable (but computationally costly).
124
128
4.2
M-step
127 Assume for moment that (6) with 16-1(6) given by (8) has a closed-form solution, yielding an estimate of B. Such an estimate would have a high variance, converging only for M + oo and The main reason that the samples would not be reused over the iterations, t. thus wasting computational resources. We consider that there is no closed-form solution of (6). and to ensure that the samples (and thus computations) are recycled over the iterations we use the stochastic approximation (SA) [5] to optimize (8). This is analogous to applying a stochastic gradient-deseee method, 8, 0 (1-)) where "/ is the step-size, satisfying the Robbins- Monro constraints, "/ E [0, 1]. Ex21 Yt = co, Eiz1 72 < oo, and V is the gradient w.r.t. O. In this way, the computations made in Q are accumulated via Og and reused over the iterations.
128
129
T + 00
130
131
132
134
135
136
The parameters Y have a different form based on a specific case of Pa, (x/=), whereas Tu is a permanent structure in (I). Therefore, without loss of generality, we split (6) into generic part and
137
138
fixed part as follows:
Th.t = Mk, + 7/6 V. (.))(f).
(9a)
M/, Mktt- 71t V., (.),
(9b)
139
where to ensure that the probabilities, (TT,t) satisfy the constraints (Section 2 we transform
140 V V. 2 via P4= log F and optimize w.r.t. V/. Then, to obtain (T), k from 14 (vk:)) we
use the softmax function, i.e. Th,t == (oftimax(P:):: = exp(ux.)) Σ exp(1I.))
141
142
Computing the gradients for all pairs of (vk, m) k would be inefficient, especially since Zi,t contains
only a small number of unique values of Z for M < K. Consequently, we compute V, Q and only for unique(z..)). We summarize the proposed approach in Algorithm
143
2
144
4.3
Proposal distribution
146
The choice of the proposal distribution has a significant impact on the speed of convergence and the computational cost of the proposed algorithm. Here, we discuss various possible choices of (<s.),.)

---

Algorithm The generalized MHSAEM algorithm
Input: obs (2i,0)=11 (x
Output: (0.)/=1
f (L....,
until convergence do
form the set (i)) by sampling (without replacement) B indices i ~ ,V)
for i E I do
set 2; as the last element of 2i,-11
for j e (1, (1..... M) do
sample 2; ~ (zz 24) sample u Uniform compute (z(z,, in
if u < aa(2,, 2,) then
set 2/, = 2, and E 23
else
set 2/,4 2
end if end for
set Zi,t
(z 2,,
end for
compute
compute 9 for (nique(zs.))
compute Th,s = softmax(v:)) for k E Z
end for
148
Optimal proposal (O). The optimal proposal distribution is q(zi(z,) = g(z)) (z((a)) This
ensures that the acceptance rate ( is always a(2,, z;) = 1. However, the need to perform K expensive evaluations of (= x4) before sampling from po (2.|1x)) is the reason we resorted to the MH sampler in the first place. We consider this only to set the upper limit admissible computational cost and to study the impact of sub-optimal proposal distribtions.
150
152
153
Uniform proposal (U). The uniform distribution on the discrete interval from 1 to K, i.c. (<((²)= =
154
((x:) = Uniform(1, K), the simplest and computationally cheapest variant of the proposal distribution. However, due to poor mixing properties, the algorithm may converge slowly for high K.
155
Tabular proposal with forgetring (TF). The key requirement to design proposal distribution is to restrict its computational complexity somewhere between that of the and O proposals. One way to satisfy this constraint
157
use the Markov chain, (zi,) to learn a transition kernel, p(z. [2), see,
158
c.g. B3) Unfortunately, this would require to store table with K entries for each (L...., N). which is very demanding even for moderate K and N. Therefore, we break the dependence in 1()==k) is the the Markov chain and define: g==//=;= = (z)) == C(as), where C(a:) II k= k
160
categorical distribution with the weights a; == (a)....... a. For C(a,) > log qa, (Zi,)) we obtain an estimate of a; at iteration as follows: Og,t = arg L(x,)
162
dad with
163
184
nhe e where is the standard basis vector (a one-hot vector) with one at kth position and zeros otherwise. This can be further rewritten into recursive form: T2s 14,1-1 + e.,. or. using the Robbins-Monee step-size, ni,t (1 C2, O ni,t11 + 7C where 1 is the vector of ones, and @ is the Hadamard product. We refer to this case simply as "table with forgetting" (TF)
185
167
168
due that it represents N K table in the memory and 7t is forgetting factor.
5
Related work
189
170
Stochastic approximation pnnectatioonmaxiiizzatiiin The application of SA to prevent the evaluation of all components in mixture models has been overlooked for a long time. The reason is that the original motivation to combine the EM algorithm with SA is to address the analytical intractability of the expected value under (=]xxx) in 4 which is, however, almost always tractable for mixture models. The intractability issue is addressed by either the Monte Carlo SAEM (MCSAEM) [12
171
172
175
175 or the Markov chain Monte Carlo SAEM (MCMCSAEM) BT. Applying the former approach to
mixture models would be inefficient, since it evaluates K joint distributions, pe(2, æ). before drawing
5

---

exacT
10.0
EM
SAEM
7.5
SSAEM
((0))))
TSAEM
5.0
MCSAEM
MHSAEM-U
2.5
MHSAEM-TF
0.0
500
1000
1500
computational time (sec)
The training log- liliikiiihodd, [(),), versus the computational time (in seconds). Here, on the x-axis, the computational time at current iteration. f, is obtained by accumulating the time from the previous iterations. corresponds to [(0rec) ). where 195 is the iteration of reaching 95% of max [(O)). The projection of o on the x-axis gives the time to reach L(ftes). ) This experiment was performed with the following settings: (D, K, N, w, B, M) (10, 100, 10k, 0.1, 200, 2, 20k),, see Section 6.1 for details. The results are averaged over five repetitions.
Figure 1:
177
M samples from Po (z)r)) Therefore, this method reduces only the computational cost of updating the sufficient statistics. This is addressed by the latter approach, where M K samples from a proposal distribution, (()[2)) is used to calculate p0(2, x) and also the sufficient statistics. However, all these
178
179
180
methods process all data at every iteration, providing only a limited advantage over the conventional
EM algorithm. Minibatch versions of these techniques have recently been proposed [27 32 1.
All the above methods commonly assume (0(2, x)) belonging the exponential family. This provides a convenient, but limiting, property which allows (6) to be computed under a closed-form solution. The main contribution of our work release this restrictive assumption by admitting that pe (2, (and thus Q) is given by possibly complex and intractable transformations.
189
184
185
Sparse and truncated variational techniques. There is only small body of methods explicitly reducing the number of evaluated components. Their common aspect is that they follow from the variational framework, where the exact posterior, pa (2( is approximated by variational posterior, (=(z)) This sparse, approximate, posterior is defined over a lower number of components, M < K, such that only the important components are selected, relying relaxation of the hard EM algorithm from taking single M = assignment [26) to taking multiple M < K assignments. The sparse SAEM (SSAEM) algorithm 24 selects the components by quick partial sorting of the posterior probabilities, pa (z|æ). Again, this requires K evaluations of pe (z, x) before the sorting, thus only reducing the amount of updated statistics. Similarly, the truncated SAEM (TSAEM) algorithm (18) selects M < K cluster-to-cluster and M < K cluster-to-datapoin minimal Euclidean distances, preventing the problem in the SSAEM algorithm. However, all these distances are evaluated for all components in pairwise manner, leading to K computaaaii complexity, which makes the saving dubious. Similarly before, these methods assume po(z, x) belong the exponential family.
187
189
190
191
193
194
185
196
197
198
199
We summarize the distinguishing features of the above discussed methods in Table
Experiments
200
To demonstrate the key features of our algorithm its low computational complexity, competitive learning performance, and generality- we use it below to train: (i) GMMs on synthetic datasets, and (ii) SPTNs (47) and (iii) mixtures of real NVP flows [48] on real datasets. All experiments have been performed on Slurm cluster equipped with Intel Xeon Scalable Gold 6146 with 384GB of RAM.
202
204
205
6.1
Gaussian mixture models
206
Consider the special case of a data- generating distribution given by (T). with the components taking the form of the multivariate Gaussian distribution, Pr, (x|=) = N(x; Ma Σz), where μı= is the mean value and Σg is the covariance matrix. The difficulty of learning GMMs heavily depends on the degree of interaction among all mixture components, hence having the ability to generate synthetic
207
209
datasets with arbitrary overlap characteristies between all pairs of components is crucial for systematic

---

EM
SAEM
SSAEM
TSAEM
MCSAFM
MHSAAMM-L (ours)
MHSAFM TF (Ours)
8
K=25
B=100
12
K=50
0 B=200
10
0
6
O K=75
O R-:00
O M
A K=100
A B=400
4 M
AE
AET
AB
D0--A
2
0
10a computational time (400)
10ah
10 computational time (sec)
10a
10 computational time (sec)
10a
Figure 2: The absolute error, AE = [C((to)) [(0)), versus the computational time (in seconds). All experiments use the following settings: where the number of components, K, (left), the batchsize, B, (middle) and the number of samples, M. (right) change for different values denoted by (t, q o, 4). At each of these points (marks), we perform an experiment as illustrated in Figure [ find L(Btan) to compute the AE, and record the time corresponding to tys- The results are averaged over five repetitions.
N,w B, M, T) (10, 100, 10k, 0.1, 200, 2, 20k).
211
evaluation of performance of learning algorithms 443 Traditional techniques usually define overlap (or separation) of components only in terms of their mean vectors and maximum cigenvalues of the covariance matrices, not accounting for their rotation and mixing weights (see [6) for a detailed treatment of the problem). We therefore a more objective measure of the clustering complexity
213
215
defined by the total probability of misclassification [41] which allows to generate data with a user-defined degree of maximum pairwise overlap, w.
216
217
Experiment settings: We generate the parameters of [ and the corresponding dataset, uniquely for a given quadruple (D, K, N, W.). Therefore, the parameters of the generative model are known and we can measure and display the convergence of the training liglliielhood, [(O)). compared to the exact 219
log- likeliiiihodd, [()), for t (I....,T) .)) We are further interested in the absolute error between the training log-likelihood at the iteration of reaching 95% of its maximum value, t95, and the exact log- likelihood, i.e. AE [((...)) L(O)
221
222
223
We also measure the computational time until reaching tys- We have used 95% of the maximum value instead of the maximum value to prevent cases, where the model oscillate around target value, making the estimate of convergence time very noisy (for example ACSAEM Figure IT.
224
225
226 Algorithms:
The GMMs belong to the exponential family of probability distributions. This allows us to find closed-form, recursive, solution of (6. relying on a Robbins-Monro type of the step-size sequence. (-o)==- 7 44). In this setting, we compare our MHSAEM algorithm with number of related methods Table Note we use the acronyms U and TF to specifyy the proposal distribution of the MHSAEM algorithm (Section 4.3) However, we do not use the O-proposal, since the MHS AEM- algorithm is equivalent to the MCSAEM algorithm. All the SA-variants Table T Jse minibatch of size B. The key quantity to reduce the number of evaluated components and/or sufficient statistics
227
228
229
230
231
232
233
in the SSAEM, TSAEM, MCSAEM and MHSAEM algorithms is collectively denoted by M (Section
5) Note that we always keep M = M in the TSAEM algorithm (see Figure and 2 for concrete 235
numbers). We use the step-size given by 7 = for L = 1, 1,..., 50 and 7t 0.05 otherwise. In this section, to counteract the issue of attaining poor local optima, we equip all algorithms with the anti-anncaling schedule (B.) starting with 3 = 0.1, reaching 82/37 = 1.2, and decreasing back to 37 1.0, see 43) for details. The initial estimates of: (i) μ4 are uniformly drawn from the unit hyper-cube, (ii) Σ are fixed to unit diagonal matrix, and (iii) Ty are uniformly drawn from the unit interval (followed by normalization.
237
239
241
Figure I shows that the EM and SAEM algorithms take the longest time to converge, attaining a poor local optima. On the other hand, the MCSAEM [ and MHSAEM (U and TF) algorithms achieve L(Orga) closest to the likelihood L(O) of the true model. Moreover, both MHSAEM algorithms reach this value in the shortest time compared to all the other methods. The SSAEM 24 and TSAEM [18] algorithms are comparable in terms of the computational time, but we investigate sensitivity of fitting the model to
Results:
242
243
245
they both provide the lowest L(0rgs). In Figure 2
7

---

247
increasing values of K, B and M by measuring the time and the likelihood again. In all the cases, the proposed MHSAEM algorithms achieve the lowest AE in the shortest time.
248
249
SSAEM and TSAEM algorithms failed to converge for M >> 2 and for K > 50 respectively. We
250
believe this is caused by selecting only M maximal probabilities in the SSAEM (or distances in the
TSAEM) algorithm (Section[5) 5 which prevents certain, but not a negligible number of, components
251
252
from being updated, thus providing only crude approximation of (2(x)) The results then suffer
from substantial variational gap to the exact log-likelihood (Figure D On the contrary, MH sampler
253
254
provides samples which consistently approximate po (2(2z) despite evaluating much lower number of components in each step.
255
256
6.2
Sum-produc transform networks
257
The sum product networks (SPNs) are a deep learning extension of finite mixture models. They can be interpreted as a mixture of trees [60], where each tree corresponds to a component. Therefore, they can be cast into the form of ( but the number of components grows exponentially with their
258
259
260
depth. In this section, we use recently proposed SPTNs which introduce additional transformation
261
nodes to provide better expressiveness than the SPNs (SPTNs effectively generalize SPNs and flow models into one large family of models).
263
Experimental settings: We use 19 real datasets from the UCI database [16 37
35
54 preprocessed
264
in the same way as in [46]. For each experiment, we randomly split the data into 64%, 16% and 20% for training, validation and testing. respectively. We calculate the average ligeeiilldood on the test set and measure again the time to reach 95% of the maximal training log-likelihood L(O..))
285
286
To evaluate various (possibly shallow and/or deep) architectures of SPTNs, we fit each dataset with all the following combinations of "uper-parameterr]] 8 E (8, 32, 128), E (2, 4, 6, 8), E (2,3, 3.4 4), where s is the number of children of each sum node, b is the number of partitions of each product node, and / is the number of layers (one layer contains sum and product nodes). The number of components of the SPTN, after its conversion into (I. is given as follows: K s Note that the maximum number of components for the investigated parameters of the SPTN is 466,443,456. To restrict ourselves only to (i) the leaf nodes given by (0, I); (ii) affine transformations fixed to the singular value decomposition, choosing the the Givens parameterization for the unitary matrices [47] and (iii) no sharing of any type of nodes (47.
268
270
272
273
reduce the space of possible architectures,
274
275
276
Algorithms: We evaluate only on the MHSAEM-U algorithm-due to its favourable computational complexity and simplicity- and compare it with the stochastic gradient-descent (SGD) algorithm, which is routinely used to train SP(T) [45] (7). In this case, SGD in each iteration performs computations over all subtrees of the network, whereas the MHSAEM-U algorithm computes with only M = subtrees, thus we should observe speed-up of the computations. In our implementation. both these methods perform optimization of their respective objective functions-thh the log-likelihood 2 for SGD and the EM objective (8) for MHSAEM-U via the use of the automatic differentiation
278
280
282
283
and the ADAM optimizer [29], using B 100 and T' 20000.
284
Results: Since each dataset might benefit from a different architecture, Table 6.2 shows the test log- likelihooo of the architectures selected according to the best likelihood measured on the validation set and the corresponding speed-up. The test log likeeihoosss reveal that the MHSAEM-U algorithm outperforms the SGD algorithm 10 out of 19 datasets, which the added stochasticity helps to escape poor local minima. The speed-up demonstrates lower computational complexity of the MHSAEM-U algorithm on 17 out of 19 datasets, which was the main goal. The magic-telescope and wine datasets show approximately 102x and 75x speed-up, respectively, while on very small datasets (pima-indian and iris), the SGD is faster due to effective implementation. In the supplementary material, we present Table 3, exhibiting the same trends on a fixed architecture.
286
not originally the goal, but
288
289
290
292
293
6.3
Mixtures of real NVP flows
285
We consider another class of mixture models II. where each component Pr, (x|²) is transformed by the flow meael-real NVP [15.. These transformations are parameterized via deep neural networks, allowing for flexible adjustment of the learning capacity of each component.
296
297
We have set hard limit train a single model to 24h, which is default Slurm cluster.
8

---

The speed-up and test log-likelihood, L for the SGD and MHSAEM-U algorithms. The test log-likelihood (higher is better) is computed for the best model, with the corresponding K, which is selected based the validation log-likelihood. The speed-up is computed as the ratio of MHSAEM-U to SGD, i.e. their time to reach 95% of the training bog-likelihood. The results are averaged over five repetitions. Then, the higher test log-likelihodd highlighted with bold blue, and and no speed-up is highlighted with red. The average rank is computed as the standard competition on each dataset (lower is better).
Table 2:
(*1224"") ranking
[1
Sum product fransform nelworks
NVP llws
SGD
MHSAEM-U
SGD
MHHAEMMU
dataset
speed-up
glea
K
K
speed-up
plo
K
breasf-cancer wisconsinnn
4.66
64
1.43
1024
0.63
128
candiotocugraply
10.55
59.52
512
31.04
1024
9.85
54.34
32
56.08
128
magic telescope
102.53
pendigits
4.89
0.88
1024
-486
16384
4.17
1.46
0.48
pina-inclians
0.37
-8.54
64
-7.62
64
1.35
-20.09
128
-16.33
128
wall- vall-following-ae wireiiinggooooooo robot
3.43
1.44
1024
-11.3
16384
22.21
-14.26
128
-17.56
128
waveform-1
4,35
-26.14
64
22,91
1024
3.72
-34.12
-33.42
4.82
-26.21
4.12
yeast
20.57
10.26
512
5,18
1024
14.49
6.61
128
9.59
128
ecoli
-5.5
2.15
-1..77
ionosphere
1.88
-20.27
64
-5,93
512
2.74
-87.01
128
-42.75
128
itis
0.23
-10.65
64
-1.49
16384
3,28
-16.34
128
-9.21
32
page blocks
12.18
12.21
512
6.84
1024
44,95
17.13
128
17,94
32
parkinsons
1.46
-2188
64
0.5
512
3.09
-566.58
128
-33,31
32
sonar
-95.39
2.52
-622.2
statlog segment
1.44
47.35
512
26.53
16384
38.49
35.84
128
42.04
32
statlog-veht vehicle
2,97
6.78
-31.34
-26.43
wine
75.42
-25.99
1024
-13,27
1024
2,05
-171.58
128
-25.57
128
rank
1.56
1.44
1.83
1.17
298
Experimental settings: We use the same experimental settings and evaluation metrics as in Section
289
We apply the mixture model on all datasets, changing the number of components as follows: (8, 32, 128). Each real NVP-based component in the mixture model has (i) the translation function parameterized via multi-layer perceptron with single hidden layer of dimension 10, using the rectified linear activation function; and (ii) the scale function parameterized via the same network except with the hyperbolic tangent activation function. We do not use the batch normalization [15] and
6.2
300
302
303
304
we stack two layers the translation-see scale transformation (we have used implementation from 20
305
Al Algorithms: The algorithms and their settings are the same as those in Section 6.2
Results:
The experimental results are presented in right part of Table 6.2 They are similar to those obtained in the previous section. terms of the test log-likelihood, the MHSAEM-U algorithm outperforms the SGD algorithm on all but three datasets, and it provides substantial speed-up on all datasets except one. The test likelihood of models with the real NVP flows is most the time worse than that of SPTNs with the affine transformations. As explained in the supplementary, this is due to the overfitting, which has been observed in [47].
307
309
311
Conclusion
312
This paper has presented a method to decrease computational complexity of fitting mixture models, including their generalizations, such as sum-product- (transforr)))) networks and mixtures flow models. The speed-up is achieved by evaluating and updating only a single component (per iteration), where the Hatropolis-Hasttin algorithm ensures sampling of components from a proper posterior. An experimental comparison all three classes of models mentioned above confirmed the theoretical expectations. The method significantly speeds-up the fitting time and, importantly, without sacrificing the quality of the fit. In fact, the likelihood was better than that of the models fitted by the EM algorithm or the SGD algorithm in more than 50% cases. We attribute this to higher stochasticity, which helps to escape from poor local minima.
314
316
318
320
322
the experiments, the proposed method has used a uniform proposal distribution the MH sampler. Despite outperforming the alternative methods, we conjecture that this limits the speed of convergence. Therefore, we believe that there still a room for improvement in the implementation. We plan to address these issues in future work.
324

---

328 8
Broader impact statement
327
The presented method decreases the computational complexity of fitting large (and deep) mixture models, which leads to five to hundred time speed-up depending size of the problem (although negative exceptions occurs). We believe this line of research, which we want to continue, to have important benefits. First, it is directly related to decrease in energy consumption and in production of
329
331
CO2 (we expect similar rates as the speedup). Second, it has a positive effect on financial aspects of deploying (and experimenting with) mixture models. Third, it decreases the hardware requirements, as in all experiments presented above the model was fitted on a single-core.
333
References
335
[1] S. Allassonnière and J. Chevallier. A new class of stochastic EM algorithms: Escaping local maxima and handling intractable sampling. Computational Statistics Data Analysis, 159:107155 159, 2021.
336
338
[2] S. Ayesha, M. K. Hanif, and R. Talib. Overview and comparative study of dimensionality reduction techniques for high dimensional data. Information Fusion, 59:44 58, 2020.
340
[3] D. S. Bai. Efficient estimation of transition probabilities a Markov chain. The Annals of Statistics, pages 1305-1317, 1975.
342
[4] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe. Variational inference: A review for statisticians. Journal of the American statistical Association, 112(518)859-877, 2017.
343
344
[5] J. Blömer and K. Bujna. Adaptive seeding for Gaussian mixture models. In Pacific-asia conference on knowledge discovery and data mining, pages 296-30969900 308. Springer, 2016.
345
[6] T. Campbell and B. Beronov. Sparse variational inference: Bayesian coresets from scratch. arXiv preprint arXiv: 190::099003399 2019.
347
77] O. Cappé and E. Moulines. On-line xxpecttionnmmxxiizaato algorithm for latent data models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), (13):593-613,, 2009.
349
351
[8] J. Chen, J. Zhu, Y. Teh, and T. Zhang. Stochastic expectation maximization with variance reduction. Advances in Neural Information Processing Systems, page 7967, 2018.
352
353
[9] N. Chopin. A sequential particle filter method for static models. Biometrika, 89(3)5539-552,
354
2002.
355
[10] G. Compiani and Y. Kitamura. Using mixtures in econometric models: A brief review and some new results. The Econometrics Journal, 1953) 2016.
356
[11] S. Dasgupta. Learning mixtures of Gaussians. In 40th Annual Symposium Foundations of Computer Science (Cat. No. 99CB37039), pages 634-644. IEEE, 1999.
358
[12] B. Delyon, M. Lavielle, E. Moulines, et al. Convergence of stochastic approximation version of the EM algorithm. The Annals of Statistics, -7(1):44-128, 1999.
360
[13] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39(1):1-22, 1977.
362
363
384
[14] J. Demšar. Statistical comparisons of classifiers over multiple data sets. The Journal of Machine Learning Research, 7:1-30, 2006.
385
366
[15] I.. Dinh, J. Sohl-Dickstein, and S. Bengio. Density estimation using real NVP. In 5th Interna- tional Conference on Learning Representations, ICLR 2017, 2017.
367
[16] D. Dua and C. Graff. UCI machine learning repository, 2017.
10

---

389 [17] D. Feldman, M. Faulkner, and A. Krause. Scalable training of mixture models via coresets. In Proceedings of the 24th International Conference on Neural Information Processing Systems, pages 2142-2150, 2011.
370
[18] D. Forster and J. Lücce.. Can clustering scale sublinearly with its clusters? A variational EM acceleration of GMMs and k-means. International Conference on Artificial Intelligence and Statistics, pages 132. PMLR, 2018
373
[19] G. Fort, E. Moulines, et al. Convergence of the Monte Carlo expectation maximization for curved exponential families. Annals of Statistics, 3240)1220-1259, 2003.
376
377
[20] J. Francū.
2020.
[21] K. Greff, S. van Steenkiste, and J. Schmidhuber. Neural expectation maximization. In Proceed- ings of the st International Conference on Neural Information Processing Systems, pages 6694-6704, 2017.
380
382
[22] R. Gribonval, A. Chatalic, N. Keriven, V. Schellekens, L. Jacques, and P. Schniter. Sketching datasets for large-scale learning (long version). arXiv preprint arXiv:2008.01833 2020.
383
384
[23] M. D. Hoffman, M. Blei, C. Wang, and J. Paisley, Stochastic variational inference. Journal of Machine Learning Research, 14(5), 2013.
385
386
[24] M. Hughes and E. B. Sudderth. Fast learning of clusters and topics via sparse posteriors. arXiv preprint 1609.0569. 2016.
387
388
[25] K. Humphreys and D. Titterington. Approximate Bayesian inference for simple mixtures, In COMPSTAT, pages 331-336, Springer, 2000.
389
390
[26] B.-H. Juang and L. Rabiner. The segmental K-means algorithm for estimating parameters of hidden Markov models. IEEE Transactions on acoustics, speech, and signal Processing, 3839)1163911441 1990.
392
[27] B. Karimi, M. Lavielle, and É. Moulines, On the convergence properties of the mini-batch EM and MCEM algorithms, 2019.
393
[28] N. Keriven, A. Bourrier, R. Gribonval, and P. Pérez. Sketching for large-scale learning of mixture models. Information and Inference: A Journal of the IMA, 7(3):47-508, 2018
396
397
[29] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. arXiv:1412.6980, 2014.
arXiv preprint
[30] A. Kucukelbir, D. Tran, R. Ranganath, A. Gelman, and D. M. Blei. Automatic differentiation variational inference. The Journal of Machine Learning Research, 18(1):430-47) 2017.
400
401
[31] E. Kuhn and M. Lavielle. Coupling a stochastic approximation version of EM with an MCMC procedure. ESAIM: Probability and Statistics, 8: 1:115-131, 2004.
[32] E. Kuhn, C. Matias, and T. Rebafka. Properties of the stochastic approximation EM algorithm with mini-batch sampling. Statistics and Computing, 3006): 2020.
404
405
[33] A. Lagrange, M. Fauvel, and M. Grizonnet. Large-scale feature selection with Gaussian mixture models for the classification of high dimensional remote sensing images. IEEE Transactions on Computational Imaging, 3(2):230-2)) 2422 2017.
406
407
408
[34] M. Lavine and M. West. A Bayesian method for classification and discrimination. Canadian Journal of Statistics, 20(4):451 461, 1992.
409
410
[35] M. Little, McSharry, S. Roberts, D. Costello, and I. Moroz. Exploiting nonlinear recurrence and fractal scaling properties for voice disorder detection. Nature Precedings, pages 2007.
[36] R. Maitra and V. Melnykov. Simulating data to study performance of finite mixture modeling and clustering algorithms. Journal of Computational and Graphical Statistics, 19(2):34-3376
413
2010.
11

---

415
[37] O. L. Mangasarian and W. H. Wolberg, Cancer diagnosis via linear programming. Technical report, University of Disconsin-Madison Department of Computer Sciences, 1990.
416
417
[38] J.-M. Marin, K. Mengersen, and C. P. Robert. Bayesian modelling and inference on mixtures of distributions. Handbook of statistics, 25:459-507, 2005.
[39] C. A. McGrory, D. C. Ahfock, J. A. Horsley, and C. L. Alston. Weighted Gibbs sampling for mixture modelling of massive datasets via coresets. Stat, 3(1):291 2014.
420
421
[40] G. J. Mel achlaan, S. X. Lee, and S. I. Rathnayake. Finite mixture models. Annual review of statistics and its application, 6:355-378, 2019.
[41] V. Melnykov, W.-C. Chen, and R. Maitra. MixSim: An R package for simulating data study performance of clustering algorithms. Journal of Statistical Software, 51(12):1, 2012.
424
425
[42] T. Monnier, T. Groueix, and M. Aubry. Deep transformation-invatte clustering, In Conference on Neural Information Processing Systems (NeurIPS 2020), 2020.
426
427
[43] L. Naim and D. Gildea. Convergence of the EM algorithm for Gaussian mixtures with unbalanced mixing coefficients. In Proceedings of the 29th International Coference on International Conference on Machine Learning, pages 1427-1431, 2012.
429
430
[44] H. D. Nguyen, F. Forbes, and G. J. McLachlan. Mini-batch learning of exponential family finite mixture models. Statistics and Computing, pages 1-18, 2020.
432
[45] R. Peharz, A. Vergari, K. Stelzner, A. Molina, X. Shao, M. Trapp, K. Kersting, and Z. Ghahra- mani. Random sum-product networks: A simple and effective approach to probabilistic deep learning. In Uncertainty in Artificial Intelligence, pages 334-344. PMLR, 2020.
433
434
435
[46] Pevný. Loda: Lightweight on-line detector of anomalies. Machine Learning, 102(2):275-304, 2016.
436
437
[47] T. Pevný, V. Smídl, M. Trapp, O. Poláček, and T. Oberhuber. Sum-product-transfioe networks:
438
Exploiting symmetries using invertible transformations. arXiv preprint arXiv:2005.01297, 2020.
[48] G. G. Pires and M. A. Figueiredo. Variational mixture of normalizing flows. arXiv preprint ar arXiv:2009.00585 2020.
441
442
[49]
A. Rau, C. Maugis M.-L.. Martin Martin-Magniette, and G. Celeux. Co-expression analysis of high-throughput transcriptome sequencing data with Poisson mixture models. Bioinformatics, 31 31(9):1420-144 1420001-122,2200 2015.
444
445
[50] R. A. Redner and H. F. Walker. Mixture densities, maximum likelihood and the EM algorithm. SIAM review, 26(2):195-25) 1984.
448
447
[51] H. Robbins and S. Monro. A stochastic approximation method. The annals of mathematical statistics, pages 400-407, 1951.
448
[52] C. Robert and G. Casella. Monte Carlo statistical methods. Springer Science & Business Media, 2013.
450
451
[53] W. Siblini, P. Kuntz, and F. Meyer. A review on dimensionality reduction for multi-label classification. IEEE Transactions on Knowledge and Data Engineering, 2019.
453
[54] J. P. Siebert. Vehicle recognition using rule based methods. 1987.
454
[55] C. Viroli and G. J. McLachlan. Deep gaussian mixture models. Statistics and Computing, 29(1):43-51, 2019.
455
456
[56] C. F. J. Wu. On the convergence properties of the EM algorithm. The Annals of statistics, pages
457
95 95-103, 1983.
458
[57] D. Wu and Ma. An effective EM algorithm for mixtures of Gaussian processes via the MCMC sampling and approximation. Neurocomputing, 331:366-374, 2019.
12

---

480
[58] W. Xiang, A. Karfoul, C. Yang. H. Shu, and R. L. B. Jeannès. An exact line search scheme to accelerate the EM algorithm: Application to Gaussian mixture models identification. Journal of computational science, 41:101073, 2020.
[59] M. Zhang, Y. Fu, K. M. Bennett, and T. Wu. Computational efficient variational Bayesian Gaus- sian mixture models via coreset. In 2016 International Conference on Computer, Information and Telecommunication Systems (CITS), pages 1-5. IEEE, 2016.
484
[60] H. Zhao, P. Poupart, and G. Gordon. A unified approach for learning the parameters of sum- product networks. In Proceedings of the 30th International Conference on Neural Information Processing Systems, pages 433-44- 2016.
467
13

---

Checklist
489
470
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's s contributions and scope? [Yes]
472
473
(b) Did you describe the limitations of your work? [Yes] Our main contribution is compu-
474
tational speedup. Cases where it was not achieved are highlighted in the experimental
475
section.
(c) Did you discuss any potential negative societal impacts of your work? [No] We do not foresee any potential negative impact.
477
478
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
479
480
2. If you are including theoretical results...
481
(a) Did you state the full set of assumptions of all theoretical results? [N/A]
(b) Did you include complete proofs of all theoretical results? [N/A]
483
3. If you ran experiments..
484
(a) Did you include the code, data, and instructions needed to reproduce the main exper- imental results (either in the supplemental material as a URL)? [Yes] The code is available in a github repository. All dataset are public from the UCI database.
486
487
(b) Did you specify all the training details (c.g., data splits, hyperparameters, how they were chosen)? [Yes] See Section
488
489
(c) Did you report error bars (e.g.. with respect to the random seed after running experi- ments multiple times)? [No] We report only average of Monte Carlo repetitions, the error bars were too small to have any visual impact in the reported logarithmic scale.
491
492
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes]
483
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
485
(a) If your work uses existing assets, did you cite the creators? [Yes] We use 20 datasets from UCI, we cite the required papers for each dataset, mostly the UCI database and few additional publications.
496
497
(b) Did you mention the license of the assets? [No] The data are publically available, comply with the requirement on citing appropriate publications.
499
500
(c) Did you include any new assets either in the supplemental material or as URL? [No]
501
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A]
502
503
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing conducted research with human subjects...
506
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
507
508
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
511
14


Reviews:
[{"review_id": "S8ARXopR14", "reviewer": "Reviewer_n8QC", "summary": "This paper proposes a stochastic EM algorithm. Specifically, the authors use a minibatch to approximate the data and use samples obtained from Metropolis-Hasting (MH) to approximate the likelihood component per datum. The proposed method is tested on a synthetic dataset and several UCI datasets.", "questions": "", "limitations": "", "rating": 4, "confidence": 3, "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": "The main contribution of this paper is the proposed likelihood component approximation using MH, as using minibatch to approximate the dataset is pretty common in ML. The empirical results of the proposed method seem promising, but I have some concerns about the methodology.\n\nIt is known that MCMC needs a burn-in phase in order to produce samples from the correct distribution. The proposed algorithm only runs M steps where M needs to be small enough to save computation. It is hard to believe the samples are from the right distribution. Moreover, the algorithm collects consecutive samples thus the samples will be highly correlated. What is the acceptance probability of MH? It seems that these will introduce bias to the approximation. Besides introducing bias, how is the variance of this approximation?\n\nIn summary, the idea is reasonable and there seems empirical performance. But the methodology is not very sound.  \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The paper does not contain this part.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "FFPprYdRALk", "reviewer": "Reviewer_XENJ", "summary": "The authors propose a sampling-based method to improve the efficiency of the E-M algorithm on mixture models. The method relies on performing the computations on a reduced number of MCMC-sampled components, rather than all of them.", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "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": "The paper is well-structured, and I did not have a problem following the exposition.\n\nThe main idea of the paper is straightforward: instead of performing computations over all mixture model components, perform an approximate computation on a smaller number of components sampled via MCMC to reduce the total complexity of E-M. While reasonable, this idea by itself does not constitute a significant contribution, and there is no further theoretical analysis or discussion. The \"TF\" proposal distribution in section 4.3 could potentially be a novel contribution, but it is only discussed very briefly (see also following comments).\n\nThe experimental section, and specifically the two non-synthetic data experiments of sections 6.2 and 6.3, seem to have some major flaws:\n- The authors do not even evaluate the previously discussed \"TF\" proposal in these experiments.\n- The experiments seem to have been run with $M = 1$, which means just choosing a single component at random. Whether this choice is uniform or not depends on whether there is a burn-in period or not, which is also not discussed in the paper. I do not see the point in discussing an MCMC method if only a single sample is used in the end.\n- I do not understand why there are such big differences in the test likelihood values in Table 2 between the two compared methods. In particular, why would the \"SGD\" method ever result in so much worse likelihoods? Is it some implicit regularization effect? Is the variance of the likelihood over different repetitions high? Why do the authors not show standard deviations, since their results are over five repetitions?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "6zXrhukATmf", "reviewer": "Reviewer_nhhz", "summary": "This paper proposes a stochastic approximation to the E-step of the EM algorithm to reduce the computational cost of fitting mixture models when the number of components $K$ is large. Specifically, the paper proposes using $M << K$ samples of the latent variables $z$ from a Metropolis-Hasting (MH) sampler to approximate the E-step. This lets us only use $M$ likelihood evaluations instead of $K$ for each observation in the E-step. The paper compares there proposed method against other EM-based methods on Gaussian mixture models on synthetic data, and sum product  transform networks and non-volume preserving flows on 19 UCI datasets. ", "questions": "", "limitations": "", "rating": 4, "confidence": 3, "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": "This paper should be rejected.\n\nOriginality: the combination of MH-sampler and EM appears new, specifically only sampling $M << K$ points to avoid evaluating the likelihood $K$ times for each point. The paper focus primarily on EM methods, but does not discuss comparisons with variational inference (both EM and VI are variational methods that maximize the ELBO).\n\nQuality: the paper is mostly correct. I have a few small concerns: \nFor Section 4.1, initializing MH sampler $z_i$ from the last sampled $z_i$ (to avoid using a burnin) seems reasonable, but for large N relative to minibatch size B, isn't it possible the last sampled $z_i$ could be very stale (if $i$ hasn't been selected in a minibatch for a long time)?\nFor Section 6.1 line 227, doesn't the GMM have an analytic closed form solution for $\\theta$ given $x,z$? Why resort to gradient based optimization?\nFor Section 6.3. I do not know if a large $\\Delta \\mathcal{L}^{test}$ is meaningful or not. Especially for large $K$. How do you interpret this practically? And for models with $K=16384$ do we expect convergence given the # iteration, $M$ and $N$? \nIn the conclusion line 316, MH when *run to convergence* ensures sampling from the proper posterior, here we explicitly are only taking $M = 2$ steps. How does this approximation affect the accuracy vs runtime tradeoff? For more complex models, where the runtime benefit of avoiding likelihood computations is larger, but this approximation is less accurate does this algorithm help or hurt?\n\nClarity: The paper's writing and organization is okay, but the paper could use additional editing to improve grammar + use more professional/less casual tone. For example: the paragraph around line 77.\n  \nSignificance: Why are we interested in fitting mixture models with large $K$? \n\nOther notes:\n* line 49: mdoels -> models\n* order of $TKND$ in $O(TKND)$ isn't consistent throughout\n* not necessary, but would be nice to include comparisons against gibbs sampling or VI to make the experiments stronger. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "N/A.", "review_text": "", "extra_scores": {"ethical_concerns": "N/A", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "38waq0hLS6f", "reviewer": "Reviewer_iZc8", "summary": "This paper proposes a new method for scalably fitting mixture models. The authors tackle the problem of fitting with a large number of mixture components. While there exist some works handling this problem, these methods are all limited to mixtures of exponential family distributions, whereas the authors’ method applies to generic mixture models (as long as the model provides a likelihood p(data point i | cluster k)). They do so by proposing a modified version of the EM algorithm: instead of exactly computing the expectations in the E step, the authors propose to use a Metropolis Hastings algorithm to approximate the required expectations. In experiments, the authors show that their method (1) outperforms pre-existing large-K methods (where such previous applications are applicable), and (2) show that their method outperforms standard stochastic gradient descent in settings where other large-K methods are not applicable. The authors’ method seems to work well, but I have a few clarification questions about the technical ideas behind the proposed algorithm and its empirical evaluation. Currently, I put a borderline score for the paper, but I’m definitely open to increasing it after discussion.", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "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": "---------------------------------------------\n\n**Update after discussion with authors and other reviewers**\n\nThanks to the authors for really engaging in the discussion period. I think it's become clear that some major revisions are needed in order to fully describe / experiment with the proposed E-step approximation. I've actually changed my score to a weak reject because of this. I do want to be clear that this is only because of the limits of the conference review cycle, where reviewers don't see an updated draft. In particular, I think a revision is likely to result in a strong accept, but the revisions are too substantial to accept without seeing them. Just to summarize, there are two things that I think remain:\n\n1. The main issue is the need for a more targeted discussion and experimentation surrounding the proposed E-step approximation. I think the directions in terms of references that the authors have brought up will strengthen things here a lot. And the proposed experiment seems like a great step in the right direction on the empirics front. Personally, I think using a continuing Markov chain between optimization iterations is pretty interesting and surprising, so I really think the paper would be both more complete and more interesting with a deeper investigation of this idea.\n\n2. This is a lot more minor, but I do want to argue a little bit more about the stopping rule used in the experiments. I agree that adding some tolerance $\\delta$ will give an extra degree of freedom in the experiments. But, since this is how the method will be used in practice, it seems worth looking into. I don't mean that it's necessary to have a whole extra set of experiments investigating the sensitivity with respect to $\\delta$. It's just good to make the experiments as close to reality as possible, and it doesn't seem like there's a strong reason to not make the experiments more realistic in this regard.\n\n\n\n---------------------------------------------\n\nI’ve listed my main comments/questions as (A), (B), (C) below.\n\n(A) The use of Monte Carlo to help with an intractable E step is not a new idea, and I think the paper needs more discussion of the related work in this area. For example, see Levine and Casella (2001) (reference below). Ideas from this literature might be helpful for the current paper – e.g. Levine and Casella (2001) recommend updating the Monte Carlo samples from the previous E step via importance sampling, rather than re-running a MCMC method. Still, I don’t know that other work has handled the issue of a high-dimensional parameter space (here, the number of mixture components K is large), which sounds like an important addition over pre-existing work.\n\n(B) I have a few technical questions about the presented algorithm.\n1. The hope of the method is that we can effectively approximate the E step by drawing $M \\ll K$ samples from the distribution $p(z_i | x_i)$. It seems that in some sense this is automatically impossible – the distribution $p(z_i | x_i)$ has support of size $K$, so we can’t, in general, approximate expectations over this distribution with only $M$ samples. So I think it’s unlikely that, in general, “the chain passes the transition period (the burn-in phase)” as the authors note. Presumably the authors’ method either relies on some regularity conditions (maybe the distribution $p(z_i | x_i)$ is mostly concentrated on a few points?) or maybe performing the M-step with an only approximate objective is fine, or maybe something else is happening. In any case, can the authors clarify with discussion and/or additional experiments what is driving the success of their method?\n2. On lines 130-138, I think the authors are proposing to somehow reuse the samples from $p(z_i | x_i)$ between iterations. But the distribution being sampled from is changing from iteration to iteration (i.e. it depends on the parameters $\\theta$ estimated from the M-step). What do the authors mean on line 131 when they talk about samples being “recycled over the iterations”?\n3. The authors parameterize the mixing weights $\\pi_k$ by a $K$ real-valued parameters $\\nu_1, \\dots, \\nu_K$, and then set $\\pi_k = e^{\\nu_k} / \\sum_{l=1}^K e^{\\nu_l}$. But this transformation leaves a degree of freedom unconstrained (i.e. scaling all $\\nu_k$ by the same constant $c$ will leave the $\\pi$’s unchanged), meaning there is no unique optimum of the log-likelihood. This seems potentially undesirable, but maybe isn’t a huge deal.\n\n(C) I also have a few questions about the empirical evaluation of the authors’ method.\n1. For the purposes of considering an algorithm terminated, the authors stop each algorithm when it reaches 95% of the maximum likelihood value that it will achieve. I have two questions about this. The first is why is this a reasonable metric? Can we be sure that the fitted parameters will not change much as we traverse the last 5% of the loss? Second, this doesn’t necessarily seem like a good way to assess convergence of practical algorithms. In particular, the only way to define this stopping criteria is to first run the algorithm to termination so that we know what the final log likelihood value is, and then go back and stop it 95% of the way there. Can the authors argue that this does not bias the results at all? For example, it could be that it takes much longer for their algorithm to traverse the final 5% of the likelihood than other algorithms.\n2. The authors choose a particular annealing schedule for use across all algorithms. How did they choose this schedule? Does it definitely not bias the results in favor of one algorithm or another?\n3. Are all the experiments done with multiple random restarts for each method? Figure 1 says that “the results are averaged over five repetitions”, but it’s not clear to me how these repetitions are done. Are these five different random initializations? Or five different runs from the same initialization? Each curve in Figure 1 looks like the run from a single random seed, so it seems to be the latter. I think there should be multiple random seeds used so that we know that, for example, TSAEM (the fastest algorithm to reach its optimum in Figure 1) isn’t just getting unluckily stuck in a local optimum in Figure 1.\n4. Why isn’t MHSAEM-TF used in the experiment in Section 6.2? The authors state that it is due to the “favorable computational complexity and simplicity” of MHSAEM-U, but it seems like in Figure 2 the U and TF versions of these algorithms have just about the same performance.\n5. Similarly to point #3, are the results in Table 2 run over multiple random restarts? If so, there should be standard errors reported so that we know how significant the differences are.\n\n-----Citations--------\n\nRichard A. Levine and George Casella. Implementations of the Monte Carlo EM Algorithm. Journal of Computational and Graphical Statistics , Sep., 2001, Vol. 10, No. 3\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors do include a broader impact statement, but it's over the nine page limit. Note that the broader impact statement is supposed to be included within the page limit of the main text.\n\nIn their paper checklist, the authors write [N/A] for whether their datasets were obtained with consent and whether their datasets contain personally identifiable information or offensive content. Their datasets were taken from the UCI repository, so presumably neither of these things is true, but that doesn't mean it isn't relevant to this paper. I think a quick appendix about the data used would resolve this issue.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}]

Author rebuttal:

Meta-review: Reviewers are intrigued by the proposed approach to learning mixture models, but voice a number of substantial concerns related to both the clarity and correctness of the technical material on Monte Carlo approximations.  The author response was considered, but for some reviewers, it raised more questions than it answered.  After some discussion, the consensus was that very-major revisions are needed to clearly explain and justify the proposed method, and that such changes require a full review at some future conference.
Decision: Reject