# Sampling Constrained Continuous Probability Distributions: A Review

Shiwei Lan  School of Mathematical and Statistical Sciences, Arizona State
University  Lulu Kang

###### Abstract

The problem of sampling constrained continuous distributions has frequently
appeared in many machine/statistical learning models. Many Monte Carlo Markov
Chain (MCMC) sampling methods have been adapted to handle different types of
constraints on the random variables. Among these methods, Hamilton Monte Carlo
(HMC) and the related approaches have shown significant advantages in terms of
computational efficiency compared to other counterparts. In this article, we
first review HMC and some extended sampling methods, and then we concretely
explain three constrained HMC-based sampling methods, reflection,
reformulation, and spherical HMC. For illustration, we apply these methods to
solve three well-known constrained sampling problems, truncated multivariate
normal distributions, Bayesian regularized regression, and nonparametric
density estimation. In this review, we also connect constrained sampling with
another similar problem in the statistical design of experiments of
constrained design space.

Keywords: constrained sampling; Hamilton Monte Carlo; Riemannian Monte Carlo;
regularized regression; truncated multivariate Gaussian.

##  1 Introduction

In many machine learning applications, it is necessary to sample from
distributions with various types of constraints. For example, the truncated
multivariate normal distribution can be difficult to sample from, especially
for high-dimensional cases. Even using the leave-one-out type of Gibbs
sampling scheme (Held and Holmes,, 2006; Kang et al.,, 2021), the algorithms
can still be computationally costly. Another common example is the regression
model with norm constraint on the parameters, $\|\bm{\beta}\|_{q}\leq C$, such
as Lasso ($l_{1}$ norm, $q=1$) (Tibshirani, 1996a, ) or bridge estimator
($l_{q}$ norm, $q\geq 0$) (Frank and Friedman, 1993b, ; Fu,, 1998). Other
examples include copula models, latent Dirichlet allocation, covariance matrix
estimation, and nonparametric density function estimation. Often, the
resulting models are intractable, and thus sampling from these constrained
distributions is a challenging task (Neal and Roberts,, 2008; Sherlock and
Roberts,, 2009; Neal et al.,, 2012; Brubaker et al.,, 2012; Pakman and
Paninski,, 2014). In this article, we give an overview of the statistical
sampling methods for constrained distributions. Specifically, we focus on the
distributions of continuous variables and provide a more detailed explanation
of methods based on Hamiltonian Monte Carlo (HMC) and some related approaches.
Sampling constrained discrete distributions (Jessen,, 1970; Jacob et al.,,
2021; Chewi et al.,, 2022) is not included in this review due to the different
nature from continuous distribution.

There are various types of boundary constraints on the parameters of many
statistical models, such as positive requirement, linear (summation)
constraint, upper bound on a vector norm of the parameters, etc. Sometimes,
the constraints can be considered as certain general manifolds, e.g. sphere,
positive definite matrices, and Stiefel manifold. Many sampling approaches
have been introduced to tackle one or several kinds of constraints. Based on
the nature of these methods, we roughly categorize them into three groups.

  1. 1.

Rejection Type. These methods simply discard samples that violate the
constraints or keep trying until the proposed sample satisfies the constraint.
Most Markov Chain Monte Carlo (MCMC) or other statistical sampling algorithms
can be easily modified to achieve this goal. For example, Lang et al., (2007)
proposed a rejection based sequential Monte Carlo for Bayesian estimation of
constrained dynamic system, and Li and Ghosh, (2015) also developed methods
based on rejection for truncated multivariate normal and student-$t$
distributions subject to linear inequality constraints. Since these approaches
do not directly address the constraints, they can be computationally
inefficient for complicated constraints and high-dimensional problems.

  2. 2.

Reflection Type. These approaches consider the boundary of the constrained
domain as an (energy) wall and make a reflection (hit-and-bounce) for the
sampler to move inward the constrained domain whenever it hits the (energy)
wall. For example, Neal, (2011) suggested modifying the standard HMC by
setting the potential energy to infinity for parameter values that violate the
constraints. Following this idea, Pakman and Paninski, (2014) proposed an
exact HMC for truncated multivariate Gaussian distributions. Betancourt,
(2011) and Olander, (2020) applied such HMC-based reflection idea to nested
sampling which requires likelihood-restricted prior sampling (Skilling,,
2006).

  3. 3.

Reformulation Type. These techniques transform the constrained sampling
problem or the constrained domain into something easier to work with.
Motivated by the constrained optimization methods, Brubaker et al., (2012)
proposed a family of HMC-based MCMC methods, which incorporated the constraint
on parameters $c(\bm{\theta})=0$ using Lagrange multipliers. Ahn and Chewi,
(2021) recently derived another optimization motivated algorithm using mirror-
Langevin dynamics. In many cases, distributions with constraints in
$\mathbb{R}^{d}$ can be transformed into distributions on manifolds. Some HMC-
related algorithms have been created to sample distributions on manifolds. For
example, Kook et al., (2022) used Riemannian manifold HMC instead of the
original HMC with a Lagrange multiplier. Byrne and Girolami, (2013) showed how
HMC methods can be designed for and applied to the distribution defined on
manifolds embedded in Euclidean space by the explicit forms for geodesics if
they exist. In particular, motivated by Byrne and Girolami, (2013), Spherical
HMC (Lan et al.,, 2014; Lan and Shahbaba,, 2016) does not require any manifold
embedding. It focuses on constraints that can be transformed into vector norms
and eventually mapped onto a hyper-sphere. Spherical HMC can be viewed as a
more efficient special case of Brubaker et al., (2012) and Kook et al.,
(2022). A related work SPInS (Chaudhry et al.,, 2021) maps a hyper-ball
containing the constrained domain inside out so that sampler can be defined on
a larger unbounded space.

In general, algorithms of the rejection type tend to be inefficient because
frequent rejected attempts cause a significant waste of computation. The
reflection type is intuitive but its efficiency usually depends on the
specific constraints. Certain constrained domains may need too many
reflections in high dimensional space and hence the excessive computational
time. The reformulation type is more intrinsic and sophisticated in
incorporating the constraints in the step of proposing new samples. More
importantly, they usually can be easily scaled to large dimensions and handle
more varieties of constraints. Therefore, we are going to focus on the
reflection and reformulation types of methods in the following review.

In recent years, Optimal Transport (OT) and the more general variational
inference methods have gained much attention from the machine learning
community and have been adapted to many statistical and machine learning
models. Although in theory variational inference methods, including OT, can
deal with distributions with any compact support regions, in numerical
implementation it is much more challenging when the support region is not the
entire $\mathbb{R}^{d}$ with $d$ being the dimension of the variables. Das and
Bhattacharya, (2020) addressed the issue of state-dependent nonlinear
equality-constrained state estimation using Bayesian filtering based on OT.
Ahn and Chewi, (2021) proposed using the mirror-Langevin algorithm, which is a
discretization of the mirror-Langevin diffusion, for constrained sampling.

Constrained sampling is naturally connected with numerical integration in
constrained domains. Quasi-Monte Carlo methods deal with numerical integration
over box constrained domains (Owen,, 2013; Leobacher and Pillichshammer,,
2014). Other quadrature methods have to be modified for non-box constraints
(Olshanskii and Safin,, 2016; Gessner et al.,, 2020; Legrain,, 2021; Saye,,
2022). Constrained sampling is also related to the statistical design of
experiments. Draguljić et al., (2012), Pratola et al., (2017), and Huang et
al., (2021) proposed different methods to generate space-filling designs,
which essentially approximate the uniform distribution, in various irregular-
constrained domains. Kang, (2019) developed algorithms to generate different
optimal designs, including a distance-criterion-based space-filling design, in
complicated design regions.

In the remaining article, we first review the necessary background on HMC and
its related methods in Section 2. Next, we review in detail three reflection
and reformulation approaches in Section 3. Three constrained sampling problems
are illustrated in Section 4, including truncated multivariate Gaussian
distributions, Bayesian regularized regression, and non-parametric density
estimation. In Section 5, we connect constrained sampling with the constrained
design of experiments and review some existing methods for constructing
different types of constrained designs. The article concludes in Section 6.

##  2 Hamiltonian Monte Carlo and Extensions

MCMC can be inherently inefficient due to its random walk nature. Different
from the Gibbs sampler and Metropolis algorithms, HMC simulates the
Hamiltonian dynamics to propose new states and reduce the local random walk
behavior, and thus moves more rapidly towards the target distribution. The
proposed states of HMC are significantly distant from the current states and
yet still have a high acceptance probability. Neal, (2011) recalled the origin
and history of HMC and detailed the HMC method and its appealing properties.
Here we briefly review the HMC algorithm and its manifold extensions such as
the Riemannian/Lagrangian Monte Carlos, geodesic Monte Carlo, etc.

###  2.1 Hamiltonian Monte Carlo

We begin with an introduction to the Hamiltonian dynamics, which is the basis
of HMC. It consists of a $d$-dimensional vector $\bm{\theta}$, called
_position_ state, and a $d$-dimensional vector $\bm{\phi}$, called _momentum_
state. For illustration, Neal, (2011) used a simple physical example, the
dynamics of a frictionless puck that slides over a surface of varying height.
The _potential energy_ of the puck, denoted by $U(\bm{\theta})$, is
proportional to the height given position $\bm{\theta}$. The _kinetic energy_
of the puck, denoted by $K(\bm{\phi})$, is equal to $|\bm{\phi}|^{2}/m$, with
$m$ the _mass_ of the puck. The total energy of the dynamic, called
_Hamiltonian function_ and denoted by $H(\bm{\theta},\bm{\phi})$, is the sum
of the potential energy and kinetic energy of the puck, i.e.,

$H(\bm{\theta},\bm{\phi})=U(\bm{\theta})+K(\bm{\phi}).$ (1)

As the puck slides over the surface, the potential energy increases as it
moves over a rising slope and the kinetic energy decreases as the velocity of
the puck, $\bm{\phi}/m$, decreases. The two energies changes in opposite
directions when the puck moves over a descending slope. Due to frictionless
assumption, the total energy remains the same as the initial state of the
system. The system of $(\bm{\theta},\bm{\phi})$ evolves following the
_Hamilton’s equations_

$\displaystyle\dot{\bm{\theta}}=\frac{\mathrm{d}\bm{\theta}}{\mathrm{d}t}$
$\displaystyle=\frac{\partial
H(\bm{\theta},\bm{\phi})}{\partial\bm{\phi}}=\nabla_{\bm{\phi}}K(\bm{\phi}),$
(2) $\displaystyle\dot{\bm{\phi}}=\frac{\mathrm{d}\bm{\phi}}{\mathrm{d}t}$
$\displaystyle=-\frac{\partial
H(\bm{\theta},\bm{\phi})}{\partial\bm{\theta}}=-\nabla_{\bm{\theta}}U(\bm{\theta}).$
(3)

HMC applies the Hamiltonian dynamics to MCMC sampling. The position vector
$\bm{\theta}$ is the random variable of interest. The potential energy is
$U(\bm{\theta})=-\log p(\bm{\theta})$, where $p(\bm{\theta})$ is the density
function of the target distribution to be sampled from. In the Bayesian
framework, $p(\bm{\theta})$ is the posterior distribution
$p(\bm{\theta}|\bm{y})$, where $\bm{y}$ represents data and $\bm{\theta}$ the
unknown parameters. The momentum vector $\bm{\phi}$ can be considered as an
auxiliary random variable that is usually assumed to follow a multivariate
normal distribution $\mathcal{N}({\bf 0},\bm{M})$, where $\bm{M}$ is a user-
specified covariance matrix (often chosen to be identity matrix $\bm{I}$),
also known as a _mass matrix_. Thus the kinetic energy becomes
$K(\bm{\phi})=-\log N(\bm{\phi}|{\bf
0},\bm{M})=\bm{\phi}^{\top}\bm{M}^{-1}\bm{\phi}/2+\text{constant}$. To
numerically solve (2) and (3) with the above potential and kinetic energies,
the _leapfrog_ method (Verlet,, 1967) is commonly used to approximate the
Hamilton’s equations by discretizing time. The HMC uses the leapfrog method to
simulate the Hamiltonian dynamics for some time horizon $\tau=L\epsilon$ to
propose new samples that are further accepted or rejected according to certain
probability as the next state. Algorithm 1 shows the HMC procedure.

Algorithm 1 Hamiltonian Monte Carlo (HMC)

1:Initialize $\bm{\theta}^{(0)}$ at current $\bm{\theta}_{t-1}$, and randomly
sample $\bm{\phi}^{(0)}\sim\mathcal{N}({\bf 0},\bm{M})$.

2:for $\ell=0$ to $L-1$ do

3: Update $\bm{\phi}$ by a half-step of $\epsilon$:
$\bm{\phi}^{(\ell+\frac{1}{2})}=\bm{\phi}^{(\ell)}+\frac{1}{2}\epsilon\nabla_{\bm{\theta}}\log
p(\bm{\theta}^{(\ell)})$.

4: Update $\bm{\theta}$ by a full-step of $\epsilon$:
$\bm{\theta}^{(\ell+1)}=\bm{\theta}^{(\ell)}+\epsilon\bm{M}^{-1}\bm{\phi}^{(\ell+\frac{1}{2})}$.

5: Update $\bm{\phi}$ by another half-step of $\epsilon$:
$\bm{\phi}^{(\ell+1)}=\bm{\phi}^{(\ell+\frac{1}{2})}+\frac{1}{2}\epsilon\nabla_{\bm{\theta}}\log
p(\bm{\theta}^{(\ell+1)})$.

6:end for

7:Set $(\bm{\theta}^{*},\bm{\phi}^{*})=(\bm{\theta}^{(L)},\bm{\phi}^{(L)})$
and compute the accept rate $r=\frac{p(\bm{\theta}^{*})N(\bm{\phi}^{*}|{\bf
0},\bm{M})}{p(\bm{\theta}_{t-1})N(\bm{\phi}_{t-1}|{\bf 0},\bm{M})}$.

8:Set $\bm{\theta}_{t}=\bm{\theta}^{*}$ with probability $\min(r,1)$ and
$\bm{\theta}_{t}=\bm{\theta}_{t-1}$ otherwise.

Following either the Hamilton’s equations (2) and (3) or the leapfrog
procedure, the intuition of HMC is straightforward as explained in Gelman et
al., (2013). When the current value of $\bm{\theta}$ is at a flat region of
$p(\bm{\theta})$, similar to the situation when the puck is on a flat surface,
the velocity of the puck, and thus its momentum $\bm{\phi}$, becomes close to
a constant. Therefore, the position $\bm{\theta}$ would move at a constant
speed exploring the flat region. If the position $\bm{\theta}$ moves to a
region with decreasing density of $p(\bm{\theta})$, which is not favorable,
then $\nabla_{\bm{\theta}}\log p(\bm{\theta})$ is negative, and thus the
momentum $\bm{\phi}$ would decrease in the direction of movement. Next, the
position $\bm{\theta}$ would move in this unfavorable direction with a reduced
velocity. The trends reverse if $\bm{\theta}$ value moves to a region with an
increasing density of $p(\bm{\theta})$. The system possesses three important
properties for the proof of ergodicity: (i) time-reversibility (going from the
end of a trajectory with the reversed momentum takes the sampler back to the
starting point); (ii) volume-preservation (volume moving along the flow
$T_{t}:(\bm{\theta},\bm{\phi})\mapsto(\bm{\theta}^{*},\bm{\phi}^{*})$ does not
change); and (iii) energy conservation (approximately under discretized
system). The exact proof of convergence of HMC can be found in Neal, (1994)
and Neal, (2011). The choice of $\epsilon$ should be sufficiently small so
that the acceptance rate is high but not too small so that the computation is
still efficient. The length of trajectory is also a crucial parameter for HMC
and it can be varied from 20 to 1000, depending on the complexity and
dimension of the problem. Trial and error can be used for setting both
$\epsilon$ and $L$. More discussion can be found in Neal, (2011).

###  2.2 Riemannian and Lagrangian Monte Carlo

Girolami and Calderhead, (2011) extended HMC to Riemannian HMC (RHMC) by
defining the Hamiltonian dynamics on a Riemannian manifold of distributions.
Compared to HMC, RHMC can exploit fully the geometric properties of the
parameter space of $\bm{\theta}$ using a position-specific mass matrix, i.e.,
$\bm{M}=\bm{G}(\bm{\theta})$, where $\bm{G}(\bm{\theta})$ is usually set as
the Fisher information matrix of $p(\bm{\theta})$. The distribution of
momentum vector $\bm{\phi}$ is $\mathcal{N}({\bf 0},\bm{G}(\bm{\theta}))$,
which is no longer independent of $\bm{\theta}$. The Hamiltonian function (1)
becomes

$H(\bm{\theta},\bm{\phi})=\psi(\bm{\theta})+K(\bm{\theta},\bm{\phi})=-\log
p(\bm{\theta})+\frac{1}{2}\log\det(\bm{G}(\bm{\theta}))+\frac{1}{2}\bm{\phi}^{\top}\bm{G}(\bm{\theta})^{-1}\bm{\phi}.$
(4)

where $\psi(\bm{\theta})=-\log
p(\bm{\theta})+\frac{1}{2}\log\det(G(\bm{\theta}))$, and
$K(\bm{\theta},\bm{\phi})=\frac{1}{2}\bm{\phi}^{\top}\bm{G}(\bm{\theta})^{-1}\bm{\phi}$.
Due to dependence of $\bm{\phi}$ on $\bm{\theta}$, the dynamics of
$\bm{\theta}$ and $\bm{\phi}$ become non-separable. The previous version of
the leapfrog is not applicable. Instead, the _generalized leapfrog_ (ISERLES,,
1986) is used, which is an implicit scheme of fixed-point iterations.

To avoid the time-consuming iterations, Lan et al., (2015) introduced the
_Lagrangian Monte Carlo_ (LMC). It uses a variable transformation approach
that changes the Hamiltonian dynamics to Lagrangian dynamics. Specifically,
let $\bm{v}=\bm{G}(\bm{\theta})^{-1}\bm{\phi}$ and its distribution is
$\bm{v}\sim\mathcal{N}({\bf 0},\bm{G}(\bm{\theta})^{-1})$. As $\bm{v}$ is the
momentum divided by mass, it can be considered as velocity intuitively. The
original Hamilton’s equations in HMC (2) and (3) become the following
Lagrangian dynamics (a.k.a. Euler-Lagrange equation):

$\displaystyle\frac{\mathrm{d}\bm{\theta}}{\mathrm{d}t}$
$\displaystyle=\bm{v},$ (5)
$\displaystyle\frac{\mathrm{d}\bm{v}}{\mathrm{d}t}$
$\displaystyle=-\bm{\eta}(\bm{\theta},\bm{v})-\bm{G}(\bm{\theta})^{-1}\nabla_{\bm{\theta}}\psi(\bm{\theta}),$
(6)

where $\bm{\eta}(\bm{\theta},\bm{v})$ is a vector whose $k$th element is
$\bm{v}^{\top}\bm{\Gamma}^{k}(\bm{\theta})\bm{v}$. Here
$\bm{\Gamma}^{k}_{i,j}(\bm{\theta}):=\frac{1}{2}\sum_{l}g^{k,l}(\partial_{i}g_{l,j}+\partial_{j}g_{i,l}-\partial_{l}g_{i,j})$
are the Christoffel symbols, where $g_{i,j}=[\bm{G}(\bm{\theta})]_{i,j}$ and
$g^{i,j}=[\bm{G}(\bm{\theta})^{-1}]_{i,j}$ and $\partial_{i}$ means partial
derivative with respect to $\theta_{i}$. Based on the Lagrangian dynamics, Lan
et al., (2015) proposed an explicit integrator, which is time reversible but
not volume preserving. This is different from the HMC and RHMC. However, one
can adjust the acceptance probability with the Jacobian determinant to satisfy
the detailed balance condition. The LMC algorithm in Lan et al., (2015) is
shown to be computationally more stable and efficient than RHMC.

Byrne and Girolami, (2013) developed these manifold HMC algorithms for a class
of problems where the geodesic equation (the system (5) and (6) without
$\bm{G}(\bm{\theta})^{-1}\nabla_{\bm{\theta}}\psi(\bm{\theta})$ term) can be
analytically solved. Lan et al., (2014); Lan and Shahbaba, (2016) proposed a
specific geodesic MC on a hyper-sphere (whose geodesic is a big circle) and
applied it to sample from distributions with constraints defined by vector
norm.

##  3 Constrained Sampling based on HMC

In this section, we explain three different constrained sampling methods which
are adapted from the original HMC. They are _Wall HMC_ , _Constrained HMC_ ,
and _Sphere HMC_ , which are among the most representative ones in the
reflection and reformulation types of algorithms.

###  3.1 Constrained HMC by Reflection

Neal, (2011) discussed a method of handling the constraint $c(\bm{\theta})\geq
0$ by modifying the original potential energy $U(\bm{\theta})$ to create a
“soft” wall:

$U_{r}(\bm{\theta})=U(\bm{\theta})+w(\bm{\theta}),\quad
w(\bm{\theta})=\begin{cases}0,&\;\text{if}\;c(\bm{\theta})\geq 0,\\\
r^{r+1}|c(\bm{\theta})|^{r},&\text{else}.\end{cases}$ (7)

Such constraint becomes a “hard” wall with infinite barrier

$\tilde{U}(\bm{\theta})=\lim_{r\to\infty}U_{r}(\bm{\theta})=\begin{cases}U(\bm{\theta}),&\;\text{if}\;c(\bm{\theta})\geq
0,\\\ \infty,&\text{else}.\end{cases}$ (8)

Suppose the sampler just hits the wall, i.e. $c(\bm{\theta})<0$. According to
(3) with $U_{r}$, we have the momentum $\bm{\phi}$ updated as

$\bm{\phi}_{t+1}=\bm{\phi}_{t}-\nabla_{\bm{\theta}}U(\bm{\theta})\Delta
t-r^{r+2}\frac{|c(\bm{\theta})|^{r}}{c(\bm{\theta})}\Delta
t\nabla_{\bm{\theta}}c(\bm{\theta})$ (9)

We can choose the time step $\Delta t=C(r)\to 0$ as $r\to\infty$ such that we
have the following perfect reflection for updating momentum (Betancourt,,
2011)

$\bm{\phi}_{t+1}=\bm{\phi}_{t}-2\langle\bm{\phi}_{t},\bm{n}\rangle\bm{n},\quad\bm{n}=\nabla_{\bm{\theta}}c(\bm{\theta})/\|\nabla_{\bm{\theta}}c(\bm{\theta})\|.$
(10)

Such reflection based HMC (named as “Wall HMC” in Lan et al., (2014)) proceeds
with $\tilde{U}(\bm{\theta})$ replacing $U(\bm{\theta})$ in Algorithm 1 where
the half-step momentum updates (lines 3 and 5) are replaced by the above
reflection (10) once hitting the wall (constraint violated). Betancourt,
(2011) and Olander, (2020) successfully applied this algorithm to nest
sampling. Pakman and Paninski, (2014) developed a more sophisticated exact HMC
algorithm for truncated multivariate Gaussian distributions based on a similar
idea of reflection.

###  3.2 Constrained HMC by Reformulation

Brubaker et al., (2012) considered HMC for a general constraint
$c(\bm{\theta})=0$ that defines a connected, differentiable submanifold of
$\mathbb{R}^{d}$, denoted as
$\mathcal{M}=\\{\bm{\theta}\in\mathbb{R}^{d}|c(\bm{\theta})=0\\}$. The
constraint determines the _tangent bundle_ of $\mathcal{M}$,
$\mathcal{TM}=\\{(\bm{\theta},\dot{\bm{\theta}})|c(\bm{\theta})=0\;\text{and}\;\frac{\partial
c}{\partial\bm{\theta}}\dot{\bm{\theta}}=0\\}$ where $\frac{\partial
c}{\partial\bm{\theta}}$ is the Jacobian of the constraints. Now we have the
new Hamiltonian as

$H(\bm{\theta},\bm{\phi},\lambda)=\hat{H}(\bm{\theta},\bm{\phi})+\lambda^{\top}c(\bm{\theta}),\quad\hat{H}(\bm{\theta},\bm{\phi})=\psi(\bm{\theta})+K(\bm{\theta},\bm{\phi})$
(11)

where $\lambda$ is the Lagrange multiplier. Then the Hamiltonian dynamics with
the above guided Hamiltonian $\hat{H}(\bm{\theta},\bm{\phi})$ become (2) and
(3) with extra an equation $c(\bm{\theta})=0$ and are defined on the
_cotangent bundle_
$\mathcal{T^{*}M}=\\{(\bm{\theta},\bm{\phi})|c(\bm{\theta})=0\;\text{and}\;\frac{\partial
c}{\partial\bm{\theta}}\frac{\partial\hat{H}}{\partial\bm{\phi}}=0\\}$.
Brubaker et al., (2012) used a consistent integrator called RATTLE (Andersen,,
1983; Leimkuhler and Reich,, 1994) to solve the constrained Hamiltonian
dynamics in their proposed constrained HMC (CHMC).

$\displaystyle\bm{\phi}^{(\ell+\frac{1}{2})}$
$\displaystyle=\bm{\phi}^{(\ell)}-\frac{\epsilon}{2}\left[\frac{\partial\hat{H}(\bm{\theta}^{(\ell)},\bm{\phi}^{(\ell+\frac{1}{2})})}{\partial\bm{\theta}}+\frac{\partial
c}{\partial\bm{\theta}}(\bm{\theta}^{(\ell)})\lambda\right]$ (12)
$\displaystyle\bm{\theta}^{(\ell+1)}$
$\displaystyle=\bm{\theta}^{(\ell)}+\frac{\epsilon}{2}\left[\frac{\partial\hat{H}(\bm{\theta}^{(\ell)},\bm{\phi}^{(\ell+\frac{1}{2})})}{\partial\bm{\phi}}+\frac{\partial\hat{H}(\bm{\theta}^{(\ell+1)},\bm{\phi}^{(\ell+\frac{1}{2})})}{\partial\bm{\phi}}\right]$
$\displaystyle 0$ $\displaystyle=c(\bm{\theta}^{(\ell+1)})$
$\displaystyle\bm{\phi}^{(\ell+1)}$
$\displaystyle=\bm{\phi}^{(\ell+\frac{1}{2})}-\frac{\epsilon}{2}\left[\frac{\partial\hat{H}(\bm{\theta}^{(\ell+1)},\bm{\phi}^{(\ell+\frac{1}{2})})}{\partial\bm{\theta}}+\frac{\partial
c}{\partial\bm{\theta}}(\bm{\theta}^{(\ell+1)})\mu\right]$ $\displaystyle 0$
$\displaystyle=\frac{\partial
c}{\partial\bm{\theta}}(\bm{\theta}^{(\ell+1)})\frac{\partial\hat{H}(\bm{\theta}^{(\ell+1)},\bm{\phi}^{(\ell+1)})}{\partial\bm{\phi}}$

where $\lambda$ and $\mu$ are Lagrange multipliers associated with the state
and momentum constraints (the third and the fifth equations of (12)). This
generalizes the leapfrog method to handle the manifold constraints. The
proposed state $\bm{\theta}^{*}$ is accepted with probability defined by $H$
in (11).

###  3.3 Spherical HMC

Lan et al., (2014) considered HMC defined on a special manifold, hyper-sphere,
denoted as
$\mathcal{S}^{d}=\\{\bm{\theta}\in\mathbb{R}^{d+1}|\|\bm{\theta}\|_{2}^{2}=\sum_{i=1}^{d+1}\theta_{i}^{2}=1\\}$.
This algorithm is particularly useful to handle a class of constraints defined
by the following vector ($\bm{\beta}\in\mathbb{R}^{d}$) $q$-norm:

$\|\bm{\beta}\|_{q}=\begin{cases}(\sum_{i=1}^{d}|\beta_{i}|^{q})^{1/q},&q\in(0,+\infty)\\\
\max_{1\leq i\leq d}|\beta_{i}|,&q=+\infty\end{cases}$ (13)

Figure 1: Transforming $q$-norm constrained domain to the unit ball. Left:
from general $q$-norm domain $\mathcal{Q}^{d}$ to unit ball
${\mathcal{B}}_{0}^{d}(1)$; Right: from the unit cube ${\mathcal{C}}^{d}$ to
the unit ball ${\mathcal{B}}_{0}^{d}(1)$.

The $q$-norm domain,
$\mathcal{Q}^{d}:=\\{\bm{\beta}\in\mathbb{R}^{d}|\|\bm{\beta}\|_{q}\leq 1\\}$,
can be transformed to the unit ball
$\mathcal{B}_{0}^{d}(1):=\\{\bm{\theta}\in\mathbb{R}^{d}|\|\bm{\theta}\|_{2}\leq
1\\}$ by either
$\beta_{i}\mapsto\theta_{i}=\mathrm{sgn}(\beta_{i})|\beta_{i}|^{q/2}$ (shown
the left panel of Figure 1) or
$\bm{\beta}\mapsto\bm{\theta}=\bm{\beta}\frac{\|\bm{\beta}\|_{\infty}}{\|\bm{\beta}\|_{2}}$
(shown in the right panel of Figure 1).

Figure 2: Transforming the unit ball $\mathcal{B}_{\bf 0}^{d}(1)$ to the
sphere $\mathcal{S}^{d}$.

To define HMC for these constrained distributions, Lan and Shahbaba, (2016)
proposed an idea of _spherical augmentation_ to further map the unit ball
$\mathcal{B}_{\bf 0}^{d}(1)$ to the hyper-sphere $\mathcal{S}^{d}$ by
appending an auxiliary variable $\theta_{d+1}$ to the original vector
$\bm{\theta}\in\mathcal{B}_{\bf 0}^{d}(1)$ such that the extended parameter
$\tilde{\bm{\theta}}=(\bm{\theta},\theta_{d+1})\in\mathcal{S}^{d}$. The lower
hemisphere ($\theta_{d+1}=-\sqrt{1-\|\bm{\theta}\|_{2}^{2}}$) is also
identified with the upper hemisphere
($\theta_{d+1}=\sqrt{1-\|\bm{\theta}\|_{2}^{2}}$) by ignoring the sign of
$\theta_{d+1}$. After collecting samples $\\{\tilde{\bm{\theta}}\\}$ using
spherical HMC defined on the sphere, $\mathcal{S}^{d}$, the last component
$\theta_{d+1}$ is discarded and the obtained samples $\\{\bm{\theta}\\}$
automatically satisfy the constraint $\|\bm{\theta}\|_{2}\leq 1$. As
illustrated in Figure 2, the boundary of the constraint, i.e.,
$\|\bm{\theta}\|_{2}=1$, corresponds to the equator on the sphere
$\mathcal{S}^{d}$. Therefore, as the sampler moves on the sphere, e.g. from
$A$ to $B$, passing across the equator from one hemisphere to the other
translates to “bouncing back” off the boundary in the original parameter
space.

On the Riemannian manifold $(\mathcal{S}^{d},{\bf G}(\bm{\theta}))$ where
${\bf G}(\bm{\theta})={\bf
I}_{d}+\bm{\theta}{\bm{\theta}}^{\top}/(1-\|\bm{\theta}\|_{2}^{2})$ is the
_canonical spherical metric_ , the tangent space at $\tilde{\bm{\theta}}$ is
defined as
$\mathcal{T_{\tilde{\bm{\theta}}}S}^{d}=\\{\tilde{\bm{v}}=(\bm{v},v_{d+1})\in\mathbb{R}^{d+1}|\tilde{\bm{\theta}}^{\top}\tilde{\bm{v}}=0\\}$.
We have the Hamiltonian (4) redefined on the tangent bundle
$\mathcal{TS}^{d}=\\{(\tilde{\bm{\theta}},\tilde{\bm{v}})|\|\tilde{\bm{\theta}}\|_{2}=1\;\text{and}\;\tilde{\bm{\theta}}^{\top}\tilde{\bm{v}}=0\\}$:

$H(\tilde{\bm{\theta}},\tilde{\bm{v}})=H^{*}(\tilde{\bm{\theta}},\tilde{\bm{v}})+\frac{1}{2}\log\det(\bm{G}(\bm{\theta})),\quad
H^{*}(\tilde{\bm{\theta}},\tilde{\bm{v}})=U(\tilde{\bm{\theta}})+K(\tilde{\bm{v}})$
(14)

where the potential energy $U(\tilde{\bm{\theta}})=U(\bm{\theta})$, and the
kinetic energy $K(\tilde{\bm{v}})=\frac{1}{2}{\bm{v}}^{\top}{\bf
G}(\bm{\theta}){\bm{v}}=\frac{1}{2}\|\tilde{\bm{v}}\|_{2}^{2}$ which is
defined for the velocity random variable
$\tilde{\bm{v}}\sim\mathcal{N}(0,\mathcal{P}(\tilde{\bm{\theta}}))$ with
$\mathcal{P}(\tilde{\bm{\theta}})={\bf
I}_{d+1}-\tilde{\bm{\theta}}\tilde{\bm{\theta}}^{\top}$ being the projection
matrix.

With such guided Hamiltonian function $H^{*}$ in (14), the Hamiltonian
dynamics can be defined on the Riemannian manifold $(\mathcal{S}^{d},{\bf
G}(\bm{\theta}))$ in terms of $(\bm{\theta},{\bf p})$, or equivalently as the
following Lagrangian dynamics in terms of
$(\tilde{\bm{\theta}},\tilde{\bm{v}})$ (Lan et al.,, 2015):

$\displaystyle\begin{aligned} &\dot{\tilde{\bm{\theta}}}&&=&&\tilde{\bm{v}}\\\
&\dot{\tilde{\bm{v}}}&&=&&-\|\tilde{\bm{v}}\|_{2}^{2}\tilde{\bm{\theta}}-\mathcal{P}(\tilde{\bm{\theta}})\nabla_{\tilde{\bm{\theta}}}U(\bm{\theta})\end{aligned}$
(15)

where the projection matrix $\mathcal{P}(\tilde{\bm{\theta}})$ maps the
directional derivative $\nabla_{\tilde{\bm{\theta}}}U(\bm{\theta})$ onto the
tangent space $\mathcal{T_{\tilde{\bm{\theta}}}S}^{d}$. The dynamics (15) can
be split into two dynamics

where the solution to (16a) only involves updating velocity and (16b) is the
geodesic equation on $\mathcal{S}^{d}$ with the solution being the big circle.
_Spherical HMC (SphHMC)_ proceeds by proposing a joint state by alternate
updates according to (16a) (lines 5 and 8 in Algorithm 2) and (16b) (lines 6-7
in Algorithm 2) and accepting or rejecting the proposal based on the
acceptance probability defined by $H$. Algorithm 2 summarizes the details of
spherical HMC (SphHMC).

Algorithm 2 Spherical HMC (SphHMC)

1:Initialize $\tilde{\bm{\theta}}^{(0)}$ at current $\tilde{\bm{\theta}}$.

2:Sample a new velocity value $\tilde{\bm{v}}^{(0)}\sim\mathcal{N}({\bf
0},{\bf I}_{d+1})$, and set
$\tilde{\bm{v}}^{(0)}\leftarrow\mathcal{P}(\tilde{\bm{\theta}}^{(0)})\tilde{\bm{v}}^{(0)}$.

3:Calculate
$H(\tilde{\bm{\theta}}^{(0)},\tilde{\bm{v}}^{(0)})=U(\bm{\theta}^{(0)})+K(\tilde{\bm{v}}^{(0)})$

4:for $\ell=0$ to $L-1$ do

5:
$\tilde{\bm{v}}^{(\ell+\frac{1}{2})}=\tilde{\bm{v}}^{(\ell)}-\frac{\epsilon}{2}\mathcal{P}(\tilde{\bm{\theta}}^{(\ell)})\nabla_{\tilde{\bm{\theta}}}U(\bm{\theta}^{(\ell)})$

6:
$\tilde{\bm{\theta}}^{(\ell+1)}=\tilde{\bm{\theta}}^{(\ell)}\cos(\|\tilde{\bm{v}}^{(\ell+\frac{1}{2})}\|\epsilon)+\frac{\tilde{\bm{v}}^{(\ell+\frac{1}{2})}}{\|\tilde{\bm{v}}^{(\ell+\frac{1}{2})}\|}\sin(\|\tilde{\bm{v}}^{(\ell+\frac{1}{2})}\|\epsilon)$

7:
$\tilde{\bm{v}}^{(\ell+\frac{1}{2})}\leftarrow-\tilde{\bm{\theta}}^{(\ell)}\|\tilde{\bm{v}}^{(\ell+\frac{1}{2})}\|\sin(\|\tilde{\bm{v}}^{(\ell+\frac{1}{2})}\|\epsilon)+\tilde{\bm{v}}^{(\ell+\frac{1}{2})}\cos(\|\tilde{\bm{v}}^{(\ell+\frac{1}{2})}\|\epsilon)$

8:
$\tilde{\bm{v}}^{(\ell+1)}=\tilde{\bm{v}}^{(\ell+\frac{1}{2})}-\frac{\epsilon}{2}\mathcal{P}(\tilde{\bm{\theta}}^{(\ell+1)})\nabla_{\tilde{\bm{\theta}}}U(\bm{\theta}^{(\ell+1)})$

9:end for

10:Calculate
$H(\tilde{\bm{\theta}}^{(L)},\tilde{\bm{v}}^{(L)})=U(\bm{\theta}^{(L)})+K(\tilde{\bm{v}}^{(L)})$

11:Calculate the acceptance probability
$\alpha=\min\\{1,\exp[-H(\tilde{\bm{\theta}}^{(L)},\tilde{\bm{v}}^{(L)})+H(\tilde{\bm{\theta}}^{(0)},\tilde{\bm{v}}^{(0)})]\\}$

12:Accept or reject the proposal according to $\alpha$ for the next state
$\tilde{\bm{\theta}}^{\prime}$

##  4 Applications of Constrained Sampling

###  4.1 Truncated Multivariate Gaussian

Figure 3: Density plots of a truncated bivariate Gaussian using exact density
function (upper leftmost) and MCMC samples from RWM, Wall HMC, exact HMC,
c-SphHMC, and s-SphHMC respectively. Solid elliptical curves always show true
unconstrained probability density contours. Dashed lines define linear
constrained domains. Colored heatmaps indicate constrained probability density
based on truth or estimation from MCMC samples.

For illustration purposes, we start with a truncated bivariate Gaussian
distribution,

$\binom{\beta_{1}}{\beta_{2}}\sim\mathcal{N}\left(\bf{0},\begin{bmatrix}1&.5\\\
.5&1\end{bmatrix}\right),\qquad 0\leq\beta_{1}\leq 5,\quad 0\leq\beta_{2}\leq
1$

This is a box type constraint with the lower and upper limits as ${\bf
l}=(0,0)$ and ${\bf u}=(5,1)$ respectively. The original rectangle domain can
be mapped to 2d unit disc $\mathcal{B}_{\bf 0}^{2}(1)$ to use c-SphHMC, or
mapped to 2d rectangle $\mathcal{R}_{\bf 0}^{2}$ where s-SphHMC can be
directly applied (Lan and Shahbaba,, 2016).

The upper leftmost panel of Figure 3 shows the heatmap based on the exact
density function, and the other panels show the corresponding heatmaps based
on MCMC samples from RWM, Wall HMC, exact HMC, c-SphHMC, and s-SphHMC
respectively. All algorithms generate probability density estimates that
visually match the true density. Table 1 compares the true mean and covariance
of the above truncated bivariate Gaussian distribution with the point
estimates using $2\times 10^{5}$ ($2\times 10^{4}$ for each of 10 repeated
experiments with different random seeds) MCMC samples in each method. Overall,
all methods estimate the mean and covariance reasonably well.

Table 1: Comparing the point estimates for the mean and covariance of a bivariate truncated Gaussian distribution using RWM, Wall HMC, exact HMC, and SphHMC. Method | Mean | Covariance
---|---|---
Truth | $\begin{bmatrix}0.7906\\\ 0.4889\end{bmatrix}$ | $\begin{bmatrix}0.3269&0.0172\\\ 0.0172&0.08\end{bmatrix}$
RWM | $\begin{bmatrix}0.7796\pm 0.0088\\\ 0.4889\pm 0.0034\end{bmatrix}$ | $\begin{bmatrix}0.3214\pm 0.009&0.0158\pm 0.001\\\ 0.0158\pm 0.001&0.0798\pm 5e-04\end{bmatrix}$
Wall HMC | $\begin{bmatrix}0.7875\pm 0.0049\\\ 0.4884\pm 8e-04\end{bmatrix}$ | $\begin{bmatrix}0.3242\pm 0.0043&0.017\pm 0.001\\\ 0.017\pm 0.001&0.08\pm 3e-04\end{bmatrix}$
exact HMC | $\begin{bmatrix}0.7909\pm 0.0025\\\ 0.4885\pm 0.001\end{bmatrix}$ | $\begin{bmatrix}0.3272\pm 0.0026&0.0174\pm 7e-04\\\ 0.0174\pm 7e-04&0.08\pm 3e-04\end{bmatrix}$
SphHMC | $\begin{bmatrix}0.79\pm 0.005\\\ 0.4864\pm 0.0016\end{bmatrix}$ | $\begin{bmatrix}0.3249\pm 0.0045&0.0172\pm 0.0012\\\ 0.0172\pm 0.0012&0.0801\pm 0.001\end{bmatrix}$

To evaluate the efficiency of the above-mentioned methods, we repeat this
experiment for higher dimensions, $D=10$, and $D=100$. As before, we set the
mean to zero and set the $(i,j)$-th element of the covariance matrix to
$\Sigma_{ij}=1/(1+|i-j|)$. Further, we impose the following constraints on the
parameters,

$0\leq\beta_{i}\leq u_{i}$

where $u_{i}$ (i.e., the upper bound) is set to 5 when $i=1$; otherwise, it is
set to $0.5$.

For each method, we obtain $10^{5}$ MCMC samples after discarding the initial
$10^{4}$ samples. We set the tuning parameters of algorithms such that their
overall acceptance rates are within a reasonable range. As shown in Table 2,
Spherical HMC algorithms are substantially more efficient than RWM and Wall
HMC. For RWM, the proposed states are rejected about $95\%$ of times due to
violation of the constraints. On average, Wall HMC bounces off the wall at
around 3.81 ($L=2$) and 6.19 ($L=5$) times per iteration for $D=10$ and
$D=100$ respectively. Exact HMC is quite efficient for relatively low
dimensional truncated Gaussian ($D=10$); however, it becomes very slow for
higher dimensions ($D=100$). In contrast, by augmenting the parameter space,
Spherical HMC algorithms handle the constraints more efficiently. Since
s-SphHMC is more suited for box-type constraints, it is substantially more
efficient than c-SphHMC in this example.

Table 2: Comparing the efficiency of RWM, Wall HMC, exact HMC, and SphHMC in
terms of sampling from truncated Gaussian distributions.

Dimension | Method | AP a | s/iter b | ESS(min,med,max) c | Min(ESS)/s d | speedup
---|---|---|---|---|---|---
| RWM | 0.62 | 5.72E-05 | (48,691,736) | 7.58 | 1.00
| Wall HMC | 0.83 | 1.19E-04 | (31904,86275,87311) | 2441.72 | 322.33
D=10 | exact HMC | 1.00 | 7.60E-05 | (1e+05,1e+05,1e+05) | 11960.29 | 1578.87
| SphHMC | 0.82 | 2.53E-04 | (62658,85570,86295) | 2253.32 | 297.46
| RWM | 0.81 | 5.45E-04 | (1,4,54) | 0.01 | 1.00
| Wall HMC | 0.74 | 2.23E-03 | (17777,52909,55713) | 72.45 | 5130.21
D=100 | exact HMC | 1.00 | 4.65E-02 | (97963,1e+05,1e+05) | 19.16 | 1356.64
| SphHMC | 0.73 | 3.45E-03 | (55667,68585,72850) | 146.75 | 10390.94

a acceptance probability
b seconds per iteration
c (minimum, median, maximum) effective sample size
d minimal ESS per second

###  4.2 Bayesian Regularized Regression

In regression analysis, overly complex models tend to overfit the data.
Regularized regression models control complexity by imposing a penalty on
model parameters. _Bridge regression_ (Frank and Friedman, 1993b, ) is a
family of regression models where the coefficients are obtained by minimizing
the residual sum of squares subject to a constraint on the magnitude of
regression coefficients as follows:

$\min_{\bm{\beta}}\sum_{i=1}^{n}(y_{i}-X_{i}\bm{\beta})^{2}\quad
subject\;to\quad\|\bm{\beta}\|_{q}\leq r$ (17)

When $q=1$, this corresponds to _Lasso_ (least absolute shrinkage and
selection operator) proposed by Tibshirani, 1996a  which allows the model to
force some of the coefficients to become exactly zero (i.e., become excluded
from the model). When $q=2$, this model is known as _ridge regression_. Bridge
regression is more flexible by allowing different $q$ norm constraints for
different effects on shrinking the magnitude of parameters (See Figure 4).

Park and Casella, (2008) and Hans, (2009) proposed _Bayesian Lasso_ by
employing a conjugate prior distribution of the form
$P(\beta)\propto\exp(-\lambda|\beta|)$. Frank and Friedman, 1993b  also
constructed a complicated prior for the Bayesian inference of bridge
regression. With spherical HMC, there is much flexibility in choosing priors
with $q$-norm constraints. We can define the following _Bayesian regularized
linear regression model_ :

$\displaystyle\bm{y}|\bm{X},\bm{\beta},\sigma^{2}_{\epsilon}$
$\displaystyle\sim\mathcal{N}(\bm{X}\bm{\beta},\sigma^{2}_{\epsilon}{\bf I})$
(18) $\displaystyle\bm{\beta}$ $\displaystyle\sim
p(\bm{\beta})\bm{1}(\|\bm{\beta}\|_{q}\leq r)$

Figure 4: Bayesian Bridge Regression by Spherical HMC: Lasso (q=1, left),
q=1.2 (middle), and q=0.8 (right).

Figure 4 compares the parameter estimates of Bayesian Lasso to the estimates
obtained from two Bridge regression models with $q=1.2$ and $q=0.8$ for the
diabetes data set (N=442, D=10) studied in Park and Casella, (2008). Truncated
Gaussian prior $\bm{\beta}\sim\mathcal{N}(0,\sigma^{2}{\bf
I})\bm{1}(\|\bm{\beta}\|_{q}\leq r)$ is considered and posterior samples are
collected using spherical HMC algorithm. Figure 4 illustrates the parameter
estimates $\hat{\bm{\beta}}$ with respect to the shrinkage factor
$s:=\|\hat{\bm{\beta}}^{\textrm{Lasso}}\|_{1}/\|\hat{\bm{\beta}}^{\textrm{OLS}}\|_{1}$
varying from 0 to 1, where $\hat{\bm{\beta}}^{\textrm{OLS}}$ denotes the
estimates obtained by ordinary least squares (OLS) regression. As expected,
tighter constraints (e.g., $q=0.8$) would lead to faster shrinkage of
regression parameters as we decrease $s$. Note, the model (18) can easily be
generalized to generalized linear models or nonparametric models such as
Gaussian process regression.

###  4.3 Non-parametric Density Estimation

In this example of non-parametric density estimation, we show how a density
function $p(x)$ can be modeled on an infinite dimensional sphere and how the
spherical HMC can be applied to the Bayesian inference to efficiently obtain
the posterior estimate.

Suppose we want to attribute a smooth density function $p(x)$ to observed data
$\\{x_{n}\\}_{n=1}^{N}$ on finite domain $\mathcal{D}\subset\mathbb{R}^{d}$.
Define the space of density functions $p(x)$ and the space of square-root
density functions, $q(x)=\sqrt{p(x)}$:

$\displaystyle\mathcal{P}:=$ $\displaystyle\left\\{p:\mathcal{D}\to\mathbb{R}\
|\ p\geq 0,\int_{\mathcal{D}}p(x)\,\mu(dx)=1\right\\}$
$\displaystyle\mathcal{Q}:=$ $\displaystyle\left\\{q:\mathcal{D}\to\mathbb{R}\
|\ \int_{\mathcal{D}}q(x)^{2}\,\mu(dx)=1\right\\}\,,$

respectively. Although the space $\mathcal{P}$ contains the functions of
interest, we instead opt to deal with the space $\mathcal{Q}$ of square-root
densities, which can be viewed as the unit sphere in the infinite-dimensional
Hilbert space $L^{2}(\mathcal{D})$. We model the square-root density with a GP
prior multiplied by the Dirac measure restricting the function to the unit
sphere:

$\displaystyle
q(\cdot)\sim\mathcal{GP}(0,\mathcal{C})\times\delta(q(\cdot)\in\mathcal{Q})\,.$
(19)

where $\mathcal{C}$ has eigen-pairs $\\{\lambda_{\ell},\phi_{\ell}(x)\\}$ such
that $\mathcal{C}\phi_{\ell}(x)=\lambda_{\ell}\phi_{\ell}(x)$ for
$\ell\in\mathbb{N}$.

Based on the Karhunen-Loève representation (Wang,, 2008) of the Gaussian
random function, we have the following expansion of $q$:

$\displaystyle q(\cdot)=\sum_{\ell=1}^{\infty}q_{\ell}\,\phi_{i}(\cdot),$
$\displaystyle q_{\ell}\stackrel{{\scriptstyle
ind}}{{\sim}}N(0,\lambda^{2}_{\ell})$ (20)

If we let $\\{\phi_{\ell}(x)\\}$ be an orthonormal basis on
$L^{(}\mathcal{D})$, then the unit sphere restriction,
$\delta(q\in\mathcal{Q})$, translates to following requirement on the infinite
sequence $q:=\\{q_{\ell}\\}$:

$q\in\mathcal{S}^{\infty}=\left\\{q\in\ell^{2}\big{|}\,\langle
q,q\rangle_{\ell^{2}}=\sum_{i=1}^{\infty}q_{i}^{2}=1\right\\}\,.$

In practice, we truncate the K-L expansion (20) at $L>0$ terms and have
$q=\\{q_{\ell}\\}_{\ell=1}^{L}\in\mathcal{S}^{L-1}$.

Now we have the prior for $q$ and the likelihood of the data
$x:=\\{x_{n}\\}_{n=1}^{N}$ given $q$ as follows

$\displaystyle\begin{aligned}
\pi(q)\propto&\delta(q\in\mathcal{S}^{L-1})\prod_{i=1}^{L}\exp\big{(}-q_{i}^{2}/(2\lambda_{i}^{2})\big{)}\,,\\\
\pi(x|q)=&\prod_{n=1}^{N}q^{2}(x_{n})=\prod_{n=1}^{N}\left|\sum_{\ell=1}^{L}q_{\ell}\,\phi_{\ell}(x_{n})\right|^{2},\end{aligned}$

We can then apply the spherical HMC to sample from the posterior
$\pi(q|x)\propto\pi(q)\pi(x|q)$ which is naturally defined on the sphere
$\mathcal{S}^{L-1}$.

Figure 5 depicts 1,000 data points (red) drawn from four different
distributions on the unit square along with the contours of the pointwise
median of 1,000 posterior draws from the model $\pi(q|x)$ as defined above.
The data in the first three plots are generated using truncated Gaussians and
mixtures of truncated Gaussians. The data for the last plot is generated by
Gaussian noise added to the uniform distribution on the circle. The model
adapts easily to multimodal and patterned data samples.

Figure 5: The contours (black) of the posterior median from 1,000 draws of the
$\chi^{2}$-process density sampler. Each posterior is conditioned on 1,000
data points (red).

##  5 Constrained Design of Experiments

The statistical design of experiments is naturally related to statistical
sampling, although they achieve different goals. Driven by practical needs,
sometimes the design variables have to satisfy certain constraints, which
makes the design space irregular, i.e., non-box.

A common type of constrained physical experiment is the mixture experiment,
which is widely used in the chemical, pharmaceutical, and food industries. The
input variables of a mixture experiment are the proportions of different
ingredients of the whole mixture. Thus the values of all the variables are
between 0 and 1 (without scaling or other transformation), and the sum of them
must be equal to 1. In Piepel et al., (2005), the experiment contains 21
ingredients including both box and linear constraints. Mixture experiments can
also have layered mixture structures which makes both design and modeling more
complicated (Kang et al.,, 2011; Shen et al.,, 2020).

Many computer experiments also involve irregular constraints, and space-
filling design methods have to be adapted to deal with these constraints.
Stinstra et al., (2003) showed a large-scale computer experiment for
television tube design that involves 23 input variables and 44 non-box
constraints. Draguljić et al., (2012) illustrated a Total Elbow Replacement
computer simulation with four input variables and two linear non-box
constraints.

There have been many design methods proposed to handle the constraints.
Draguljić et al., (2012) developed a column-wise construction method for
space-filling design. In Pratola et al., (2017), a dense candidate set is
generated and then the unfeasible candidate points are removed and the optimal
design points are selected from the feasible ones. Maximin distance and the
(robust) IMSPE criteria are used for the space-filling design. Huang et al.,
(2021) proposed a constrained minimum energy design method for constructing
space-filling designs in any non-regular bounded space, and its key idea is to
apply the minimum energy design (a deterministic sampling algorithm) on the
target distribution using the probabilistic constraints proposed in
sequentially constrained Monte Carlo.

The aforementioned methods are only for constrained space-filling designs,
which can be considered as sampling from uniform distributions with
constraints. Kang, (2019) proposed a generic method to construct optimal
designs for an irregular constrained design space. A stochastic coordinate-
exchange (SCE) method is developed. In each iteration of coordinate exchange,
the multi-dimensional constraints are projected into the dimension of the
coordinate to be exchanged (or improved) with the other coordinates fixed at
the current values. Therefore, the multi-dimensional constraints are reduced
to one-dimensional box constraints. The generic method can be adapted for
different design criteria, including the $D$\- and linear-optimal design for
physical experiments, and the $\phi_{p}$-optimal space-filling design for
computer experiments. Figure 6 shows three optimal designs constructed by the
SCE method for different kinds of constraints.

Figure 6: Different constrained designs returned by SCE in Kang, (2019): (a)
$D$-optimal design; (b) $I$-optimal design; and (c) space-filling design for
the glacier area in the case study of Pratola et al., (2017).

##  6 Conclusion

In this article, we review in detail the Hamilton Monte Carlo sampling and its
variants, and more importantly, how they are modified to overcome various
constraints. Based on the nature of the algorithms, we categorize them into
three groups, rejection, reflection, and reformulation. Specifically, we
explain three constrained HMC-based sampling algorithms, Wall HMC, Constrained
HMC, and Sphere HMC. Wall HMC is a reflection-type algorithm and the other two
belong to the reformulation group. Three important applications of constrained
sampling algorithms are illustrated. They are truncated multivariate Gaussian,
Bayesian regularized regression, and non-parametric density estimation.

Constrained sampling is an important problem and has broad applications in
many statistical/machine learning models. Some models that are typically not
considered to be constrained sampling problems can be solved by constrained
sampling algorithms, such as regularized regression/classification and density
estimation. With the rapid development of AI and big data technologies, future
research on this topic is most likely in the direction of creating fast
algorithms for high-dimensional distributions with complex constraints, their
theoretical foundations, and software implementations.

## Funding Information

Shiwei Lan’s work is supported by NSF grant DMS-2134256. Lulu Kang’s work is
supported by NSF grants DMS-1916467 and DMS-2153029.

## Acknowledgments

We thank the editor, the associate editor, and reviewers in advance for
reviewing this paper. We look forward to their comments and suggestions.

## References

  * Ahn and Chewi, (2021) Ahn, K. and Chewi, S. (2021).  Efficient constrained sampling via the mirror-langevin algorithm.  Advances in Neural Information Processing Systems, 34:28405–28418.
  * Andersen, (1983) Andersen, H. C. (1983).  Rattle: A “velocity” version of the shake algorithm for molecular dynamics calculations.  Journal of Computational Physics, 52(1):24–34.
  * Betancourt, (2011) Betancourt, M. (2011).  Nested sampling with constrained hamiltonian monte carlo.  In AIP Conference Proceedings, volume 1305, pages 165–172. American Institute of Physics.
  * Brubaker et al., (2012) Brubaker, M. A., Salzmann, M., and Urtasun, R. (2012).  A family of mcmc methods on implicitly defined manifolds.  In Lawrence, N. D. and Girolami, M. A., editors, Proceedings of the Fifteenth International Conference on Artificial Intelligence and Statistics (AISTATS-12), volume 22, pages 161–172.
  * Byrne and Girolami, (2013) Byrne, S. and Girolami, M. (2013).  Geodesic monte carlo on embedded manifolds.  Scandinavian Journal of Statistics, 40(4):825–845.
  * Chaudhry et al., (2021) Chaudhry, S., Lautzenheiser, D., and Ghosh, K. (2021).  An efficient scheme for sampling in constrained domains.  10.48550/ARXIV.2110.10840.
  * Chewi et al., (2022) Chewi, S., Gerber, P. R., Lu, C., Le Gouic, T., and Rigollet, P. (2022).  Rejection sampling from shape-constrained distributions in sublinear time.  In International Conference on Artificial Intelligence and Statistics, pages 2249–2265. PMLR.
  * Das and Bhattacharya, (2020) Das, N. and Bhattacharya, R. (2020).  Optimal transport based filtering with nonlinear state equality constraints.  IFAC-PapersOnLine, 53(2):2373–2378.
  * Draguljić et al., (2012) Draguljić, D., Santner, T. J., and Dean, A. M. (2012).  Noncollapsing space-filling designs for bounded nonrectangular regions.  Technometrics, 54(2):169–178.
  * (10) Frank, I. E. and Friedman, J. H. (1993a).  A Statistical View of Some Chemometrics Regression Tools.  Technometrics, 35(2):109–135.
  * (11) Frank, L. E. and Friedman, J. H. (1993b).  A statistical view of some chemometrics regression tools.  Technometrics, 35(2):109–135.
  * Fu, (1998) Fu, W. J. (1998).  Penalized regressions: the bridge versus the lasso.  Journal of computational and graphical statistics, 7(3):397–416.
  * Gelman et al., (2013) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., and Rubin, D. B. (2013).  Bayesian data analysis.  Chapman and Hall/CRC, New York.
  * Gessner et al., (2020) Gessner, A., Kanjilal, O., and Hennig, P. (2020).  Integrals over gaussians under linear domain constraints.  In International Conference on Artificial Intelligence and Statistics, pages 2764–2774. PMLR.
  * Girolami and Calderhead, (2011) Girolami, M. and Calderhead, B. (2011).  Riemann manifold Langevin and Hamiltonian Monte Carlo methods.  Journal of the Royal Statistical Society, Series B, (with discussion) 73(2):123–214.
  * Hans, (2009) Hans, C. (2009).  Bayesian lasso regression.  Biometrika, 96(4):835–845.
  * Held and Holmes, (2006) Held, L. and Holmes, C. C. (2006).  Bayesian auxiliary variable models for binary and multinomial regression.  Bayesian analysis, 1(1):145–168.
  * Huang et al., (2021) Huang, C., Joseph, V. R., and Ray, D. M. (2021).  Constrained minimum energy designs.  Statistics and Computing, 31(6):1–15.
  * ISERLES, (1986) ISERLES, A. (1986).  Generalized Leapfrog Methods.  IMA Journal of Numerical Analysis, 6(4):381–392.
  * Jacob et al., (2021) Jacob, P. E., Gong, R., Edlefsen, P. T., and Dempster, A. P. (2021).  A gibbs sampler for a class of random convex polytopes.  Journal of the American Statistical Association, 116(535):1181–1192.
  * Jessen, (1970) Jessen, R. J. (1970).  Probability sampling with marginal constraints.  Journal of the American Statistical Association, 65(330):776–796.
  * Kang, (2019) Kang, L. (2019).  Stochastic coordinate-exchange optimal designs with complex constraints.  Quality Engineering, 31(3):401–416.
  * Kang et al., (2011) Kang, L., Roshan Joseph, V., and Brenneman, W. A. (2011).  Design and modeling strategies for mixture-of-mixtures experiments.  Technometrics, 53(2):125–136.
  * Kang et al., (2021) Kang, X., Ranganathan, S., Kang, L., Gohlke, J., and Deng, X. (2021).  Bayesian auxiliary variable model for birth records data with qualitative and quantitative responses.  Journal of Statistical Computation and Simulation, 91(16):3283–3303.
  * Kook et al., (2022) Kook, Y., Lee, Y. T., Shen, R., and Vempala, S. S. (2022).  Sampling with riemannian hamiltonian monte carlo in a constrained space.  10.48550/ARXIV.2202.01908.
  * Lan and Shahbaba, (2016) Lan, S. and Shahbaba, B. (2016).  Algorithmic Advances in Riemannian Geometry and Applications, chapter 2-Sampling Constrained Probability Distributions Using Spherical Augmentation, pages 25–71.  Advances in Computer Vision and Pattern Recognition. Springer International Publishing, 1 edition.
  * Lan et al., (2015) Lan, S., Stathopoulos, V., Shahbaba, B., and Girolami, M. (2015).  Markov chain monte carlo from lagrangian dynamics.  Journal of Computational and Graphical Statistics, 24(2):357–378.
  * Lan et al., (2014) Lan, S., Zhou, B., and Shahbaba, B. (2014).  Spherical hamiltonian monte carlo for constrained target distributions.  In Eric P. Xing, T. J., editor, The 31st International Conference on Machine Learning, volume 32, pages 629–637, Beijing.
  * Lang et al., (2007) Lang, L., shiang Chen, W., Bakshi, B. R., Goel, P. K., and Ungarala, S. (2007).  Bayesian estimation via sequential monte carlo sampling—constrained dynamic systems.  Automatica, 43(9):1615–1622.
  * Legrain, (2021) Legrain, G. (2021).  Non-negative moment fitting quadrature rules for fictitious domain methods.  Computers & Mathematics with Applications, 99:270–291.
  * Leimkuhler and Reich, (1994) Leimkuhler, B. and Reich, S. (1994).  Symplectic integration of constrained hamiltonian systems.  Mathematics of Computation, 63(208):589–589.
  * Leobacher and Pillichshammer, (2014) Leobacher, G. and Pillichshammer, F. (2014).  Introduction to quasi-Monte Carlo integration and applications.  Springer.
  * Li and Ghosh, (2015) Li, Y. and Ghosh, S. K. (2015).  Efficient sampling methods for truncated multivariate normal and student-t distributions subject to linear inequality constraints.  Journal of Statistical Theory and Practice, 9(4):712–732.
  * Neal and Roberts, (2008) Neal, P. and Roberts, G. O. (2008).  Optimal scaling for random walk metropolis on spherically constrained target densities.  Methodology and Computing in Applied Probability, Vol.10(No.2):277–297.
  * Neal et al., (2012) Neal, P., Roberts, G. O., and Yuen, W. K. (2012).  Optimal scaling of random walk metropolis algorithms with discontinuous target densities.  Annals of Applied Probability, Volume 22(Number 5):1880–1927.
  * Neal, (1994) Neal, R. M. (1994).  An improved acceptance procedure for the hybrid monte carlo algorithm.  Journal of Computational Physics, 111(1):194–203.
  * Neal, (2011) Neal, R. M. (2011).  MCMC using Hamiltonian dynamics.  In Brooks, S., Gelman, A., Jones, G., and Meng, X. L., editors, Handbook of Markov Chain Monte Carlo, pages 113–162. Chapman and Hall/CRC.
  * Olander, (2020) Olander, J. (2020).  Constrained space mcmc methods for nested sampling bayesian computations.  Master’s thesis, CHALMERS UNIVERSITY OF TECHNOLOGY, Gothenburg, Sweden.
  * Olshanskii and Safin, (2016) Olshanskii, M. A. and Safin, D. (2016).  Numerical integration over implicitly defined domains for higher order unfitted finite element methods.  Lobachevskii Journal of Mathematics, 37(5):582–596.
  * Owen, (2013) Owen, A. B. (2013).  Monte Carlo theory, methods and examples.
  * Pakman and Paninski, (2014) Pakman, A. and Paninski, L. (2014).  Exact hamiltonian monte carlo for truncated multivariate gaussians.  Journal of Computational and Graphical Statistics, 23(2):518–542.
  * Park and Casella, (2008) Park, T. and Casella, G. (2008).  The bayesian lasso.  Journal of the American Statistical Association, 103(482):681–686.
  * Piepel et al., (2005) Piepel, G., Cooley, S., and Jones, B. (2005).  Construction of a 21-component layered mixture experiment design using a new mixture coordinate-exchange algorithm.  Quality Engineering, 17(4):579–594.
  * Pratola et al., (2017) Pratola, M. T., Harari, O., Bingham, D., and Flowers, G. E. (2017).  Design and analysis of experiments on nonconvex regions.  Technometrics, 59(1):36–47.
  * Saye, (2022) Saye, R. I. (2022).  High-order quadrature on multi-component domains implicitly defined by multivariate polynomials.  Journal of Computational Physics, 448:110720.
  * Shen et al., (2020) Shen, S., Kang, L., and Deng, X. (2020).  Additive heredity model for the analysis of mixture-of-mixtures experiments.  Technometrics, 62(2):265–276.
  * Sherlock and Roberts, (2009) Sherlock, C. and Roberts, G. O. (2009).  Optimal scaling of the random walk metropolis on elliptically symmetric unimodal targets.  Bernoulli, Vol.15(No.3):774–798.
  * Skilling, (2006) Skilling, J. (2006).  Nested sampling for general Bayesian computation.  Bayesian Analysis, 1(4):833 – 859.
  * Stinstra et al., (2003) Stinstra, E., den Hertog, D., Stehouwer, P., and Vestjens, A. (2003).  Constrained maximin designs for computer experiments.  Technometrics, 45(4):340–346.
  * (50) Tibshirani, R. (1996a).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288.
  * (51) Tibshirani, R. (1996b).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society, Series B, 58(1):267–288.
  * Verlet, (1967) Verlet, L. (1967).  Computer ”Experiments” on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules.  Phys. Rev., 159(1):98–103.
  * Wang, (2008) Wang, L. (2008).  Karhunen-Loeve expansions and their applications. PhD thesis, London School of Economics and Political Science (United Kingdom).