  We develop sampling algorithms to fit Bayesian hierarchical models, the
computational complexity of which scales linearly with the number of
observations and the number of parameters in the model. We focus on crossed
random effect and nested multilevel models, which are used ubiquitously in
applied sciences. The posterior dependence in both classes is sparse: in
crossed random effects models it resembles a random graph, whereas in nested
multilevel models it is tree-structured. For each class we identify a framework
for scalable computation, building on previous work. Methods for crossed models
are based on extensions of appropriately designed collapsed Gibbs samplers,
where we introduce the idea of local centering; while methods for nested models
are based on sparse linear algebra and data augmentation. We provide a
theoretical analysis of the proposed algorithms in some simplified settings,
including a comparison with previously proposed methodologies and an
average-case analysis based on random graph theory. Numerical experiments,
including two challenging real data analyses on predicting electoral results
and real estate prices, compare with off-the-shelf Hamiltonian Monte Carlo,
displaying drastic improvement in performance.
