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-Hassinns,, 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
39
41
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
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 k 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. .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)m;), 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 ( -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 po (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 [I] 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-Hastting 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 xxpecttionnmmxxiizaaio 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": "t-abNQ-97Ve", "reviewer": "Reviewer_4qCw", "summary": "The paper suggests a method for efficiently evaluating mixture models. The main idea is to use Expectation-Maximization (EM) and Metropolis-Hastings (MH) to evaluate only a small number of mixture components, reducing the computational cost by consequence.\nThe theoretical contributions of the manuscript are two-folded. First, the authors proposed an EM-based algorithm for training, which uses MH for sampling from a smaller set of mixture components. Second, the manuscript generalizes the EM-based algorithm by making it compatible with gradient-based optimization methods. The practical contributions include experiments with Gaussian mixture models (GMMs), sum-product networks (SPNs), and mixtures of real-valued flows. Results show a benchmark speed-up of 100x.", "questions": "", "limitations": "", "rating": 6, "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": "This work is original. As mentioned in the \"Related Work\" section, the main accomplishment of this work is to release a limiting restriction of the previous method. The restriction involves assuming an exponential family for the joint distribution. The helpful Table 1 highlights the theoretical extension proposed in this paper. This extension is non-trivial, and thus, deserving of merit.\n\nThe manuscript proposed methods are sound, to the best of my knowledge. I found the theoretical discussion straightforward to follow and well analyzed. Specifically, Section 4.3 brings a valuable debate on different possible distributions that can be applied in Algorithm 1. These distributions yield various runtime complexities. However, the authors can improve sections 4.2 and 4.3 by writing formal derivations and complexity notation. The discussion here can be loose at times. For example, the claim in \"[...] due to poor mixing properties, the algorithm may converge slowly for high K\" is unclear.\n\nAlthough the paper's writing is sharp and well-structured, the authors could improve the presentation. First, the motivation for evaluating fewer mixture components is not well-founded in real-world needs. The Introduction section cites papers with models that allow for a large number of mixture components. However, I think the paper's motivation could be strengthened by showing real-world applications which demand a considerably large amount of mixture components.\nAnother (minor) presentation issue is the lack of proper discussion on SPNs. Section 6 exploits these models characteristics during experiments, but the manuscript does not provide a formal definition or introduction, as it does for GMMs.\n\nThe manuscript's significance is not clear. In this paper, there are theoretical and practical contributions. The theoretical work is mostly derived from previous methods. There are no new Theorem or extra Proofs.  Although - as a reader - I welcome the simplicity of the theoretical development, the practical contributions have then a higher weight in this work's relevance.\nUnfortunately, it is not clear to me the relevance of the results in the experiments of Section 6. I appreciate the thorough experiments and the impressive 100x average speed-up. However, the manuscript is not clear on the relevance of these speed-ups. For example, Figures 1 and 2 show some runtime improvements in seconds. It is not clear if the minutes saved here are relevant in specific applications or circumstances.\n\nMinor typos: \nSection 1: MH is not defined.\nSection 6.2: SPTNs are not defined or cited.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The discussion around limitations is minor. In the Conclusion section, there is a sentence conjecturing a limit in the speed of convergence. It would be interesting to discuss the possible disadvantages of sampling the mixture components. Could that limit some applications?\n\nSection 8 discusses broader societal impact. In summary, the idea is that the speed-up provided by this paper's method has a beneficial impact on the financial and environmental aspects of a project.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "kEAcHEiSxus", "reviewer": "Reviewer_weQ2", "summary": "The paper presents a modified version of the EM algorithm that allows for faster fitting of mixture models. The approach is based on the idea of evaluating and updating a single mixture component per iteration. Furthermore, they present a generalized version that extends beyond the standard case of the exponential family. \n\nThe authors then follow with an empirical demonstration of the benefits: up to two orders of magnitude faster training while in some cases increasing the likelihood compared to the standard SGD approach. \nThe evaluation is carried on two architectures, SPTNs and mixture of real NVP flows.\n", "questions": "", "limitations": "", "rating": 7, "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-written and technically sound. The authors present an approach that is original and significant especially for the tractable probability models community. However, it can be applied to other sub-fields as demonstrated by the mixture of real NVP flows experiment. \n\nThe paper first presents the problem formulation and follows with the EM algorithm. The paper presents a generalization beyond the exponential family using both an MCMC approach and a stochastic approximation for the updates. For the sampler, they pick the Metropolis-Hastings and follow a rejection sample approach based on an uniform distribution.\n\nThe authors explore the gap between the uniform sampling and an optimal one, noticing potential room for improvements.\nAll in all, the paper is solid and well-founded. \n\nThe authors cover the related works section appropriately, although for the particular case of mini-batch EM in SPNs, a reference is missing to [1]. And this is my only significant question for this paper. Would it make sense to compare to the standard mini-batch EM approach presented in [1]? It would be very nice, as it allows to compare against other papers in the same category although with a reduced scope.  \nIndeed, this ignores one of the advantages of the approach, the generalization beyond the exponentials.\n\nMore precisely, an experiment comparing with M=K, but B < N for the exponential family with closed-form solutions for eq. 8 could be interesting.\n\nNevertheless, I believe the paper is a solid contribution.\n\n\n[1] Peharz, Robert, et al. \"Einsum networks: Fast and scalable learning of tractable probabilistic circuits.\" International Conference on Machine Learning. PMLR, 2020. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors address the impact.", "review_text": "", "extra_scores": {"ethical_concerns": "No Ethical Concerns.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "Oc2qMDFFhH1", "reviewer": "Reviewer_UbDq", "summary": "This work addresses the challenge of applying the EM algorithm to mixture models with a large number of components (K). The key algorithmic idea is using a Metropolis-Hastings sampler for approximating the E-step, which avoids evaluating the posterior probability for K times. ", "questions": "", "limitations": "", "rating": 3, "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": "* Although an important problem is addressed, the proposed approach is very elementary and cannot be justified as an original contribution.\n\n* Using MCMC for estimating E-step in EM algorithms is not new and there is a huge literature on this. To provide useful insights, the authors either need to analyze the bias introduced by using a finite-step MH chain or provide novel solution to correct it.\n\n", "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": "1"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The goal of this paper is to speed up the EM algorithm for mixtures with a large number of mixture components. You might hope to do this by simply replacing the E step with a stochastic sample, but this isn't feasible since all components must be processed to draw one sample. The central idea of this paper is to instead use persistent markov chains to draw approximate samples from this distribution in each iteration.

Reviewers generally felt the paper was clearly written. There were two basic concerns. First, there was a concern that the existing literature on Monte Carlo EM was not sufficiently discussed. The authors suggest that these algorithms typically have a different motivation, but readers would still benefit from better contextualization.

Second, there were concern about the significance of the paper. The algorithm is a relatively simple variant of existing works (albeit perhaps overlooked and with different motivation) and no new theory is presented. Thus, the significance ultimately boils down to the experimental evidence.

Regarding guarantees for the method, it would be better to either provide better guarantees, or to be more explicit that there are no such guarantees. One concern is that the of the method is not clear. At one point in the discussion, the authors suggest that convergence is guaranteed by Thm 1 of [31]. However, this is not so simple. That theorem requires many technical assumptions, which aren't obviously satisfied. In particular, assumption (SAEM3') assumes that the Markov chains is uniformly ergodic in terms of constants K and ρ that can be bounded over all possible Markov chains. Is that satisfied here? Satisfied or not, the authors should at least give some discussion of the issue if this is intended to be used as theoretical evidence for the convergence of the method. (At a later point, the authors explicitly state that there are no guarantees.)

Of course, it is fine to have a method that relies exclusively on experimental evidence. However, some reviewers felt that the current experiments are not totally convincing that the speedups given are significant. There were several constructive suggestions for other experiments.
Decision: Reject