Hierarchical Bayesian models often capture distributions over a very large number of distinct atoms. The need for this arises when organizing huge amount of unsupervised data, for instance, features extracted using deep convnets can be exploited to organize abundant unlabeled images. Inference for hierarchical Bayesian models in such cases can be rather nontrivial, leading to approximate approaches. In this work, we propose a sampler based on Cover Trees that is exact and that has guaranteed runtime logarithmic in the number of atoms and is polynomial in the inherent dimensionality of the underlying parameter space. In other words, the algorithm is as fast as search over a hierarchical data structure and we demonstrate the effectiveness on both synthetic and real datasets, consisting of over 100 million images.