Title: VARIATIONAL INFORMATION BOTTLENECK FOR UN-SUPERVISED CLUSTERING: DEEP GAUSSIAN MIX-TURE EMBEDDING

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.

1. INTRODUCTION:

Clustering consists in partitioning a given data set into various groups (clusters) based on some similarity metric, such as Euclidean distance, L 1 norm, L 2 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) . The K-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.

Encoder 

f θ Decoder g φ U ∼ c π c N (u; µ c , Σ c ) X X

2. PROBLEM DEFINITION AND MODEL:

Consider a dataset that is composed of N samples {x i } N i=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 = x i 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 P C|X (•|x i ) for all elements x i 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., P C|U = Q C|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 P X|C P U|X Figure 2: Inference Network C U X Q U|C Q X|U Figure 3: Generative Network 2.1 INFERENCE NETWORK MODEL

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., P X|C ∼ N (x; μc , Σc ).

3. Encoder maps x to a latent representation u according to P U|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)

2.2. GENERATIVE NETWORK MODEL:

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 Q C .

2.. The representation:

u is generated from the c-th component, i.e., Q U|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 Q X|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)

3. PROPOSED METHOD:

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.

3.1. BRIEF REVIEW OF VARIATIONAL INFORMATION BOTTLENECK FOR UNSUPERVISED:



LEARNING:

As described in Chapter 2, the stochastic encoder P U|X maps the observed data x to a representation u. Similarly, the stochastic decoder Q X|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 L s (P) = I(C; U) -sI(X; U) , (3) where s ≥ 0 designates the Lagrange multiplier and for convenience P denotes the conditional distribution P U|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

Ls (P) : = -H(X|U) -s[H(U) -H(U|X)] (4) = E P X E P U|X [log P X|U + s log P U -s log P U|X ] .

(5)

For a variational distribution Q U on U (instead of the unknown P U ) and a variational stochastic decoder Q X|U (instead of the unknown optimal decoder P X|U ), let

EQUATION

)

Lemma 1. For given P, we have

L VB s (P, Q) ≤ Ls (P)

, for all Q . In addition, there exists a unique Q that achieves the maximum max Q L VB s (P, Q) = Ls (P), and is given by

Q * X|U = P X|U , Q * U = P U .

Using Lemma 1, maximization of equation 4 can be written in term of the variational IB cost as follows

EQUATION

)

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

EQUATION

Let, for an arbitrary non-negative parameter s, L VaDE s be a generalization of the ELBO bound in equation 8 of (Jiang et al., 2017) given by

L VaDE s := E P X E P U|X [log Q X|U ] -sD KL (P C|X Q C ) -sE P C|X [D KL (P U|X Q U|C )] . (9)

Investigating the RHS of equation 9, we get

EQUATION

)

Thus, by the non-negativity of relative entropy it is clear that L VaDE s is a lower bound on L VB s (P, Q). Also, if variational distribution Q is such that the conditional marginal Q C|U is equal to P C|X the bound is tight since the relative entropy term is zero in this case.

3.2. PROPOSED ALGORITHM: VIB-GMM:

In order to compute equation 7, we parametrize the distributions P U|X and Q X|U using DNNs. For instance, let the stochastic encoder P U|X be a DNN f θ and the stochastic decoder Q X|U be a DNN g φ . That is

EQUATION

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 ψ :

EQUATION

Using the parametrizations above, the optimization of equation 7 can be rewritten as

EQUATION

where the cost function L NN s (θ, φ, ψ) given by

EQUATION

Then, for a given observations of N samples, i.e., {x i } N i=1 , equation 13 can be approximated in terms of an empirical cost as follows

EQUATION

where L emp s,i (θ, φ, ψ) is the empirical cost for the i-th observation x i , and given by

EQUATION

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 P Z , i.e., P Z = N (0, I), then transforming the samples using some function fθ : X × Z → U, i.e., u = fθ (x, z). Thus,

E P θ (Ui|Xi) [log Q φ (X i |U i )] = 1 M M m=1 log q(x i |u i,m ), u i,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

EQUATION

)

In particular, in the specific case in which the covariance matrices are diagonal, i.e., Σ θ,i := diag({σ 2 θ,i,j } nu j=1 ) and Σ c := diag({σ 2 c,j } nu j=1 ), with n u denoting the latent space dimension, equation 17 can be computed as follows

EQUATION

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 σ 2 c,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., Q C|U = P C|X . Hence, the probability that the observed data x i belongs to the c-th cluster is computed as follows

EQUATION

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 Compute m samples u i,j = fθ (x i , z j )

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 P C|U = Q C|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. 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 D KL (P θ (U|X) Q ψ (U)) is nearly zero.

3.3. EFFECT OF THE HYPER-PARAMETER:

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.

4. EXPERIMENTS:



4.1. DESCRIPTION OF USED DATASETS:

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 n x 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.

4.2. NETWORK SETTINGS AND OTHER PARAMETERS:

We use the following network architecture: the encoder is modeled with NNs with 3 hidden layers with dimensions n x -500-500-2000-J, where n x is the input dimension and n u is the dimension of the latent space. The decoder consists of NNs with dimensions n u -2000 -500 -500n x . 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 n u = 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. 

4.3. CLUSTERING ACCURACY:

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 [s min , s max ], starting from the smaller valuer s 1 and gradually increasing the value of s every n epoch epochs. For the MNIST dataset, we set (s min , s max , n epoch ) = (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 VIB-GMM 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 = s 1 and we increase its value gradually every n epoch = 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).

4.4. VISUALIZATION ON THE LATENT SPACE:

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) .

Hence, we have the following relation

Ls (P) -L VB s (P, Q) = E P X [D KL (P X|U Q X|U )] + sD KL (P U Q U ) ≥ 0

where equality holds under equalities Q X|U = P X|U and Q U = P U . We note that s ≥ 0.

B THE PROOF OF ALTERNATIVE EXPRESSION L VADE s:

Here we show how we obtained equation 10.

To do so, p(u, c, x) log p(u|x)p(c|x) q(u)q(c|u) du dx

L VaDE s = E P X E P U|X [log Q X|U ] -sD KL (P U|X Q U )

= E P X E P U|X [log Q X|U ] -s ux c

p(u, c, x) log p(c|x) q(c) p(u|x) q(u|c) du dx 

= E P X E P U|X [log Q X|U ] -

C KL DIVERGENCE BETWEEN MULTIVARIATE GAUSSIAN DISTRIBUTIONS:

The KL divergence between two multivariate Gaussian distributions P 1 ∼ N (µ 1 , Σ 1 ) and

EQUATION

)

For the case in which Σ 1 and Σ 2 covariance matrices are diagonal, i.e., Σ 1 := diag({σ 2 1,j } J j=1 ) and Σ 2 := diag({σ 2 2,j } J j=1 ), equation 21 boils down to the following

D KL (P 1 P 2 ) = 1 2 J j=1

(µ 1,jµ 2,j ) 2 σ 2 2,j

EQUATION

D KL DIVERGENCE BETWEEN GAUSSIAN MIXTURE MODELS:

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 .

In this paper, we are interested, in particular, M = 1. Hence, equation 23 simplifies to

EQUATION

where D KL (• •) is the KL divergence between single component multivariate Gaussian distribution, defined as in equation 21.

:



Figures and Tables:

Figure 1: Variational Information Bottleneck with Gaussian Mixtures.

input: Dataset D := {x i } n i=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 {x i } b i=1 from D. 6:Draw m random i.i.d samples {z j } m j=1 from P Z .7:

Figure 4: Accuracy vs. number of epochs for the STL-10 dataset.

Figure 5: Information plane for the STL-10 dataset.

E P X E P U|X [log Q X|U ] -sD KL (P C|X Q C ) -sE P C|X [D KL (P U|X Q U|C )] (c) = L VB s (P, Q) -sE P X E P U|X D KL (P C|X Q C|U ) ,where (a) and (b) follow due to the Markov chain C --X --U; (c) follows from the definition of L VB s (P, Q) in equation 6.

c (x).The KL divergence between two Gaussian mixtures f an g can be approximated as followsD vKL (f g) := M m=1 ω m log m ∈M\m ω m exp (-D KL (f m f m )) C c=1 π c exp (-D KL (f m g c ))

values are taken from VaDE (Jiang et al., 2017) ‡ values are taken from DEC (Xie et al., 2016) Comparison of clustering accuracy of various algorithms.

Table content: <table><tr><td/><td colspan="3">MNIST STL-10 REUTERS10K</td></tr><tr><td>GMM</td><td>50.4</td><td>77.1</td><td>53.74</td></tr><tr><td>DEC</td><td>84.3  ‡</td><td>80.6  †</td><td>72.17  ‡</td></tr><tr><td>VaDE</td><td>94.5  †</td><td>84.3</td><td>79.8  †</td></tr><tr><td>VIB-GMM</td><td>96.2</td><td>91.6</td><td>80.4</td></tr><tr><td/><td/><td/><td/></tr></table>

-sE P U|X D KL (P C|X Q C|U ) = E P X E P U|X [log Q X|U ]s = E P X E P U|X [log Q X|U ]s

Table content: <table><tr><td/><td/><td/><td/><td/><td/><td>x</td><td colspan="2">p(x)</td><td>u</td><td>p(u|x) log</td><td>p(u|x) q(u)</td><td>du dx</td></tr><tr><td>-s</td><td>x</td><td>p(x)</td><td>u</td><td>p(u|x)</td><td>c</td><td colspan="4">p(c|x) log</td><td>p(c|x) q(c|u)</td><td>du dx</td></tr><tr><td>(a)</td><td/><td/><td/><td/><td/><td/><td>ux</td><td colspan="2">p(x)p(u|x) log</td><td>p(u|x) q(u)</td><td>du dx</td></tr><tr><td>-s</td><td/><td>ux c</td><td colspan="7">p(x)p(u|c, x)p(c|x) log</td><td>p(c|x) q(c|u)</td><td>du dx</td></tr><tr><td>= E P</td><td/><td/><td/><td/><td/><td/><td/><td/></tr></table>

Peer review (iclr2020):
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. The idea of using a latent mixture of Gaussian*s to variationally encode high-dimensional data is not new. The author(s) appropriately cite VaDE (Jiang, 2017) and DEC (Xie, 2016), *Deep Unsupervised Clustering with Gaussian Mixture Variational Autoencoders* (Dilokthanakul, 2017). However, the author(s) unfortunately did not mention *Approximate Inference for Deep Latent Gaussian Mixtures* (Nalisnick, 2016). Nalisnick et al. consider the exact generative process as the proposed by the author(s), but with the addition of a Dirichlet prior on the Categorical distribution. Nalisnick et al. fit their model with a VAE and circumvent Dirichlet intractability by positing a Kumarswamy stick-breaking variational posterior (the resulting distribution is on the simplex). They achieve 91.58% accuracy, but do so after fitting a KNN to the latent space. New techniques such as *Implicit Reparameterization Gradients* (Figurnov, 2018) and *Pathwise Derivatives Beyond the Reparameterization Trick* (Jankowiak, 2018) allow direct use of a Dirichlet posterior in VAEs. These methods are respectively implemented in the TensorFlow and PyTorch APIs. My point here is that there are many ways to fit this pre-existing model. From my review of these other works we have, for *best run* on MNIST, that author(s) > VaDE > Nalisnick > Dilokthanakul > DEC. Thus, the author(s) are SoTA for this particular generative process for their best run. It would be nice to see their standard deviation to gauge how statistically significant their results are. However, I am dubious of the SoTA claim as I detail later. The author(s) derivation was sound, but a bit confusing for me. In particular, I found keeping track of P*s and Q*s very burdensome after reading sections 2.1 and 2.2. In my experience, Q is typically reserved for a variational distribution that approximates some intractable P distribution, while P is used to describe the generative process (i.e. likelihood) and other exact distributions. Furthermore, one typically introduces the generative process first using P distributions. Once, I got past this confusion everything else made sense. I might suggest introducing the generative process first and with P distributions instead of Q*s. The variational model can follow with the Q distributions as the author(s) have it. Equations 4, 5, and 6 all exactly match the unsupervised information bottleneck objective (Alemi, 2017)--see appendix B. I am therefore confident in those equations. I carefully checked their derivation of the VADe comparison (equation 10 and appendix B). Their derivation is straightforward. The principle trick is using the MC assumption C->X->U for P distributions to claim p(u|x,c) = p(u|x). Equations 12-16 all follow naturally from equation 11. If equation 17 is indeed an approximation or bound approximation, I would suggest not using the equals sign. Instead, consider another appropriate operator or rename Dkl to indicate it is the approximated version (just as in equation 24). If the author(s) like my suggestion regarding generative vs variational nomenclature, I would also change the second sentence of page 6 to something like *We use our variational Qc|u distribution to compute assignments.* Thereafter, I would drop the star indicator for optimal parameters. These parameters are not necessarily optimal given the non-convexity of the DNN. Replace with, *after convergence, we ...* Similarly, drop optimal from line 2 of algorithm 1. Circling back to the experiments, the author(s) use reported values from DEC and VaDE. Those works compute cluster assignments using a KNN classifier on the latent space. This paper however uses the arg max of equation 19. I much prefer this article*s method, but for comparison purposes, the author(s) should similarly use a KNN classifier on their latent space to compute accuracy in the same manner. The use of KNN in these other works allows them to consider a number of latent clusters larger than the number of true classes (e.g. 20 clusters for MNIST). I like that the author(s) stick to 10 clusters for MNIST, but for comparison I would have liked to see a KNN generated accuracy alongside their equation 19 based accuracy. It looks like the author(s) implemented VADe for STL-10 based on table 1. If so, it seems they could easily implement equation 19 for their STL-10 value for VADe. If not, please correct this. Not to belabor further, but I would really like to see a table 2 that reports both equation 19 and KNN accuracies when available. Namely, the author(s) should report both values for their model and can leave equation 19 accuracies blank for reported values. Lastly, I always raise an eyebrow when I see tSNE latent space representations. STL-10 was used to generate figure 6, where one needs dim(u) >> 2. However, MNIST can be well reconstructed using just 2 latent dimensions. In this case, tSNE is unnecessary. The author(s) state *Colors are used to distinguish between clusters.* This statement is unclear as to whether the author(s) are using the class label or learned latent cluster. If it is the former, figure 6 makes sense in that it shows misclassifications (i.e. red dots in the green cluster). However, if it is the latter then I am concerned tSNE is doing something weird. To summarize, I enjoyed the paper. My only concern is novelty. Being the first to pair an existing model with an existing method, in my eyes, does not necessarily meet the ICLR bar. The author(s) seemingly achieve SoTA, but without KNN-based accuracies for their model it is hard to compare to cited works. Having these KNN results and error bars would strengthen their case substantially to: achieving SoTA for an existing model by being the first to pair it with an existing method.