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.
x +ε∇ log ϕ(x) x +ε∇ log ϕ(x)
Figure 1: Our proposed smLHR sampler preconditions the direction using local curvature and clips the natural gradient step to prevent it from leaving the feasible region.
Hit-&-Run (HR, Smith, 1984;Bélisle et al., 1993) is a well-established Markov chain Monte Carlo (MCMC) method for sampling from probability densities defined on constrained domains. HR is applied in diverse domains like operations research (Tervonen et al., 2013), cosmology (Lubini et al., 2013), systems biology (Herrmann et al., 2019;Theorell et al., 2017) and ecology (Gellner et al., 2023), where modelling principles give rise to linearly constrained domains. Moreover, HR has been used as a building block in advanced sampling algorithms (Theorell and Nöh, 2019;Ashton et al., 2022;Łatuszyński and Rudolf, 2024). Its main strength for constrained sampling is the algorithm’s guarantee to only generate feasible samples. Lovász and Vempala (2006) prove HR mixes rapidly for uniform sampling of bounded and convex linearly constrained domains, referred to as convex polytopes.
However, real-world inverse problems, e.g. those arising from Bayesian inference (Theorell et al., 2017;Borah Slater et al., 2023), lead to more complex non-uniform densities. These densities present additional sampling challenges such as strong nonlinear correlations between variables, or concentration of measure in small regions of the constrained domain. For sampling densities defined on unconstrained domains, the incorporation of higher-order information of the density has shown great success in dealing with these challenges. Prominent examples are gradient-guided samplers such as the Metropolis-adjusted Langevin algorithm (MALA, Rossky et al., 1978) and Hamiltonian Monte Carlo (HMC, Duane et al., 1987), or curvature-based preconditioning as in Riemannian Manifold MALA (mMALA) and Riemannian HMC (RHMC, Girolami and Calderhead, 2011). In addition, recent advances in automatic differentiation have greatly benefited the calculation of the higher-order information (Baydin et al., 2018).
To leverage higher-order information for sampling densities defined on convexly constrained domains, alternatives to HR have been proposed. (Kook et al., 2022) introduce a constrained version of Riemannian HMC for sampling log-concave densities, but test it only for uniform problems. Recently, (Srinivasan et al., 2024) proposed the first-order Metropolis-adjusted Preconditioned Langevin algorithm (MAPLA). Built upon Girolami and Calderhead (2011), MAPLA preconditions a Gaussian proposal density with a metric induced by the log-barrier function (Fiacco and McCormick, 1990) to incorporate information about the constraints instead of the curvature information of the domainconstrained density. Although MAPLA aims to guide the sampling such that proposed samples respect the constraints, its proposal density remains Gaussian, thus admitting some positive probability of proposing infeasible samples. This raises the question whether higher-order algorithms that only propose feasible samples improve the sampling efficiency for non-uniform constrained sampling problems.
Contribution In this work, we propose feasible higher-order samplers, which have positive proposal density just exactly on the convex-constrained domain. To this end, we combine the HR proposal mechanism with existing firstand second-order samplers like MALA and the simplified mMALA (Girolami and Calderhead, 2011). As a result, we develop three new samplers Langevin Hit-&-Run (LHR), simplified manifold Hit-&-Run (smHR), and simplified manifold Langevin Hit-&-Run (smLHR). Our approach is mainly driven by the observation that a Gaussian random variable can be sampled in a HR-like fashion by decomposing it into direction and magnitude. We provide theoretical analysis to prove convergence of our samplers to the desired target densities. We numerically evaluate our algorithms on a systematically constructed benchmark consisting of 2240 problems, which combine parametrizable polytopes and probability densities of varying geometry or ill condition, as well as on real-world examples from the field of Bayesian 13 C metabolic flux analysis ( 13 C-MFA, Theorell et al., 2017). Our results show that our combination of both the incorporation of higher-order information and the constraining of the proposal density by the HR proposal mechanism improve sampling efficiency.
We consider the problem of MCMC sampling from smooth probability densities π(x) with support restricted to non-empty, convex-constrained domains
defined by m linear constraints Ax ≤ b, i.e. A ∈ R m×d , b ∈ R m . If P is bounded, we refer to it as a polytope.
The Metropolis-Hastings (MH, Metropolis et al., 1953;Hastings, 1970) algorithm is a widely applied method to draw samples from a target density π : R d → R, even if access is only provided to some unnormalized density φ(x) = Z π(x), where Z is the normalization constant. It works by iterativel
This content is AI-processed based on open access ArXiv data.