Higher-Order Hit-&-Run Samplers for Linearly Constrained Densities

Markov chain Monte Carlo (MCMC) sampling of densities restricted to linearly constrained domains is an important task arising in Bayesian treatment of inverse problems in the natural sciences. While efficient algorithms for uniform polytope sampling exist, much less work has dealt with more complex constrained densities. In particular, gradient information as used in unconstrained MCMC is not necessarily helpful in the constrained case, where the gradient may push the proposal’s density out of the polytope. In this work, we propose a novel constrained sampling algorithm, which combines strengths of higher-order information, like the target’s log-density’s gradients and curvature, with the Hit-&-Run proposal, a simple mechanism which guarantees the generation of feasible proposals, fulfilling the linear constraints. Our extensive experiments demonstrate improved sampling efficiency on complex constrained densities over various constrained and unconstrained samplers.

💡 Research Summary

The paper addresses the problem of sampling probability densities that are confined to linearly constrained domains—a situation that frequently arises in Bayesian inverse problems, structural optimization, and physics‑based simulations. While uniform sampling of convex polytopes can be performed efficiently with classic Hit‑and‑Run (HR) algorithms, these methods ignore the shape of the target distribution and therefore converge slowly in high dimensions. Conversely, modern unconstrained MCMC techniques such as the Metropolis‑Adjusted Langevin Algorithm (MALA) or Hamiltonian Monte Carlo (HMC) exploit first‑ and second‑order information (gradients and curvature) of the log‑density to propose moves that are well‑aligned with the target, but they break down when a proposal is forced outside the feasible region, leading to high rejection rates.

To bridge this gap, the authors propose a novel sampler called Higher‑Order Hit‑and‑Run (HO‑HR). The algorithm combines the feasibility guarantees of HR with the efficiency gains of higher‑order information. At each iteration the following steps are performed:

1. Higher‑order information extraction – Compute the gradient ∇log π(x) and the Hessian ∇²log π(x) of the target log‑density at the current state x.
2. Projection onto the constraint subspace – Decompose the gradient into components parallel and orthogonal to the linear constraints defined by matrix A. The orthogonal component is discarded, leaving a tangent direction d that respects the constraints. The Hessian is used to adapt the local step‑size, effectively scaling the proposal according to the curvature of the target in the feasible subspace.
3. Closed‑form feasible interval determination – For the line x + t d, solve the inequalities A(x + t d) ≤ b analytically to obtain the admissible interval