Title: HIERARCHICALLY CLUSTERED REPRESENTATION LEARNING

Abstract: The joint optimization of representation learning and clustering in the embedding space has experienced a breakthrough in recent years. In spite of the advance, clustering with representation learning has been limited to flat-level categories, which oftentimes involves cohesive clustering with a focus on instance relations. To overcome the limitations of flat clustering, we introduce hierarchically clustered representation learning (HCRL), which simultaneously optimizes representation learning and hierarchical clustering in the embedding space. Specifically, we place a nonparametric Bayesian prior on embeddings to handle dynamic mixture hierarchies under the variational autoencoder framework, and to adopt the generative process of a hierarchical-versioned Gaussian mixture model. Compared with a few prior works focusing on unifying representation learning and hierarchical clustering, HCRL is the first model to consider a generation of deep embeddings from every component of the hierarchy, not just leaf components. This generation process enables more meaningful separations and mergers of clusters via branches in a hierarchy. In addition to obtaining hierarchically clustered embeddings, we can reconstruct data by the various abstraction levels, infer the intrinsic hierarchical structure, and learn the level-proportion features. We conducted evaluations with image and text domains, and our quantitative analyses showed competent likelihoods and the best accuracies compared with the baselines.

1. INTRODUCTION:

Clustering is one of the most traditional and frequently used machine learning tasks. Clustering models are designed to represent intrinsic data structures, such as latent Dirichlet allocation (Blei et al., 2003) . The recent development of representation learning has contributed to generalizing model feature engineering, which also enhances data representation (Bengio et al., 2013) . Therefore, representation learning has been merged into the clustering models, e.g., variational deep embedding (VaDE) (Jiang et al., 2017) . Besides merging representation learning and clustering, another critical line of research is structuring the clustering result, e.g., hierarchical clustering. This paper introduces a unified model enabling nonparametric Bayesian hierarchical clustering with neural-network-based representation learning.

Autoencoder (Rumelhart et al., 1985) is a typical neural network for unsupervised representation learning and achieves a non-linear mapping from a high-dimensional input space to a lowdimensional embedding space by minimizing reconstruction errors. To turn the low-dimensional embeddings into random variables, a variational autoencoder (VAE) (Kingma & Welling, 2014) places a Gaussian prior on the embeddings. The autoencoder, whether it is probabilistic or not, has a limitation in reflecting the intrinsic hierarchical structure of data. For instance, VAE assuming a single Gaussian prior needs to be expanded to suggest an elaborate clustering structure.

Due to the limitations of modeling the cluster structure with autoencoders, prior works combine the autoencoder and the clustering algorithm. While some early cases pipeline just two models, e.g., Huang et al. (2014) , a typical merging approach is to model an additional loss, such as a clustering loss, in the autoencoders (Xie et al., 2016; Guo et al., 2017; Yang et al., 2017; Nalisnick et al., 2016; Chu & Cai, 2017; Jiang et al., 2017) . These suggestions exhibit gains from unifying the encoding and the clustering, yet they remain at the parametric and flat-structured clustering. A more recent development releases the previous constraints by using the nonparametric Bayesian approach. Figure 1 : Example of hierarchically clustered embeddings on MNIST with three levels of hierarchy, the reconstructed digits from the hierarchical Gaussian mixture components, and the extracted level proportion features. We marked the mean of a Gaussian mixture component with the colored square, and the digit written inside the square refers to the unique index of the mixture component.

For example, the infinite mixture of VAEs (IMVAE) (Abbasnejad et al., 2017) explores the infinite space for VAE mixtures by looking for an adequate embedding space through sampling, such as the Chinese restaurant process (CRP). Whereas IMVAE remains at the flat-structured clustering, VAEnested CRP (VAE-nCRP) (Goyal et al., 2017) captures a more complex structure, i.e., a hierarchical structure of the data, by adopting the nested Chinese restaurant process (nCRP) prior (Griffiths et al., 2004) into the cluster assignment of the Gaussian mixture model.

This paper proposes hierarchically clustered representation learning (HCRL) that is a joint model of 1) nonparametric Bayesian hierarchical clustering, and 2) representation learning with neural networks. HCRL extends a previous work on merging flat clustering and representation learning, i.e., VaDE, by incorporating inter-cluster relation modelings. Unlike a previous work of VAE-nCRP, HCRL learns the full spectrum of hierarchical clusterings, such as the level assignment and the level proportion of generating a component hierarchy. These level assignments and proportions were not modeled in VAE-nCRP, so each data instance cannot be analyzed from the perspective of generalization and specialization in a hierarchy. On the contrary, by adding level assignment and proportion modeling, a data instance can be generated from an internal component of the hierarchy, which is limited to the leaf component in VAE-nCRP. Hierarchical mixture density estimation (Vasconcelos & Lippman, 1999) , where all internal and leaf components are directly modeled to generate data, is a flexible framework for hierarchical mixture modeling, such as hierarchical topic modeling (Mimno et al., 2007; Griffiths et al., 2004) , with regard to the learning of the internal components.

Specifically, HCRL jointly optimizes soft-divisive hierarchical clustering in an embedding space from VAE via two mechanisms. First, HCRL includes a hierarchical-versioned Gaussian mixture model (HGMM) with a mixture of hierarchically organized Gaussian distributions. Then, HCRL sets the prior of embeddings by adopting the generative processes of HGMM. Second, to handle a dynamic hierarchy structure dealing with the clusters of unequal sizes, we explore the infinite hierarchy space by exploiting an nCRP prior. These mechanisms are fused as a unified objective function; this is done rather than concatenating the two distinct models of clustering and autoencoding. The quantitative evaluations focus on density estimation quality and hierarchical clustering accuracy, which shows that HCRL has competent likelihoods and the best accuracies compared with the baselines. When we observe our results qualitatively, we visualize 1) the hierarchical clusterings, 2) the embeddings under the hierarchy modeling, and 3) the reconstructed images from each Gaussian mixture component, as shown in Figure 1 . These experiments were conducted by crossing the data domains of texts and images, so our benchmark datasets include MNIST, CIFAR-100, RCV1 v2, and 20Newsgroups.

2. PRELIMINARIES:



2.1. VARIATIONAL DEEP EMBEDDING:

Figure 2 presents a graphical representation and a neural architecture of VaDE (Jiang et al., 2017) .

The model parameters of κ, µ 1:K , and σ 2 1:K , which are a proportion, means, and covariances of et al., 1999) . VaDE uses the Gaussian mixture model (GMM) as the prior, whereas VAE assumes a single standard Gaussian distribution on embeddings. Following the generative process of GMM, VaDE assumes that 1) the embedding draws a cluster assignment, and 2) the embedding is generated from the selected Gaussian mixture component.

VaDE uses an amortized inference as VAE, with a generative and inference networks; L(x) in Equation 1 denotes the evidence lower bound (ELBO), which is the lower bound on the log likelihood. It should be noted that VaDE merges the ELBO of VAE with the likelihood of GMM.

log p(x) ≥ L(x) = E q log p(c, z, x) q(c, z|x) = E q log K c=1 κ c N (z|µ c , σ 2 c I J ) p(c|z)N (z| µ, σ 2 I J ) + log p(x|z) (1)

2.2 VARIATIONAL AUTOENCODER NESTED CHINESE RESTAURANT PROCESS VAE-nCRP uses the nonparametric Bayesian prior for learning tree-based hierarchies, the nCRP (Griffiths et al., 2004) , so the representation could be hierarchically organized. The nCRP prior defines the distributions over children components for each parent component, recursively in a topdown way. The variational inference of the nCRP can be formalized by the nested stick-breaking construction (Wang & Blei, 2009) , which is also kept in the VAE setting. The distribution over paths on the hierarchy is defined as being proportional to the product of weights corresponding to the nodes lying in each path. The weight, π i , for the i-th node follows the Griffiths-Engen-McCloskey (GEM) distribution (Pitman et al., 2002) , where π i is constructed as

π i = v i i-1 j=1 (1 -v j ), v i ∼ Beta(1, γ

) by a stick-breaking process. Since the nCRP provides the ELBO with the nested stick-breaking process, VAE-nCRP has a unified ELBO of VAE and the nCRP in Equation 2.

L(x) = E q log p(v) q(v|x) + log p(α par(p) |α * )p(α p |α par(p) , σ 2 N ) q(α p , α par(p) |x) (3.1) p(ζ|v) q(ζ|x) p(z|α p , ζ, σ 2 D ) q(z|x) (3.2) + log p(x|z) (2)

Given the ELBO of VAE-nCRP, we recognized a number of potential improvements. First, term (3.1) is for modeling the hierarchical relationship among clusters, i.e., each child is generated from its parent. VAE-nCRP trade-off is the direct dependency modeling among clusters against the meanfield variational approach. This modeling may reveal that the higher clusters in the hierarchy are more difficult to train. Second, in term (3.2), leaf mixture components generate embeddings, which implies that only leaf clusters have direct summarization ability for sub-populations. Additionally, in term (3.2), variance parameter σ 2 D is modeled as the hyperparameter shared by all clusters. In other words, only with J-dimensional parameters, α, for the leaf mixture components, the local density modeling without variance parameters has a critical disadvantage. For all of these weaknesses, we were able to compensate with the level proportion modeling and HGMM prior. The level assignment generated from the level proportion allows a data instance to select among all mixture components. We do not need direct dependency modeling between the parents and their children because all internal mixture components also generate embeddings.

3. METHODOLOGY:



3.1. GENERATIVE PROCESS:

The generative process of HCRL resembles the generative process of hierarchical clusterings, such as the hierarchical latent Dirichlet allocation (Griffiths et al., 2004) . In detail, the generative process departs from selecting a path ζ, from the nCRP prior (phase 1). Then, we sample a level proportion (phase 2) and a level, l (phase 3), from the sampled level proportion to find the mixture component in the path, and this component of ζ l provides the Gaussian distribution for the latent representation (phase 4). Finally, the latent representation is exploited to generate an observed datapoint (phase 5). The below formulas are the generative process with its density functions. In addition, Figure 3 illustrates a graphical representation corresponding to the described generative process. The generative process also presents our formalization of corresponding prior distributions, denoted as p(•), and variational distributions, denoted as q(•), by generation phases. The variational distributions are used in our inference methods called mean-field variational inference (MFVI) (Jordan et al., 1999) as detailed in Section 3.3.

1. Choose a path ζ ∼ nCRP(ζ|γ) • p(ζ) = L l=1 π 1,ζ2,...,ζ l where π 1,ζ2,...,ζ l = l l =1 {v 1,ζ2,...,ζ l ( ζ l -1 j=1 (1 -v 1,ζ2,...,j ))}, q(ζ|x) ∝ S ζ ζ∈child(ζ) S ζ 2. Choose a level proportion η ∼ Dirichlet(η|α) • p(η) = Dir(η|α), q φη (η|x) = Dirichlet(η| α) ≈ LogisticNormal(η| µ η , σ 2 η I L ) where [ µ η ; log σ 2 η ] = g φη (x), α l = 1 σ 2 η l (1 -2 L + e -µη l L 2 l e -µη l ) 3. Choose a level l ∼ Multinomial(l|η) • p(l) = Multinomial(η), q(l|x) = Multinomial(l|ω)

where

ω l ∝ exp ψ( α l ) -ψ( α 0 ) + ζ S ζ J j=1 -1 2 log(2πσ 2 ζ l ,j ) - σ 2 z j 2σ 2 ζ l ,j - ( µz j -µ ζ l ,j ) 2 2σ 2 ζ l ,j 4. Choose a latent representation z ∼ N (z|µ ζ l , σ 2 ζ l I J ) • p(z) = ζ,l p(ζ|γ) • η l • N (z|µ ζ l , σ 2 ζ l I J ), q φz (z|x) = N (z| µ z , σ 2 z I J ) where [ µ z ; log σ 2 z ] = g φz (x) 5. Choose an observed datapoint x ∼ N x|µ x , σ 2 x I D where [µ x ; log σ 2 x ] = f θ (z) 2 3.2 NEURAL ARCHITECTURE

The neural architecture of HCRL consists of two probabilistic encoders on z and η, and one probabilistic decoder on z as shown in the right part of Figure 3 . This unbalanced architecture originates from our modeling assumption of p(x|z), not p(x|z, η). The reconstruction design of x depending on the two stochastic variables of z and η may lead to a large variance of the reconstruction on x.

Additionally, we cannot guarantee that both z and η contribute to the the reconstruction on x (Chen et al., 2016) . Although the decoding structure of η is not included explicitly in the neural network architecture of HCRL, we provide the formalization of p(η|z) in Table 1 according to our generative assumptions. We call this reconstruction process, which is inherently a generative process of the traditional probabilistic graphical model (PGM), PGM reconstruction (see the decoding neural network part of Figure 3 ).

Table 1 : Encoding and decoding structure on z and η in HCRL. (s) indicates the s-th sample.

Encoding Decoding z z ∼ q φz (z|x), z (s) = g φz ( (s) , x) x ∼ p θ (x|z), x (s) = f θ (z (s) ) η η ∼ q φη (η|x), η (s) = g φη ( (s) , x) p(x|η) ∝ v,z ζ,l p(x|z)p(z|ζ, l)p(l|η)p(ζ|v)p(v)

3.3. MEAN-FIELD VARIATIONAL INFERENCE:

The formal specification can be a factorized probabilistic model as Equation 3, where Φ = {v, ζ, η, l, z} denotes the set of latent variables, and M T denotes the set of all nodes in tree T .

p(Φ, x) = j / ∈M T p(v j |γ) i∈M T p(v i |γ) N n=1 p(ζ n |v)p(η n |α)p(l n |η n )p(z n |ζ n , l n )p θ (x n |z n ) (3)

The proportion and assignment on the mixture components for the n-th data instance are modeled by ζ n as a path assignment; η n as a level proportion; and l n as a level assignment. v is a Beta draw used in the stick-breaking construction. The latent variables are inferred through MFVI, and therefore we assume the variational distributions are as Equation 4:

q(Φ|x) = j / ∈M T p(v j |γ) i∈M T q(v i |a i , b i ) N n=1 q(ζ n |x n )q φη (η n |x n )q(l n |ω n , x n )q φz (z n |x n ) (4)

where q φη (η n |x n ) and q φz (z n |x n ) should be noted because these two variational distributions follow the amortized inference of VAE. q(ζ|x) ∝ S ζ ζ∈child(ζ) S ζ is the variational distribution over path ζ, where child(ζ) means the set of all full paths that are not in T but include ζ as a sub path. Because we specified both generative and variational distributions, we define the ELBO of HCRL, L = E q log p(Φ,x) q(Φ|x) , in Equation 5. Appendix F enumerates the full derivation in detail. We report that the Laplace approximation with the logistic normal distribution is applied to model the prior, α, of the level proportion, η. We choose a conjugate prior of a multinomial, so p(η n |α) follows the Dirichlet distribution. To configure the inference network on the Dirichlet prior, the Laplace approximation is used (MacKay, 1998; Srivastava & Sutton, 2017; Hennig et al., 2012) .

L(x) = E q log p(v) q(v|x) + log p(η) q(η|x) + log ζ,l p(ζ|v) q(ζ|x) p(l|η) q(l|x) p(z|µ ζ l , σ 2 ζ l ) q(z|x) + log p(x|z) (5)

3.4. TRAINING ALGORITHM OF CLUSTERING HIERARCHY:

This model is formalized according to the stick-breaking process scheme. Unlike the CRP, the stickbreaking process does not represent the direct sampling of the mixture component at the data instance level. Therefore, it is necessary to devise a heuristic algorithm for operations, such as GROW, PRUNE, and MERGE, to refine the hierarchy structure. Appendix C provides details about each operation, together with the overall training algorithm of HCRL. In the below description, an inner path and a full path refer to the path ending with an internal node and a leaf node, respectively.

• GROW expands the hierarchy by creating a new branch under the heavily weighted internal node.

Compared with the work of Wang & Blei (2009) , we modified GROW to first sample a path, ζ * , proportional to n q(ζ n = ζ * ), and then to grow the path if the sampled path is an inner path.

• PRUNE cuts a randomly sampled minor full path, ζ * , satisfying n q(ζn=ζ * ) n,ζ q(ζn=ζ) < δ, where δ is the pre-defined threshold. If the removed leaf node of the full path is the last child of the parent node, we also recursively remove the parent node.

• MERGE combines two full paths, ζ (i) and ζ (j) , with similar posterior probabilities, measured by

J(ζ (i) , ζ (j) ) = q i q T j /|q i ||q j |, where q i = [q(ζ 1 = ζ (i) ), • • • , q(ζ N = ζ (i) )].

4. EXPERIMENTS:



4.1. DATASETS AND BASELINES:

Datasets: We used various hierarchically organized benchmark datasets as well as MNIST.

• MNIST (LeCun et al., 1998) : 28x28x1 handwritten image data, with 60,000 train images and 10,000 test images. We reshaped the data to 784-d in one dimension. • CIFAR-100 (Krizhevsky & Hinton, 2009) : 32x32x3 colored images with 20 coarse and 100 fine classes. We used 3,072-d flattened data with 50,000 training and 10,000 testing. • RCV1 v2 (Lewis et al., 2004) : The preprocessed text of the Reuters Corpus Volume. We preprocessed the text by selecting the top 2,000 tf-idf words. We used the hierarchical labels up to the 4-level, and the multi-labeled documents were removed. The final preprocessed corpus consists of 11,370 training and 10,000 testing documents randomly sampled from the original test corpus. • 20Newsgroups (Lang, 1995) : The benchmark text data extracted from 20 newsgroups, consisting 11,314 training and 7,532 testing documents. We also labeled by 4-level following the annotated hierarchical structure. We preprocessed the data through the same process as that of RCV1 v2.

Baselines: We completed our evaluation in two aspects: 1) optimizing the density estimation, and 2) clustering the hierarchical categories. First, we evaluated HCRL from the density estimation perspective by comparing it with diverse flat clustered representation learning models, and VAE-nCRP. Second, we tested HCRL from the accuracy perspective by comparing it with multiple divisive hierarchical clusterings. The below is the list of baselines. We also added the two-stage pipeline approaches, where we trained features from VaDE first and then applied the hierarchical clusterings.

We reused the open source codes3 provided by the authors for several baselines, such as IDEC, DCN, VAE-nCRP, and SSC-OMP. (Yang et al., 2017) : optimizes the K-means-related cost defined on the embedding space. We used the open source code provided by the authors. 5. Infinite Mixture of Variational Autoencoders (IMVAE) (Abbasnejad et al., 2017) : searches for the infinite embedding space by using a Bayesian nonparametric prior. 6. Variational Autoencoder -nested Chinese Restaurant Process (VAE-nCRP) (Goyal et al., 2017) : We used the open source code provided by the authors. 7. Hierarchical K-means (HKM) (Nister & Stewenius, 2006) : performs K-means (Lloyd, 1982) recursive in a top-down way. 8. Mixture of Hierarchical Gaussians (MOHG) (Vasconcelos & Lippman, 1999 ): infers the level-specific mixture of Gaussians. 9. Recursive Gaussian Mixture Model (RGMM): runs GMM recursively in a top-down manner. 10. Recursive Scalable Sparse Subspace Clustering by Orthogonal Matching Pursuit (RSS-COMP): performs SSC-OMP (You et al., 2016) recursively for hierarchical clustering. SSC-OMP is a well-known methods for image clustering, and we used the open source code.

4.2. QUANTITATIVE ANALYSIS:

We used two measures to evaluate the learned representations in terms of the density estimations: 1) negative log likelihood (NLL), and 2) reconstruction errors (REs). Autoencoder models, such as IDEC and DCN, were tested only for the REs. The NLL is estimated with 100 samples. VaDE generally performed better than VAE did, whereas other flat clustered representation learning models tended to be slightly different for each dataset. HCRL showed overall competent performance and better results with a deeper hierarchy of level four than of level three, which implies that capturing the deeper hierarchical structure is likely to be useful for the density estimation.

Additionally, we evaluated hierarchical clustering accuracies by following Xie et al. (2016) , except for MNIST that is flat structured. Table 3 points out that HCRL has significantly better microaveraged F-scores compared with every baseline. HCRL is able to reproduce the ground truth hierarchical structure of the data, and this trend is consistent when HCRL compared with the pipelined model, such as VaDE with a clustering model. The result of the comparisons with the clustering models, such as HKM, MOHG, RGMM, and RSSCOMP, is interesting because it experimentally proves that the joint optimization of hierarchical clustering in the embedding space improves hierarchical clustering accuracies. HCRL also presented better hierarchical accuracies than VAE-nCRP. We conjecture the reasons for the modeling aspect of VAE-nCRP: 1) the simplified prior modeling on the variance of the mixture component as just constants, and 2) the non-flexible learning of the internal components. 

4.3. QUALITATIVE ANALYSIS:



MNIST::

In Figure 1 , the digits {4, 7, 9} and the digits {3, 8} are grouped together with a clear hierarchy, which was consistent between HCRL and VaDE. Also, some digits {0, 4, 2} in a round form are grouped, together, in HCRL. In addition, among the reconstructed digits from the hierarchical mixture components, the digits generated from the root have blended shapes from 0 to 9, which is natural considering the root position.

CIFAR-100: Figure 4 shows the hierarchical clustering results on CIFAR-100. Given that there were no semantic inputs from the data, the color was dominantly reflected in the clustering criteria. However, if one observes the second hierarchy, the scene images of the same sub-hierarchy are semantically consistent, although the background colors are slightly different. 20Newsgroups: Figure 6 shows the example sub-hierarchies on 20Newsgroups. We enumerated topic words from documents with top-five likelihoods for each cluster, and we filtered the words by tf-idf values. We observe relatively more general contents in the internal clusters than in the leaf clusters of each internal cluster. In this paper, we have introduced a hierarchically clustered representation learning framework for the hierarchical mixture density estimation on deep embeddings. HCRL aims at encoding the relations among clusters as well as among instances to preserve the internal hierarchical structure of data. The main differentiated features of HCRL are 1) the crucial assumption regarding the internal mixture components for having the ability to generate data directly, and 2) the unbalanced autoencoding neural architecture for the level proportion modeling as the encoding structure, and the probabilistic model as the decoding structure. From the modeling and the evaluation, we found that HCRL enables the improvements due to the high flexibility modeling compared with the baselines.

A SYNTHETIC DEMO:

We created a synthetic dataset that has a hierarchical structure and is sampled from the 50dimensional Gaussian distributions, presented in Figure 7 . The hierarchy, which has a branch factor of two and a depth of four, has a total of eight leaf clusters. Figure 7a shows the raw synthetic dataset in the input space of R 50 , and after running HCRL, we plot the hierarchically clustered embeddings in the latent space in Figure 7b . In addition to the embeddings, we also present a confidence ellipse with dashed lines for each learned Gaussian mixture component. Because the root component is involved in generating all of the data, it forms a large ellipse, while the leaf component summarizes the local density, so the small ellipse is learned. 

B EXPERIMENTAL SETTINGS:

We conducted experiments for all autoencoder-based models with a neural architecture whose encoder network was set as fully connected layers with dimensions D-2000-2000-500-J for z, and D-10-10-L for η, and the decoder network is a mirror of the encoder network for z. The hyperparameters of HCRL given by users, γ and α, was set to 1.0, and a vector of all entries 1 sized of L, respectively. We used the Adam optimizer (Kingma & Ba, 2014) with an init learning rate of 0.001 for MNIST dataset and 0.0001 for other datasets. Meanwhile, VAE-nCRP is targeted for grouped data. For experiments with our non-grouped datasets, we treated the group instance as a group instance having a single data instance. For parametric hierarchical clustering models, we gave the branch factor as the input parameter, [1, 20, 5], [1, 4, 7, 9] , and [1, 6, 4, 3], for CIFAR-100, RCV1 v2, and 20Newsgroups, respectively. For VaDE, we set the number of clusters to the number of leaf clusters; 100 for CIFAR-100, 252 for RCV1 v2, and 72 for 20Newsgroups.

C ALGORITHMS C.1 TRAINING ALGORITHM:

Algorithm 1 summarizes the overall algorithm for HCRL. The tree-based hierarhcy T is defined as (N, P), where N and P denote a set of nodes and paths, respectively. We refer to the node at level l lying on path ζ, as N(ζ 1:l ) ∈ N. The defined paths, P, consist of full paths (ending at a leaf node), P full , and inner paths (ending at an internal node), P inner , as a union set.

Algorithm 1 selects an operation out of three operations: GROW, PRUNE, and MERGE. The GROW algorithm is executed for every specific iteration period, t G . After ellapsing t b iterations since performing the GROW operation, we begin to check whether the PRUNE or MERGE operation should be performed. We prioritize the PRUNE operation first, and if the condition of performing PRUNE is not satisfied, we check for the MERGE operation next. After performing any operation, we initialize n b to 0, which is for locking the changed hierarchy during minimum t b iterations to be fitted to the training data. φ z , φ η , θ ← Update the network weight parameters using gradients ∇ φz,φη,θ L(x)

T (E) , φ z , φ η , θ, ω, {a i , b i , µ i , σ 2 i } i∈M T (E) 1: µ ζ 1:L , σ 2

6: {a i , b i , µ i , σ 2 i } i∈M T (e-1) ← Update node-specific params. using gradients ∇ a,b,µ,σ 2 L(x) 7:

Update other variational parameters using gradients ∇L(x)

8:

if mod(e, t grow ) = 0 then 9:

T (e) , Q ← GROW(T (e-1) , Q, s grow , Q max ) // See Algorithm 2 10:

end if 11:

if T (e) = T (e-1) and n lock ≥ t lock then 12:

T (e) , Q ← PRUNE(T (e-1) , Q, δ prune ) // See Algorithm 3 13:

if T (e) = T (e-1) then T (e) , Q ← MERGE(T (e-1) , Q, δ merge , Q max ) // See Algorithm 4 14:

end if 15:

if T (e) = T (e-1) then n lock ← 0 else n lock ← n lock + 1 16: end for Q ← φ // Temporary set of changed paths in this epoch 5: )) where ∼ N (0, n g I J )

if ζ * ∈ P inner and ζ * / ∈ Q s.t. Q ∈ Q then 6: l 0 ← |ζ * | 7: for l = l 0 , • • • , L -1 do 8: j 0 ← Maximum index

11:

Q ← Q ∪ {ζ * 1:l +1 } 12:

if l < L -1 then T ← (N, P) The PRUNE operation cuts a minor path, which is sampled according to

N n=1 q(ζ n = ζ) among the full paths satisfying N n=1 q(ζ n = ζ) < δ,

where δ is the pre-defined threshold parameter. If the removed leaf node of the full path is the last child of the parent node, we also recursively remove the parent node as shown in the upper case of Figure 10 .

Algorithm 3 PRUNE Operation 1: function PRUNE(T, Q, δ prune ) 2: J(ζ) ← N n=1 q(ζ n = ζ), ζ ∈ P // Calculate the measure 3: Ω ← {ζ | ζ ∈ P full , J(ζ) ζ J(ζ ) < δ prune } 4: Randomly sample a full path ζ * ∼ Ω 5: if |P full | > 1 and ζ * / ∈ Q s.t. Q ∈ Q then 6: P full ← P full \{ζ * 1:L } 7: for l = L -1, • • • , 1 do 8: N(ζ * 1:l +1 ) ← φ 9:

if l < L -1 then The MERGE operation combines two full paths with similar posterior probabilities, measured by Algorithm 4 MERGE Operation

J(ζ (i) , ζ (j) ) = q i q T j /|q i ||q j |, where q i = [q(ζ 1 = ζ (i) ), ..., q(ζ N = ζ (i) )]. We

1: function MERGE(T, Q, δ merge , Q max ) 2: J(ζ (i) , ζ (j) ) ← qiq T j |qi||qj | s.t. q i = [q(ζ 1 = ζ (i) ), ..., q(ζ N = ζ (i) )] // Calculate the measure 3: Ω ← {(ζ (i) , ζ (j) ) | J(ζ (i) , ζ (j) ) ≥ δ merge , {ζ (i) , ζ (j) } ⊂ P full } 4:

Randomly sample a pair of paths (ζ

EQUATION

8:

EQUATION

EQUATION

EQUATION

w (1) +w (2) 13: j 0 ← Maximum index for the child node whose parent path is ζ * 1:l 14:

ζ * 1:l +1 ← ζ * 1:l , j 0 + 1 15: N(ζ * 1:l +1 ) ← N (µ * , diag(σ 2 * )) 16: N(ζ (1) 1:l +1 ) ← φ , N(ζ (2) 1:l +1 ) ← φ 17:

if l < L -1 then 18:

EQUATION

1:l +1 } 19:

else 20:

EQUATION

1:L } 21: end if 22: Q ← Q ∪ {ζ * 1:l +1 } 23:

end for 24:

end if

25: enqueue Q to Q 26: while Q max < |Q| do dequeue Q 27:

P ← P full ∪ P inner 28:

T ← (N, P)

29::

return T, Q 30: end function

D NOTATIONS:

The following Table 4 lists the notations used throughout this paper. The input / latent dimensionality g * (x)

A encoder network parametrized by * , whose input is x f θ (z)

A decoder network parametrized by θ, whose input is z θ

The variational parameters and weights of the decoder network f θ μz, σ2

z

The variational mean and variance for Gaussian distribution q φz (z|x) µx, σ 2

x The prior parameters, mean and variance, for Gaussian distribution p θ (x|z)

VaDE & VAE-nCRP φ

The variational parameters and weights of the encoder network g φ μ, σ2

The variational mean and variance for Gaussian distribution q φ (z|x)

VaDE & HCRL N The number of datapoints xn=1,••• ,N n-th observed datapoint zn=1,••• ,N n-th latent representation corresponding to xn VAE-nCRP & HCRL

L:

The height of the tree-based hierarchy

VaDE:



K:

The number of (finite

) clusters cn=1,••• ,N

The cluster assignment of zn, ∈ {1, ..., K} κ

The prior parameter for multinomial distribution p(c) µc, σ 2 c The prior parameters, mean and variance, for Gaussian distribution of c-th cluster, p(z)

VAE-nCRP:



M:

The number of sequences The variational parameters, for Beta distribution q(vmp|xm) ζmn

Nm=1,••• ,M The number datapoints in m-th sequence xm,n=1,••• ,N n-th

The path assignment of zmn S * mn The variational parameter for multinomial distribution q(ζmn|xmn) α par(p)

The J-dimensional parameter vector for the parent node of p α *

The prior parameter for Gaussian distribution p(αp) for the root node µ par(p) , σ 2 par(p)

The variational mean and variance for Gaussian distribution q(α par(p) |x) αp

The J-dimensional parameter vector for node p σ 2

N:

The prior parameter, variance, for Gaussian distribution p(αp|α par(p) ) µp, σ 2 p The variational mean and variance for Gaussian distribution q(αp|x) σ 2

D:

The prior parameter, variance, for Gaussian distribution p(zmn|ζmn, αp)

HCRL φz:

The variational parameters and weights of the encoder network g φz φη

The variational parameters and weights of the encoder network g φη μz, σ2

z

The variational mean and variance for Gaussian distribution q φz (z|x) μη, σ2

η

The variational mean and variance for logistic normal distribution q φη (η|x) α

The variational parameter for Dirichlet distribution q φη (η|x) vi

The Beta draws for the tree-based stick-breaking construction of node i γ

The prior parameter for Beta distribution p(vi) ai, bi

The variational parameters, for Beta distribution q(vi|x) ζn

The path assignment of zn Sn

The variational parameter for multinomial distribution q(ζn|xn) ηn

The level proportion of zn α

The prior parameter for Dirichlet distribution p(ηn) ln

The level assignment of zn, ∈ {1, ..., L} ωn

The variational parameter for multinomial distribution q(ln|xn) µi, σ 2 i The prior parameters, mean and variance, for Gaussian distribution of node i, p(zn|ζn, ηn)

E GENERATIVE AND INFERENCE MODEL FOR HCRL:

HCRL assumes the generative process as described in Section 3.1. Section E.1 describes the joint probability distribution, and Section E.2 presents the corresponding variational distributions. We adopt the much notation-related conventions from Wang & Blei (2009) , especially on paths.

E.1 GENERATIVE MODEL p θ (v, ζ, η, l, z, x) = p(v|γ) N n=1 p(ζ n |v)p(η n |α)p(l n |η n )p(z n |ζ n , l n , µ 1:∞ , σ 2 1:∞ )p θ (x n |z n ) = j / ∈M T p(v j |γ) i∈M T p(v i |γ) N n=1 p(ζ n |v)p(η n |α)p(l n |η n )p(z n |ζ n , l n )p θ (x n |z n ) • M T : Set of all nodes in truncated tree T • For j / ∈ M T , p(v j |γ) = Beta(v j |1, γ) • For i ∈ M T , p(v i |γ) = Beta(v i |1, γ) • p(ζ n = [1, ζ 2 , ..., ζ L ]|v) p(ζ n = [1, ζ 2 , ..., ζ L ]|v) = L l=1 π 1,ζ2,...,ζ l = L l=1 π 1,ζ2,...,ζ l-1 v 1,ζ2,...,ζ l ζ l -1 j=1 (1 -v 1,ζ2,...,j ) = L l=1 l l =1 v 1,ζ2,...,ζ l ζ l -1 j=1 (1 -v 1,ζ2,...,j ) -π 1,ζ2,...,ζ l = l l =1 {v 1,ζ2,...,ζ l ( ζ l -1 j=1 (1 -v 1,ζ2,...,j ))} • p(η n |α) = Dirichlet(η n |α) • p(l n |η n ) = Multinomial(η n ) • p(z n |ζ n = ζ, l n = l) = N (z n |µ ζ l , σ 2 ζ l I J ) • p θ (x n |z n ) :

Probabilistic decoding of x n parametrized by θ, whose input is z n • Tree-based stick-breaking construction -We will denote all Beta draws as v, each of which is an independent draw from Beta(v|1, γ)

(except for root v 1 = 1) * v i ∼ Beta(v i |1, γ) -The root nodes stick length: π 1 = v 1 ≡ 1 -Stick length at second level: π 1i = π 1 v 1i i-1 j=1 (1 -v 1j ), ∞ i=1 π 1i = π 1 = 1 -For the segment π 1k , the stick lengths of its children are π 1ki = π 1k v 1ki i-1 j=1 (1 -v 1kj ), for i = 1, 2, ..., ∞, ∞ i=1 π 1ki = π 1k E.2 INFERENCE MODEL

As VAE, we infer the random variables via the mean-field approximation, where the variational distribution, q φη,φz (v, ζ, η, l, z|x) , approximates the intractable posterior. We model the variational distributions as follows:

q φη,φz (v, ζ, η, l, z|x) = q(v|a, b, x) N n=1 q(ζ n |x n )q φη (η n |x n )q(l n |ω n , x n )q φz (z n |x n ) = j / ∈M T p(v j |γ) i∈M T q(v i |a i , b i ) N n=1 q(ζ n |x n )q φη (η n |x n )q(l n |ω n , x n )q φz (z n |x n ) • For j / ∈ M T , p(v j |γ) = Beta(v j |1, γ) • For i ∈ M T , q(v i |a i , b i ) ∝ v ai-1 i (1 -v i ) bi-1 = Beta(v i |a i , b i ) -a i = 1 + (L -l i + 1) N n=1 ζ l 0 +1 ,...,ζ L q(ζ n = [1, ζ 2 , ..., ζ l0 , ζ l0+1 , ..., ζ L ]) -b i = γ + (L -l i + 1) N n=1 j,ζ l 0 +1 ,...,ζ L :j>ζ l 0 q(ζ n = [1, ζ 2 , ..., ζ l0-1 , j, ζ l0+1 , ..., ζ L ]) * l i : The level of the mixture component i • q(ζ n |x n ) ∝ S nζ ζ∈child(ζ) S nζ

ζ: a path in the truncated tree T , either an inner path (a path ending at an internal node) or a full path (a path ending at a leaf node)

child(ζ): the set of all full paths that are not in T but include ζ as a sub path * As a special case, if ζ is a full path, child(ζ) just contains itself -In the case of a full path, 

EQUATION

S nζ = ζ∈child(ζ) S nζ = ζ∈child(ζ) exp E q L l=1 l l =1 log v 1,ζ2,...,ζ l + ζ l -1 j=1 log(1 -v 1,ζ2,...,j ) + log N (z n |µ ζ nl , σ 2 ζ nl I J ) = exp (E q [ L l=l0+1 log N (z n |µ ζ nl , σ 2 ζ nl I J )] + L l=l0+1 log(L -l + 1) + (L -l 0 )(ψ(1) -ψ(1 + γ)))} (1 -exp{ψ(γ) -ψ(1 + γ)}) L-l0 exp l0 l=1 (L -l + 1) E q [log v 1,ζ2,...,ζ l ] + ζ l -1 j=1 E q [log(1 -v 1,ζ2,...,j )] + log N (z n |µ ζ nl , σ 2 ζ nl I J ) exp j0 j=1 E q (L -l 0 ) log(1 -v 1,ζ2,...,ζ l 0 ,j ) -q φη (η n |x n ) = Dirichlet(η n | α n ) * α nl = 1 σ 2 η nl (1 -2 L + e -µη nl L 2 L l =1 e -µη nl ) -q(l n |ω n , x n ) = Multinomial(l n |ω n ) * ω nl ∝ exp{ψ( α nl ) -ψ( α n0 ) + ζ S nζ ( J j=1 -1 2 log(2πσ 2 ζ nl ,j ) - σ 2 z nj 2σ 2 ζ nl ,j - ( µz nj -µ ζ nl ,j ) 2 2σ 2 ζ nl ,j )} • α n0 = L i=1 α ni * Derivation for ω nl L ω nl = l ω nl (ψ( α nl ) -ψ( α n0 )) - l ω nl log ω nl + ζ S nζ l ω nl J j=1 - 1 2 log(2πσ 2 ζ nl ,j ) - σ 2 znj 2σ 2 ζ nl ,j - ( µ znj -µ ζ nl ,j ) 2 2σ 2 ζ nl ,j = ω nl (ψ( α nl ) -ψ( α n0 )) + ζ S nζ ω nl J j=1 - 1 2 log(2πσ 2 ζ nl ,j ) - σ 2 znj 2σ 2 ζ nl ,j - ( µ znj -µ ζ nl ,j ) 2 2σ 2 ζ nl ,j -ω nl log ω nl + λ( l ω nl -1) ∂L ω nl ∂ω nl = (ψ( α nl ) -ψ( α n0 )) + ζ S nζ J j=1 - 1 2 log(2πσ 2 ζ nl ,j ) - σ 2 znj 2σ 2 ζ nl ,j - ( µ znj -µ ζ nl ,j ) 2 2σ 2 ζ nl ,j -(log ω nl + 1) + λ ω nl = exp ψ( α nl ) -ψ( α n0 ) + ζ S nζ J j=1 - 1 2 log(2πσ 2 ζ nl ,j ) - σ 2 znj 2σ 2 ζ nl ,j - ( µ znj -µ ζ nl ,j ) 2 2σ 2 ζ nl ,j -1 + λ ω nl ∝ exp ψ( α nl ) -ψ( α n0 ) + ζ S nζ J j=1 - 1 2 log(2πσ 2 ζ nl ,j ) - σ 2 znj 2σ 2 ζ nl ,j - ( µ znj -µ ζ nl ,j ) 2 2σ 2 ζ nl ,j -q φz (z n |x n ) = N (z n | µ zn , σ 2 zn I J )

F EVIDENCE LOWER BOUND:

In this section, we present the detailed derivation of the ELBO in Equation 6, which is the objective function for learning HCRL.

EQUATION

F.1 DETAILED DERIVATION FOR ELBO

The followings are additional notations used for the detailed derivation:

• ψ : The digamma function (α i -1)(ψ( α ni ) -ψ( α n0 ))

• α n0 = L i=1 α ni , α 0 = L i=1 α i E q [log p(v i |γ)] = v z η l ζ j / ∈M T p(v j = v ) i∈M T q(v i = v ) N n=1 q(ζ n = ζ )q(η n = η ) q(l n = l )q(z n = z ) log p(v i = v )dη dz dv = v j / ∈M T p(v j = v ) i∈M T q(v i = v ) log p(v i = v )dv = v q(v i = v )

E q [log p(l n |η n )] = v z η l ζ j / ∈M T p(v j = v ) i∈M T q(v i = v )q(ζ n = ζ)q(η n = η )

q(l n = l )q(z n = z ) log p(l n = l |η n )dηdz dv x nd ) -

(x nd -µ (s)

x nd

) 2 2σ (s) 2 x nd if x n is real-valued

E q [log q(v i |a i , b i )] = v z η l ζ j / ∈M T p(v j = v ) i∈M T q(v i = v )q(ζ n = ζ)q(η n = η )

q(l n = l )q(z n = z) log q(v i = v )dη dz dv

= v j / ∈M T p(v j = v ) i∈M T q(v i = v ) log q(v i = v )dv = v q(v i = v ) log q(v i = v )dv = v Beta(v |a i , b i ) • log Beta(v |a i , b i )dv = log Γ(a i + b i ) -log Γ(a i ) -log Γ(b i ) + (a i -1)ψ(a i ) + (b i -1)ψ(b i ) +(-a i -b i + 2)ψ(a i + b i ) E q [log q(ζ n |x n )] = v z η l ζ j / ∈M T p(v j = v ) i∈M T q(v i = v )q(ζ n = ζ)q(η n = η )

q(l n = l )q(z n = z) log q(ζ n = ζ)dη dz dv

= ζ q(ζ n = ζ) log q(ζ n = ζ) = ζ S nζ ζ S nζ log S nζ ζ S nζ E q [log q φη (η n |x n )] = v z η l ζ j / ∈M T p(v j = v ) i∈M T q(v i = v )q(ζ n = ζ)q(η n = η )

q(l n = l )q(z n = z) log q(η n = η )dη dz dv ( α ni -1)(ψ( α ni ) -ψ( α n0 ))

E q [log q(l n |ω n , x n )] = v z η l ζ j / ∈M T p(v j = v ) i∈M T q(v i = v )q(ζ n = ζ)q(η n = η )

q(l n = l )q(z n = z) log q(l n = l )dη dz dv

= l q(l n = l ) log q(l n = l ) = l Mult(l |ω n ) log Mult(l |ω n ) = l ω nl • log ω nl E q [log q φz (z n |x n )] = v z η l ζ j / ∈M T p(v j = v ) i∈M T q(v i = v )q(ζ n = ζ)q(η n = η )

q(l n = l )q(z n = z) log q(z n = z)dη dz dv (1 + log σ 2 znj )

:

Appendix D enumerates the symbols in this paper.

We introduce the sample distribution for the real-valued data instances, and Appendix F provides the binary case as well, which we use for MNIST.

https://github.com/XifengGuo/IDEC (IDEC); https://github.com/boyangumn/DCN (DCN); https://github.com/prasoongoyal/bnp-vae (VAE-nCRP); http://vision.jhu.edu/code/ (SSC-OMP)

Figures and Tables:

Figure 2: Graphical representation of VaDE (Jiang et al., 2017) (left), VAE-nCRP (Goyal et al., 2017) (center), and neural architecture of both models (right). In the graphical representation, the white/shaded circles represent latent/observed variables. The black dots indicate hyper or variational parameters. The solid lines represent a generative model, and dashed lines represent a variational approximation. A rectangle box means a repetition for the number of times denoted by the bottom right of the box.

Figure 3: A simple depiction (left) of the key notations, where each numbered circle refers to the corresponding Gaussian mixture component. The graphical representation (center) and the neural architecture (right) of our proposed model, HCRL. The neural architecture of HCRL consists of two probabilistic encoder networks, g φη and g φz , and one probabilistic decoder network, f θ .

Figure 4: Example extracted sub-hierarchies on CIFAR-100 RCV1 v2: Figure 5 shows the embedding of RCV1 v2. VAE and VaDE show no hierarchy, and close sub-hierarchies are distantly embedded. VAE-nCRP guides the internal mixture components to be agglomerated at the center, and the cause of agglomeration is the generative process of VAE-nCRP, where the parameter of the internal components are inferred without direct information from data. HCRL shows a clear separation between the sub-hierarchy without the agglomeration.

Figure 5: Comparison of embeddings on RCV1 v2, plotted using t-SNE (Maaten & Hinton, 2008). We mark the mean of a mixture component with a numbered square, colored in {red} for VaDE, {red (root), green (internal), blue (leaf)} for VAE-nCRP and HCRL. The first-level sub-hierarchies are indicated with four colors.

Figure 6: Example extracted sub-hierarchies on 20Newsgroups

(a) Raw synthetic data ∈ R 50 (b) Hierarchically clustered embeddings ∈ R 2

Figure 7: Synthetic data in the input space of R 50 (left), which is visualized via t-SNE (Maaten& Hinton, 2008), and hierarchically clustered embeddings in the latent space of R 2 (right). We additionally show a 95% confidence ellipse with a dashed line for each Gaussian mixture component.

Initialize L Gaussian mixture components 2: T (0) ← Initialize the tree-based hierarchy having a single path with µ ζ 1:L , σ 2 ζ 1:L 3: n lock ← 0 // for counting the number of epochs, where the hierarchy has not changed 4: for each epoch e = 1, • • • , E do 5:

Figure 9: The illustration of GROW operation

for the child node whose parent path is ζ * 1

Figure 10: The illustration of PRUNE operation

Figure 11: The illustration of MERGE operation

merged two Gaussian components by following Ueda et al. (1999). The specific meaning of combining the two paths is merging the paired two Gaussian distributions lying on the two paths by level, if the two Gaussian distribtions are different. The estimation of merged Gaussian parameters, µ and σ, is the weighted summation of two subject Gaussian parameters. The propbability of the node at level l lying on a path ζ given x, p(ζ l |x), is proportional to n {q(l n = l) • ζ∈Λ q(ζ n = ζ)}, where Λ = {ζ |ζ l = ζ l and ζ ∈ P full }.

observed datapoint in m-th sequence zm,n=1,••• ,N n-th latent representation corresponding to xmn vmp The Beta draws of m-th sequence on node p, for the tree-based stickbreaking construction γ *The prior parameter for Beta distribution p(vmp) γ

In the case of an inner path, ζ [1, ζ 2 , ..., ζ l0 ] ⊂ M T * child(ζ) {[ζ, ζ l0+1 , ..., ζ L ] : ζ l0+1 > j 0 } * j 0 : maximum index for the child node whose parent path is ζ

log p(v i = v )dv = v Beta(v |a i , b i ) • log Beta(v |1, γ)dv = log Γ(1 + γ) -log Γ(1) -log Γ(γ) + (1 -1)ψ(a i ) + (γ -1)ψ(b i ) +(-1 -γ + 2)ψ(a i + b i ) = log Γ(1 + γ) -log Γ(γ) + (γ -1)ψ(b i ) + (1 -γ)ψ(a i + b i ) = log Γ(1 + γ) -log Γ(γ) + (γ -1)(ψ(b i ) -ψ(a i + b i )) = log(γΓ(γ)) -log Γ(γ) + (γ -1)(ψ(b i ) -ψ(a i + b i )) = log γ + log Γ(γ)) -log Γ(γ) + (γ -1)(ψ(b i ) -ψ(a i + b i )) = log γ + (γ -1)(ψ(b i ) -ψ(a i + b i )) E q [log p(ζ n |v)] = exp{log Γ(L -l 0 + 1) + (L -l 0 )(ψ(1) -ψ(1 + γ))} (1 -exp{ψ(γ) -ψ(1 + γ)}) l + 1) E q [log v 1,ζ2,...,ζ l ] + ζ l -1 j=1 E q [log(1 -v 1,ζ2,...,j )] × exp E q (L -l 0 ) j0 j=1 log(1 -v 1,ζ2,...,ζ l 0 ,j ) E q [log p(η n |α)] i = v )q(ζ n = ζ )q(η n = η ) q(l n = l )q(z n = z) log p(η n = η |α)dη dz dv = η q(η n = η )log p(η n = η |α)dη = η Dir(η | α n ) • log Dir(η |α)dη = log Γ(α 0 ) -

n = η )q(l n = l) log p(l n = l|η n = η )dη = l q(l n = l ) η q(η n = η ) log Mult(l |η )dη = l ω nl η Dir(η | α n ) log η l dη = l ω nl (ψ( α nl ) -ψ( α n0 )) E q [log p(z n |ζ n , l n )] = i = v )q(ζ n = ζ)q(η n = η ) q(l n = l )q(z n = z) log p(z n = z|ζ n = ζ, l n = l )dη dz dv = z l ζ q(ζ n = ζ)q(l n = l )q(z n = z) log p(z n = z|ζ n = ζ, l n = l )dz = l ζ q(ζ n = ζ)q(l n = l ) z q(z n = z) log p(z n = z|ζ n = ζ, l n = l )dz = l ζ q(ζ n = ζ )q(l n = l ) z N (z| µ z , σ 2 z I J ) • log N (z|µ ζ nl , -µ ζ nl ,j ) 2 2σ 2 ζ nl ,j E q [log p θ (x n |z n )] = v z η l ζ j / ∈M T p(v j = v ) i∈M T q(v i = v )q(ζ n = ζ)q(η n = η ) q(l n = l )q(z n = z) log p θ (x n |z n = z )dη dz dv = z q(z n = z) log p θ (x n |z n = z )dz ≈ 1 S S s=1 log p θ (x n |z (s) n ), where z (i,s) = µ(i)  x + σ (i)x (s) and (s) ∼ N (0, I J )

n = η ) log q(η n = η )dη = η Dir(η | α n ) • log Dir(η | α n )dη = log Γ( α n0 ) -

n = z) log q(z n = z)dη dz = z q(z n = z) log q(z n = z)dz = z N (z| µ z , σ 2 z I J ) • log N (z| µ z , σ 2 z I

Table2indicates that HCRL is best in the NLL and is competent in the REs which means that the hierarchically clustered embeddings preserve the intrinsic raw data structure. Test set performance of the negative log likelihood (NLL) and the reconstruction errors (REs). Replicated ten times, and the best in bold. P † < 0.05 (Student's t-test). Model-L# means that the model trained with the #-depth hierarchy.

Table content: <table><tr><td/><td colspan="2">MNIST</td><td colspan="2">CIFAR-100</td><td colspan="2">RCV1 v2</td><td colspan="2">20Newsgroups</td></tr><tr><td>Model</td><td>NLL</td><td>REs</td><td>NLL</td><td>REs</td><td>NLL</td><td>REs</td><td>NLL</td><td>REs</td></tr><tr><td>VAE</td><td colspan="2">230.71 10.46</td><td>1960.06</td><td>57.54</td><td colspan="3">2559.46 1434.59 2735.80</td><td>1788.22</td></tr><tr><td>VaDE</td><td colspan="2">217.20 10.35</td><td>1921.85</td><td>53.60</td><td colspan="3">2558.32 1426.38 2733.46</td><td>1782.86</td></tr><tr><td>IDEC</td><td colspan="2">N/A 12.75</td><td>N/A</td><td>64.09</td><td colspan="2">N/A 1376.26</td><td colspan="2">N/A 1660.61  †</td></tr><tr><td>DCN</td><td colspan="2">N/A 11.30</td><td>N/A</td><td>44.26</td><td colspan="2">N/A 1361.98</td><td>N/A</td><td>1691.17</td></tr><tr><td>IMVAE</td><td colspan="2">296.57 10.69</td><td colspan="2">1992.83 40.45  †</td><td colspan="3">2566.01 1387.02 2722.81</td><td>1718.08</td></tr><tr><td>VAE-nCRP-L3</td><td colspan="2">718.78 32.67</td><td colspan="2">2969.62 198.66</td><td colspan="3">2642.88 1538.42 2712.28</td><td>1680.56</td></tr><tr><td>VAE-nCRP-L4</td><td colspan="2">721.00 32.53</td><td colspan="2">2950.73 198.97</td><td colspan="3">2646.48 1542.81 2713.58</td><td>1680.71</td></tr><tr><td>HCRL-L3</td><td colspan="2">203.24  † 8.70  †</td><td>1843.40  †</td><td colspan="2">50.44 2554.50  †</td><td colspan="2">1395.05 2726.75</td><td>1828.71</td></tr><tr><td>HCRL-L4</td><td colspan="3">203.91  † 8.16  † 1849.13  †</td><td colspan="4">50.47 2535.43  † 1353.34 2702.88</td><td>1711.30</td></tr></table>

Hierarchical clustering accuracies with F-scores, on CIFAR-100 with a depth of three, RCV1 v2 with a depth of four, and 20Newsgroups with a depth of four. Replicated ten times, and a confidence interval with 95%. Best in bold.

Table content: <table><tr><td>Model</td><td>CIFAR-100</td><td>RCV1 v2 20Newsgroups</td></tr><tr><td>HKM</td><td colspan="2">0.1620±0.0077 0.2564±0.0679 0.4088±0.0426</td></tr><tr><td>MOHG</td><td colspan="2">0.0846±0.0378 0.1026±0.0135 0.0402±0.0119</td></tr><tr><td>RGMM</td><td colspan="2">0.1686±0.0115 0.2743±0.0521 0.4351±0.0369</td></tr><tr><td>RSSCOMP</td><td colspan="2">0.1461±0.0228 0.2657±0.0545 0.2953±0.0474</td></tr><tr><td>VAE-nCRP</td><td colspan="2">0.2011±0.0076 0.4128±0.0242 0.5584±0.0267</td></tr><tr><td>VaDE+HKM</td><td colspan="2">0.1637±0.0116 0.3308±0.0664 0.4850±0.0558</td></tr><tr><td>VaDE+MOHG</td><td colspan="2">0.1659±0.0155 0.4227±0.0927 0.4915±0.0713</td></tr><tr><td>VaDE+RGMM</td><td colspan="2">0.1806±0.0132 0.3858±0.0615 0.4095±0.0651</td></tr><tr><td colspan="3">VaDE+RSSCOMP 0.1923±0.0211 0.2718±0.0444 0.2905±0.0431</td></tr><tr><td>HCRL</td><td colspan="2">0.2245±0.0137 0.4553±0.0295 0.6008±0.0973</td></tr></table>

Table of symbols

Table content: <table><tr><td>Models</td><td>Symbol</td><td>Definition</td></tr><tr><td/><td>x/x</td><td>An observed / reconstructed datapoint</td></tr><tr><td/><td>z</td><td>A latent representation</td></tr><tr><td/><td>D/J</td><td/></tr><tr><td>All</td><td/><td/></tr></table>

Peer review (iclr2019):
 The paper presents a novel hierarchical clustering method over an embedding space. In the presented approach, both the embedding space and the hierarchical clustering are simultaneously learnt. The hierarchical clustering algorithm aims to recover complex clustering hierarchies which cannot be captured by previously proposed methods. The paper address a relevant problem, which is of great interest for extracting knowledge from data. In general, the quality of the paper is high. The presented approach is based on a sound formalization of hierarchical clustering and deep generative models. The paper is easy to follow in spite of the technical difficulty. The experimental evaluation is really extensive. It compares against many state-of-the-art methods. And the results are promising from both a quantitative and qualitative point view. The only issue with this paper is its degree of novelty, which is narrow. The proposed method adapt a previously presented hierarchical clustering method in the "standard space" (Griffiths et al., 2004) to an embedding space defined by a variational autoencoder model. The inference algorithm builds on standard techniques of deep generative models and, also, on previously proposed methods (Wand and Blei, 2003) for dealing with the complex hierarchical priors involved in this kind of models.