# Monte Carlo on manifolds in high dimensions

Kerun Xu1 and Miranda Holmes-Cerfon2

(1Courant Institute of Mathematical Sciences, New York University
2 Department of Mathematics, University of British Columbia and
Courant Institute of Mathematical Sciences, New York University.
1984 Mathematics Rd, Vancouver, BC, V6T 1Z2<EMAIL_ADDRESS>
)

###### Abstract

We introduce an efficient numerical implementation of a Markov Chain Monte
Carlo method to sample a probability distribution on a manifold (introduced
theoretically in Zappa, Holmes-Cerfon, Goodman (2018) [48]), where the
manifold is defined by the level set of constraint functions, and the
probability distribution may involve the pseudodeterminant of the Jacobian of
the constraints, as arises in physical sampling problems. The algorithm is
easy to implement and scales well to problems with thousands of dimensions and
with complex sets of constraints provided their Jacobian retains sparsity. The
algorithm uses direct linear algebra and requires a single matrix
factorization per proposal point, which enhances its efficiency over
previously proposed methods but becomes the computational bottleneck of the
algorithm in high dimensions. We test the algorithm on several examples
inspired by soft-matter physics and materials science to study its complexity
and properties.

##  1 Introduction

Equality constraints occur widely in physical and statistical sampling
problems. In a physical system, such as a collection of interacting particles,
constraints may be used to remove stiff degrees of freedom [42, 18, 19, 32].
For example, a stiff spring between particles may be replaced by a fixed
distance between them, which allows taking a longer step on the submanifold
that preserves this distance constraint. In statistics, constraints or
approximate constraints may arise as a condition on the set of parameters in a
parametric or Bayesian model [27, 10, 37, 17, 16], and they may accelerate
training of neural networks [31]. In algebraic geometry, constraints may be
used to study properties of algebraic varieties, such as surface or boundary
integrals [29], metric properties [20], or topological properties [4]. In all
such systems, equality constraints reduce the intrinsic dimension of the space
being sampled, so the probability measure of interest lives on a submanifold
of the ambient space.

Sampling a measure supported on a submanifold requires special consideration,
and a number of ideas and algorithms have been developed to simulate or sample
from systems subject to holonomic constraints. These go back to the Shake [42]
and Rattle [1] algorithms from molecular dynamics, which are still widely used
to this day, and have been extended to a variety of other systems such as
Langevin samplers [44, 7, 32, 30, 49, 43], Hamiltonian samplers [5, 33],
geodesic samplers [6], and Monte Carlo samplers [48]. Yet constraints don’t
come for free: since constraints are typically nonlinear functions of a
system’s variables, imposing these constraints requires solving a nonlinear
system of equations for the Lagrange multipliers at each step of the
algorithm, except for special manifolds for which an explicit parameterization
is available [10, 46]. Therefore, applying these methods to high-dimensional
systems requires carefully considering the numerical aspects of the
algorithms. Such consideration has been applied to algorithms in molecular
dynamics [2, 28, 47, 45, 15, 11], which may consider systems with $O(1000+)$
degrees of freedom, but has been less applied to Monte Carlo algorithms, where
a slightly different set of issues arise.

Our goal in this paper is to describe an efficient numerical implementation of
a Markov Chain Monte Carlo (MCMC) sampling method with constraints, which
allows sampling problems with $O(1000)$ dimensions while remaining easy to
implement using standard linear algebra packages. We directly build on [48],
which introduced an MCMC algorithm on manifolds but for which the
implementation described in that reference could not scale to high dimensions.
We focus on physically-motivated problems inspired by soft-matter physics and
materials science [34, 25, 3] where the constraints are approximations for
stiff forces, because these have an additional challenge associated with the
presence of an “effective” measure on the constraint surface, which can be
hard to calculate in high dimensions (and appears to be often neglected in
molecular dynamics simulations). However, our algorithm can be applied to
systems without this effective measure simply by not including it in the
Metropolis ratio. Apart from this effective measure, we mainly concentrate on
problems with no additional force field, since we wish to focus on the
numerical issues associated with handling constraints, and not on problems
such as evaluating forces that are common to both constrained and
unconstrained sampling problems. The implementation we describe requires a
single matrix factorization per proposal point, which we demonstrate is the
main computational bottleneck of the algorithm for complex constraints. The
algorithm is feasible for problem sizes for which this factorization is
possible.

A key idea that allows for just a single matrix factorization is to solve the
nonlinear constraint equations using a quasi-Newton method, by replacing the
exact Jacobian of these equations with a symmetric approximation that doesn’t
change with the iterations of a single Newton-like solve. Such an idea was
first introduced in [2] for molecular dynamics, and subsequently explored in
other molecular dynamics algorithms, such as Ref. [45], which used this idea
to develop a parallel implementation of SHAKE, Ref. [15], which used this idea
along with a preconditioner to speed up SHAKE, and Ref. [11], which proposed
explicitly inverting an approximate Jacobian and using this same inverse
throughout the simulation. Notably, while this idea in the context of
molecular dynamics leads to a typically modest decrease in total computational
time, generally 2-10% and at most 50% in the examples considered in [2], in
our examples it leads to a huge decrease in computational time, a minimum of
50% but up to a 98% decrease. We believe this is because the examples
considered in molecular dynamics have relatively simple constraints and more
complex force fields which are expensive to evaluate, whereas we consider
examples with complex constraints and simple force fields so our primary
computation is in handling the constraints.

Here is an overview of the paper. Section 2 introduces the setup, by
explaining the precise problem formulation, including how an effective measure
arises from delta-functions in a physical setting where constraints are
approximations for stiff forces. It also gives an overview of the sampling
algorithm. Section 3 describes the numerical implementation of the algorithm.
Section 4 applies the algorithm to four examples and analyzes its properties,
such as computational time and complexity, the nature of the rejection
statistics, the proportion of time occupied by matrix factorizations. It also
studies the effective diffusivity as a function of the average acceptance
ratio, finding that the optimal acceptance ratio depends on dimension for
examples with flatter probability distributions. Section 5 summarizes the
results and puts them in a broader context. Code which implements the examples
in this paper can be found at [23].

##  2 Setup and overview of the Monte Carlo algorithm

We consider a system described by $n$ variables, $x\in\mathbb{R}^{n}$, which
is subject to $m$ constraints of the form $q_{i}(x)=0$ $(i=1,\ldots,m)$, where
$q_{i}:\mathbb{R}^{n}\to\mathbb{R}$ are continuously differentiable functions,
and $m<n$. The system is therefore constrained to a set of the form

$M=\\{x\in\mathbb{R}^{n}:q_{i}(x)=0,i=1,\ldots,m\\}.$ (1)

We write $q(x)=(q_{1}(x),\ldots,q_{m}(x))^{T}$ to denote the vector-valued
constraint function.

The geometric nature of $M$ is in part determined by the Jacobian of the
constraints. Let

$Q_{x}=\begin{pmatrix}\nabla q_{1}(x)&\cdots&\nabla q_{m}(x)\end{pmatrix}$ (2)

be the $n\times m$ matrix whose columns are the gradients of the constraints
(taking gradients to be column vectors). This is the transpose of the Jacobian
$\nabla q(x)=\frac{\partial q}{\partial x}\Big{|}_{x}$; we refer to it simply
as the Jacobian. We make the following assumption:

$Q_{x}$ has full rank everywhere on $M$. (3)

This assumption implies, by the inverse function theorem, that $M$ is a
differentiable manifold, and its dimension is $d=n-m$. The tangent space to
$M$ at $x$, call it $T_{x}$, is the null space of $Q^{T}_{x}$, and the normal
space to $M$ at $x$, call it $N_{x}$, is spanned by the columns of $Q_{x}$.

We wish to sample a probability measure on $M$. We consider a measure that
arises naturally in physical modelling problems, which has the form

$\rho(dx)=Z^{-1}f(x)\prod_{i=1}^{m}\delta(q_{i}(x))dx.$ (4)

Here $\delta$ is a Dirac delta-function, $f:\mathbb{R}^{n}\to\mathbb{R}$ is
integrable in that sense that
$\int_{\mathbb{R}^{n}}|f(x)|\prod_{i=1}^{m}\delta(q_{i}(x))dx<\infty$, and $Z$
is a normalization constant, which is typically not known.

This measure arises from the Boltzmann distribution when there are stiff
forces that constrain a system to a neighbourhood of $M$. To see why, consider
a contribution to a system’s internal energy of the form
$\sum_{i=1}^{m}\frac{1}{2}kq_{i}^{2}(x)$, where $k$ is a parameter. This gives
a contribution to the partition function of the form
$e^{-\sum_{i=1}^{m}\frac{1}{2}kq_{i}^{2}(x)/k_{B}T}$, where $T$ is temperature
and $k_{B}$ is Boltzmann’s constant. If either $k\to\infty$ (stiff forces) or
equivalently $T\to 0$ (low temperature), this contribution approaches a
multiple of $\prod_{i=1}^{m}\delta(q_{i}(x))$ [12, 24]. For example, if
particles $x_{1},x_{2}\in\mathbb{R}^{d}$ are bound by a spring with rest
length $l>0$ and Hookean energy $\frac{1}{2}k(|x_{1}-x_{2}|-l)^{2}$, then as
$k\to\infty$ the Boltzmann distribution approaches a multiple of
$\delta(|x_{1}-x_{2}|-l)$. The stiff spring now acts as a distance constraint
$|x_{1}-x_{2}|-l=0$.

This distance constraint is however a _soft_ constraint – the measure it
induces depends on the particular function (or particular physical forces)
used to impose the constraint, via this function’s appearance in the argument
of the $\delta$-function. This measure is notably different from the _hard_
distance constraint one would obtain by simply restricting the Boltzmann
distribution to the manifold $M$ [13, 26, 21], which gives rise to a
restricted Boltzmann distribution proportional to the surface measure on $M$.

Our sampling algorithm is most easily formulated in terms of the natural
surface measure (Hausdorff measure) on $M$, which we call $\mu$. We may
convert (4) to a density with respect to $\mu$ using the coarea formula [36],
which is written formally as $\delta(q_{i}(x))dx=|Q_{x}|^{-1}\mu(dx)$. Here
$|Q_{x}|$ is the pseudodeterminant of $Q_{x}$, equal to the product of its
singular values or equivalently to $|Q^{T}_{x}Q_{x}|^{1/2}$, where $|\cdot|$
is the usual determinant of a square matrix. Therefore, the measure (4) we
wish to sample, written in terms of the surface measure, is

$\rho(dx)=Z^{-1}f(x)|Q_{x}|^{-1}\mu(dx).$ (5)

If one is interested in “hard” constraints, which don’t arise from delta-
functions but rather are directly imposed by restriction, then one can omit
the factor $|Q_{x}|^{-1}$ in the density above.

We will sample (5) using a Markov Chain Monte Carlo algorithm that was
proposed in [48]. The algorithm generates a Markov chain $X_{0},X_{1},\ldots$
whose stationary distribution is (5). Here is a high-level overview of the
algorithm.

Suppose $X_{k}=x\in M$. Repeat the following steps.

  1. 1.

Generate a random step in the tangent space, $v_{x}\in T_{x}$. Let the density
of this step with respect to the Euclidean measure on $T_{x}$ be
$p_{x}(v_{x})$.

  2. 2.

Project to manifold: solve $q(y)=0$ for the proposal $y\in M$. We require
$y=x+v_{x}+w_{x}$, where $w_{x}\in N_{x}$: the projection is normal to the
tangent space at the original point $x$. If the numerical solver fails to find
a solution, reject the proposal, and set $X_{k+1}=x$; repeat from Step 1.

  3. 3.

Solve for the reverse tangent step $v_{y}\in T_{y}$. For the reverse move,
from $y\to x$, we require $x=y+v_{y}+w_{y}$ where $w_{y}\in N_{y}$.

  4. 4.

Metropolis-Hastings rejection step. The acceptance probability is given by

$a_{xy}=\min\left(\frac{p_{y}(v_{y})f(y)|Q_{y}|^{-1}}{p_{x}(v_{x})f(x)|Q_{x}|^{-1}},1\right).$
(6)

One can omit the factors $|Q_{x}|^{-1},|Q_{y}|^{-1}$, if one is interested in
hard constraints. Use the acceptance probability above in a Metropolis-
Hastings rejection step: let $U\sim\mbox{Unif}([0,1])$. If $U>a_{xy}$, reject
the move, and set $X_{k+1}=x$; repeat from Step 1.

  5. 5.

Reverse check: determine whether the reverse move is feasible. Attempt to
solve for $w^{\prime}\in N_{y}$ such that $x=y+v_{y}+w^{\prime}$. If the
solver fails, or if it finds a solution such that
$x^{\prime}=y+v_{y}+w^{\prime}\neq x$, then reject the proposal, and set
$X_{k+1}=x$; repeat from Step 1.

  6. 6.

Accept proposal. If we made it here, the proposal should be accepted: set
$X_{k+1}=y$, and copy over Jacobian and Cholesky factorization,
$Q_{x}\leftarrow Q_{y},L_{x}\leftarrow L_{y}$. Repeat from Step 1.

The stationary distribution of this chain was shown in [48] to be (5).

##  3 Numerical implementation of the algorithm

This section describes an efficient numerical implementation of the algorithm
above to sample from (5). There are two primary numerical challenges. One is
in Step 2, the projection step, which requires solving a nonlinear system of
equations. The second is in Step 5, evaluating the acceptance probability,
which requires computing the pseudodeterminants $|Q_{x}|,|Q_{y}|$. We will
describe a method which efficiently addresses both of these issues. The method
will require only a single matrix factorization at each iteration of the
algorithm. This matrix factorization will turn out to be the primary
contribution to the computational time in high dimensions.

The matrix factorization required is the Cholesky decomposition of
$Q^{T}_{x}Q_{x}$. Let this decomposition be

$L_{x}L_{x}^{T}=Q^{T}_{x}Q_{x},$ (7)

where $L_{x}$ is a lower triangular matrix. This decomposition exists because
$Q^{T}_{x}Q_{x}$ is positive definite, since its null space is trivial by
assumption (3).

Here is how each of the steps 1-5 are implemented, in turn. Pseudocode of the
algorithm is shown in Algorithm 1 and 2.

Algorithm 1 ManifoldSampler: Given $x=X_{k}$, generate the next point
$X_{k+1}$ in the Markov chain

1:Parameters: $\sigma,\texttt{xtol}$

2:procedure ManifoldSampler($x=X_{n},Q_{x},L_{x}$)

3: Choose a random step $v_{x}$ in tangent space:

4: Generate $\xi_{i}\sim N(0,1)$ for $i=1,\ldots,n$ $\triangleright$
$\xi=(\xi_{i})_{i=1}^{n}$ is vector of i.i.d. normals

5: Solve $L_{x}L_{x}^{T}z=Q_{x}^{T}\xi$ for $z$ $\triangleright$ use forward-
backward substitution

6: $v_{x}\leftarrow\sigma\cdot(\xi-Q_{x}z)$ $\triangleright$ $v_{x}$ is
isotropic Gaussian in tangent space

7: Project to manifold: solve for proposal $y\in M$

8: $y,\texttt{newtonflag}\leftarrow$ Project($x+v_{x},Q_{x},L_{x}$)
$\triangleright$ See Algorithm 2

9: if newtonflag == fail then $\triangleright$ Projection failed to converge

10: Reject proposal: $X_{k+1}=X_{k}$. Return.

11: Find reverse tangent step $v_{y}$:

12: $Q_{y}\leftarrow\textsc{Jacobian}(q,y)$ $\triangleright$ Jacobian $\nabla
q(y)^{T}$ provided by user

13: $L_{y}\leftarrow\textsc{Cholesky}(Q_{y}^{T}Q_{y})$ $\triangleright$
Cholesky decomposition provided by user/package

14: $r\leftarrow x-y$

15: Solve $L_{y}L_{y}^{T}z=Q_{y}^{T}r$ for $z$

16: $v_{y}\leftarrow r-Q_{y}z$

17: Metropolis-Hastings Step:

18: $\texttt{vdiff}\leftarrow(|v_{y}|^{2}-|v_{x}|^{2})/(2\sigma^{2})$

19: $\texttt{udiff}\leftarrow\log(f)(y)-\log(f)(x)$ $\triangleright$ Function
$\log(f)$ provided by user

20: $\texttt{qdet}\leftarrow\prod_{i}l_{i}$, $l_{i}=(L_{x})_{ii}/(L_{y})_{ii}$
$\triangleright$ Set qdet=1 if you have hard constraints

21: $a\leftarrow\exp(-\texttt{vdiff})*\exp(-\texttt{udiff})*\texttt{qdet}$
$\triangleright$ Acceptance probability

22: if $U\sim\text{Unif}([0,1])>a$ then

23: Reject proposal: $X_{k+1}=X_{k}$. Return.

24: Reverse Projection:

25: $x^{\prime},\texttt{newtonflag}\leftarrow$
Project($y+v_{y},Q_{y},L_{y},v_{y}$)

26: if $\texttt{newtonflag}==\texttt{fail}$ then $\triangleright$ Reverse
projection failed to converge

27: Reject proposal: $X_{k+1}=X_{k}$. Return.

28: if $|x^{\prime}-x|>\texttt{xtol}$ then $\triangleright$ Projection
converged, but to wrong point

29: Reject proposal: $X_{k+1}=X_{k}$. Return.

30: Accept Proposal:

31: $x\leftarrow y$, $Q_{x}\leftarrow Q_{y}$, $L_{x}\leftarrow L_{y}$
$\triangleright$ Copy matrices/factorizations by exchanging pointers

32:end procedure

Algorithm 2 Project to manifold (symmetric Newton): Solve $q(y)=0$ s.t.
$y=z+Qa$

1:Parameters: $\texttt{tol},\eta,\texttt{MaxIter}$

2:procedure Project[SymmetricNewton]($z,Q,L$)

3: $a\leftarrow 0$-vector $\triangleright$ Initial guess (must be
deterministic)

4: $y\leftarrow z+Q*a$, $\texttt{qerr}\leftarrow\infty$ $\triangleright$
Initialization

5: for $\texttt{niter}=0:\texttt{MaxIter}-1$ do

6: $\texttt{qerr0}=\texttt{qerr}$ $\triangleright$ Save error from previous
step

7: $\texttt{qval}\leftarrow q(y)$ $\triangleright$ Evaluate constraints

8: $\texttt{qerr}\leftarrow|q(y)|_{\infty}$ $\triangleright$ Could use other
norms

9: if qerr < tol then $\triangleright$ We found a solution

10: return $y$, Success

11: if $\texttt{qerr}>\eta*\texttt{qerr0}$ then $\triangleright$ Error isn’t
decreasing fast enough; no solution

12: return $y$, Fail

13: Solve $LL^{t}\delta a=\texttt{qval}$ for $\delta a$ $\triangleright$
Forward-backward substitution

14: $a\leftarrow a+\delta a$, $y\leftarrow z+Q*a$ $\triangleright$ Next guess

15: end for

16: return $y$, Fail $\triangleright$ We didn’t find a solution to sufficient
accuracy

17:end procedure

Algorithm 3 Project to manifold (traditional Newton): Solve $q(y)=0$ s.t.
$y=z+Qa$

1:Parameters: $\texttt{tol},\eta,\texttt{MaxIter}$

2:procedure Project[TraditionalNewton]($z,Q$)

3: $a\leftarrow 0$-vector $\triangleright$ Initial guess (must be
deterministic)

4: $y\leftarrow z+Q*a$, $\texttt{qerr}\leftarrow\infty$ $\triangleright$
Initialization

5: for $\texttt{niter}=0:\texttt{MaxIter}-1$ do

6: $\texttt{qerr0}=\texttt{qerr}$ $\triangleright$ Save error from previous
step

7: $\texttt{qval}\leftarrow q(y)$ $\triangleright$ Evaluate constraints

8: $\texttt{qerr}\leftarrow|q(y)|_{\infty}$ $\triangleright$ Could use other
norms

9: if qerr < tol then $\triangleright$ We found a solution

10: return $y$, Success

11: if $\texttt{qerr}>\eta*\texttt{qerr0}$ then $\triangleright$ Error isn’t
decreasing fast enough; no solution

12: return $y$, Fail

13: $Q_{y}\leftarrow\textsc{Jacobian}(q,y)$

14: $L,U\leftarrow\textsc{LU}(Q_{y}^{T}Q)$ $\triangleright$ LU decomposition
provided by user/package

15: Solve $LU\delta a=\texttt{qval}$ for $\delta a$ $\triangleright$ Forward-
backward substitution

16: $a\leftarrow a+\delta a$, $y\leftarrow z+Q*a$ $\triangleright$ Next guess

17: end for

18: return $y$, Fail $\triangleright$ We didn’t find a solution to sufficient
accuracy

19:end procedure

###  3.1 Generate a random step in the tangent space $v_{x}\in T_{x}$

We generate an isotropic Gaussian in the tangent space. Let
$\xi\in\mathbb{R}^{n}$ be a vector of i.i.d. standard normals. We find $v_{x}$
by projecting $\xi$ orthogonally onto the tangent space, and then scaling by a
parameter $\sigma>0$ which controls the overall size of the step. The
orthogonal projection onto the tangent space is
$(I-Q_{x}(Q_{x}Q_{x}^{T})^{-1}Q_{x}^{T})\xi$. We implement this projection
without matrix inversion, by solving

$L_{x}L_{x}^{T}z=Q^{T}_{x}\xi\qquad\text{for $z$},$ (8)

using the Cholesky decomposition (7) to solve the system of equations by
forward-backward substitution. Then we set

$v_{x}=\sigma(\xi-Q_{x}z).$ (9)

The density of this tangent space proposal is

$p_{x}(v_{x})=\frac{1}{(2\pi\sigma^{2})^{d/2}}e^{-\frac{|v_{x}|^{2}}{2\sigma^{2}}}.$
(10)

In [48] it was suggested to choose $v_{x}$ by first computing an orthonormal
basis of the tangent space $T_{x}$ using the QR decomposition of $Q_{x}$ (from
the last $d$ columns of the Q-matrix), storing the basis in the columns of a
matrix $E_{x}$, and then setting $v_{x}=\sigma E_{x}\xi^{\prime}$ where
$\xi^{\prime}\in\mathbb{R}^{d}$ is a vector of i.i.d. standard normals. In
high dimensions, a similar approach could be implemented using a sparse QR
decomposition of $Q_{x}$, which represents the Q-matrix as a product of
Householder transformations or Givens rotations [9], thereby producing a basis
of the normal space $N_{x}$ which can be used for projection. However, we
expect such an approach to be less efficient than the one proposed above,
partly because the QR decomposition is more expensive to compute than the
Cholesky decomposition, and partly because we will require the Cholesky
decomposition in later steps of the algorithm in places where it cannot be
substituted with the QR decomposition.

We remark that constrained molecular dynamics simulations (e.g. [2]) generate
an initial proposal in the full, ambient space, equivalent to our random
vector $\xi$, and then project to $M$ to find $y=x+\sigma\xi+w$, thereby
avoiding the need to calculate the tangent step. It is not possible to
Metropolize such a proposal without additional processing, because for a given
$x,y$, there is a positive-dimensional space of possible $\xi,w$, hence, we
cannot determine the probability density of $\xi$ which is required to
evaluate the acceptance probability.

###  3.2 Project to manifold: Solve for the proposal $y\in M$

The next step is to solve for $y$ such that $q(y)=0$, with $y=x+v_{x}+w_{x}$
and $w_{x}\in N_{x}$ is the vector to be determined. Observe that we can write
any $w\in N_{x}$ as $w=\sum_{i=1}^{m}a_{i}\nabla q_{i}(x)$ for a unique choice
of $a=(a_{1},\ldots,a_{m})^{T}\in\mathbb{R}^{m}$, by assumption (3).
Therefore, we must solve the system of equations

$q_{i}\left(x+v_{x}+\sum_{i=1}^{m}a_{i}\nabla q_{i}(x)\right)=0,\qquad
i=1,\ldots,m,$ (11)

for the unknown $a$.

We propose to solve this system using a quasi-Newton method that replaces the
exact Jacobian with an approximate Jacobian for which we already have a
factorization; we call this symmetric Newton, following [2]. We compare this
method with Newton’s method, which we refer to as traditional Newton. One
contribution of this paper is to demonstrate that symmetric Newton
_significantly_ increases the overall efficiency of the algorithm compared to
traditional Newton.

Let us first describe traditional Newton’s method to solve (11). (Pseudocode
is found in Algorithm 3.) The method generates a sequence of approximations
$a^{(0)},a^{(1)},a^{(2)},\ldots$ with corresponding estimates of $y$ computed
as $y^{(n)}=x+v_{x}+Q_{x}a^{(n)}$. The Jacobian of (11) is
$Q_{y^{(n)}}^{T}Q_{x}$, so the iterates are updated as
$a^{(n+1)}=a^{(n)}+\delta a^{(n)}$, where $\delta a^{(n)}$ solves

$Q_{y^{(n)}}^{T}Q_{x}\delta a^{(n)}=-q(y^{(n)}).$ (12)

We solve the above linear system of equations using an LU decomposition of
$Q_{y^{(n)}}^{T}Q_{x}$, and initial condition $a^{(0)}=0$, the zero vector.
The method iterates until a stopping criterion is satisfied; we discuss the
stopping criterion at the end of this section. Notice that at each iteration,
we must evaluate $Q_{y^{(n)}}$ and compute an LU decomposition of
$Q_{y^{(n)}}^{T}Q_{x}$, both of which are expensive operations in high
dimensions.

Our second method, symmetric Newton, replaces the exact Jacobian
$Q_{y^{(n)}}^{T}Q_{x}$ with an approximation $Q_{x}^{T}Q_{x}$. (Pseudocode is
found in Algorithm 2.) We already have its Cholesky factorization (7) which we
used to generate a tangent step in Section 3.1, so there is no further
Jacobian evaluation nor matrix factorization required. Therefore, symmetric
Newton computes updates by solving

$Q_{x}^{T}Q_{x}\delta a^{(n)}=-q(y^{(n)}),$ (13)

using the Cholesky decomposition (7) to solve this linear system of equations
by forward-backward substitution.

We expect symmetric Newton to converge to the solution more slowly than
traditional Newton, since it should converge linearly to the true solution
whereas traditional Newton converges quadratically [38]. It could also fail to
converge in some cases even when Newton’s method converges, leading to an
increase in rejections. However, we also expect an increase in efficiency
gained by avoiding a Jacobian evaluation and matrix factorization at each
iteration of the solver. It is not immediately clear how this tradeoff will
play out; we will determine this via examples.

We turn to discussing the stopping criterion, which is the same for both
Newton methods. The stopping criterion is important because it is inefficient
to spend a lot of computational effort searching for a solution, when we could
instead give up easily and try again in a different random direction.
Furthermore, there is sometimes no solution to (11), and indeed the sampler
can be more efficient when operated in a regime where there is frequently no
solution [33]; it is important not to waste time searching for a solution that
doesn’t exist.

We have found the following stopping criterion to be effective. Terminate
Newton’s method when any one of the following is satisfied:

  1. (i)

$\texttt{err}_{n}<\texttt{tol}$, where tol is a tolerance parameter
determining how close we are to the manifold, and
$\texttt{err}_{n}=|q(y^{(n)})|_{\infty}=\max_{i}|q_{i}(y^{(n)})|$ is a measure
of error. When this happens we return a _success_ : we have found a solution
sufficiently close to $M$.

  2. (ii)

$\texttt{err}_{n}>\eta\texttt{err}_{n-1}$, where $\eta$ is a parameter
governing how much the error must shrink at each step (typically $\eta<1$). If
this criterion is satisfied, we return a _fail_.

  3. (iii)

$n\geq\texttt{MaxIter}$, where MaxIter is a parameter controlling the maximum
number of iterations. If this criterion is satisfied, we return a _fail_.
Because of the previous stopping criterion, we rarely (if ever) reach this
step.

We explored other stopping criteria that would allow us to spend more time
looking for a solution. For example, we explored a method where if case (ii)
occurs, then instead of stopping, we tried shrinking the Newton step size by
half, i.e. we set $\delta a^{(n)}\to\delta a^{(n)}/2$, and re-computed the
error, and re-evaluated this stopping criterion; repeating the shrinking step
several times if necessary. We found this increased the number of acceptances
by a negligible amount, so we dropped this extra method.

###  3.3 Solve for the reverse tangent step $v_{y}\in T_{y}$

The reverse move proposes $x=y+v_{y}+w_{y}$ with $v_{y}\in T_{y}$ and
$w_{y}\in N_{y}$. Since, $v_{y},w_{y}$ live in orthogonal subspaces, they may
each be determined by projecting $x-y$ orthogonally onto the appropriate
subspace. Therefore, we find $v_{y}$ using the same steps as in (8), (9):

$\text{Solve }\quad L_{y}L_{y}^{T}z^{\prime}=Q_{y}^{T}(x-y)\quad\text{ for
$z^{\prime}$;}\qquad\text{set}\quad v_{y}=(x-y)-Q_{y}z^{\prime}.$ (14)

To do this we must first evaluate $Q_{y}$ and compute the Cholesky
decomposition (7) of $Q_{y}Q_{y}^{T}$. In practice, at each iteration of the
Monte-Carlo algorithm, _this is the only step where a Jacobian evaluation and
a Cholesky decomposition is computed_ – we have the Cholesky decomposition of
$Q_{x}Q_{x}^{T}$ because we computed it when $x$ was first proposed.

###  3.4 Metropolis-Hastings rejection step

We test whether the acceptance criterion is satisfied _before_ we perform the
reverse check (which is also necessary for the acceptance criterion to hold),
because it is more efficient in this order: if we reject a move because of the
Metropolis criterion, which compares $f,v_{(\cdot)},|Q_{(\cdot)}|^{-1}$, then
we do not need to proceed to a reverse projection, which is computationally
much more expensive.

Evaluating the acceptance probability in high dimensions requires a few
remarks on implementation. The main challenge is that the terms in (6) become
increasingly tiny as the dimension is increased, resulting in numerical
underflow. The solution is generally to work with the logarithms of the
quantities involved. That is, to compute $f(y)/f(x)$ we provide a function
which returns $\log f(x)$, then exponentiate the difference $\log f(y)-\log
f(x)$. For the ratio
$p_{y}(v_{y})/p_{x}(v_{x})=e^{-|v_{y}|^{2}/2\sigma^{2}}/e^{-|v_{x}|^{2}/2\sigma^{2}}$
(see (10)) we first compute $-(|v_{y}|^{2}-|v_{x}|^{2})/2\sigma^{2}$, and then
exponentiate.

The final term to evaluate is $|Q_{y}|^{-1}/|Q_{x}|^{-1}$. These
pseudodeterminants may be expressed in terms of the Cholesky factorizations
(7) as

$|Q_{x}|=|Q_{x}Q_{x}^{T}|^{1/2}=\prod_{i=1}^{m}(L_{x})_{ii},\qquad|Q_{y}|=|Q_{y}Q_{y}^{T}|^{1/2}=\prod_{i=1}^{m}(L_{y})_{ii}.$
(15)

In high-dimensional problems where these determinants result in numerical
underflow, we evaluate the contribution to the acceptance probability by
computing the ratio of individual diagonal elements of the Cholesky matrices,
and then multiplying these ratios together:

$r_{i}=\frac{(L_{x})_{ii}}{(L_{y})_{ii}},\qquad\quad\frac{|Q_{y}|^{-1}}{|Q_{x}|^{-1}}=\prod_{i=1}^{m}r_{i}.$
(16)

This is effective as the corresponding diagonal elements $(L_{x})_{ii}$,
$(L_{y})_{ii}$, tend to be similar to each other.

We could alternatively compute
$\log|Q_{(\cdot)}|=\sum_{i}\log(L_{(\cdot)})_{ii}$ and then compute the ratio
of determinants as
$|Q_{y}|^{-1}/|Q_{x}|^{-1}=\exp\left(\log|Q_{x}|-\log|Q_{y}|\right)$. We
verified that this gives the same ratio of determinants and avoids numerical
underflow, however we prefer (16) as it avoids computing logs and
exponentials.

###  3.5 Reverse check: determine whether the reverse move is feasible

The nonlinear solver in Section 3.2 does not always produce a solution to
(11), even if one exists. Therefore, it is possible that the solver produces a
proposal $y\in M$, but then when asked to solve for $x^{\prime}\in M$ such
that $x^{\prime}=y+v_{y}+w^{\prime}$, it fails to find a solution. If this
happens we must reject the proposal, to satisfy detailed balance. (See [48]
for an enhanced discussion of this issue.)

We check this possibility by attempting to solve (11) for $x^{\prime}$,
starting from $y+v_{y}$, using exactly the same solver as for the forward
move: either traditional Newton, or symmetric Newton. If the solver fails, or
if it finds a solution $x^{\prime}\neq x$ (in practice, such that
$|x^{\prime}-x|>\texttt{xtol}$, where xtol is a tolerance parameter), we
reject the proposal.

This additional nonlinear solve increases the computational time. Whether or
not it this additional time is worth it, depends on the problem. Ref. [33]
found that in their Hamiltonian sampler, for the low-dimensional examples they
considered, their sampler was more efficient when operated in a regime where
there was a significant fraction (up to 15%) of moves rejected during this
step. In 3 out orfour 4 examples, we find a negligible percentage of moves are
rejected during this step (see Section 4.3), except occasionally in the
lowest-dimensional versions of each example. Therefore, in many cases, this
extra projection check could be dropped, at the expense of introducing a
slight bias in the measure being sampled.

##  4 Numerical Examples and Results

We apply our algorithm to four different examples, described in Section 4.1.
We aim to demonstrate the algorithm is effective in systems with $O(1000)$
variables, in the sense that it allows generating $O(10^{6})$ points or more
on reasonable timescales, e.g. no more than a timescale of days. We make no
claims about the _convergence_ of a Monte Carlo method on these examples
(which would require analyzing for example the correlation time of statistics
of interest), as this depends sensitively on each example and introduces
considerations beyond the numerical ones considered in this paper. We did
verify that the sampler produces points with the correct distribution, on some
of the examples considered in [48, 22] for which a nontrivial analytic
marginal distribution is known.

We further wish to compare the two projection methods, and show that the
symmetric Newton projection method significantly increases the efficiency of
the sampler compared to traditional Newton, by a factor of 3-23 depending on
the example. In some examples, this makes an otherwise intractable problem
tractable. We show that a primary reason for the increase in efficiency is it
reduces the number of matrix factorizations.

We make a few remarks on implementation. We implement our algorithm in C++. We
use Eigen for linear algebra operations, an open-source library for linear
algebra computations, that is comparable to other high-performance linear
algebra libraries. We consider problems where $Q_{x}$ is sparse. If it isn’t,
then constraint-based methods will be inefficient. We use AMD ordering to
compute a Cholesky desompotion, and COLAMD ordering for LU decompositions,
which we found empirically to be slightly faster than AMD (Eigen does not have
COLAMD ordering for Cholesky decompositions). The algorithm is compiled and
run on a 2.2GHz Intel Xeon processor. These implementation choices will affect
the timings reported, but are not expected to affect their qualitative
behaviour; furthermore we provide heuristics for estimating timings for other
implementations.

We chose parameters tol = $10^{-5}$, $\eta=0.95$, MaxIter = 100 for the
projection step, and $\texttt{xtol}=\texttt{tol}*n*10$. The maximum number of
iterations was rarely (if ever) reached, because of the other termination
criteria.

Code which implements all of the examples below and several other simpler ones
can be found at [23].

###  4.1 Examples

We test our algorithm on four examples. Each has a parameter, $n$, controlling
the size of the problem; in most examples $n$ represents a number of particles
or vertices (this is different from how we used the variable $n$ in the
previous sections). The number of variables $n_{\rm vars}$ and the number of
constraints $m$ are each proportional to $n$. The examples explore different
regimes for the fraction of constraints $m/n_{\rm vars}$, both small to large,
and different patterns for the constraints, leading to different sparsity
patterns in $Q_{x}$. In all but one example we set $f(x)=1$ to focus purely on
the geometrical aspects of the sampler. Problems where $f$ is significant and
dominate the computational time were explored in [2].

Figure 1: Illustrations of selected examples. Left: Polymer (Example 4.1.1)
with $n=12$ particles, with fixed endpoints shown in red. Middle: Lattice
(Example 4.1.2) with $n=16$ particles. Right: Ngon with extra random edges
(Example 4.1.4) with $n=12$ particles.

####  4.1.1 Polymer with fixed endpoints

We consider a three-dimensional polymer with $n$ vertices, bound by bars of
length 1, whose endpoints are bound with bars of length 1 to fixed vertices at
$a_{0}=(0,0,0)$ and $a_{1}=(n/2,0,0)$; see Figure 1. This models a polymer of
length $n+2$ with fixed endpoints, whose end-to-end length is fixed at roughly
half the total length of the polymer. We remark that a polymer without
endpoint constraints is straightforward to parameterize and sample directly,
while adding the endpoint constraints makes sampling significantly more
challenging and requires specialized methods [39]. We have

$n_{\rm vars}=3n,\qquad m=n+1.$

In this example the number of constraints is reasonably small, $m/n_{\rm
vars}\approx 1/3$. For point particles in three dimensions, this is the
smallest scaling for the number of distances constraints that allows for a
connected set of vertices.

The system is described by the positions of the vertices,
$x=(x_{1},\ldots,x_{n})^{T}\in\mathbb{R}^{3n}$ where $x_{i}\in\mathbb{R}^{3}$
is the position of the $i$th vertex. The bars are implemented with constraints
of the form

$q_{k}(x)=|x_{k}-x_{k+1}|^{2}-1,\qquad k=1,\ldots,m-1,$ (17)

and for the endpoints we have $q_{m}(x)=|x_{1}-a_{0}|^{2}-1$,
$q_{m+1}(x)=|x_{n}-a_{1}|^{2}-1$. We set $f(x)=1$.

We remark that although it might be more physical to sample a constraint of
the form $q_{i}(x)=|x_{i}-x_{i+1}|-1$, since this is the constraint that would
arise if the bar were a stiff spring, one can verify that $|Q_{x}|$ varies by
a constant between the two cases, so it doesn’t alter the acceptance
probability. Therefore, we choose the constraint whose gradient is less
expensive to evaluate.

####  4.1.2 Lattice

We consider a two-dimensional square lattice with $n$ vertices bound by bars
of unit length, modelling a crystalline material, see Figure 1. We must choose
$n=s^{2}$ where $s$ in an integer, representing the number of vertices on each
side of the lattice. We have

$n_{\rm vars}=2n,\qquad m=2n-2\sqrt{n}.$

This example has close to the maximum number of constraints, $m/n_{\rm
vars}=1-1/\sqrt{n}\approx 1$.

The system is described by $x=(x_{1},\ldots,x_{n})^{T}\in\mathbb{R}^{2n}$
where $x_{i}\in\mathbb{R}^{2}$ is the position of the $i$th vertex. The
constraints have the form

$q_{k}(x)=|x_{i_{k}}-x_{j_{k}}|^{2}-1,$ (18)

where $k$ indexes the constraints and where $(i_{k},j_{k})$ give the vertices
associated with the $k$th edge.

We found that without imposing an additional energy function on the lattice,
it collapsed in our simulations to a configuration where one line of vertices
was nearly on top of another. The sampler was unable to successfully make a
proposal $y$ from this configuration, so it remained in that configuration
thereafter. Such a configuration is geometrically “singular”, with
$|Q_{x}|\approx 0$, violating our assumption (3).

Therefore we impose an energy on the lattice to prevent it from collapsing.
This is the only example we consider with a nontrivial energy, so this energy
also serves to explore a regime where the energy function is nontrivial. We
found the following function to be effective, allowing the lattice to deform
slightly yet not collapse:

$f(x)=\exp(-k\sum_{i}(d_{i}-\sqrt{n}/\sqrt{2})^{2}).$ (19)

Here $d_{i}$ is the length of a diagonal for a small unit block and $k=5$ is a
parameter. The sum indexed by $i$ runs over all the diagonals in both
directions of the lattice.

####  4.1.3 Orthonormal matrix

We consider the special orthogonal group of $n=s\times s$ orthonormal matrices
with positive determinant, $SO(s)=\\{A\in\mathbb{R}^{s\times
s}:A^{T}A=I,|A|=1\\}$. Each matrix can be reshaped into a vector with
$n=s^{2}$ entries. The set $SO(s)$ is characterized by a condition on the dot
products of its rows, which form our constraints:

$A_{i,\cdot}^{T}A_{i,\cdot}-1=0,\qquad A_{i,\cdot}^{T}A_{j,\cdot}=0\quad(i\neq
j).$

Here $A_{i,\cdot}$ is the $i$th row of $A$. We have

$n_{\rm vars}=n,\qquad m=\frac{n+\sqrt{n}}{2}.$

Here the number of constraints is $m/n_{\rm vars}\approx 1/2$.

We wish to sample points in proportion to the natural surface measure (the
Haar measure) on this space. It can be shown that $|Q_{x}|$ is constant on
$SO(s)$ for the constraints as written above, so we keep this factor in our
sampling as it makes no difference to the acceptance probability. We also note
that one can show that $|v_{x}|=|v_{y}|$ for all $x,y\in SO(s)$. Therefore,
with $f(x)=1$, the acceptance probability should be identically $1$, and the
only way to reject a proposal should be if the projection step fails to
converge.

We remark that our sampler is not the best way to sample from the uniform
distribution on $SO(s)$ [35]. We include this example merely to test our
sampler with an additional type of constraint structure.

####  4.1.4 N-gon with random vertices and edges

Our next example considers random graphs, chosen to have a fairly large number
of constraints, and such that there is no structure in the sparsity pattern in
$Q_{x}$. We consider $n$ vertices in three-dimensional space, which are
initially equally spaced along a circle in the $xy$-plane such that the
distance between neighbouring points on the circle is 1. The vertices are then
connected by $n$ edges to form a polygon, and then an additional $n$ edges are
chosen randomly from the remaining pairs. All edges represent distance
constraints. The vertices are then perturbed: the vertical ($z$) coordinates
are set to be independent Gaussian random variables with standard deviation
0.5, while the horizontal radial coordinate is multiplied by an independent
random variable that is uniform on $[0.6,1]$. (These perturbations allow the
graph to have some flexibility; the particular choices of perturbations are
not important.) The lengths of all edges are then computed from this
configuration. See Figure 1 for one such random graph.

For this system

$n_{\rm vars}=3n,\qquad m=2n,$

and hence $m/n_{\rm vars}=2/3$. Because each realization of the example is
different, we typically consider 3 random examples for each value of $n$. We
always use $f(x)=1$.

###  4.2 Comparing projection methods: Traditional Newton versus Symmetric
Newton

We are interested in comparing the efficiency of the algorithm for each type
of projection algorithm, traditional Newton and symmetric Newton. Each method
has a parameter $\sigma$ governing the size of the tangent step, and any
measure of efficiency will be sensitive to this parameter. Comparing the
methods therefore requires choosing a value of $\sigma$ for each method that
allows for a fair comparison. Ideally we should choose $\sigma$ to minimize
the computational time required to generate an “effectively independent”
sample point. However, because measuring the level of convergence of a Markov
chain is computationally not straightforward [8, 40], a common practice is to
choose a target average acceptance ratio $a$ and tune the step size parameter
$\sigma$ to achieve this ratio. Theory on idealized problems suggests that
$a\approx 0.234$ leads to optimal correlation times in high dimensions [14,
41].

We follow this same practice, and determine for each example and each method
the value of $\sigma$ which leads to a desired average acceptance ratio $a$;
call these values $\sigma^{(N)}_{a},\sigma^{(S)}_{a}$ for traditional Newton
and symmetric Newton respectively. For most of our tests we choose $a=0.25$;
subsequent tests (Section 4.5) show that changing the value of $a$ doesn’t
alter our qualitative observations.

Figure 2: Estimated values of $\sigma_{0.25}$, the step-size parameter that
leads to average acceptance probability $a=0.25$, as a function of system size
parameter $n$, for our four examples. The data for the N-gon is an average
over 3 independent examples at each $n$, with $\sigma_{0.25}$ determined
separately for each example, and is shown with one-standard-deviation error
bars.

Figure 2 shows the estimated values of
$\sigma^{(N)}_{0.25},\sigma^{(S)}_{0.25}$ for each example at different system
sizes $n$. These were estimated using a bisection algorithm that repeats until
the value of $a$ estimated using $10^{5}$ sampled points is within
$1.25/\sqrt{10^{5}}$ of the target value $a$. The values of
$\sigma^{(N)}_{0.25},\sigma^{(S)}_{0.25}$ are generally very close. A slight
exception is the matrix (Example 4.1.3), where $\sigma^{(N)}_{0.25}$ is about
8-10% larger than $\sigma^{(S)}_{0.52}$ (the difference decreases slightly
with $n$). We would expect this difference to be associated with a
$(\sigma^{(N)}_{0.25}/\sigma^{(S)}_{0.25})^{2}\approx 16-20$% increase in
overall efficiency for traditional Newton, assuming diffusive dynamics. As a
side remark, the step sizes decrease with $n$ as $\sim
n^{-0.2},n^{-0.5},n^{-0.3}$ for the Polymer, Lattice, and Matrix respectively.
They remain roughly constant for the Ngon, presumably because the geometry of
the example changes with dimension.

Figure 3: Time (in seconds) to sample $10^{6}$ points for each example, using step-size parameter $\sigma_{0.25}$, for each of traditional Newton and symmetric Newton. The ratio (traditional:symmetric) is shown in grey on the right axes. Data for the Ngon is an average over 3 independent examples; the one-standard deviation error bars are smaller than the markers.  | Polygon | Lattice | Matrix | Ngon
---|---|---|---|---
Newton | 0.98 | 1.27 | 2.24 | 2.35
Symmetric | 0.99 | 1.21 | 2.16 | 2.08
LU | 0.97 | 1.36 | 2.35 | 2.48
Cholesky | 0.97 | 1.56 | 2.79 | 2.75
Table 1: Complexity of each method or matrix factorization, reported as the
power $\alpha$ in the scaling law $\text{time}\sim O(n^{\alpha})$, obtained by
using least-squares to fit our empirical timing data (minus the 3 smallest
values of $n$) on a log-log scale.

We then measured the total computational time to generate $10^{6}$ points for
each example at several values of $n$, using the estimated step sizes
$\sigma^{(N)}_{0.25},\sigma^{(S)}_{0.25}$. Figure 3 shows that symmetric
Newton is _uniformly_ faster than traditional Newton, by a factor that depends
on the example: for the Polymer it is roughly 2.5-3 times faster; for the
Lattice it is roughly 5 times faster; for the Matrix it is roughly 7-9 times
faster (though accounting for the slightly greater step size for traditional
Newton, the actual speedup is about a factor of 6-8). The Ngon shows the
greater speedup, and the strongest dependence on $n$, increasing from an
average of about 8 times faster for $n=12$ to more than 23 times faster for
$n=750$. Notably, this increase changes the order of magnitude of the
computations: at $n=750$, symmetric Newton took just over $6$ _hours_ while
traditional Newton took just over $6$ _days_.

The computational complexities of traditional and symmetric Newton are similar
(Table 1 ), and depend on the example, being smallest for the Polymer, roughly
$O(n)$, and largest for the Matrix and Ngon, roughly $O(n^{2.1})$. The
complexity for traditional Newton is slightly higher than for symmetric
Newton. The complexities are similar to the complexities of performing sparse
matrix factorizations for each example (Section 4.4).

We remark that we verified that the actual measured acceptance ratio $\hat{a}$
for the $10^{6}$ points was close to the desired one, in the range
$\hat{a}\in[0.24,0.26]$. For the Ngon the actual acceptance ratio sometimes
varied more than this. We suspect this could be because the random geometry
caused the configuration space of the N-gon to contain some regions allowing
for different (usually larger) step sizes, but that reaching these regions was
a rare event requiring a long burn-in time. When $\hat{a}$ was too different
from $a$, we re-ran the bisection algorithm using more sample points and a
longer burn-in time.

###  4.3 Rejection statistics and projection diagnostics

Figure 4:  Rejection statistics for each example and each method (traditional
Newton / symmetric Newton), as a function of $n$. Proposals may be rejected
because the forward projection fails (Project), because the Metropolis-
Hastings criterion is not satisfied (Metropolis), or because the reverse
projection fails or produces a point different from the starting point
(Reverse). The Polymer, Lattice, and Ngon have essentially 0 Reverse
rejections; for these examples traditional and symmetric Newton have virtually
identical statistics of each type. The Lattice has 0 Project rejections. Data
for the Ngon is an average over 3 independent examples; the one-standard
deviation error bars are smaller than the markers.

We next examine the rejection statistics and projection diagnostics, to see if
there are notable differences between the methods and examples. Figure 4 shows
the fraction of each type of rejection, for each method: either the projection
step failed to converge, or it was rejected in the Metropolis-Hastings step,
or the reverse projection failed to converge. In most examples, at smaller
$n$, rejections are mainly due to projections, while at larger $n$, Metropolis
rejections become more important. The Lattice was an exception, with only
Metropolis rejections and no Projection rejections; this is because $\sigma$
was small to accommodate a strongly varying energy.

It is notable that the Polymer, Lattice, and Ngon, have identical rejection
statistics for traditional Newton and symmetric Newton. The Matrix, on the
other hand, has slightly different rejection statistics for traditional Newton
and symmetric Newton – traditional Newton has mostly projection rejections,
with a small number of Metropolis rejections, while symmetric Newton has fewer
rejections in the forward projection step, and some rejections in the reverse
projection step (an amount that increases with $n$), with virtually no
Metropolis rejections. Therefore, symmetric Newton in this example achieves
the same overall acceptance rate $a=0.25$, by distributing its rejections
between forward and reverse projections.

In theory, the Matrix should have no rejections in the Metropolis step, as was
discussed in Section 4.1.3. The reason there are Metropolis rejections in
practice for traditional Newton, is that we are not generating samples which
satisfy the constraints exactly, so the tangent steps $v_{x},v_{y}$ do not
have exactly the same norm. We verified that decreasing the tolerance tol in
Newton’s method decreases the number of Metropolis rejections for traditional
Newton. It is interesting that symmetric Newton does not have a significant
number of Metropolis rejections – we hypothesize this is because proposals
where $|v_{x}|,|v_{y}|$ are more different are also more likely to not be
found by the sloppier symmetric Newton projection.

In all examples except the Matrix with symmetric Newton, the fraction of
rejections in the reverse projection step was tiny – less than 0.3% and
decreasing with $n$ for nearly all examples except the Ngon, for which the
fraction was 2.7%, 1.5%, 0.5% for $n=12,24,48$ respectively, and negligible
thereafter. Therefore, it appears that the reverse projection step could be
skipped in these examples without significantly affecting the stationary
distribution, particularly at large $n$.

Figure 5: Average number of iterations $n_{\rm iter}$ required for the
projection method to converge, for traditional Newton and symmetric Newton,
for different examples. The ratio symmetric:traditional is shown in grey on
the right. Also shown are the average values of $n_{\rm iter}$ when the
projection succeeds, and when it fails. Data for the Ngon is an average over 3
independent examples; the one-standard deviation error bars are generally
smaller than the markers.

Since the rejection statistics are nearly identical for traditional and
symmetric Newton for most examples, the algorithms are performing the same
average amount of computation, except during the projection step. Let’s study
how the algorithms compare during this step.

In each projection step, each algorithm performs a number $n_{\rm iter}$ of
Newton or quasi-Newton iterations, during which there are various costly
linear algebra operations. Figure 5 shows the average $n_{\rm iter}$ for each
method and example. Symmetric Newton generally requires more iterations than
traditional Newton, but notably, not that many more: typically a factor of 3-5
times more, depending on the example. An exception is the Lattice at large $n$
where $n_{\rm iter}\approx 1$ for both methods, presumably because $\sigma$ is
so small. Figure 5 also shows the average values of $n_{\rm iter}$ conditional
on a projection succeeding and conditional on it failing. For traditional
Newton these conditional values are similar, with a very slightly higher value
for successful projections. For symmetric Newton, which conditional value is
higher depends on $n$ and the example; though the difference is quite striking
for the Ngon, where unsuccessful projections require far fewer iterations.
This suggests our termination criterion for Newton’s method is effective; we
are not wasting time with projections that will not succeed.

###  4.4 Matrix factorizations explain why symmetric Newton is faster than
traditional Newton

Figure 6: Times to generate 1 point with each method, $T_{{\rm new}},T_{{\rm sym}}$, and estimated times based on matrix factorizations, $T_{{\rm new}}^{{\rm est}},T_{{\rm sym}}^{{\rm est}}$ from (20).  | Newton | Symmetric
---|---|---
| small $n$ | large $n$ | small $n$ | large $n$
Polymer | 56% ($n{=5}$) | 40% (n=10,240) | 4.6% ($n{=5}$) | 3.9% ($n{=}10,240$)
Lattice | 68% ($n{=}36$) | 99% ($n{=}1225$) | 15% ($n{=}36$) | 43% ($n{=}1225$)
Matrix | 40% ($n{=}36$) | 83% ($n{=}$1225) | 4.8% ($n{=}$36) | 39% ($n{=}1225$)
Ngon | 56% ($n{=}24$) | 104% ($n{=}740$) | 5.1% ($n{=}24$) | 57% ($n{=}750$)
Table 2: Percentage of total computational time explained by matrix
factorizations, $T_{{\rm new}}^{{\rm est}}/T_{{\rm new}}$, $T_{{\rm
sym}}^{{\rm est}}/T_{{\rm sym}}$, for small $n$ and large $n$, for each
example. We chose the smallest value of $n$ for the Lattice, Matrix, Ngon such
that the ratios increased monotonically as $n$ increased beyond that value.
(The smallest values of $n$ generally had higher ratios.)

What sets the overall computational time of the algorithm? Let’s count the
number of linear algebra operations involving matrices that contribute to the
computational time (those involving only vectors are expected to be
insignificant).

  * •

Proposing a tangent step $v_{x}$: 2 matrix-vector multiplications, one
forward-backward substitution.

  * •

Projecting to obtain $y$: $n_{{\rm iter}}$ times all of the following: one
matrix-vector multiplication, one forward-backward substitution, and for
traditional Newton, one LU decomposition and one Jacobian evaluation.

If $y$ is succcessfully proposed:

  * •

Computing $v_{y}$: one matrix-matrix multiplication, one Cholesky
decomposition, 2 matrix-vector multiplications, one forward-backward
substitution.

  * •

If $y$ is not rejected by the Metropolis step: one reverse projection,
involving the same steps as the forward projection.

We expect the matrix factorizations to be the most expensive step, so we
estimate the amount they contribute to the total computational time. Let
$T_{{\rm new}},T_{{\rm sym}}$ be the average computational time to generate 1
new point for each of traditional Newton and symmetric Newton respectively,
determined empirically in Section 4.2, and let $T_{{\rm new}}^{{\rm
est}},T_{{\rm sym}}^{{\rm est}}$ be estimates obtained by counting the time
associated with matrix factorizations. These estimates are constructed from
the times $T_{{\rm LU}},T_{{\rm chol}}$ required to perform one LU
factorization and one Cholesky factorization respectively, multiplied by the
average number of times each such factorization occurs. In each step of the
Markov chain, symmetric Newton requires an average of approximately $1-r_{{\rm
sym}}(1-a)$ Cholesky factorizations, where $r_{{\rm sym}}$ ($r_{{\rm new}}$
for traditional Newton) is the fraction of rejections in the forward
projection step, so that $r_{{\rm sym}}(1-a)$ is the total fraction of
proposals rejected in the projection step, hence which don’t proceed to a
Cholesky decomposition of a proposal $y$. Traditional Newton requires
approximately $1-r_{{\rm new}}(1-a)$ Cholesky factorizations, and
approximately $(1+a)n_{\rm iter}$ LU factorizations, where $an_{{\rm
iter}}=0.25n_{{\rm iter}}$ is added to account for the reverse projection
step, which occurs approximately once per point accepted. Therefore, we may
estimate the computational times as

$T_{{\rm sym}}^{{\rm est}}=(1-0.75r_{{\rm sym}})T_{{\rm chol}},\qquad T_{{\rm
new}}^{{\rm est}}=1.25n_{{\rm iter}}T_{{\rm LU}}+(1-0.75r_{{\rm new}})T_{{\rm
chol}}.$ (20)

We constructed these estimates by first estimating $T_{{\rm chol}},T_{{\rm
LU}}$, which we did for each example by constructing a matrix $A$ with the
same sparsity pattern as $Q_{x}^{T}Q_{x}$, and repeatedly performing either a
Cholesky decomposition or a LU decomposition of this matrix. We measured the
average time to perform $10^{4}$ factorizations of each (and repeated for
several random matrices to ensure the results were independent of the
particular entries, which they were).

Figure 6 shows the actual and estimated timings, $T_{{\rm sym}},T_{{\rm new}}$
and $T_{{\rm sym}}^{{\rm est}},T_{{\rm new}}^{{\rm est}}$ respectively, and
Table 2 shows the percentages explained by matrix factorizations,
$T_{(\cdot)}^{{\rm est}}/T_{(\cdot)}$, for small and large $n$. For all
examples except the Polymer, the estimates (20) are significantly smaller than
the measured values for small $n$, but their agreement increases with $n$,
until at large $n$ they explain most of the timing – 80-100% for traditional
Newton, and 40-60% for symmetric Newton. For the Polymer, the percentage
explained by the estimates actually decreases with $n$, though slowly, and is
approximately 40-50% for traditional Newton, and 4% for symmetric Newton.

We believe the Polymer behaves differently because it has so few constraints,
with such a simple sparsity structure ($Q_{x}$ has 9 nonzero elements per
column, and in $Q_{x}^{T}Q_{x}$ the nonzero elements are almost all
concentrated near the diagonal), that performing matrix factorizations is
roughly comparable to the other linear algebra operations. There are several
reasons why $T_{{\rm sym}}^{{\rm est}}>T_{{\rm sym}}$ for the Ngon with
traditional Newton at large $n$ – perhaps the 3 realizations behave
differently and our estimate ignores such correlations in statistics; perhaps
the slight difference in $n_{{\rm iter}}$ for successful and failed
projections, which we haven’t accounted for in (20), is playing a role.

For the Lattice, Matrix, and Ngon, our estimates (20) show that at large $n$,
the primary contribution to the total time of the algorithm are the matrix
factorizations. This is why symmetric Newton is so much faster than
traditional Newton – because it significantly reduces the number of matrix
factorizations. This agreement also allows one to generalize our findings to
other implementations of the algorithm, e.g. different machine architectures,
sparsity orderings, linear algebra routines, etc – we expect matrix
factorizations to still dominate given other choices, so one can estimate the
total time of the algorithm using the time it takes to perform a matrix
factorization, as in (20). (This model also uses $r_{{\rm sym}},r_{{\rm
new}}$, which can be dropped without changing the estimates significantly.)
Furthermore, the model shows that improving the factorizations, e.g. by
performing them in parallel, or by avoiding them altogether (an issue we
discuss in Section 5), would have the greatest effect on timing.

It is notable however that for symmetric Newton, matrix factorizations
represent about 50% of the timing at large $n$ (though this percentage would
be higher if we explored even larger $n$), and much less for smaller $n$,
meaning that we have made matrix factorizations equivalent to the other linear
algebra operations for the values of $n$ explored. This suggests that even if
we were able to avoid a matrix factorization, e.g. by using iterative methods
to solve systems of equations instead of the direct methods we proposed in
this paper, the effect on timing would not be significant, except for larger
values of $n$ than we have explored here.

###  4.5 Effective diffusivity as a function of $a$

Figure 7: Estimated diffusivity $a\sigma_{a}^{2}/\tau_{c}(a)$ as a function of
average acceptance probability $a$, for several values of $n$. Traditional
Newton is shown in shades of blue/circle markers; symmetric Newton is shown in
shades of red/diamond markers. One sample was used for the Ngon.

We repeated the experiments above with different values of $a$ to test the
robustness of our observations, and found they remained qualitatively the same
for all $a$. Quantitatively, the computational times changed, generally
increasing weakly with $a$ by a factor of at most 2, except for the Lattice
for which it was nearly constant. However, the ratio of times between
traditional Newton and symmetric Newton remained roughly the same. The
rejection statistics were also roughly similar, though with a smaller
proportion of Projection rejections and a greater proportion of Metropolis
rejections at larger values of $a$. This shift in rejection statistics could
underlie the slight increase in computational time, since a decreasing
proportional of forward projection rejections means an increase in the average
number of matrix factorizations, since there are a greater fraction of
successfully proposed $y$s.

Which value of $a$ makes the method most efficient? We tested this by
computing the effective diffusivity $D_{\rm eff}(a)$, assuming steps are small
enough that diffusivity is estimated as for a random walk in a flat space by
the averaged step size squared, divided by the mean computational time per
step:

$D_{\rm eff}(a)=a\sigma_{a}^{2}/T(a)$ (21)

Figure 7 shows that the value of $a$ that maximizes $D_{\rm eff}(a)$ generally
depends on $n$, as well as on the presence of an ambient energy function. For
the three examples with no ambient energy function – the Polymer, Matrix, and
Ngon – $D_{\rm eff}(a)$ is maximized at the smallest values of $a$ when $n$ is
small, gradually transitioning to a maximum at large values of $a$ for larger
$n$. For the Lattice, which has an ambient energy function, $D_{\rm eff}(a)$
is generally maximized near $a=0.25$; this is consistent with the theoretical
predictions of $a=0.23$ (a value we didn’t test) [14].

We may speculate why large $a$ is more efficient for large $n$, for the
Polymer, Matrix, and Ngon. For these examples, at large $n$, $\sigma_{a}$
decreases only weakly with $a$. A small decrease in $\sigma$ therefore gives a
large increase in $a$; this is mostly because the number of Metropolis
rejections decreases – apparently the difference in $|v_{x}|,|v_{y}|$ is felt
more strongly at large $n$. For these values of $n$ it is advantageous to use
a large acceptance ratio so the effective step size $a\sigma_{a}^{2}$ is
larger; the slight increase in computational time and slight decrease in
$\sigma$ do not appear to outweigh the gain in effective diffusivity (21) due
to increasing $a$. This issue does not apply to the Lattice, where
$\sigma_{a}$ decreases more rapidly with $a$ at all values of $n$. In this
example it is the energy function which determines the step size, and not the
limits set by asking for approximate reversibility in the geometry.

It is also notable that $D_{\rm eff}(a)$ does not vary strongly with $a$ – the
most it varied over the range $a\in[0.15,0.7]$ for the Polymer, Matrix, Ngon
was by a factor of 2 (a factor of 3 for the Ngon with small $n$), and by a
factor of about 4 for the Lattice. Therefore, the choice of $a$ does not seem
particularly important within reasonable limits.

##  5 Discussion

We have introduced a numerical algorithm to implement a Markov Chain Monte
Carlo method to sample probability distributions defined on manifolds, where
the manifolds are defined by the level sets of constraint functions. The
algorithm builds on the method proposed in [48] by proposing a numerically
efficient implementation of this method, which is easy to implement, and which
scales well to problems with thousands of dimensions. The algorithm requires
one matrix factorization, a Cholesky factorization of a large sparse matrix,
per step, and it uses this factorization both to directly solve linear systems
of equations, and to compute the determinant of a Jacobian matrix ($Q_{x}$)
which is required in physical systems where the constraints are approximations
for stiff forces. A key step in making the method maximally efficient is to
use a quasi-Newton method to solve the nonlinear system of equations defining
the manifold, by replacing the exact Jacobian of the equations with a
symmetric approximation that doesn’t change with iteration, an idea that was
first explored in the molecular dynamics community [2].

We applied the algorithm to several examples with varying constraint
structures, most of which were inspired by physics and materials science.
Three out of four examples had no ambient energy function, which allowed us to
focus purely on the geometrical aspects of the sampler. The computational
complexity of the algorithm varies with the example, and is mainly determined
by the complexity of performing a sparse matrix factorization. This in turn
depends on the complexity of the constraints, with simpler constraints (such
as for a polymer with fixed endpoints) having the smallest complexity, and
dense constraints with no structure (such as for a random graph) having the
highest complexity. The algorithm becomes infeasible when performing the
requisite Cholesky factorization becomes infeasible.

We estimated the fraction of time taken by the Cholesky factorization and
found it was negligible for small system sizes, though increasing to about 50%
for the largest system sizes we considered (provided they had sufficiently
complex constraints); presumably this fraction would continue to increase for
even larger systems. Therefore, for large system sizes a parallel
implementation of the requisite linear algebra operations should be
beneficial, though we haven’t explored parallelization here.

We estimated the effective diffusivity of the sampler (measured in units of
space squared divided by computational time) as a function of the average
acceptance probability $a$, and found that for problems with no ambient
energy, the sampler was most efficient at small $a$ in low dimensions and at
large $a$ in high dimensions. However, as the effective diffusivity did not
vary by more than a factor of 2-3 with $a$, we conclude that the particular
choice of $a$ is unimportant.

For very large systems one might wish to use a method that doesn’t require a
matrix factorization. In fact, _the only step in our algorithm for which the
factorization is required_ is computing the ratio of pseudodeterminants
$|Q_{y}|^{-1}/|Q_{x}|^{-1}$ – we know of no iterative way to compute this
ratio. Other steps which used the factorization could in principle be done
using iterative methods, such as conjugate gradient to solve equations
(8),(14),(13), or GM-RES to solve (12). Such iterative methods require only
matrix-vector multiplications, and hence may be applied to extremely high-
dimensional problems, however direct methods are generally preferable to
iterative methods when available.

Therefore, it is ultimately computing the ratio $|Q_{y}|^{-1}/|Q_{x}|^{-1}$
which appears to be the limiting factor in determining the sizes of the
systems one may sample using the method we have introduced, provided one is
interested in systems which require this factor. (In molecular dynamics this
factor appears to be omitted, e.g. [2], however if a system is quite floppy
this should lead to nonnegligible errors in statistics [26].) Interestingly,
one may avoid computing determinants by using a flow-based method to solve
SDEs on manifolds [49], though it seems impossible to metropolize such a
method. This flow-based method operates by taking an unconstrained step in the
ambient space $x\to x+\xi$, and then projecting to the manifold by solving the
flow $\frac{d\phi}{dt}=-\nabla F(\phi)$, $\phi(0)=x+\xi$, where
$F(x)=|q(x)|^{2}$; this projection implicitly accounts for the determinant
factor. While solving the ODEs for the flow is expected to be slow, it would
be interesting to test whether this is still faster than computing a matrix
factorization in high dimensions.

We anticipate the collection of numerical techniques developed here can be
extended to other sampling and simulations algorithms on manifolds. We
anticipate in particular that they will be useful when extended to a Monte
Carlo method that samples a collection of manifolds (a stratification), which
is necessary for example to simulate sticky particles where bonds are not
constant but can break and form [22]. Most ideas used here should apply,
though new ones will be needed to efficiently handle the Jacobian factors
associated with adding and subtracting bonds, which are not directly obtained
from $Q_{x}$, or more generally in jumping to manifolds of lower or higher
dimension.

## Acknowledgements

MHC would like to thank George Stadler and Jonathan Goodman for many helpful
conversations, and Benedict Leimkuhler for introducing us to [2]. KX and MHC
acknowledge support from Research Training Group in Modeling and Simulation
funded by the National Science Foundation via grant RTG/DMS–1646339. MHC
acknowledges support from NSF-DMS-2111163, NSERC-GR026736, and the Alfred P.
Sloan Foundation.

## References

  * Andersen [1983] Hans C Andersen.  Rattle: A “velocity” version of the shake algorithm for molecular dynamics calculations.  _Curr Opin Struct Biol._ , 52:24 – 34, 00 1983.  URL https://www.sciencedirect.com/science/article/pii/0021999183900141.
  * Barth et al. [1995] Eric Barth, Krzysztof Kuczera, Benedict Leimkuhler, and Robert D. Skeel.  Algorithms for constrained molecular dynamics.  _Journal of Computational Chemistry_ , 16(10):1192–1209, 00 1995.  ISSN 1096-987X.  doi: 10.1002/jcc.540161003.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.540161003.
  * Bertoldi et al. [2017] Katia Bertoldi, Vincenzo Vitelli, Johan Christensen, and Martin van Hecke.  Flexible mechanical metamaterials.  _Nature Reviews Materials_ , 2(11):17066, 2017\.  doi: 10.1038/natrevmats.2017.66.
  * Breiding and Marigliano [2020] Paul Breiding and Orlando Marigliano.  Random Points on an Algebraic Manifold.  _SIAM Journal on Mathematics of Data Science_ , 2(3):683–704, 2020.  doi: 10.1137/19m1271178.
  * Brubaker et al. [2012] Marcus A Brubaker, Mathieu Salzmann, and Raquel Urtasun.  A Family of MCMC Methods on Implicitly Defined Manifolds.  _Proceedings of the 15th International Conference on Artificial Intelligence and Statistics_ , pages 1 – 12, 04 2012.
  * Byrne and Girolami [2013] S Byrne and M Girolami.  Geodesic Monte Carlo on embedded manifolds.  _Scandinavian Journal of Statistics_ , 00 2013.  doi: 10.1111/sjos.12036.  URL http://onlinelibrary.wiley.com/doi/10.1111/sjos.12036/full.
  * Ciccotti et al. [2007] Giovanni Ciccotti, Tony Lelièvre, and Eric Vanden-Eijnden.  Projection of diffusions on submanifolds: Application to mean force computation.  _Communications on Pure and Applied Mathematics_ , 61(3):371 – 408, 00 2007.  doi: 10.1002/cpa.20210.  URL http://doi.wiley.com/10.1002/cpa.20210.
  * Cowles and Carlin [1996] Mary Kathryn Cowles and Bradley P. Carlin.  Markov Chain Monte Carlo Convergence Diagnostics: A Comparative Review.  _Journal of the American Statistical Association_ , 91(434):883–904, 1996.  ISSN 0162-1459.  doi: 10.1080/01621459.1996.10476956.
  * Davis [2006] Timothy Davis.  _Direct methods for sparse linear systems_.  2006\.
  * Diaconis et al. [2013] P Diaconis, S Holmes, and Mehrdad Shahshahani.  Sampling from a manifold.  _Advances in Modern Statistical Theory and Applications_ , 10:102 – 125, 00 2013.  doi: 10.1214/12-imscoll1006.  URL http://projecteuclid.org/euclid.imsc/1379942050.
  * Eastman and Pande [2010] Peter Eastman and Vijay S. Pande.  Constant Constraint Matrix Approximation: A Robust, Parallelizable Constraint Method for Molecular Simulations.  _Journal of Chemical Theory and Computation_ , 6(2):434–437, 00 2010.  ISSN 1549-9618.  doi: 10.1021/ct900463w.  URL https://pubs.acs.org/doi/abs/10.1021/ct900463w.
  * Fatkullin et al. [2010] Ibrahim Fatkullin, Gregor Kovacic, and Eric Vanden-Eijnden.  Reduced Dynamics of Stochastically Perturbed Gradient Flows.  _Communications in Mathematical Sciences_ , 8(2):439 – 461, 06 2010.  URL http://gateway.webofknowledge.com/gateway/Gateway.cgi?GWVersion=2&SrcAuth=mekentosj&SrcApp=Papers&DestLinkType=FullRecord&DestApp=WOS&KeyUT=000272842200008.
  * Fixman [1978] Marshall Fixman.  Simulation of polymer dynamics. I. General theory.  _The Journal of Chemical Physics_ , 69(4):1527–1537, 1978.  ISSN 0021-9606.  doi: 10.1063/1.436725.
  * Gelman et al. [1996] A. Gelman, G. O. Roberts, and W. R. Gilks.  Efficient Metropolis Jumping Rules.  _Bayesian Statistics_ , 5:599–607, 1996.
  * Gonnet [2007] Pedro Gonnet.  P-SHAKE: A quadratically convergent SHAKE in O(n2).  _Journal of Computational Physics_ , 220(2):740–750, 2007.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2006.05.032.
  * Graham et al. [2022] Matthew M. Graham, Alexandre H. Thiery, and Alexandros Beskos.  Manifold Markov chain Monte Carlo methods for Bayesian inference in diffusion models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 84(4):1229–1256, 2022.  ISSN 1369-7412.  doi: 10.1111/rssb.12497.
  * Griffiths and Hernández-Lobato [2019] Ryan-Rhys Griffiths and José Miguel Hernández-Lobato.  Constrained Bayesian optimization for automatic chemical design using variational autoencoders.  _Chemical Science_ , 11(2):577–586, 2019.  ISSN 2041-6520.  doi: 10.1039/c9sc04026a.
  * Gunsteren [1980] W F van Gunsteren.  Constrained dynamics of flexible molecules.  _Molecular Physics_ , 40(4):1015 – 1019, 07 1980\.  ISSN 0026-8976.  doi: 10.1080/00268978000102101.  URL http://www.tandfonline.com/doi/abs/10.1080/00268978000102101.
  * Hartmann [2007] Carsten Hartmann.  An Ergodic Sampling Scheme for Constrained Hamiltonian Systems with Applications to Molecular Dynamics.  _Journal of Statistical Physics_ , 130(4):687–711, 2007.  ISSN 0022-4715.  doi: 10.1007/s10955-007-9470-2.
  * Headrick and Nassar [2013] M Headrick and A Nassar.  Energy functionals for Calabi-Yau metrics.  _Journal of Physics: Conference Series_ , 462(1):012019, 2013.  ISSN 1742-6596.  doi: 10.1088/1742-6596/462/1/012019.
  * Hinch [1994] E J Hinch.  Brownian motion with stiff bonds and rigid constraints.  _Journal of Fluid Mechanics_ , 271:219 – 234, 00 1994.  ISSN 0022-1120.  doi: 10.1017/s0022112094001746.  URL http://dx.doi.org/10.1017/S0022112094001746.
  * Holmes-Cerfon [2020] Miranda Holmes-Cerfon.  Simulating sticky particles: A Monte Carlo method to sample a stratification.  _Journal of Chemical Physics_ , 153(16):164112, 10 2020.  doi: 10.1063/5.0019550.  URL http://aip.scitation.org/doi/10.1063/5.0019550.
  * Holmes-Cerfon [2023] Miranda Holmes-Cerfon.  ManifoldSampler, 2023.  URL https://github.com/mirandaholmes/ManifoldSampler/tree/main.
  * Holmes-Cerfon et al. [2013] Miranda Holmes-Cerfon, Steven J Gortler, and Michael P Brenner.  A geometrical approach to computing free-energy landscapes from short-ranged potentials.  _Proceedings of the National Academy of Sciences_ , 110(1):E5 – E14, 00 2013.  doi: 10.1073/pnas.1211720110.  URL http://www.pnas.org/cgi/doi/10.1073/pnas.1211720110.
  * Holmes-Cerfon [2016] Miranda C Holmes-Cerfon.  Enumerating Rigid Sphere Packings.  _SIAM Review_ , 58(2):229 – 244, 05 2016.  ISSN 0036-1445.  doi: 10.1137/140982337.  URL http://epubs.siam.org/doi/10.1137/140982337.
  * Kampen [1981] N G van Kampen.  Statistical mechanics of trimers.  _Applied Scientific Research. An International Journal on Thermal, Mechanical, and Electromagnetic Phenomena in Continua_ , 37(1-2):67 – 75, 00 1981.  doi: 10.1007/bf00382618.  URL http://link.springer.com/10.1007/BF00382618.
  * Klugkist and Hoijtink [2007] Irene Klugkist and Herbert Hoijtink.  The Bayes factor for inequality and about equality constrained models.  _Computational Statistics & Data Analysis_, 51(12):6367–6379, 2007.  ISSN 0167-9473.  doi: 10.1016/j.csda.2007.01.024.
  * Kräutler et al. [2001] Vincent Kräutler, Wilfred F. van Gunsteren, and Philippe H. Hünenberger.  A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations.  _Journal of Computational Chemistry_ , 22(5):501–508, 2001.  ISSN 0192-8651.  doi: 10.1002/1096-987x(20010415)22:5<501::aid-jcc1021>3.0.co;2-v.
  * Lasserre and Magron [2020] Jean B. Lasserre and Victor Magron.  Computing the Hausdorff Boundary Measure of Semialgebraic Sets.  _SIAM Journal on Applied Algebra and Geometry_ , 4(3):441–469, 2020.  doi: 10.1137/20m1314392.
  * Leimkuhler and Matthews [2016] Benedict Leimkuhler and Charles Matthews.  Efficient molecular dynamics using geodesic integration and solvent–solute splitting.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 472(2189):20160138 – 22, 05 2016.  ISSN 1364-5021.  doi: 10.1098/rspa.2016.0138.  URL https://royalsocietypublishing.org/doi/10.1098/rspa.2016.0138.
  * Leimkuhler et al. [2020] Benedict Leimkuhler, Timothée Pouchon, Tiffany Vlaar, and Amos Storkey.  Constraint-Based Regularization of Neural Networks.  _arXiv_ , 2020.  doi: 10.48550/arxiv.2006.10114.
  * Lelièvre et al. [2012] Tony Lelièvre, Mathias Rousset, and Gabriel Stoltz.  Langevin dynamics with constraints and computation of free energy differences.  _Mathematics of Computation_ , 81(280):2071 – 2125, 00 2012.  doi: 10.1090/s0025-5718-2012-02594-4.  URL http://www.ams.org/jourcgi/jour-getitem?pii=S0025-5718-2012-02594-4.
  * Lelièvre et al. [2019] Tony Lelièvre, Mathias Rousset, and Gabriel Stoltz.  Hybrid Monte Carlo methods for sampling probability measures on submanifolds.  _Numerische Mathematik_ , 143(2):379 – 421, 2019\.  doi: 10.1007/s00211-019-01056-4.  URL http://link.springer.com/10.1007/s00211-019-01056-4.
  * Manoharan [2015] V N Manoharan.  Colloidal matter: Packing, geometry, and entropy.  _Science_ , 349(6251):1253751 – 1253751, 08 2015\.  ISSN 0036-8075.  doi: 10.1126/science.1253751.  URL http://www.sciencemag.org/cgi/doi/10.1126/science.1253751.
  * Mezzadri [2007] Francesco Mezzadri.  How to generate random matrices from the classical compact groups.  _Notices of the American Mathematical Society_ , 54(5):592–604, 2007.
  * Morgan [2016] Morgan.  _Geometric measure theory: a beginner’s guide_.  5th edition, 2016.
  * Mulder [2014] Joris Mulder.  Prior adjusted default Bayes factors for testing (in)equality constrained hypotheses.  _Computational Statistics & Data Analysis_, 71:448–463, 2014.  ISSN 0167-9473.  doi: 10.1016/j.csda.2013.07.017.
  * Nocedal and Wright [2006] J Nocedal and S Wright.  _Numerical optimization_.  Springer. Springer, 2nd edition, 00 2006.
  * Oyarzún and Mognetti [2018] Bernardo Oyarzún and Bortolo Matteo Mognetti.  Efficient sampling of reversible cross-linking polymers: Self-assembly of single-chain polymeric nanoparticles.  _Journal of Chemical Physics_ , 148(11):114110 – 15, 03 2018.  ISSN 0021-9606.  doi: 10.1063/1.5020158.  URL http://aip.scitation.org/doi/10.1063/1.5020158.
  * Qin, Galin Jones and Qian [2022] Qin, Galin Jones and Qian.  Markov Chain Monte Carlo in Practice.  _Annual Review of Statistics and Its Application_ , 9:557–578, 2022.
  * Roberts et al. [1997] G. O. Roberts, A. Gelman, and W. R. Gilks.  Weak convergence and optimal scaling of random walk Metropolis algorithms.  _The Annals of Applied Probability_ , 7(1), 1997.  ISSN 1050-5164.  doi: 10.1214/aoap/1034625254.
  * Ryckaert et al. [1977] J P Ryckaert, G Ciccotti, and HJC Berendsen.  Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes.  _Journal of Computational Physics_ , 00 1977.  URL http://www.sciencedirect.com/science/article/pii/0021999177900985.
  * Sharma and Zhang [2021] Upanshu Sharma and Wei Zhang.  NonReversible Sampling Schemes on Submanifolds.  _SIAM Journal on Numerical Analysis_ , 59(6):2989–3031, 2021.  ISSN 0036-1429.  doi: 10.1137/20m1378752.
  * Vanden-Eijnden and Ciccotti [2006] Eric Vanden-Eijnden and Giovanni Ciccotti.  Second-order integrators for Langevin equations with holonomic constraints.  _Chemical Physics Letters_ , 429(1-3):310 – 316, 09 2006.  doi: 10.1016/j.cplett.2006.07.086.  URL http://linkinghub.elsevier.com/retrieve/pii/S0009261406011092.
  * Weinbach and Elber [2005] Yael Weinbach and Ron Elber.  Revisiting and parallelizing SHAKE.  _Journal of Computational Physics_ , 209(1):193 – 206, 10 2005.  doi: 10.1016/j.jcp.2005.03.015.  URL https://linkinghub.elsevier.com/retrieve/pii/S0021999105001828.
  * Yang et al. [2022] Jun Yang, Krzysztof Łatuszyński, and Gareth O Roberts.  Stereographic Markov Chain Monte Carlo.  _arXiv_ , 2022.  doi: 10.48550/arxiv.2205.12112.
  * Yoneya [2001] Makoto Yoneya.  A Generalized Non-iterative Matrix Method for Constraint Molecular Dynamics Simulations.  _Journal of Computational Physics_ , 172(1):188 – 197, 09 2001.  doi: 10.1006/jcph.2001.6819.  URL https://linkinghub.elsevier.com/retrieve/pii/S0021999101968196.
  * Zappa et al. [2018] Emilio Zappa, Miranda Holmes-Cerfon, and Jonathan Goodman.  Monte Carlo on Manifolds: Sampling Densities and Integrating Functions.  _Communications on Pure and Applied Mathematics_ , 71(12):2609 – 2647, 00 2018.  doi: 10.1002/cpa.21783.  URL http://doi.wiley.com/10.1002/cpa.21783.
  * Zhang [2020] Wei Zhang.  Ergodic SDEs on submanifolds and related numerical sampling schemes.  _ESAIM: Mathematical Modelling and Numerical Analysis_ , 54(2):391–430, 2020.  ISSN 0764-583X.  doi: 10.1051/m2an/2019071.