High-throughput sequencing technology brings unprecedented opportunities to quantitatively explore human gut microbiome and its relation to diseases. Microbiome data are compositional, sparse, and heterogeneous which pose serious challenges for statistical modeling. We propose a Bayesian multinomial matrix factorization model to infer overlapping clusters on both microbes and human hosts. The proposed method represents the observed zero-inflated count matrix as Dirichlet-multinomial mixtures on which latent cluster structures are built hierarchically. Under the Bayesian framework, the number of clusters is automatically determined and prior information based on a taxanomic rank tree of the microbes is incorporated into the matrix factorization which greatly improves the interpretability of the findings. We demonstrate the utility of the proposed approach using extensive simulations. Application to an inflammatory bowel disease microbiome dataset reveals interesting results, some of which contain known bacteria that are related to the disease and agree with biological knowledge.