PaperlayerView exact revision on arXiv

Source-linked AI summary

Hamiltonian Monte Carlo for Hierarchical Models

M. J. Betancourt, Mark Girolami

arXiv:1312.0906v1stat.ME

TL;DR

Hierarchical models create funnel-shaped posteriors, correlations, and curvature variations that limit common samplers. This paper applies Hamiltonian Monte Carlo, including careful parameterizations and curvature-aware metrics, to enable efficient inference as model complexity grows.

  • Problem

    Hierarchical models exhibit funnel-shaped densities and correlated targets that make Random Walk Metropolis, Gibbs sampling, and some parameter expansions inefficient or narrowly applicable.

  • Method

    The paper uses Hamiltonian Monte Carlo with Euclidean metrics, non-centered parameterizations, and Riemannian metrics such as SoftAbs to explore complex hierarchical targets.

  • Results

    Hamiltonian Monte Carlo is at least an order of magnitude more efficient than Random Walk Metropolis and Gibbs sampling in the reported comparison, with advantages that scale with model dimensionality.

  • Takeaways & Limitations

    Hamiltonian Monte Carlo can provide practical inference for complex hierarchical models whose local distributions and correlations impede conventional methods.

  • Takeaways & Limitations

    Hamiltonian Monte Carlo is harder to implement because it requires Hamiltonian integration and derivatives of the target distribution.

Abstract

from arXiv · show

Hierarchical modeling provides a framework for modeling the complex interactions typical of problems in applied statistics. By capturing these relationships, however, hierarchical models also introduce distinctive pathologies that quickly limit the efficiency of most common methods of in- ference. In this paper we explore the use of Hamiltonian Monte Carlo for hierarchical models and demonstrate how the algorithm can overcome those pathologies in practical applications.

I. HIERARCHICAL MODELS

Hierarchical models organize parameters into exchangeable groups that share information while preserving group-specific context. This structure creates funnel-shaped posteriors with position-dependent correlations and sharply varying curvature.

  • Hierarchical models organize parameters into exchangeable groups with shared parent information, allowing groups to learn from one another without losing unique context.
  • A one-level hierarchy distinguishes global parameters φ from local parameters θ, with data D represented through observations y_i and uncertainties σ_i.
  • Sparse data can produce a funnel density containing a high-density, low-volume region beside a low-density, high-volume region with comparable probability mass.
  • The funnel’s curvature varies strongly with position, requiring samplers to handle dramatic local deviations from isotropy to explore the posterior fully.
  • Position-dependent correlations between v and each θ_i prevent a single global rotation or rescaling from simplifying the distribution.

II. COMMON IMPLEMENTATIONS OF HIERARCHICAL MODELS

Common samplers struggle with the local correlations induced by hierarchical structure, while reparameterizations can improve geometry under restricted conditions. These efficient implementations therefore remain limited in scope.

  • A. Naïve Implementations: Random Walk Metropolis and Gibbs sampling explore highly correlated targets through small conditional-variance transitions, producing slow random walks.
  • A. Naïve Implementations: As hierarchies gain groups and levels, their worsening correlations make naïve MCMC implementations quickly impractical.
  • B. Efficient Implementations: Global covariance corrections cannot adequately address hierarchical models’ local correlations, so efficient implementations must exploit conditional-independence structure explicitly.
  • B. Efficient Implementations: Non-centered parameterizations shift dependence from latent parameters toward data and can improve Random Walk Metropolis and Gibbs sampling when data are weak.
  • B. Efficient Implementations: The funnel example shows Gibbs and Random Walk Metropolis failing to explore the marginal distribution of v after 2500 iterations.
  • B. Efficient Implementations: Parameter expansion and non-centered parameterizations apply mainly when analytic conditional results or generalized linear dependence are available.

III. HAMILTONIAN MONTE CARLO FOR HIERARCHICAL MODELS

Hamiltonian Monte Carlo uses differential-geometric dynamics and gradients to generate transitions spanning the target’s marginal variation. This enables efficient exploration of correlated hierarchical posteriors, though momentum choices can still create challenges.

  • Hamiltonian Monte Carlo uses differential geometry to generate transitions spanning full marginal variance, avoiding the random-walk behavior of Metropolis and Gibbs samplers.
  • The algorithm augments target parameters q with auxiliary momenta p, defines a Hamiltonian, and evolves the joint system through Hamilton’s equations.
  • Gradients guide transitions through high-probability regions, while integration time must adapt to the target’s shape without breaking detailed balance.
  • Hamiltonian trajectories retain efficient exploration as targets become correlated, including the largely local correlations typical of hierarchical models.

A. Euclidean Hamiltonian Monte Carlo

Euclidean Hamiltonian Monte Carlo uses q-independent Gaussian momenta and a metric that globally rotates and rescales the target, but hierarchical models expose characteristic length-scale and density-variation weaknesses.

  • Euclidean Hamiltonian Monte Carlo uses Gaussian momenta independent of q, producing a quadratic kinetic energy.
  • In the 50 + 1-dimensional funnel, reducing the integrator step size raises average acceptance above 0.651 while sampler output converges to a consistent distribution.
  • The metric Σ induces a global rotation and rescaling of the target distribution and is often taken to be the identity.
  • Euclidean Hamiltonian Monte Carlo has two weaknesses accentuated in hierarchical models: a characteristic length scale and limited density variations.

1. Characteristic Length Scale

Numerical integration introduces a characteristic step-size scale, and acceptance-based tuning can miss instability caused by spatially varying curvature. Robust adaptation therefore requires checking stabilized inferences and rare divergent transitions.

  • Symplectic numerical integration introduces a characteristic length scale through the time-discretization step size ǫ.
  • Stable integration requires tuning the step size to match curvature, because acceptance-based optimality criteria can ignore integrator instability.
  • Acceptance-based adaptation can become overly aggressive near the chain’s starting neighborhood, potentially biasing resulting inferences.
  • A robust strategy weakens adaptation until inferences stabilize and divergent transitions become rare, then selects the step size.
  • A well-chosen metric reduces the condition number of M averaged over the target distribution, improving Euclidean Hamiltonian Monte Carlo efficiency.

2. Limited Density Variations

Conservation of the Hamiltonian limits potential-energy changes by available kinetic-energy variation. Hierarchical correlations can therefore prevent Euclidean HMC transitions from spanning the full marginal variation, while non-centered parameterizations may restore performance.

  • Hamiltonian conservation requires changes in potential energy to be compensated by opposite changes in kinetic energy.
  • Because Euclidean HMC kinetic energy varies by only half the target dimensionality in expectation, each transition has limited density variation.
  • Hierarchical correlations induce huge density variations that can prevent transitions from spanning full marginal variation, causing random-walk behavior and sharply reduced efficiency.
  • Non-centered parameterizations can reduce hierarchical density variations and drastically increase Euclidean HMC performance, although efficacy depends on relative data strength.

B. Riemannian Hamiltonian Monte Carlo

Riemannian Hamiltonian Monte Carlo uses a position-dependent metric to correct local correlations and a log determinant term that can absorb and release energy. These features are intended to improve exploration in hierarchical models with difficult local geometry.

  • Riemannian Hamiltonian Monte Carlo allows the covariance, and hence the metric, to vary with position.
  • The dynamic metric provides local corrections for position-dependent correlations and can reduce integrator step-size sensitivity.
  • The log determinant term acts as an energy reservoir, potentially allowing much larger variations in potential energy.
  • The metric is chosen to approximate local curvature, but the raw Hessian may fail to be positive-definite and therefore requires manipulation.
  • Euclidean HMC with a unit metric can reduce to random-walk behavior when hierarchical trajectories are limited to moderate potential-energy variations.

IV. EXAMPLE

The example compares sampling algorithms and parameterizations on hierarchical funnel-like models. Euclidean HMC is substantially more efficient than competing samplers, while Riemannian HMC with SoftAbs produces nearly independent transitions in the funnel.

  • The log determinant term allows kinetic energy to vary sufficiently that the potential energy becomes essentially unconstrained in practice.
  • Riemannian HMC with the SoftAbs metric admits transitions spanning the entirety of the funnel distribution, yielding nearly independent transitions and much smaller autocorrelations.
  • The experiments compare Random Walk Metropolis, Metropolis-within-Gibbs, and Euclidean HMC across centered and non-centered parameterizations.The step sizes were tuned for consistent estimates, and both parameterizations were included despite expected inferiority of the centered form.
  • The comparison procedure optimized each algorithm only while maintaining mutually consistent estimates.
  • Euclidean HMC is at least an order of magnitude more efficient than Random Walk Metropolis and Gibbs under the centered parameterization.
  • The advantage of Hamiltonian Monte Carlo scales with increasing model dimensionality and can make it the only practical solution for the most complex models.

V. CONCLUSION

The paper concludes that Hamiltonian Monte Carlo uses local curvature to explore complex hierarchical models efficiently. Carefully parameterized Euclidean HMC and SoftAbs Riemannian HMC can scale with hierarchy size and local-distribution complexity, although implementation is more difficult.

  • Hamiltonian Monte Carlo uses local curvature to provide efficient exploration for complex hierarchical models.
  • Hamiltonian Monte Carlo is harder to implement because it requires Hamiltonian integration and derivatives of the target distribution.

Appendix A: Stan Models

The appendix provides Stan code for the hierarchical models used in the experiments, including parameter declarations, priors, generated data, and transformed quantities.

  • The Stan models declare latent parameters, group-level variables, observations, and scale vectors.
  • The funnel-style model assigns v a normal(0, 3) prior and theta a normal distribution with scale exp(v/2).
  • The generated one-way normal data use mu = 8, tau = 3, alpha = 10, and N = 800.
  • The models preserve generated mu and tau in print variables for reporting.
  • For each of N observations, theta[i] is sampled around mu with scale tau, sigma[i] is set to alpha, and y[i] is sampled around theta[i] with scale sigma[i].

Configuration

The experiments use hierarchical-model implementations with centered and non-centered parameterizations, alongside HMC sampling configurations for the eight-schools example.

  • The sampling commands use specified random seeds, the big_schools.dat data file, and separate CSV output files.
  • The centered model samples mu, tau, theta, and y using normal and Cauchy distributions.
  • The centered eight-schools run uses 100,000 warmup iterations, 1,000,000 samples, thinning of 1,000, adapt delta 0.999, and NUTS with max_depth 20.
  • The non-centered model samples mu, tau, var_theta, and y, then defines theta[j] as tau * var_theta[j] + mu.
  • The baseline non-centered run uses 5,000 warmup iterations, 100,000 samples, and adapt delta 0.99.
  • The non-centered E-HMC run uses 5,000 warmup iterations, 50,000 samples, and adapt delta 0.8.
Loading 1312.0906v1…