Year: 2019

Abstract:
In this paper, we develop an unsupervised generative clustering framework that combines variational information bottleneck and the Gaussian Mixture Model. Specifically, in our approach we use the variational information bottleneck method and model the latent space as a mixture of Gaussians. We derive a bound on the cost function of our model that generalizes the evidence lower bound (ELBO); and provide a variational inference type algorithm that allows to compute it. In the algorithm, the coders’ mappings are parametrized using neural networks and the bound is approximated by Markov sampling and optimized with stochastic gradient descent. Numerical results on real datasets are provided to support the efficiency of our method.

text: In this paper, we develop an unsupervised generative clustering framework that combines variational information bottleneck and the Gaussian Mixture Model. Specifically, in our approach we use the variational information bottleneck method and model the latent space as a mixture of Gaussians. We derive a bound on the cost function of our model that generalizes the evidence lower bound (ELBO); and provide a variational inference type algorithm that allows to compute it. In the algorithm, the coders’ mappings are parametrized using neural networks and the bound is approximated by Markov sampling and optimized with stochastic gradient descent. Numerical results on real datasets are provided to support the efficiency of our method.

heading: 1 INTRODUCTION
text: Clustering consists in partitioning a given data set into various groups (clusters) based on some similarity metric, such as Euclidean distance, L1 norm, L2 norm, L∞ norm, the popular logarithmic loss measure or others. The principle is that each cluster should contain elements of the data that are closer to each other than to any other element outside that cluster, in the sense of the defined similarity measure. If the joint distribution of the clusters and data is not known, one should operate blindly in doing so, i.e., using only the data elements at hand; and the approach is called unsupervised clustering. Unsupervised clustering is perhaps one of the most important tasks of unsupervised machine learning algorithms nowadays, due to a variety of application needs and connections with other problems.
Examples of unsupervised clustering algorithms include the so-popular K-means (Hartigan & Wong, 1979) and expectation maximization (EM) (Dempster et al., 1977). TheK-means algorithm partitions the data in a manner that the Euclidean distance among the members of each cluster is minimized. With the EM algorithm, the underlying assumption is that the data comprises a mixture of Gaussian samples, namely a Gaussian Mixture Model (GMM); and one estimates the parameters of each component of the GMM while simultaneously associating each data sample to one of those components. Although they offer some advantages in the context of clustering, these algorithms suffer from some strong limitations. For example, it is well known that the K-means is highly sensitive to both the order of the data and scaling; and the obtained accuracy depends strongly on the initial seeds (in addition to that it does not predict the number of clusters or K-value). The EM algorithm suffers mainly from low convergence, especially for high dimensional data.
Recently, a new approach has emerged that seeks to perform inference on a transformed domain (generally referred to as latent space), not the data itself. The rationale is that because the latent space often has fewer dimensions it is more convenient computationally to perform inference (clustering) on it rather than on the high dimensional data directly. A key aspect then is how to design a latent space that is amenable to accurate low-complex unsupervised clustering, i.e., one that preserves only those features of the observed high dimensional data that are useful for clustering while removing out all redundant or non-relevant information. Along this line of work, we can mention (Ding & He, 2004) which utilizes Principal Component Analysis (PCA) (Wold et al., 1987) for dimensionality reduction followed by K-means for clustering the obtained reduced dimension data; or (Roweis, 1997) which uses a combination of PCA and the EM algorithm. Other works that use alternatives
for the linear PCA include Kernel PCA (Hofmann et al., 2008), which employs PCA in a non-linear fashion to maximize variance in the data.
The usage of deep neural networks (DNN) for unsupervised clustering of high dimensional data on a lower dimensional latent space has attracted considerable attention, especially with the advent of autoencoder (AE) learning and the development of powerful tools to train them using standard backpropagation techniques (Kingma & Welling, 2014; Rezende et al., 2014). Advanced forms include Variational autoencoders (VAE) (Kingma & Welling, 2014; Rezende et al., 2014) which are generative variants of AE that regularize the structure of the latent space and the more general Variational Information Bottleneck (VIB) of (Alemi et al., 2017) which is a technique that is based on the Information Bottleneck method (Tishby et al., 1999) and seeks a better trade-off between accuracy and regularization than VAE via the introduction of a Lagrange-type parameter s which controls that trade-off and whose optimization is similar to deterministic annealing (Slonim, 2002) or stochastic relaxation.
In this paper, we develop an unsupervised generative clustering framework that combines VIB and the Gaussian Mixture Model. Specifically, in our approach we use the variational information bottleneck method and model the latent space as a mixture of Gaussians. The encoder and decoder of the model are parametrized using neural networks (NN). The cost-function is calculated approximatively by Markov sampling and optimized with stochastic gradient descent. Furthermore, the application of our algorithm to the unsupervised clustering of various datasets, including the MNIST (Lecun et al., 1998), REUTERS (Lewis et al., 2004) and STL-10 (Coates et al., 2011), allows a better clustering accuracy than previous state of the art algorithms. For instance, we show that our algorithm performs better than the variational deep embedding (VaDE) algorithm of (Jiang et al., 2017) which is based on VAE and performs clustering by maximizes the ELBO and can be seen as a specific case of our algorithm (Section 3.1). Our algorithm also generalizes the VIB of (Alemi et al., 2017) which models the latent space as an isotropic Gaussian which is generally not expressive enough for the purpose of unsupervised clustering. Other related works, but which are of lesser relevance to the contribution of this paper, are the deep embedded clustering (DEC) of (Xie et al., 2016), the improved deep embedded clustering (IDEC) of (Guo et al., 2017) and (Dilokthanakul et al., 2017). For a detailed survey of clustering with deep learning, the readers may refer to (Min et al., 2018).
To the best of our knowledge, our algorithm performs the best in terms of clustering accuracy by using deep neural networks without any prior knowledge regarding the labels (except the usual assumption regarding the number of the classes) compared to the state-of-the-art algorithms of this category. In order to achieve the aforementioned accuracy, i) we derive a cost-function that contains the IB hyper parameter s that controls the trade-off between over-fit and generalization of the model and we used an approximation of KL divergence that avoid assumptions which do not hold in the beginning of the learning process and lead to convergence issues; ii) evaluate the hyper-parameter s by following an annealing approach that improves both the convergence and the accuracy of the proposed algorithm.

heading: 2 PROBLEM DEFINITION AND MODEL
text: Consider a dataset that is composed of N samples {xi}Ni=1 which we wish to partition into |C| ≥ 1 clusters. Let C = {1, . . . , |C|} be the set of all possible clusters; and C designate a categorical random variable that lies in C and stands for the index of the actual cluster. If X is a random variable that models elements of the dataset, given X = xi induces a probability distribution on C which the learner should learn. Thus, mathematically the problem is that of estimating the values of the
unknown conditional probability PC|X(·|xi) for all elements xi of the dataset. The estimates are sometimes referred to as the assignment probabilities.
As mentioned previously, we use the VIB framework and model the latent space as a GMM. The resulting model is depicted in Figure 1, where the parameters πc, µc, Σc, for all values of c ∈ C, are to be optimized jointly with those of the employed NNs as instantiation of the coders. Also, the assignment probabilities are estimated based on the values of latent space vector instead of the observation themselves, i.e., PC|U = QC|U. In the rest of this section, we elaborate on the inference and generative network models for our method, which are illustrated below.
C X U PX|C PU|X
Figure 2: Inference Network
C U X QU|C QX|U
Figure 3: Generative Network

heading: 2.1 INFERENCE NETWORK MODEL
text: We assume that an observed data x is generated from a GMM with |C| components. Then, the latent representation u is inferred according the following procedure:
1. One of the components of the GMM is chosen according to a categorical variable C.
2. The data x is generated from the c-th competent of the GMM, i.e., PX|C ∼ N (x; µ̃c, Σ̃c). 3. Encoder maps x to a latent representation u according to PU|X ∼ N (µθ,Σθ).
3.1. The encoder is modeled with a DNN fθ which maps x to the parameters of a Gaussian distribution, i.e., [µθ,Σθ] = fθ(x). 3.2. The representation u is sampled from N (µθ,Σθ).
For the inference network, shown in Figure 2, the following Markov chain holds
C − −X− −U . (1)

heading: 2.2 GENERATIVE NETWORK MODEL
text: Since encoder extracts useful representations of the dataset and we assume that the dataset is generated from a GMM, we model our latent space also with a mixture of Gaussians. To do so, the categorical variable C is embedded with the latent variable U. The reconstruction of the dataset is generated according to the following procedure:
1. One of the components of the GMM is chosen according to a categorical variable C, with a prior distribution QC .
2. The representation u is generated from the c-th component, i.e., QU|C ∼ N (u;µc,Σc). 3. The decoder maps the latent representation u to x̂ which is the reconstruction of the source
x by using the mapping QX|U.
3.1. The decoder is modeled with a DNN gφ, that maps u to the estimate x̂, i.e., [x̂] = gφ(u).
For the generative network, shown in Figure 3, the following Markov chain holds
C − −U− −X . (2)

heading: 3 PROPOSED METHOD
text: In this section we present our clustering method. First, we provide a general cost function for the problem of the unsupervised clustering that we study here based on the variational IB framework; and we show that it generalizes the ELBO bound developed in (Jiang et al., 2017). We then parametrize our model using NNs whose parameters are optimized jointly with those of the GMM. Furthermore, we discuss the influence of the hyper-parameter s that controls optimal trade-offs between accuracy and regularization.

heading: 3.1 BRIEF REVIEW OF VARIATIONAL INFORMATION BOTTLENECK FOR UNSUPERVISED LEARNING
text: As described in Chapter 2, the stochastic encoder PU|X maps the observed data x to a representation u. Similarly, the stochastic decoder QX|U assigns an estimate x̂ of x based on the vector u. As per the IB method (Tishby et al., 1999) a suitable representation U should strike the right balance between capturing all information about the categorical variable C that is contained in the observation X and using the most concise representation for it. This leads to maximizing the following Lagrange problem Ls(P) = I(C; U)− sI(X; U) , (3) where s ≥ 0 designates the Lagrange multiplier and for convenience P denotes the conditional distribution PU|X.
Instead of equation 3 which is not always computable in our unsupervised clustering setting, we use a modified version of it (so-called unsupervised IB objective (Alemi et al., 2017)) given by
L̃s(P) : = −H(X|U)− s[H(U)−H(U|X)] (4) = EPX [ EPU|X [logPX|U + s logPU − s logPU|X] ] . (5)
For a variational distribution QU on U (instead of the unknown PU) and a variational stochastic decoder QX|U (instead of the unknown optimal decoder PX|U), let Q := {QX|U, QU}. Also, let
LVBs (P,Q) := EPX [ EPU|X [logQX|U]− sDKL(PU|X‖QU) ] . (6)
Lemma 1. For given P, we have
LVBs (P,Q) ≤ L̃s(P), for all Q . In addition, there exists a unique Q that achieves the maximum maxQ LVBs (P,Q) = L̃s(P), and is given by
Q∗X|U = PX|U , Q ∗ U = PU .
Using Lemma 1, maximization of equation 4 can be written in term of the variational IB cost as follows
max P L′s(P) = max P max Q LVBs (P,Q) . (7)
Remark 1. As we already mentioned in the beginning of this chapter, the related work (Jiang et al., 2017) performs unsupervised clustering by combining VAE with GMM. Specifically, it maximizes the following ELBO bound
LVaDE1 := EPX [ EPU|X [logQX|U]−DKL(PC|X‖QC)− EPC|X [DKL(PU|X‖QU|C)] ] . (8)
Let, for an arbitrary non-negative parameter s, LVaDEs be a generalization of the ELBO bound in equation 8 of (Jiang et al., 2017) given by
LVaDEs := EPX [ EPU|X [logQX|U]− sDKL(PC|X‖QC)− sEPC|X [DKL(PU|X‖QU|C)] ] . (9)
Investigating the RHS of equation 9, we get
LVBs (P,Q) = LVaDEs + sEPX [ EPU|X [DKL(PC|X‖QC|U)] ] . (10)
Thus, by the non-negativity of relative entropy it is clear that LVaDEs is a lower bound on LVBs (P,Q). Also, if variational distribution Q is such that the conditional marginal QC|U is equal to PC|X the bound is tight since the relative entropy term is zero in this case.

heading: 3.2 PROPOSED ALGORITHM: VIB-GMM
text: In order to compute equation 7, we parametrize the distributions PU|X and QX|U using DNNs. For instance, let the stochastic encoder PU|X be a DNN fθ and the stochastic decoder QX|U be a DNN gφ. That is
Pθ(u|x) = N (u;µθ,Σθ) , where [µθ,Σθ] = fθ(x) , Qφ(x|u) = gφ(u) = [x̂] ,
(11)
where θ and φ are the weight and bias parameters of the DNNs. Furthermore, the latent space is modeled as a GMM with |C| components with parameters ψ := {πc,µc,Σc}|C|c=1, i.e.,
Qψ(u) = ∑
c
πc N (u;µc,Σc) . (12)
Using the parametrizations above, the optimization of equation 7 can be rewritten as
max θ,φ,ψ
LNNs (θ, φ, ψ) (13)
where the cost function LNNs (θ, φ, ψ) given by
LNNs (θ, φ, ψ) := EPX [ EPθ(U|X)[logQφ(X|U)]− sDKL(Pθ(U|X)‖Qψ(U)) ] . (14)
Then, for a given observations of N samples, i.e., {xi}Ni=1, equation 13 can be approximated in terms of an empirical cost as follows
max θ,φ,ψ
1
n
n∑
i=1
Lemps,i (θ, φ, ψ) , (15)
where Lemps,i (θ, φ, ψ) is the empirical cost for the i-th observation xi, and given by
Lemps,i (θ, φ, ψ) = EPθ(Ui|Xi)[logQφ(Xi|Ui)]− sDKL(Pθ(Ui|Xi)‖Qψ(Ui)) . (16)
Furthermore, the first term of the RHS of equation 16 can be computed using Monte Carlo sampling and the re-parametrization trick (Kingma & Welling, 2014). In particular, Pθ(u|x) can be sampled by first sampling a random variable Z with distribution PZ, i.e., PZ = N (0, I), then transforming the samples using some function f̃θ : X × Z → U , i.e., u = f̃θ(x, z). Thus,
EPθ(Ui|Xi)[logQφ(Xi|Ui)] = 1
M
M∑
m=1
log q(xi|ui,m), ui,m = µθ,i+Σ 1 2 θ,i· m, m ∼ N (0, I) ,
where M is the number of samples for the Monte Carlo sampling step.
The second term of the RHS of equation 16 is the KL divergence between a single component multivariate Gaussian and a Gaussian Mixture Model with |C| components. An exact closed-form solution for the calculation of this term does not exist. However, a variational lower bound approximation (Hershey & Olsen, 2007) of it can be obtained as
DKL(Pθ(Ui|Xi)‖Qψ(Ui)) = − log |C|∑
c=1
πc exp (−DKL(N (µθ,i,Σθ,i)‖N (µc,Σc)) . (17)
In particular, in the specific case in which the covariance matrices are diagonal, i.e., Σθ,i := diag({σ2θ,i,j}nuj=1) and Σc := diag({σ2c,j}nuj=1), with nu denoting the latent space dimension, equation 17 can be computed as follows
DKL(Pθ(Ui|Xi)‖Qψ(Ui))
= − log |C|∑
c=1
πc exp
( − 1
2
nu∑
j=1
[ (µθ,i,j − µc,j)2 σ2c,j + log σ2c,j σ2θ,i,j − 1 + σ2θ,i,j σ2c,j ]) , (18)
where µθ,i,j and σ2θ,i,j are the mean and variance of the i-th representation in the j-th dimension of the latent space. Furthermore, µc,j and σ2c,j represent the mean and variance of the c-th component of the GMM in the j-th dimension of the latent space.
Finally, we train NNs to maximize the cost function equation 14 over the parameters θ, φ, as well as those ψ of the GMM. For the training step, we use the ADAM optimization tool (Kingma & Ba, 2015). The training procedure is detailed in Algorithm 1.
Once our model is trained, we assign the given dataset into the clusters. As mentioned in Section 2, we do the assignment from the latent representations, i.e., QC|U = PC|X. Hence, the probability that the observed data xi belongs to the c-th cluster is computed as follows
p(c|xi) = q(c|ui) = qψ?(c)qψ?(ui|c)
qψ?(ui) = π?cN (ui;µ?c ,Σ?c)∑ c π ? cN (ui;µ?c ,Σ∗c) , (19)
where ? indicates optimal values of the parameters as found at the end of the training phase. Finally, the right cluster is picked based on the largest assignment probability value.
Algorithm 1 VIB-GMM algorithm for unsupervised learning
1: input: Dataset D := {xi}ni=1, parameter s ≥ 0. 2: output: Optimal DNN weights θ?, φ? and GMM parameters ψ? = {π?c , µ?c , Σ?c}|C|c=1. 3: initialization Initialize θ, φ, ψ. 4: repeat 5: Randomly select b mini-batch samples {xi}bi=1 from D. 6: Draw m random i.i.d samples {zj}mj=1 from PZ. 7: Compute m samples ui,j = f̃θ(xi, zj) 8: For the selected mini-batch, compute gradients of the empirical cost equation 15. 9: Update θ, φ, ψ using the estimated gradient (e.g. with SGD or ADAM).
10: until convergence of θ, φ, ψ.
Remark 2. It is worth to mention that with the use of the KL approximation in equation 17, our algorithm does not use the assumption PC|U = QC|U (not as in Jiang et al. (2017)), which does not hold in the beginning of the training phase and leads to convergence issues. This assumption is only used in the final assignment after the training phase is over.

heading: 3.3 EFFECT OF THE HYPER-PARAMETER
text: Algorithm 2 Annealing Algorithm Pseudo-Code
input: Dataset D := {xi}ni=1, hyper-parameter interval [smin, smax]. output: Optimal DNN weights θ?, φ?, GMM parameters ψ? = {π?c , µ?c , Σ?c}|C|c=1, assignment probability PC|X. initialization Initialize θ, φ, ψ. repeat
Apply VIB-GMM algorithm. Update ψ, θ, φ. Update s, e.g., s = (1 + s)sold.
until s does not exceed smax. As we already mentioned, the hyper-parameter s controls the trade-off between the relevance of the representation U and its complexity. As it can be seen from equation 14 for small values of s, it is the cross-entropy term that dominates, i.e., the algorithm trains the parameters so as to reproduce X as accurate as possible. For large values of s, however, it is most important for the NN to produce an encoded version of X whose distribution matches the prior distribution of the latent space, i.e., the term DKL(Pθ(U|X)‖Qψ(U)) is nearly zero. In the beginning of the training process, the GMM components are randomly selected; and so starting with a large value of the hyper-parameter s is likely to steer the solution towards an irrelevant prior. Hence, for the tunning of the hyper-parameter s in practice it is more efficient to start with a small value of s and gradually increase it with the number of epochs. This has the advantage to avoid possible local minimas, an aspect that is reminiscent of deterministic annealing (Slonim, 2002), where s plays the role of the temperature parameter. The experiments that will be reported in the next section show that proceeding in the above described manner for the selection of the parameter s helps getting better accuracy results and better robustness to the initialization (i.e., no need for a strong pretraining). A pseudo-code for annealing is given in Algorithm 2. We note that tuning s is very critical, such that the step size s in update of s should be chosen carefully, otherwise phase transitions might be skipped that would cause a bad ACC score.

heading: 4 EXPERIMENTS

heading: 4.1 DESCRIPTION OF USED DATASETS
text: In our empirical experiments, we apply our algorithm to the clustering of the following datasets.
MNIST: A dataset of gray-scale images of 70000 handwritten digits of dimensions 28× 28 pixel. STL-10: A dataset of color images collected from 10 categories. Each category consists of 1300 images of size of 96× 96 (pixels) ×3 (rgb code). Hence, the original input dimension nx is 27648. For this dataset, we use a pretrained convolutional NN model, i.e., ResNet-50 (He et al., 2016) to reduce the dimensionality of the input. This preprocessing reduces the input dimension to 2048. Then, our algorithm and other baselines are used for clustering.
REUTERS10K: A dataset that is composed of 810000 English stories labeled with a category tree. As in (Xie et al., 2016), 4 root categories (corporate/industrial, government/social, markets, economics) are selected as labels and all documents with multiple labels are discarded. Then, tf-idf features are computed on the 2000 most frequently occurring words. Finally, 10000 samples are taken randomly, which are referred to as REUTERS10K dataset.

heading: 4.2 NETWORK SETTINGS AND OTHER PARAMETERS
text: We use the following network architecture: the encoder is modeled with NNs with 3 hidden layers with dimensions nx−500−500−2000−J , where nx is the input dimension and nu is the dimension of the latent space. The decoder consists of NNs with dimensions nu − 2000− 500− 500− nx. All layers are fully connected. For comparison purposes, we chose the architecture of the hidden layers as well as the dimension of the latent space nu = 10 to coincide with those made for the DEC algorithm of (Xie et al., 2016) and the VaDE algorithm of (Jiang et al., 2017). All except the last layers of the encoder and decoder are activated with ReLU function. For the last (i.e., latent) layer of the encoder we use a linear activation; and for the last (i.e., output) layer of the decoder we use sigmoid function for MNIST and linear activation for the remaining datasets. The batch size is 100 and the variational bound equation 15 is maximized by the Adam optimizer of (Kingma & Ba, 2015). The learning rate is initialized with 0.002 and decreased gradually every 20 epochs with a decay rate of 0.9 until it reaches a small value (0.0005 is our experiments). The reconstruction loss is calculated by using the cross-entropy criterion for MNIST and mean squared error function for the other datasets.

heading: 4.3 CLUSTERING ACCURACY
text: We evaluate the performance of our algorithm in terms of the so-called unsupervised clustering accuracy (ACC), which is a widely used metric in the context of unsupervised learning (Min et al., 2018). For comparison purposes, we also present those of algorithms from previous art.
For each of the aforementioned datasets, we run our VIB-GMM algorithm for various values of the hyper-parameter s inside an interval [smin, smax], starting from the smaller valuer s1 and gradually increasing the value of s every nepoch epochs. For the MNIST dataset, we set (smin, smax, nepoch) = (1, 5, 500); and for the STL-10 dataset and the REUTERS10K datset we choose these parameters to be (1, 20, 500) and (1, 5, 100), respectively. The obtained ACC accuracy results are reported in the Table 1 from which it can be seen that our algorithm outperforms significantly the DEC algorithm of (Xie et al., 2016) as well as the VaDE algorithm of (Jiang et al., 2017) and GMM on the same datsets. Important to note, for the MNIST dataset the reported ACC accuracy of 96.2% using our VIBGMM algorithm is obtained as the best case run out of ten times run all with random initializations. For instance, we do not use any pretrained values for the initialization of our algorithm in sharp contrast with the VaDE of (Jiang et al., 2017) and the DEC of (Xie et al., 2016). For the STL-10 dataset, none of the compared algorithms use a pretrained network except the intimal ResNet-50 for dimensionality reduction. For REUTERS10K, we used the same pretrain parameters as DEC and VaDE. Figure 4 depicts the evolution of the ACC accuracy with iterations (number of epochs) for the four compared algorithms.
Figure 5 shows the evolution of the reconstruction loss of our VIB-GMM algorithm for the STL-10 dataset, as a function of simultaneously varying values of the hyper-parameter s and the number of epochs (recall that, as per-the described methodology, we start with s = s1 and we increase its value gradually every nepoch = 500 epochs). As it can be seen from the figure, the few first epochs are spent
almost entirely on reducing the reconstruction loss (i.e., a fitting phase) and most of the remaining epochs are spent in making the found representation more concise (i.e., smaller KL-divergence). This is reminiscent of the two-phase (fitting v.s. compression) that was observed for supervised learning using VIB in (Schwartz-Ziv & Tishby, 2017).

heading: 4.4 VISUALIZATION ON THE LATENT SPACE
text: In this section, we investigate the evolution of the unsupervised clustering of the STL-10 dataset on the latent space using our VIB-GMM algorithm. For this purpose, we find it convenient to visualize the latent space through application of the t-SNE algorithm of (van der Maaten & Hinton, 2008) in order to generate meaningful representations in a two-dimensional space. Figure 6 shows 4000 randomly chosen latent representations before the start of the training process and respectively after 1, 5 and 500 epochs. The shown points (with a · marker in the figure) represent latent representations of data samples whose labels are identical. Colors are used to distinguish between clusters. Crosses (with an x marker in the figure) correspond to the centroids of the clusters. More specifically, Figure 6-(a) shows the initial latent space before the training process. If the clustering is performed on the initial representations it allows ACC accuracy of as small as 10%, i.e., as bad as a random assignment. Figure 6-(b) shows the latent space after one epoch, from which a partition of some of the points starts to be already visible. With five epochs, that partitioning is significantly sharper and the associated clusters can be recognized easily. Observe, however, that the cluster centers seem still not to have converged. With 500 epochs, the ACC accuracy of our algorithm reaches %91.6 and the clusters and their centroids are neater as visible from Figure 6-(d).

heading: A THE PROOF OF LEMMA 1
text: First, we expand L̃s(P) as follows L̃s(P) =−H(X|U)− sI(X; U)
=−H(X|U)− s[H(U)−H(U|X)]
=
∫∫
ux
p(u,x) log p(x|u) du dx
+ s
∫
u
p(u) log p(u) du− s ∫∫
ux
p(u,x) log p(u|x) du dx.
Then, LVBs (P,Q) is defined as follows
LVBs (P,Q) := ∫∫
ux
p(u,x) log q(x|u) du dx
+ s
∫
u
p(u) log q(u) du− s ∫∫
ux
p(u,x) log p(u|x) du dx. (20)
Hence, we have the following relation
L̃s(P)− LVBs (P,Q) = EPX [DKL(PX|U‖QX|U)] + sDKL(PU‖QU) ≥ 0 where equality holds under equalities QX|U = PX|U and QU = PU. We note that s ≥ 0.
B THE PROOF OF ALTERNATIVE EXPRESSION LVADEs Here we show how we obtained equation 10.
To do so,
LVaDEs = EPX [ EPU|X [logQX|U]− sDKL(PU|X‖QU)− sEPU|X [ DKL(PC|X‖QC|U) ]
= EPX [ EPU|X [logQX|U] ] − s
∫
x
p(x)
∫
u
p(u|x) log p(u|x) q(u) du dx
− s ∫
x
p(x)
∫
u
p(u|x) ∑
c
p(c|x) log p(c|x) q(c|u) du dx
(a) = EPX [ EPU|X [logQX|U] ] − s
∫∫
ux
p(x)p(u|x) log p(u|x) q(u) du dx
− s ∫∫
ux
∑
c
p(x)p(u|c,x)p(c|x) log p(c|x) q(c|u) du dx
= EPX [ EPU|X [logQX|U] ] − s
∫∫
ux
∑
c
p(u, c,x) log p(u|x)p(c|x) q(u)q(c|u) du dx
= EPX [ EPU|X [logQX|U] ] − s
∫∫
ux
∑
c
p(u, c,x) log p(c|x) q(c) p(u|x) q(u|c) du dx
= EPX [ EPU|X [logQX|U] ] − s
∫
x
∑
c
p(c,x) log p(c|x) q(c) dx
− s ∫∫
ux
∑
c
p(x)p(c|x)p(u|c,x) log p(u|x) q(u|c) du dx
(b) = EPX [ EPU|X [logQX|U]− sDKL(PC|X‖QC)− sEPC|X [DKL(PU|X‖QU|C)] ]
(c) = LVBs (P,Q)− sEPX [ EPU|X [ DKL(PC|X‖QC|U) ]] ,
where (a) and (b) follow due to the Markov chain C − −X− −U; (c) follows from the definition of LVBs (P,Q) in equation 6.

heading: C KL DIVERGENCE BETWEEN MULTIVARIATE GAUSSIAN DISTRIBUTIONS
text: The KL divergence between two multivariate Gaussian distributions P1 ∼ N (µ1,Σ1) and P2 ∼ N (µ2,Σ2) in RJ is
DKL(P1‖P2) = 1
2
( (µ1 −µ2)TΣ−12 (µ1 −µ2) + log |Σ2| − log |Σ1| − J + tr(Σ−12 Σ1) ) . (21)
For the case in which Σ1 and Σ2 covariance matrices are diagonal, i.e., Σ1 := diag({σ21,j}Jj=1) and Σ2 := diag({σ22,j}Jj=1), equation 21 boils down to the following
DKL(P1‖P2) = 1
2
( J∑
j=1
(µ1,j − µ2,j)2 σ22,j + log σ22,j σ21,j − 1 + σ21,j σ22,j
) . (22)

heading: D KL DIVERGENCE BETWEEN GAUSSIAN MIXTURE MODELS
text: An exact close form for the calculation of the KL divergence between two Gaussian mixture models does not exist. In this paper, we use a variational lower bound approximation for calculations of KL between two Gaussian mixture models. Let f and g be GMMs and the marginal densities of x under f and g are
f(x) =
M∑
m=1
ωmN (x;µfm,Σfm) = M∑
m=1
ωmfm(x)
g(x) =
C∑
C=1
πcN (x;µgc ,Σgc) = C∑
c=1
πcgc(x).
The KL divergence between two Gaussian mixtures f an g can be approximated as follows
DvKL(f‖g) := M∑
m=1
ωm log ∑ m′∈M\m ωm′ exp (−DKL(fm‖fm′))∑C
c=1 πc exp (−DKL(fm‖gc)) . (23)
In this paper, we are interested, in particular, M = 1. Hence, equation 23 simplifies to
DvKL(f‖g) = − log C∑
c=1
πc exp (−DKL(f‖gc)) , (24)
where DKL(·‖·) is the KL divergence between single component multivariate Gaussian distribution, defined as in equation 21.

Review summaries:
This paper considers the autoencoder model combining the usual information bottleneck and the Gaussian mixture model (GMM). Using an approximation to deal with GMMs, the authors derive a bound on the cost function generalizing the ELBO. The performance of the proposed method is tested on three benchmark datasets and compared with existing methods combining VAE with GMM.

The author(s) posit a Mixture of Gaussian's prior for a compressed latent space representation of high-dimensional data (e.g. images and documents). They propose fitting this model using the Variational Information Bottleneck paradigm and explicate its derivation and tie it to the variational objective used by similar models. They empirically showcase their model and optimization methodology on the MNIST, STL-10, and Reuters10k benchmarks.