Source-linked AI summary
Exact Hamiltonian Monte Carlo for Truncated Multivariate Gaussians
Ari Pakman, Liam Paninski
TL;DR
Efficient sampling from truncated multivariate Gaussians is difficult when constraints induce strong correlations and high-dimensional exploration. The paper presents an exactly integrable HMC sampler for linear and quadratic constraints, extends it to piecewise-quadratic log-densities, and reports faster, more efficient sampling than Gibbs in its examples. Its scope includes runtime dependence on constraint geometry and assumptions such as positive-definite precision matrices.
Problem
Sampling truncated multivariate Gaussians can be inefficient because Gibbs sampling scales linearly with dimension and may converge slowly under strongly correlated constraints.
Method
The paper uses Hamiltonian Monte Carlo with exact Hamiltonian integration, elastic reflections at constraint walls, and extensions for piecewise-quadratic log-densities.
Results
The sampler mixes faster and is more efficient than Gibbs sampling, with a reported two- to three-order-of-magnitude difference in CPU-time-normalized comparisons for two Probit variables.
Takeaways & Limitations
The algorithm supports parallelizable implementations and can exploit special covariance structure to reduce runtime, while also covering the Bayesian Lasso through a simple extension.
Takeaways & Limitations
The method assumes a positive-definite precision matrix, and the paper reports that its adaptive travel-time strategy is unstable, recommending T = π/2 instead.
Abstract
from arXiv · showhide
We present a Hamiltonian Monte Carlo algorithm to sample from multivariate Gaussian distributions in which the target space is constrained by linear and quadratic inequalities or products thereof. The Hamiltonian equations of motion can be integrated exactly and there are no parameters to tune. The algorithm mixes faster and is more efficient than Gibbs sampling. The runtime depends on the number and shape of the constraints but the algorithm is highly parallelizable. In many cases, we can exploit special structure in the covariance matrices of the untruncated Gaussian to further speed up the runtime. A simple extension of the algorithm permits sampling from distributions whose log-density is piecewise quadratic, as in the "Bayesian Lasso" model.
1 Introduction
The paper targets efficient sampling from truncated multivariate Gaussians, whose constraints arise in many statistical and machine-learning models. It introduces an exact HMC sampler for linear and quadratic constraints, extending it to piecewise-quadratic log-densities.
- Motivation: Truncated multivariate Gaussians arise when data or parameter spaces impose linear or quadratic constraints.The constraints may be products of linear and quadratic polynomials.
- Motivation: These distributions underlie models including Probit, Tobit, dichotomized Gaussian, stochastic integrate-and-fire, and Bayesian isotonic regression.They also appear in Bayesian bridge models expressed as mixtures of Bartlett-Fejer kernels.
- Limitations of existing methods: Gibbs sampling reduces the problem to one-dimensional truncated Gaussians but can scale linearly with dimension and mix slowly under strongly correlated constraints.Hit-and-run can improve runtime yet shares the slow-convergence problem under strong constraint-induced correlations.
- Proposed approach: The proposed alternative uses HMC for constraints given by linear or quadratic functions or products thereof.HMC introduces Gaussian momentum variables and generally avoids random walks.
- Proposed approach: Exact integration of the Hamiltonian dynamics yields fast mixing, no tuning parameters, and always-accepted Metropolis steps through exact energy conservation.Constraints are enforced as hard walls with elastic particle reflections.
- Implementation and scope: Runtime depends on constraint shape and location, while wall-bounce computations are parallelizable and may enable fast implementations.Special structure in the matrix M or its inverse can further accelerate the sampler, and the method is implemented in the R package “tmg.”
- Extension: A simple extension applies the truncated-Gaussian HMC idea to distributions with piecewise-quadratic log-densities, including the Bayesian Lasso.Earlier work had also considered exactly solvable Hamiltonian equations for integrable Hamiltonians and approximating target distributions.
2 The Sampling Algorithm
The algorithm samples truncated Gaussian targets by exact Hamiltonian trajectories, reflecting elastically from constraint walls. Its efficiency depends on travel time, wall geometry, and computational structure, while extensions cover polynomial constraints.
- Hamiltonian Monte Carlo: Quadratic Gaussian forms are transformed to canonical form by a linear change of variables, which preserves the Hamiltonian and enables simple trajectory solutions.The reflected velocity then supplies the initial condition for continuing the trajectory.
- Hamiltonian Monte Carlo: The sampler alternates Gaussian momentum draws with deterministic Hamiltonian evolution for a fixed travel time T, producing the next position in the Markov chain.The joint Hamiltonian dynamics preserve the target distribution through energy and phase-space volume conservation.
- Wall reflections: Constraint boundaries are handled as hard walls: the particle finds the earliest reachable boundary, then reverses the velocity component perpendicular to that wall.For quadratic constraints, boundary times arise from an algebraic equation; products of linear or quadratic factors use the earliest zero among their factors.
- Runtime: The runtime scales linearly with the number of constraints for basic evaluations, while wall-hit calculations depend on travel time and constraint geometry.General matrix-vector operations cost O(md) and are highly parallelizable; sparsity can reduce this cost to O(d), and covariance structure can accelerate O(d^2) transformations.
- Travel time and efficiency: Travel time strongly affects efficiency: super-efficient sampling with ESF > 1 is possible, but the ESF is unstable, making T = π/2 the safe choice in practice.The paper adopts T = π/2 for subsequent examples after noting that wall shape and location also influence efficiency.
3 Examples
The examples apply exact HMC to Probit, Tobit, Brownian bridge, positive spline, and quantized Gaussian-process models, exploiting covariance structure to accelerate sampling. In the Probit example, HMC is substantially more efficient than Gibbs sampling.
- Probit and Tobit Models: The Probit posterior becomes a truncated multivariate Gaussian over regression coefficients and latent variables, with constraints yiwi ≥ 0.Only the signs of latent variables are observed, and retaining sampled β values yields posterior samples.
- Probit and Tobit Models: O(N): exploiting the Probit covariance factor Z^T reduces each matrix action from O(N^2) to O(N).The factorization is enabled by the model’s special covariance structure.
- Probit and Tobit Models: Two to three orders of magnitude: HMC ESS/CPU exceeds Gibbs ESS/CPU consistently in the 803-dimensional Probit example.The comparison summarizes 10 runs of 6,000 samples after 2,000 burn-in samples.
- Brownian bridge: O(T): the Brownian bridge has tridiagonal precision, making its Cholesky factor bidiagonal and reducing inverse-factor actions from O(T^2).The sampler targets the (T − 1)-dimensional space of possible paths.
- Bayesian reconstruction of quantized stationary Gaussian processes: O(N log N): embedding the stationary Gaussian-process covariance in a circulant matrix enables FFT multiplication for quantized-process reconstruction.The first N elements of the resulting vector have the desired covariance, and the method extends to higher-dimensional functions such as images.
4 The Bayesian Lasso
The paper extends its exact Hamiltonian Monte Carlo techniques to piecewise-quadratic log-densities, using the Bayesian Lasso as the detailed example. The method is comparable in efficiency to existing Bayesian Lasso samplers and can incorporate additional coefficient constraints.
- The techniques extend to multivariate distributions with piecewise-quadratic log-density and linear or elliptical boundaries between regions.
- The Bayesian Lasso uses a λ-dependent sparsening penalty on regression coefficients through its prior.
- Hamiltonian Monte Carlo is applied by constructing a Hamiltonian for the Bayesian Lasso posterior and evolving the associated particle trajectory.
- At coefficient-boundary crossings, the algorithm updates the piecewise dynamics while preserving continuity of β and its velocity.
- The algorithm has efficiency comparable to other Bayesian Lasso samplers and readily combines with additional coefficient constraints.
- For discontinuous piecewise log-densities, velocity jumps across region boundaries so that total energy remains conserved.