Inferring Hierarchical Mixture Structures: A Bayesian Nonparametric Approach

TL;DR

This paper introduces a new Bayesian nonparametric method combining nCRP and HDP for hierarchical clustering, effectively handling complex latent mixture data structures.

Contribution

It presents a novel approach that integrates nCRP and HDP to address hierarchical clustering of complex mixture data, a problem not previously tackled.

Findings

Achieves superior clustering accuracy on multiple datasets.

Effectively models complex latent mixture structures.

Outperforms existing Bayesian hierarchical clustering methods.

Abstract

This paper focuses on the problem of hierarchical non-overlapping clustering of a dataset. In such a clustering, each data item is associated with exactly one leaf node and each internal node is associated with all the data items stored in the sub-tree beneath it, so that each level of the hierarchy corresponds to a partition of the dataset. We develop a novel Bayesian nonparametric method combining the nested Chinese Restaurant Process (nCRP) and the Hierarchical Dirichlet Process (HDP). Compared with other existing Bayesian approaches, our solution tackles data with complex latent mixture features which has not been previously explored in the literature. We discuss the details of the model and the inference procedure. Furthermore, experiments on three datasets show that our method achieves solid empirical results in comparison with existing algorithms.

Equations

$$p(c_{n+1}=k\mid {\mathbf{c}}_{1:n})=\{\begin{array}{ll}\frac{N_k}{n+\alpha }& \text{mbox}existingk\\ \frac{\alpha }{n+\alpha }& \text{mbox}newk\end{array}{\theta }_{n+1}^{\ast }\mid {{\mathbf{\theta }}^{\ast }}_{1:n}$$

$$p(\{z_n^{\prime }{\}}_{n=1}^{m_z}\mid \alpha )=\frac{{\alpha }^{m_z}\Gamma (\alpha )}{\Gamma (N_z+\alpha )}\prod _{n=1}^{m_z}\Gamma \left(N_{z_n^{\prime }}\right)$$

$$p(\mathbf{V}\mid \alpha )=\Gamma (\alpha {)}^{|{\mathcal{I}}_{\mathcal{H}}|}{\prod }_{z\in {\mathcal{I}}_{\mathcal{H}}}\frac{{\alpha }^{m_z}}{\Gamma (N_z+\alpha )}{\prod }_{(z,z^{\prime })\in E_{\mathcal{H}}}\Gamma (N_{z^{\prime }})$$

$$p\left({\mathbf{x}}_n\mid v_{nL},\mathbf{\theta },{\mathbf{\beta }}_{v_{nL}}\right)={\sum }_{k=1}^K{\beta }_{v_{nL}k}f({\mathbf{x}}_n;{\theta }_k)+{\beta }_{v_{nL}}^{\ast }{f}^{\ast }({\mathbf{x}}_n)$$

$$p\left(c_n=k\mid {\mathbf{x}}_n,v_{nL}=z,\mathbf{B},\mathbf{\theta }\right)\propto \{\begin{array}{ll}{\beta }_{zk}f({\mathbf{x}}_n;{\theta }_k)& \text{mbox}existingk\\ {\beta }_z^{\ast }{f}^{\ast }({\mathbf{x}}_n)& \text{mbox}newk.\end{array}$$

$$p\left(v_{nL}=z^{\prime }\mid v_{n(L-1)}=z\right)p\left(c_n=k\mid {\mathbf{\beta }}_z,v_{nL}=z^{\prime },v_{n(L-1)}=z\right).$$

$$\int p(c_n=k\mid {\mathbf{\beta }}_{z^{\prime }})dp({\mathbf{\beta }}_{z^{\prime }}\mid {\mathbf{\beta }}_z)$$

$$=\frac{\Gamma \left({\sum }_{j=1}^K\gamma {\beta }_{zj}\right)}{{\prod }_{j=1}^K\Gamma (\gamma {\beta }_{zj})}\int \prod _{j=1}^K{\beta }_{z^{\prime }j}^{\text{mathds}1\{c_n=k\}+\gamma {\beta }_{zk}-1}d{\mathbf{\beta }}_{z^{\prime }}$$

$$=\frac{\Gamma \left({\sum }_{j=1}^K\gamma {\beta }_{zj}\right)}{{\prod }_{j=1}^K\Gamma (\gamma {\beta }_{zj})}\frac{{\prod }_{j=1}^K\Gamma (\text{mathds}1\{c_n=k\}+\gamma {\beta }_{zj})}{\Gamma (1+{\sum }_{j=1}^K\gamma {\beta }_{zj})}={\beta }_{zk}$$

$$\text{mbox}root{z}_0$$

$$\forall z^{\prime }:(z,z^{\prime })\in E_{\mathcal{H}}$$

$${\tilde{\mathbf{\mu }}}_k={\tilde{\mathbf{\Sigma }}}_k({\mathbf{\Sigma }}_0^{-1}{\mathbf{\mu }}_0+N_k{\mathbf{\Sigma }}^{-1}{\bar{\mathbf{x}}}_k){\tilde{\mathbf{\Sigma }}}_k=({\mathbf{\Sigma }}_0^{-1}+N_k{\mathbf{\Sigma }}^{-1}{)}^{-1}.$$

$$\mathcal{A}=\min \left\{1,\frac{p({\mathbf{x}}_n\mid {\beta }_{v_{nL}^{\prime }},v_{nL}^{\prime },\mathbf{\theta })p({\mathbf{V}}^{\prime })}{p({\mathbf{x}}_n\mid {\beta }_{v_{nL}},v_{nL},\mathbf{\theta })p(\mathbf{V})}\frac{\mathrm{nCRP}({\mathbf{v}}_n;{\mathbf{V}}^{\prime }\setminus \{{\mathbf{v}}_n^{\prime }\})}{\mathrm{nCRP}({\mathbf{v}}_n^{\prime };\mathbf{V}\setminus \{{\mathbf{v}}_n\})}\right\}$$

Full text

1 Introduction

We study the problem of Hierarchical Clustering (HC) via a Bayesian Nonparametric (BNP) modelling perspective. A BNP proposes a generative model for the observed data whose dimension is not fixed, but rather is learned from the data. In the case of HC, this allows the structure of the hierarchy to be inferred along with its clusters. Considering a model for generating the data, a BNP for HC typically associates each node in the hierarchy with a particular choice of parameters. Then the model consists of two components: 1) the random process for generating a path through the hierarchy from the root node to a leaf; 2) a parent-to-child transition kernel that models how a child node’s parameters are related to those of its parent. The generative model posits that each observation firstly selects a path randomly, and is then sampled through the distribution associated with the leaf node of the path.

An early example of the parametric generative probabilistic model is the Gaussian tree generative process [6], such that each non-root node is sampled from a Gaussian distribution with the mean of the parent node and predefined level-wise covariances. In general, for the first component, there are some well-known random processes, e.g. the nested Chinese Restaurant Process (nCRP) [4, 1], the Dirichlet Diffusion Tree (DTT) [21], the Pitman-Yor Diffusion Tree (PYDT) [15], which is a generalisation of the DTT, and the tree-structured stick-breaking construction (TSSB) [1], which generalises the (nCRP). It is proved in [15], that the nCRP and PYDT (and so also the DTT) are asymptotically equivalent—thus, this work focuses on the second component and selects the nCRP for the first component, given its simplicity. Existing BNP methods mostly utilise Gaussian-to-Gaussian (G2G) diffusion kernels for the node parameters [21, 6, 1, 15]. Within such a setting, e.g., if a parent node’s distribution is (denoted by $\sim$) $\operatorname{\mathrm{Normal}}(\bm{\mu}_{p},\mathbf{\Sigma}_{p})$, its child has the node parameters $\bm{\mu}_{c}\sim\operatorname{\mathrm{Normal}}(\bm{\mu}_{p},\mathbf{\Sigma}_{p})$ and is associated with the data distribution $\operatorname{\mathrm{Normal}}(\bm{\mu}_{c},\mathbf{\Sigma}_{c})$ where $\mathbf{\Sigma}_{p}$ and $\mathbf{\Sigma}_{c}$ can be predefined hyperparameters.

G2G kernels are easy to apply but somehow lack the ability to handle more complex patterns. In our work, we investigate the Hierarchical Dirichlet Process Mixture Model (HDPMM) for the parent-to-node transition. In this setting, each node is actually a mixture distribution and the node parameters maintain the mixture weights for a global book of components. Then, the parent-to-node transition follows a Hierarchical Dirichlet Process (HDP) [28]. As a simple example, in Fig. 1, a single level of clustering is applied to the shown data points, using G2G (Fig. 1a) and HDPMM (Fig. 1b). The HDPMM fits a mixture model to the purple nodes, which can be further refined in lower levels of the tree. Note that, the HDP [29] is discussed in [1], for applications based on Latent Dirichlet Allocations. However, our setting HDPMM is more suited to general clustering.

1.0.1 Related work

Our focus is on statistical methods, despite that there are many non-statistical methods for HC e.g. [30, 19, 16, 5], to name just a few. Statistical extensions to Agglomerative clustering (AC) have been proposed in [26, 12, 11, 18] but these are not generative models. On the other hand, Teh et al. [27] has proposed applying the Kingman’s coalescent as a prior to the HC, which is similar in spirit to DDT and PYDT, however it is a backward generative process.

It is worth mentioning that the tree generation and transition kernels that we exploit have been used in a number of other contexts, beyond HC. For example, Paisley et al. [22] discuss the nested Hierarchical Dirichlet Process (nHDP) for topic modelling which, similar to the method we present here, uses the nCRP to navigate through the hierarchy, but associates nodes in the tree with topics. The nHDP first generates a global topic tree using the nCRP where each node relates to one atom (topic) drawn from the base distribution. In summary, nHDP constitutes multiple trees with one global tree and many local trees, whereas our construction has a single tree capturing the hierarchical structure of the data. Ahmed et al. [2] also proposed a model that is very similar to the nHDP, but appealed to different inference procedures.

2 Preliminary

Dirichlet Process

We briefly describe the CRP and the stick-breaking process which are two forms of the Dirichlet Process (DP).

In the CRP, we imagine a Chinese restaurant consisting of an infinite number of tables, each with sufficient capacity to seat an infinite number of customers. A customer enters the restaurant and picks one table at which to sit. The $n^{th}$ customer picks a table based on the previous customers’ choices. That is, assuming $c_{n}$ is the table assignment label for customer $n$ and $N_{k}$ is the number of customers at table $k$, one obtains

[TABLE]

where $\bm{c}_{1:n}=\{c_{1},\ldots,c_{n}\}$, likewise for $\bm{\theta^{*}}_{1:n}$. The right hand side indicates how the parameter $\theta^{*}_{n+1}$ is drawn given the previous parameters, where each $\theta_{k}$ is sampled from a base measure $H$, and $\theta_{1},\ldots,\theta_{K}$ are the unique values among $\theta^{*}_{1},\ldots,\theta^{*}_{n}$.

Denoting a distribution by $G$, the stick-breaking process can depicted by $G=\sum_{k=1}^{\infty}\beta_{k}\delta_{\theta_{k}}$, $\{\theta_{k}\}_{k=1}^{\infty}\sim H$ and $\bm{\beta}\sim\operatorname{\mathrm{GEM}}(\alpha)$. Also, $\operatorname{\mathrm{GEM}}$ (named after Griffiths, Engen and McCloskey) known as a stick-breaking process, is analogous to iteratively breaking a portion from the remaining stick which has the initial length $1$. In particular, we write $\bm{\beta}\sim\operatorname{\mathrm{GEM}}(\alpha)$ when $u_{k}\sim\operatorname{\mathrm{Beta}}(1,\alpha)$, $\beta_{1}=u_{1}$, and $\beta_{k}=u_{k}\prod_{l=1}^{k-1}(1-u_{l})$.

Nested CRP

In the nCRP [4], customers arrive at a restaurant and choose a table according to the CRP, but at each chosen table, there is a card leading to another restaurant, which the customer visits the next day, again using the CRP. Each restaurant is associated with only a single card. After $L$ days, the customer has visited $L$ restaurants, by choosing a particular path in an infinitely branching hierarchy of restaurants.

Hierarchical Dirichlet Process Mixture Model

When the number of mixture components is infinite, we connect components along a path in the hierarchy using a HDP. The 1-level HDP [29] connects a set of DPs, $G_{j}$, to a common base DP, $G_{0}$. It can be simply written as $G_{j}\sim\operatorname{\mathrm{DP}}(\gamma,G_{0})$ and $G_{0}\sim\operatorname{\mathrm{DP}}(\gamma_{0},H)$. It has several equivalent representations while we will focus on the following form: $\bm{\beta}_{0}\sim\operatorname{\mathrm{GEM}}(\gamma_{0})$, $\bm{\beta}_{j}\sim\operatorname{\mathrm{DP}}(\gamma,\bm{\beta}_{0})$, and $\theta_{k}\sim H$ to obtain $G_{0}=\sum_{k=1}^{\infty}\beta_{0k}\delta_{\theta_{k}}$ and then $G_{j}=\sum_{k}\beta_{jk}\delta_{\theta_{k}}$. Hence $G_{j}$ has the same components as $G_{0}$ but with different mixing proportions. It may be shown that $\bm{\beta}_{j}$ can be sampled by firstly drawing $u_{jk}\sim\operatorname{\mathrm{Beta}}(\gamma\beta_{k},\gamma(1-\sum_{\ell=1}^{k}\beta_{\ell}))$ and then $\beta_{jk}=u_{jk}\prod_{\ell=1}^{k-1}(1-u_{j\ell})$. Considering each data in $x$ belongs to one of the mixture models and in particular $\bm{x}_{n}$ belongs to the group associated with $G_{j}$, the 1-level HDPMM is completed by having $\bm{x}_{n}\sim F(\theta_{c_{n}})$ where $c_{n}\sim\operatorname{\mathrm{Categorical}}(\bm{\beta}_{j})$. This process can be extended to multiple levels, by defining another level of DPs with $G_{j}$ as a base distribution and so on to higher levels. In fact, we can build a hierarchy where, for any length $L$ path in the hierarchy, the nodes in the path correspond to an $L$-level HDP.

3 Generative Process

We call our model BHMC, for Bayesian Hierarchical Mixture Clustering and illustrate it using Fig. 2.

Let $z$ be the label for a certain node in the tree. We can denote the probability to choose the first child under $z$ by $w_{z1}$. Also, let us denote the mixing proportion for a node $z$ by $\bm{\beta}_{z}$, and denote the global component assignment for the $n^{th}$ data item, $\bm{x}_{n}$, by $c_{n}$. To draw $\bm{x}_{n}$, the hierarchy is traversed to a leaf node, $z$, at which $c_{n}$ is drawn from mixing proportions $\bm{\beta}_{z}$. A path through the hierarchy is denoted by a vector e.g. $\bm{v}=\{z_{0},z_{1},z_{4}\}$. In a finite setting, with a fixed number of branches, such a path would be generated by sampling from weights $\bm{w}_{z_{0}}\sim\operatorname{\mathrm{Dir}}(\alpha)$, at the first level, then $\bm{w}_{z_{1}}\sim\operatorname{\mathrm{Dir}}(\alpha)$ at the second level and so on. The nCRP enables the path to be sampled from infinitely branched nodes. Mixing proportions, $\bm{\beta}_{z}$, along a path are connected via a multilevel HDP. Hence, if the probability of a mixture component goes to zero at any node in the tree, it will remain zero for any descendant nodes. Moreover, the smaller $\gamma$ is, the sparser the resulting distribution drawn from the HDP [20].

The generative process with an infinite configuration is shown in Algorithm 1. For a finite setting, the mixing proportions are on a finite set and Line 1 should be changed to $\bm{\beta}_{z_{0}}\sim\operatorname{\mathrm{Dir}}(\gamma_{0}/K,\dots,\gamma_{0}/K)$. Correspondingly, Line 9 has to be changed to $\bm{\beta}_{z^{\prime}}\sim\operatorname{\mathrm{Dir}}(\gamma\bm{\beta}_{z})$ according to the preliminary section.

3.1 Properties

Let us denote the tree by $\mathcal{H}=(Z_{\mathcal{H}},E_{\mathcal{H}})$. First, the node set is $Z_{\mathcal{H}}=\{z_{0},\dots,z_{M-1}\}$ where $M$ is the number of nodes. Then, $E_{\mathcal{H}}$ is the set of edges in the tree, where $(z,z^{\prime})\in E_{\mathcal{H}}$ means that there exists some path $\bm{v}_{n}$, such that $\bm{x}_{n}$ moves from $z$ to $z^{\prime}$ in the path $\bm{v}_{n}$. Hence, we infer the following variables: $\mathbf{B}=\{\bm{\beta}_{z_{0}},\dots,\bm{\beta}_{z_{M-1}}\}$, the mixing proportions for the components in the node $z$; $\bm{\theta}=\{\theta_{1},\dots,\theta_{K}\}$, component parameters; $\mathbf{V}=\{\bm{v}_{1},\dots,\bm{v}_{N}\}$ where each $\bm{v}_{n}=\{v_{n0},v_{n1},\ldots,v_{nL}\}$ is the ordered set of nodes in $Z_{\mathcal{H}}$ corresponding to the path of $\bm{x}_{n}$; $\bm{c}=\{c_{1},\dots,c_{N}\}$, the component label for all the observations. Denoting $\Phi=\{{\gamma_{0}},\gamma,\alpha,H,L\}$, we focus on the marginal prior $p(\mathbf{V},\bm{\theta},\mathbf{B}\mid\Phi)$ and obtain $p(\mathbf{V},\bm{\theta},\mathbf{B}\mid\Phi)=p(\mathbf{B}\mid\mathbf{V},\gamma_{0},\gamma)p(\mathbf{V}\mid\alpha,L)p(\bm{\theta}\mid H)$.

The first term is $p(\mathbf{B}\mid\mathbf{V},\gamma_{0},\gamma)=p\left(\bm{\beta}_{0}\mid{\gamma_{0}}\right)\prod_{(z,z^{\prime})\in E_{\mathcal{H}}}p(\bm{\beta}_{z^{\prime}}\mid\gamma,\bm{\beta}_{z})$. To expand the second term, we first denote by $m_{z}$ the number of children of $z$. Hence, the CRP probability for a set of clusters $\{z^{\prime}_{n}\}_{n=1}^{m_{z}}$ under the same parent $z$ [8, 3]:

[TABLE]

where $N_{z}$ is the number of observations in $z$. With this equation, one can observe the exchangeability of the order of arriving customers—the probability of obtaining such a partition is not dependent on the order. The tree $\mathcal{H}$ is constructed via $\mathbf{V}$ with the empty nodes all removed. For such a tree, the above result can be extended to

[TABLE]

where $\mathcal{I}_{\mathcal{H}}$ is the set of internal nodes in $\mathcal{H}$. Next, we obtain $p(\bm{\theta}\mid H)=\prod_{k=1}^{K}p(\theta_{k}\mid H)$. The likelihood for a single observation is

[TABLE]

where $\beta^{*}_{z}$ denotes the probability of drawing a new component which is always the last element in the vector $\bm{\beta}_{z}$. Here, $f(\cdot)$ is the corresponding density function for the distribution $F$ and we obtain $f(\bm{x};\theta)\equiv p(\bm{x}\mid\theta)$. Furthermore, $f^{*}(\bm{x})=\int p(\bm{x}\mid\theta)dp(\theta\mid H)$. The above presentation is similar to the Polyá urn construction of DP [23], which trims the infinite setting of DP to a finite configuration. Then, one can see $p(\mathbf{X}\mid\mathbf{V},\bm{\theta},\mathbf{B})=\prod_{n}p\left(\bm{x}_{n}\mid v_{nL},\bm{\theta},\bm{\beta}_{v_{nL}}\right)$. Finally, the unnormalised posterior is $p(\mathbf{V},\bm{\theta},\mathbf{B}\mid\mathbf{X},\Phi)\propto p(\mathbf{X}\mid\mathbf{V},\bm{\theta},\mathbf{B})p(\mathbf{V},\bm{\theta},\mathbf{B}\mid\Phi)$.

4 Inference

We appeal to Markov Chain Monte Carlo (MCMC) for inferring the model. One crucial property for facilitating the sampling procedure is exchangeability, which, as noted in the previous section, follows from the model’s connection to the CRP.

Sampling $\mathbf{V}$ and $\bm{c}$

Following [4], we sample a path for a data index as a complete variable using nCRP and decide to preserve the change based on a Metropolis-Hastings (MH) step.

Recall that the component will only be drawn at the leaf. Thus,

[TABLE]

Sampling the set $\bm{c}$ is necessary for updating the set of mixtures at each node i.e., $\mathbf{B}$. Our MH scheme applies a partially collapsed Gibbs step. Following the principles in [7], our algorithm first samples $\mathbf{V}$ with $\bm{c}$ being collapsed out, and subsequently updates $\bm{c}$ based on $\mathbf{V}$.

Sampling $\mathbf{B}$

It is straightforward to decide the sampling for a leaf node. For non-root nodes, we would like to find out $p\left(c_{n}=k\mid\mathbf{B}\setminus\{\mathbf{B}_{\bar{z}}\},v_{n(L-1)}=z\right)$ where $\mathbf{B}_{\bar{z}}$ is the set of mixing proportions in the sub-tree rooted at $z$ excluding $\bm{\beta}_{z}$. We write $p\left(c_{n}=k\mid\mathbf{B}\setminus\{\mathbf{B}_{\bar{z}}\},v_{n(L-1)}=z\right)$ to be

[TABLE]

Logically, $p(c_{n}=k\mid\mathbf{B},\bm{v}_{n})=p\left(c_{n}=k\mid\{\bm{\beta}_{v_{n\ell}}\}_{\ell=0}^{L}\right)$ which is then $p\left(c_{n}=k\mid\bm{\beta}_{v_{nL}}\right)$. We also derive $p\left(c_{n}=k\mid\bm{\beta}_{z},v_{nL}=z^{\prime},v_{n(L-1)}=z\right)$ to be

[TABLE]

given that $\Gamma(x+1)=x\Gamma(x)$ holds when $x$ is any complex number except the non-positive integers. Therefore, $p\left(c_{n}=k\mid\mathbf{B}\setminus\{\mathbf{B}_{\bar{z}}\},v_{n(L-1)}=z\right)=\beta_{zk}$. This indicates that, marginalising out the subtree rooted at $z$, the component assignment is thought to be drawn from $\bm{\beta}_{z}$, which can be seen through induction. Therefore, it allows us to conduct the size-biased permutation,

[TABLE]

where $\tilde{N}_{z^{\prime}k}=N_{z^{\prime}k}+\gamma\beta_{zk}$. Eq. 3 employs a Polyá urn posterior construction of the DP to preserve the exchangeability when carrying out the size-biased permutation [23]. This step in our inference enables an analytical form of node parameter update.

Sampling $\bm{\theta}$

Even though there are many options for $H$, we choose a Gaussian distribution in this paper. With $\theta_{k}\coloneqq\bm{\mu}_{k}$, we define $f(\bm{x};\bm{\mu}_{k})=\operatorname{\mathrm{Normal}}(\bm{x};\bm{\mu}_{k},\mathbf{\Sigma})$ and $H=\operatorname{\mathrm{Normal}}(\bm{\mu}_{0},\mathbf{\Sigma}_{0})$, where covariance matrices $\mathbf{\Sigma}$ and $\mathbf{\Sigma}_{0}$ are known and fixed. Collapsing out the unused terms, we can write $p(\theta_{k}\mid\mathbf{X},\bm{c})\sim p(\theta_{k}\mid H)\prod_{n:c_{n}=k}p(\bm{x}_{n}\mid\theta_{k})\,.$ Given the conjugacy of a Gaussian prior with a Gaussian of known covariance, by considering $\bar{\bm{x}}_{k}=\frac{1}{N_{k}}\sum_{n}\mathds{1}\{c_{n}=k\}\bm{x}_{n}$ and $N_{K}=\sum_{n}\mathds{1}\{c_{n}=k\}$, we obtain $\bm{\mu}_{k}\mid\mathbf{X},\bm{c}\sim\operatorname{\mathrm{Normal}}(\tilde{\bm{\mu}}_{k},\tilde{\mathbf{\Sigma}}_{k})$ where

[TABLE]

4.1 Algorithmic Procedure

In practice, one useful step of the inference is to truncate the infinite setting of $\bm{\beta}_{z_{0}}$ to a finite setting. Referring back to Eq. 1, one can have a threshold such that the sampling of $\bm{\beta}_{z_{0}}$ terminates when the remaining length of the stick $\beta_{z_{0}}^{*}$ is shorter than that threshold. Once a new component is initialised, each node will update its $\bm{\beta}$ by one more stick-breaking step. That is, for the root node, it samples one $u$ from $\operatorname{\mathrm{Beta}}(1,\gamma_{0})$, and assigns $\beta_{{z_{0}}(K+1)}=\beta_{z_{0}}^{*}u$ and the remaining stick length $1-\sum_{k=1}^{K+1}\beta_{{z_{0}}k}$ as a new $\beta_{z_{0}}^{*}$. For a non-root node $z^{\prime}$ inheriting from node $z$, we apply the results from the preliminary section such that $u$ is sampled from $\operatorname{\mathrm{Beta}}(\gamma\beta_{z(K+1)},\gamma\beta_{z}^{*})$ and then update in the same manner as the root node.

Algorithm 2 depicts the procedure. The MH scheme samples a proposal path and the corresponding $\bm{\beta}$’s if new nodes are initialised. MH considers an acceptance variable $\mathcal{A}$ such that $\mathcal{A}=\min\left\{1,\frac{\operatorname{\mathcal{P}}^{\prime}q(\mathbf{V},\mathbf{B}\mid\mathbf{V}^{\prime},\mathbf{B}^{\prime})}{\operatorname{\mathcal{P}}q(\mathbf{V}^{\prime},\mathbf{B}^{\prime}\mid\mathbf{V},\mathbf{B})}\right\}$. Here, $\operatorname{\mathcal{P}}$ and $\operatorname{\mathcal{P}}^{\prime}$ are the posteriors at the current and proposed states, respectively. Then, $q$ is the proposal for sampling $\mathbf{V}^{\prime}$ and $\mathbf{B}^{\prime}$, which in our case is the nCRP and HDP. At each iteration for a certain data index $n$, $\mathbf{V}$ changes to $\mathbf{V}^{\prime}$ by replacing $\bm{v}_{n}$ with $\bm{v}^{\prime}_{n}$. Thus, in this example, $q(\mathbf{V}^{\prime}\mid\mathbf{V})=\operatorname{\mathrm{nCRP}}(\bm{v}^{\prime}_{n};\mathbf{V}\setminus\{\bm{v}_{n}\})$ and vice versa. The term $q(\mathbf{B}\mid\mathbf{V})/q(\mathbf{B}^{\prime}\mid\mathbf{V}^{\prime})$ is cancelled out by the terms $p(\mathbf{B}^{\prime}\mid\mathbf{V}^{\prime})/p(\mathbf{B}\mid\mathbf{V})$ in the posterior, as $q(\mathbf{B}\mid\mathbf{V})$ and $p(\mathbf{B}\mid\mathbf{V})$ are identical. Apart from that, $\bm{\theta}$ will be updated only after all the paths are decided, and hence gain no changes. Therefore, for a specific $\bm{x}_{n}$, we have that

[TABLE]

given that the likelihood for $\mathbf{X}\setminus\{\bm{x}_{n}\}$ remains unaltered. After the paths $\mathbf{V}$ for all the observations are sampled, the process updates $\mathbf{B}$ and $\bm{\theta}$ using the manner discussed above.

4.2 Time Complexity

Contrary to first impression, at each iteration, the algorithm complexity is only $N^{2}$ for the worst case, and log-linear for the expected case.

Assume that the maximum number of children is $M_{z_{\ell}}^{*}$ at the level $\ell$, then the cost of sampling a path for a single observation with the nCRP is $O(\sum_{\ell=1}^{L}M^{*}_{z_{\ell}})$. After that, sampling $c_{n}$ is carried out with time $O(Kg(D))$ where $g(D)$ is the time for computing the Gaussian likelihood of $D$-dimensional data. In regard to the global variables, $\mathbf{B}$ will be updated with every node and thus achieves $O(K\lvert Z_{\mathcal{H}}\rvert)$. Lastly, $\bm{\theta}$ will be updated with time $O(Kg_{s}(D))$ where $g_{s}(D)$ is the complexity of sampling the Gaussian mean. In addition, we notice $\sum_{\ell=1}^{L}M^{*}_{z_{\ell}}=O(\lvert Z_{\mathcal{H}}\rvert)$ and $g(D)=O(g_{s}(D))$. Overall, for one iteration, we can summarise the complexity by $\textstyle N\sum_{\ell=1}^{L}M^{*}_{z_{\ell}}+K\lvert Z_{\mathcal{H}}\rvert+NKg_{s}(D)=O((N+K)\lvert Z_{\mathcal{H}}\rvert+NKg_{s}(D))\,.$ Since the number of components will not exceed the data size, $K\leq N$ holds. With respect to $\lvert Z_{\mathcal{H}}\rvert$, at each level, it is no more than $N$. The extreme case is that each datum is a node and extends to $L$ levels. It follows that $\lvert Z_{\mathcal{H}}\rvert\leq NL$. Hence, $O((N+K)\lvert Z_{\mathcal{H}}\rvert+NKg_{s}(D))=O\left(N^{2}(L+g_{s}(D))\right)$.

The expected number of DP components, considering $\gamma_{0}$, is $\gamma_{0}\log N$ for sufficiently large $N$ [20]. Likewise for the first level in the nCRP, which is then $\alpha\log N$. Let $\operatorname{\mathbb{E}}[\cdot]$ denote the expectation over all random $\mathcal{H}$ drawn from the BHMC. If $N$ is sufficiently large, we can have $\sum_{\ell=1}^{L}\operatorname{\mathbb{E}}[M_{z_{\ell}}]=O(L\alpha\log N)$, as the expected number of nodes will be no greater than $\alpha\log N$. However, $\operatorname{\mathbb{E}}[\lvert Z_{\mathcal{H}}\rvert]$ is hard to decide but is known to be $\leq\min(NL,(\alpha\log N)^{L})$. If $N\gg 0$ and $\alpha>0$, then $NL\leq(\alpha\log N)^{L}$. We can obtain $N(\sum_{\ell=1}^{L}\operatorname{\mathbb{E}}[M_{z}]+\operatorname{\mathbb{E}}[K]g_{s}(D))=O(N\log N(\alpha L+\gamma_{0}g_{s}(D)))$ and $\operatorname{\mathbb{E}}[K\lvert Z_{\mathcal{H}}\rvert]=O(\gamma_{0}\log N\times NL)$. Combining the two expressions, we derive the upper bound $O(N\log N(L+g_{s}(D)))$ for the average case.

5 Experiments

For qualitative analysis, we use the datasets: Animals [14] and MNIST-fashion [31]. Following this, we carry out a quantitative analysis for the Amazon text data [10] since it contains the cluster labels of all the items at multiple levels.

We fix the parameters for $H$ as $\bm{\mu}_{0}=\mathbf{0}$ and $\mathbf{\Sigma}_{0}=\mathbf{I}$. Additionally, we set the covariance matrix for $F$ as $\mathbf{\Sigma}=\sigma^{2}\mathbf{I}$. We set the number of levels $L$ intuitively based on the data. For more complex cases, one can consider the theoretical results of the effective length of the nCRP [25].

5.1 Convergence Analysis

We examine the convergence of the algorithm using Animals. Fig. 4 shows the unnormalised log likelihood for five individual simulations. The simulations quickly reach a certain satisfactory level. Apart from that, the fluctuation shows that the algorithm keeps searching the solution space over the iterations.

5.2 Results

Animals

This dataset contains 102 binary features, e.g. “has 6 legs”, “lives in water”, “bad temper”, etc. Observing the heat-map of the empirical covariance of the data, there are not many influential features. Hence, we employ Principal Component Analysis (PCA) [20] to reduce it to a seven-dimensional feature space. The MCMC burns $500$ runs and then reports the one with the greatest complete data likelihood $p(\mathbf{X},\bm{c},\mathbf{V}\mid\mathbf{B},\bm{\theta})$ amongst the following $5,000$ draws (following [1]).

We apply the hyperparameters: $\alpha=0.4,\gamma_{0}=1,\gamma=0.5,\sigma^{2}=1,L=4$. Fig. 4 shows a rather intuitive hierarchical structure. From the left to the right, there are insects, (potentially) aggressive mammals, herbivores, water animals, and birds.

MNIST-fashion

This data is a collection of fashion images. Each image is represented as a $28\times 28$ vector of grayscale pixel values. For better visualisation, we sample 100 samples evenly from each class. PCA transforms the data to 22 dimensions via the asymptotic root mean square optimal threshold [9] for keeping the singular values. Using the same criterion as for Animal, we output two hierarchies with two sets of hyperparameters. We set $\alpha=0.5,\gamma_{0}=1.5,\gamma=2,\sigma^{2}=1,L=5$. This follows precisely the same running settings as for Animals.

Fig. 5 reflects a property of the CRP which is “the rich get richer”. As it is grayscale data, in addition to the shape of the items, other factors affecting the clustering might be, e.g., the foreground/background colour area, the percentage of non-black colours in the image, the darkness/lightness of the item, etc. The hierarchical structure forms a hierarchy with high purity per level. Some mis-labelled items are expected in a clustering task.

Amazon

We uniformly down-sample the indices in the fashion category of Amazon and reserve $2,303$ entries from the data, which contain textual information about items such as their titles and descriptions111http://jmcauley.ucsd.edu/data/amazon: the data is available upon request.. The data is preprocessed via the method in [9] and only $190$ features are kept, formed as the Term Frequency and Inverse Document Frequency (TF-IDF) of the items’ titles and descriptions.

In this study, we adopt a Bayesian approach to dealing with the hyperparameters. Due to a compromise to the runtime efficiency, we run BHMC once with $500$ runs for burn-in, during which a MH-based hyperparameter sampling procedure is performed with hyperpriors $\alpha\sim\operatorname{\mathrm{Ga}}(3,1),\gamma_{0}\sim\operatorname{\mathrm{Ga}}(3.5,1),\gamma\sim\operatorname{\mathrm{Unif}}(0,1.5)$ and $\sigma^{2}\sim\operatorname{\mathrm{Unif}}(0,0.1)$. We limit $\sigma^{2}$ as mentioned since the feature values in the data are all less than $1$ and some are far less. In addition, $L$ is fixed to be $7$ (presuming that we have little information about the real number of levels). The tree with the maximal complete likelihood in the subsequent $2,500$ draws is reported.

We use the evaluation methodology in [17] to compare the clustering results against the ground-truth labels level by level. We compare 6 levels, which is the maximum branch length of the items in the ground truth. When extracting the labels from the trees (either for the algorithm outputs or the ground truth), items that are on a path of length less than 6 are extended to level 6 by assigning the same cluster label as that in their last level to the remaining levels. This is to keep the consistency of the number of items for computing metrics at each level.

For the comparison, we first consider the gold standard AC with Ward distance [30]. We adopt the existing implementations for PERCH and PYDT222https://github.com/iesl/xcluster and https://github.com/davidaknowles/pydt. For PYDT which is sensitive to the hyperparameters, we applied the authors’ implemented hyperparameter optimisation to gather the hyperparameters prior to running the repeated simulations. At each level, we consider four different evaluation metrics, namely, the purity, the normalised mutual information (NMI), the adjusted rand index (ARI), and the F-Measure [13]. Level 1 is the level for the root node.

Fig. 6 depicts that our method achieves clearly better scores with respect to purity and NMI. In the figure, PYDT-SL and PYDT-MH correspond to the slice and the MH sampling solutions, respectively. As the tree approaches to a lower level, our method also achieves a better performance in F-measure. For ARI, despite that PERCH performs the best, all numerical values are exceedingly close to [math]. However, some theoretical work of [24] suggests that ARI is more preferred in the scenario that the data contains big and equal-sized clusters. This is opposed to our ground truth which is highly unbalanced among the clusters at each level. BHMC does show the potential to perform well according to certain traditional metrics.

6 Conclusion

This paper has discussed a new perspective for Bayesian nonparametric HC. Our model, BHMC, develops an infinitely branching hierarchy of mixture parameters, that are linked along paths in the hierarchy through a multilevel HDP. A nested CRP is used to select a path in the hierarchy and mixture components are drawn from the mixture distribution in the leaf node of the selected path. The evaluation shows that BHMC is able to provide good hierarchical clustering results on three real-world datasets with different types of characteristics (i.e., binary, visual and textual) and clearly performs better than other methods with respect to purity and NMI on the Amazon dataset with ground truth, which shows the promising potential of the model.

Acknowledgements

We thank the reviewers for the helpful feedback. This research has been supported by SFI under the grant SFI/12/RC/2289_P2.