Bayesian methods, powered by Markov Chain Monte Carlo estimates of posterior densities, have become a cornerstone of geophysical inverse theory. These methods have special relevance to the deep Earth, where data are sparse and uncertainties are large. We present a strategy for efficiently solving hierarchical Bayesian geophysical inverse problems for fixed parametrizations using Hamiltonian Monte Carlo sampling, and highlight an effective methodology for determining optimal parametrizations from a set of candidates by using efficient approximations to leave-one-out cross-validation for model complexity. To illustrate these methods, we use a case study of differential traveltime tomography of the lowermost mantle, using short period P-wave data carefully selected to minimize the contributions of the upper mantle and inner core. The resulting tomographic image of the lowermost mantle has a relatively weak degree 2—instead there is substantial heterogeneity at all low spherical harmonic degrees less than 15. This result further reinforces the dichotomy in the lowermost mantle between relatively simple degree 2 dominated long-period S-wave tomographic models, and more complex short-period P-wave tomographic models.