# Monte Carlo with Soft Constraints:
the Surface Augmented Sampler

BRANDO MAGNANI  Courant Institute of Mathematical Sciences, New York
University, 251 Mercer Street, New York, NY 10012, USA<EMAIL_ADDRESS>
Abstract

> We describe an MCMC method for sampling distributions with soft constraints,
> which are constraints that are almost but not exactly satisfied. We sample a
> total distribution that is a convex combination of the target ”soft”
> distribution with the nearby ”hard” distribution supported on the constraint
> surface. Hard distribution moves lead to performance that is uniform in the
> softness parameter. On/Off moves related to the Holmes-Cerfon Stratification
> Sampler enable sampling the target soft distribution. Computational
> experiments verify that performance is uniform in the soft constraints
> limit.

## 1 Introduction

Many physical and engineering systems involve motion with soft constraints,
that is, with constraints that are only approximately satisfied. This type of
constraints arise in the systems’ low temperature limit, and correspond to a
situation in which motion occurs (with high probability) within a small
neighborhood of a surface that is given as a set of equality constraints. On
the other hand, hard constraints are those which are imposed directly on the
equations of motion of the system.

When simulating particles interacting with bonds that can break and form, one
often treats the bond lengths as fixed (as hard constraints) to speed up
sampling from the equilibrium distribution. This ad hoc assumption is
oftentimes not physically realistic. For example, in molecular dynamics bonds
between particles are never fixed at one length, but rather vibrate around
that length (an example of soft constraints). The effects of this vibration
are not well studied because the available sampling methods either can’t
handle soft constraints or, if they can, they require a step size that is much
smaller than the scale at which it is possible to observe interesting
rearrangement of the system.

This paper will introduce the lgorithm, a Markov chain Monte Carlo method that
efficiently samples distributions which arise from the imposition of soft
constraints to a given physical system. We are dealing with the the problem of
sampling from distributions that are concentrated around a smooth
$d$-dimensional surface embedded in a $d_{a}$-dimensional ambient space.
Suppose we have $m$ smooth functions
$q_{i}:\mathbb{R}^{d_{a}}\rightarrow{\mathbb{R}}$, $i=1,...,m$. Together,
these define a constraint function
$q:\mathbb{R}^{d_{a}}\rightarrow\mathbb{R}^{m}$ and a hard constraint surface

$\displaystyle S=\big{\\{}x\in\mathbb{R}^{d_{a}}\ :\ q(x)=0\big{\\}}.$ (1.1)

Throughout the rest of the paper, we will use $\nabla q(x)$ to denote the
matrix whose colums are the gradients $\\{\nabla q_{i}(x)\\}_{i=1}^{m}$. This
matrix corresponds to the transpose of the Jacobian of the full constraint
function $q(x)$. We will assume that $\nabla q(x)$ has full column rank $m$
everywhere in $\mathbb{R}^{d_{a}}$. The implicit function theorem then implies
that the dimension of $S$ is $d=d_{a}-m$ and that the tangent space
$\mathcal{T}_{x}\equiv\mathcal{T}_{x}S$ is well-defined for every $x\in S$. In
this case, the vectors $\\{\nabla q_{i}(x)\\}_{i=1}^{m}$ form a basis for the
normal space $\mathcal{N}_{x}\equiv\mathcal{T}_{x}S^{\perp}$. We will also use
$\sigma(dx)$ to denote the $d$-dimensional Hausdorff measure on $S$.

We are interested in sampling from the distribution

$\displaystyle\pi_{\varepsilon}(dx)=\frac{1}{Z_{\varepsilon}}\exp\Big{(}-\frac{U(x)}{2\varepsilon^{2}}\Big{)}dx,\quad
U(x)=|q(x)|^{2}=\sum_{i=1}^{m}q_{i}(x)^{2},$ (1.2)

where $\varepsilon$ is a fixed parameter and
$Z_{\varepsilon}=\int_{\mathbb{R}^{d_{a}}}\exp(-\frac{U(x)}{2\varepsilon^{2}})dx$.
Note that $\varepsilon$ determines the degree to which the distribution is
concentrated near $S$. We want to be able to efficiently sample (1.2) in the
soft constraints limit, that is, as the temperature (or equivalently,
$\varepsilon$) is taken to be arbitrarily small, in which case the
distribution is highly concentrated near $S$. The problem of sampling in the
soft constraints limit is generally difficult. A typical sample $X$ from (1.2)
will have $\text{dist}(X,S)=O(\varepsilon)$. Therefore, any standard method
would require an $O(\varepsilon)$ step size which, in turn, would result in
highly correlated samples. Physically, this means that the step size is much
smaller than the scale at which interesting rearrangements of a given system
can be observed. Our final goal is to construct a sampler that is effective in
the sense that it produces chains whose auto-correlation time stays bounded as
$\varepsilon\rightarrow 0$.

We note here that, at low temperatures, $\pi_{\varepsilon}(dx)$ behaves
approximately as a degenerate Gaussian distribution. If we Taylor expand
$q(x)$ around a point $x_{s}\in S$, then for (say) $|x-x_{s}|\leq\varepsilon$
we have

$\displaystyle U(x)\approx(x-x_{s})^{T}\nabla q(x)\nabla q(x)^{T}(x-x_{s}).$
(1.3)

The degeneracy comes from the fact the rank of $\nabla q(x)\nabla q(x)^{T}$ is
$m<d_{a}$. A related result is that in the limit as $\varepsilon\rightarrow 0$
we have

$\displaystyle\int_{\mathbb{R}^{d_{a}}}\varphi(x)\pi_{\varepsilon}(dx)\longrightarrow\int_{S}\varphi(x)\pi_{s}(dx),$
(1.4)

for any smooth test function $\varphi$ (see appendix A), where

$\displaystyle\pi_{s}(dx)\propto\det(\nabla q(x)^{T}\nabla
q(x))^{-1/2}\sigma(dx).$ (1.5)

The presence of the gradient factor in (1.5) is what distinguishes hard
constraints from soft constraints. While hard constraints are imposed directly
on the equations of motion of a given system, soft constraints arise as we
take the temperature to be arbitrarily small. The fact that
$\pi_{s}\neq\sigma$ shows that the statistical properties of soft constraints
are fundamentally different from those of hard constraints.

The an be seen as an extension to either the ZHG Surface Sampler [4] or the
Stratification Sampler [1]. In addition to surface moves, our method sometimes
proposes jumps to points that are off the constraint surface $S$. This move is
similar to the ”gain” move in [1], with the difference that our method
attempts to drop all the constraints at once, whereas a ”gain” move drops only
one of the constraints at a time. If the current state of the chain is a point
that does not lie on the constraint surface $S$, our method proposes to either
take an isotropic Gaussian Metropolis step in $\mathbb{R}^{d_{a}}$ or to move
back to $S$. A move back to $S$ is similar to the ”lose” move in [1], the
difference being that our method adds all the constraints at once, while a
”lose” move only adds one constraint at a time. Overall, our sampler
randomizes between four different move types. The idea is that one of them,
the Surface Sampler move, will have an $O(1)$ step size as
$\varepsilon\rightarrow 0$. This should substantially reduce the auto-
correlation within the sample.

In order to implement these ideas, we need to work with an augmented
distribution that allows us to combine surface moves along $S$ with moves that
sample points in $S^{c}$. We denote the augmented distribution by $\rho$ and
define it as

$\displaystyle\rho(dx)\propto\pi_{\varepsilon}(dx)+\pi_{s}(dx).$ (1.6)

The ill preserve distribution (1.6). This means that, at equilibrium, samples
in $S^{c}$ will be distributed according to $\pi_{\varepsilon}$, while samples
in $S$ will be distributed according to $\pi_{s}$. We decided to ”weight” the
hard constraint part of $\rho$ by $\pi_{s}$ not only because it’s the
physically correct measure in the soft constraints limit, but also because
this will help us achieve $100\%$ acceptance probabilities in case $S$ is a
flat manifold (more on this in Section 3).

Acknowledgments. Thank you to my advisor Jonathan Goodman and to Miranda
Holmes-Cerfon for the helpful discussions. Also thank you to Georg Stadler for
reading earlier drafts of this work.

## 2 Algorithm

In this section we give an overview of our algorithm, the or generating a
Markov Chain $X_{1},X_{2},...\in\mathbb{R}^{d_{a}}$ with stationary
distribution $\rho$ as in (1.6). We will rewrite distribution (1.6) in a
slightly more abstract form that will (hopefully) simplify the forthcoming
description. For this, we need to introduce the following notation.

Notation. Let $\Omega_{1}=S^{c}=\mathbb{R}^{d_{a}}\setminus S$,
$\Omega_{2}=S$, and $\Omega=\Omega_{1}\cup\Omega_{2}=\mathbb{R}^{d_{a}}$.
Additionally, let $\mu_{1}(dx)=dx$, the Lebesgue measure on
$\mathbb{R}^{d_{a}}$, and $\mu_{2}(dx)=\sigma(dx)$, the $d$-dimensional
Hausdorff measure on $S$. Lastly, let
$f_{1}(x)=k_{1}\exp(-\frac{U(x)}{2\varepsilon^{2}})$ and
$f_{2}(x)=k_{2}\det(\nabla q(x)^{T}\nabla q(x))^{-1/2}$, where $k_{1},k_{2}$
are constants that will be determined in Section 3.

One observation is that the surface set $S$ has Lebesgue measure zero, so we
may think of the two measures $\mu_{1}$ and $\mu_{2}$ as being mutually
singular on the whole space $\Omega$. We can now rewrite (1.6) as

$\displaystyle\rho(dx)$
$\displaystyle=\frac{1}{Z}\Big{[}f_{1}(x)\mu_{1}(dx)+f_{2}(x)\mu_{2}(dx)\Big{]},$
(2.1)

where
$Z=\int_{\Omega_{1}}f_{1}(x)\mu_{1}(dx)+\int_{\Omega_{2}}f_{2}(x)\mu_{2}(dx)$.

The $k$-th state of the Markov chain generated by the s represented as a pair
$(X_{k},i_{k})$, where $X_{k}\in\Omega_{i_{k}}$ and $i_{k}\in\\{1,2\\}$.
Suppose that the $k$-th state is $(X_{k},i_{k})=(x,i)$. The next point
$(X_{k+1},i_{k+1})$ is generated by first constructing a proposal $(y,j)$ and
then accepting or rejecting it by Metropolis-Hastings.

The proposal is constructed via the following two steps:

  1. (1)

Choose a new label $j$ with probability $\lambda_{ij}(x)$;

  2. (2)

Chooose a new point $y\in\Omega_{j}$ with probability density $h_{ij}(x,y)$
with respect to $\mu_{j}(dy)$.

In step (1), $\lambda_{ij}(x)$ represents the probability of choosing label
$j$ given the current state $(x,i)$. Note that we must have
$\lambda_{i1}(x)+\lambda_{i2}(x)=1$ for any point $x\in\Omega_{i}$ and label
$i\in\\{1,2\\}$. The function $h_{ij}(x,y)$ in step (2) is the probability
density of proposing point $y\in\Omega_{j}$ given the that the current state
is $(x,i)$ and that we chose a new label $j$.

The proposal $(y,j)$ is accepted with probability $a_{ij}(x,y)$ defined by

$\displaystyle a_{ij}(x,y)=\min\Bigg{\\{}1,\
\frac{f_{j}(y)\lambda_{ji}(y)h_{ji}(y,x)}{f_{i}(x)\lambda_{ij}(x)h_{ij}(x,y)}\Bigg{\\}}.$
(2.2)

The algorithm generates $U\sim\text{Unif}([0,1])$ and if $U<a_{ij}(x,y)$ then
we accept the proposal and set $(X_{k+1},i_{k+1})=(y,j)$, otherwise we reject
it and set $(X_{k+1},i_{k+1})=(x,i)$.

We will show in Section 5 that the Metropolis-Hastings acceptance probability
leads to a Markov transition kernel that is self-adjoint in
$L^{2}_{\rho}(\Omega)$, which implies that the chain preserves the target
distribution $\rho$. Hence, we expect that for reasonable choices of proposal
densities the Markov chain is ergodic. If $X_{k}\sim\rho_{k}$, this means that
for any initial condition we have $\rho_{k}\rightarrow\rho$ as
$k\rightarrow\infty$.

Important Remarks: Some additional considerations are required to implement
the he first is that when $i\neq j$ the proposal densities in (2.2) are taken
with respect to different reference measures, therefore one has to account for
the Jacobian of the maps between $\Omega_{i}$ and $\Omega_{j}$. Additionally,
the proposal has to be constructed so that there is a unique map between the
current state $(x,i)$ and the proposal $(y,j)$. Last but not least, we want to
find proposal densities that produce an average acceptance probability that is
neither too high nor too low. This is non-trivial when $i\neq j$, that is,
when the sampler proposes to add or drop the constraints. The corresponding
moves need to be constructed so that the numerator and the denominator in the
Metropolis-Hastings ratio in (2.2) are asymptotically equivalent in the soft
constraints limit ($\varepsilon\downarrow 0$). If not, we would have either
$a_{12}$ or $a_{21}$ vanishing uniformly as $\varepsilon\downarrow 0$, which
in turn would make the ery inefficient111Very inefficient because in this case
transitioning from the $\pi_{s}$-regime to the $\pi_{\varepsilon}$-regime (or
viceversa) becomes a rare event. If the sampler gets stuck in the
$\pi_{s}$-regime for a long time, then many of the generated samples will be
discarded since we are only interested in draws from $\pi_{\varepsilon}$. On
the other hand, if the sampler gets stuck in the $\pi_{\varepsilon}$-regime,
we will retain the samples but these will be highly correlated since most of
them will be generated via isotropic Gaussian Metropolis with an
$O(\varepsilon)$ step size (see the Soft move in Section 2.2). As we will see
in Section 3, one way to achieve such asymptotic equivalence is to construct
moves that produce a $100\%$ acceptance probability when $S$ is a flat
manifold defined by linear constraints (such as a line or a plane). When the
surface is not flat, the smoothness of the constraint function $q$ then
guarantees that $S$ is locally well-approximated by a plane (its tangent
space), so that we can expect $a_{12}(x,y),a_{21}(x,y)\nearrow 1$ as
$\varepsilon\downarrow 0$. This will allow us to control the average
acceptance probability independently of the size of $\varepsilon$.

### 2.1 Label proposal and move type

We will distinguish between four different move types. The new label proposal
$j$ and the current label $i$ will determine the move type as follows (we
recall here that $\Omega_{1}=S^{c}$ and $\Omega_{2}=S$):

  * •

Hard move: We choose $j=i=2$. That is, $x,y\in\Omega_{2}$.
This means that we take a step along the hard constraint $S$. Clearly, this
move keeps all the constraints active.

  * •

Off move: $i=2$ and choose $j=1$. That is, $x\in\Omega_{2}$ and
$y\in\Omega_{1}$.
This means that we have current state $x\in S$ and propose a point $y\in
S^{c}$. This corresponds to a jump ”off” the hard constraint, which is
equivalent to ”dropping” all the constraints at once.

  * •

On move: $i=1$ and choose $j=2$. That is, $x\in\Omega_{1}$ and
$y\in\Omega_{2}$.
This means that we have current state $x\in S$ and propose a point $y\in
S^{c}$. This corresponds to a jump back ”on” the hard constraint, which is
equivalent to ”adding” all the constraints at once.

  * •

Soft move: We choose $j=i=1$. That is, $x,y\in\Omega_{1}$.
This amounts to a step in $S^{c}$.

A practical way to propose a new label is the following. We first decide the
overall probabilities $\lambda_{11},\lambda_{12},\lambda_{21},\lambda_{22}\geq
0$ such that $\lambda_{11}+\lambda_{12}=1$ and $\lambda_{21}+\lambda_{22}=1$.
Suppose that the current state is $(x,i)$. At this point, the algorithm draws
$U\sim\text{Unif}([0,1])$ and

  * •

if $i=1$ and $U<\lambda_{11}$, then we choose a Soft move. Else, if
$U>\lambda_{11}$ we choose an On move.

  * •

if $i=2$ and $U<\lambda_{22}$, then we choose a Hard move. Else, if
$U>\lambda_{22}$ we choose an Off move.

Note that this is equivalent to setting
$\lambda_{ij}(x)=\lambda_{ij}\mathds{1}_{\Omega_{i}}(x)$ for all pairs
$i,j\in\\{1,2\\}$. This implies that $\sum_{i,j}\lambda_{ij}(x)=1$ for any
$x\in\Omega$.

### 2.2 Point proposal

Once we have chosen a label $j$, we must propose a point $y\in\Omega_{j}$. We
now present a proposal mechanism for each of the four move types.

As we will see, three of the four moves depend on taking steps in the tangent
space of the hard constraint surface $S$. We use $\mathcal{T}_{x}$ to denote
the tangent space of $S$ at a point $x\in S$. Recall that $\mathcal{T}_{x}$ is
the $d$-dimensional space of vectors $v\in\mathbb{R}^{d_{a}}$ such that
$\nabla q(x)^{T}v=0$. This represents the linear approximation to $S$ near
$x$. The normal space of $S$ at point $x$ is denoted by
$\mathcal{N}_{x}=\mathcal{T}_{x}^{\perp}$ and it is defined as the
$m$-dimensional space of vectors $v\in\mathbb{R}^{d_{a}}$ which are orthogonal
to $\mathcal{T}_{x}$.

We will use $T_{x}$ to denote an orthonormal matrix such that
$\text{Col}(T_{x})=\mathcal{T}_{x}$. The matrix $T_{x}$ can be computed by
taking the last $d$ columns in the QR decomposition of $\nabla q(x)$.

#### Hard move: $x,\ y\in\Omega_{2}.$

To generate a point $y\in S$ we use the ZHG Surface Sampler [4]. We start by
taking a tangential step $x\mapsto x+v$ with $v\in\mathcal{T}_{x}$. We
generate $v$ by sampling a degenerate isotropic Gaussian defined on the linear
subspace $\mathcal{T}_{x}$ with variance $\sigma_{\text{hrd}}^{2}$. One way to
sample $v$ is to take $v=T_{x}R$, where $R\in\mathbb{R}^{d}$ is a random
vector with i.i.d. entries $R_{i}\sim N(0,\sigma_{\text{hrd}}^{2})$, for
$i=1,...,d$. The density for $v$ is given by

$\displaystyle
p_{\text{hrd}}(v)=\frac{1}{(2\pi)^{d/2}\sigma_{\text{hrd}}^{d}}e^{-\frac{|v|^{2}}{2\sigma_{\text{hrd}}^{2}}}.$
(2.3)

Since the tangent step usually brings us off the surface, we need to project
back onto $S$. Given $x$ and $v$, the projection step finds a
$w\in\mathcal{N}_{x}$ such that $x+v+w\in S$. This is done by first finding a
vector $a\in\mathbb{R}^{m}$ such that

$\displaystyle q\big{(}x+v+\nabla q(x)a\big{)}=0,$ (2.4)

and then setting $w=\nabla q(x)a$. The point proposal will be $y=x+v+w$.

In our code, we numerically solve equation (2.3) by standard Newton’s method.
There are two cases when we must reject the proposal $(y,2)$ so that detailed
balance is preserved. First, if Newton’s method for (2.3) does not converge,
we reject $(y,2)$ immediately and set $(X_{k+1},i_{k+1})=(x,2)$. Second, we
must check whether the reverse move is possible, that is, whether we can
propose point $x$ starting from $y$. The reverse move is constructed as
$x=y+v^{\prime}+w^{\prime}$ where $v^{\prime}\in\mathcal{T}_{y}$ and
$w^{\prime}\in\mathcal{N}_{y}$. Since we already know both $x$ and $y$, we can
always find $v^{\prime}$ and $w^{\prime}$ by projecting $x-y$ onto
$\mathcal{T}_{y}$ and $\mathcal{N}_{y}=\mathcal{T}_{y}^{\perp}$. Then, we must
check whether the Newton solver would find $x$ starting from $y+v^{\prime}$.
If not, we reject $(y,2)$ and set $(X_{k+1},i_{k+1})=(x,2)$.

The density for the point proposal $y$ is given by
$h_{22}(x,y)=p_{\text{hrd}}(x,v)J(x,y)$, where $J(x,y)$ is the inverse of the
determinant of the Jacobian of the projection map $v\rightarrow y$. The
Jacobian factor accounts for how the area of a small patch near $x$ is
distorted upon mapping it to $y$. One can linearize the projection map to find
that $J(x,y)=\det(T_{x}^{T}T_{y})$ (see Section 2.5 for details on how to
compute the Jacobian). The point proposal density for the Hard move is thus
given by

$\displaystyle\begin{split}h_{22}(x,y)&=p_{\text{hrd}}(v)\det(T_{x}^{T}T_{y})\\\
\\\
&=\frac{1}{(2\pi)^{d/2}\sigma_{\text{hrd}}^{d}}e^{-\frac{|v|^{2}}{2\sigma_{\text{hrd}}^{2}}}\det(T_{x}^{T}T_{y}).\end{split}$
(2.5)

Note that in practice we do not need to compute $\det(T_{x}^{T}T_{y})$. This
is because the Jacobian for the reverse move is
$J(y,x)=\det(T_{y}^{T}T_{x})=\det(T_{x}^{T}T_{y})=J(x,y)$, so that the
Jacobian factors simplify in the Metropolis-Hastings ratio.

#### Off move: $x\in\Omega_{2},\ y\in\Omega_{1}.$

Given $x\in S$, an Off move generates a point proposal $y\in S^{c}$. One
possibility would be to take an isotropic Gaussian step in
$\mathbb{R}^{d_{a}}$. As we will see below, this becomes unfeasible in
practice because it is impossible to make such proposal symmetric with the On
move. By contrast, we take a step $v$ that has different sizes in the
directions normal to and tangential to the boundary of $S$. The point proposal
will then be $y=x+v$. For the step in the normal direction, we first need the
transpose of the pseudo-inverse of $\nabla q(x)$, which we denote by $N_{x}$.
If $\nabla q(x)=U\Sigma V^{T}$ is the singular value decomposition for $\nabla
q(x)$, then $N_{x}=U\Sigma^{-1}V^{T}$, where $\Sigma^{-1}$ is obtained by
taking the reciprocals of the non-zero diagonal elements in $\Sigma$. We
construct $v$ by combining a random step in the directions spanned by $N_{x}$
and a random step in the directions spanned by $T_{x}$, so that

$\displaystyle\begin{split}v=N_{x}R_{n}+T_{x}R_{t},\end{split}$ (2.6)

where $R_{n}\in\mathbb{R}^{m}$ is a random vector with i.i.d. entries
$R_{n,i}\sim N(0,\sigma_{\text{prp}}^{2})$ and $R_{t}\in\mathbb{R}^{d}$ is a
random vector with i.i.d. entries $R_{t,i}\sim N(0,\sigma_{\text{tan}}^{2})$.

Note that $\text{col}(N_{x})=\mathcal{N}_{x}$, however $N_{x}$ does not have
orthonormal columns. As we will see in Section 3, this particular choice of
basis matrix for the step in the normal direction plays a key role in
determining a $100\%$ acceptance probability for the Off move in case $S$ is a
flat manifold. The point proposal density for the Off move is denoted by

$\displaystyle
h_{21}(x,y)=p_{\text{off}}^{\perp}(v_{n})p_{\text{off}}^{\|}(v_{t}),$

where $p_{\text{off}}^{\perp}(v_{n})$ is the density for the normal step
$V_{n}=N_{x}R_{n}$, and $p_{\text{off}}^{\|}(v_{t})$ is the density for the
tangent step $V_{t}=T_{x}R_{t}$.

Observe that $V_{n}$ is a degenerate Gaussian vector with mean zero and
covariance matrix

$\begin{split}C&=\sigma_{\text{prp}}^{2}N_{x}N_{x}^{T}\\\
&=\sigma_{\text{prp}}^{2}U\Sigma^{-1}\Sigma^{-T}U^{T}\end{split}$ (2.7)

(this is degenerate since the rank$(N_{x}N_{x}^{T})$ = rank$(\nabla q(x)\nabla
q(x)^{T})=m<d$). Therefore, $V_{n}=N_{x}R_{n}$ has a density

$\displaystyle
p_{\text{off}}^{\perp}(v_{n})=\frac{1}{\sqrt{(2\pi)^{m}\det^{*}(C)}}e^{-\frac{1}{2}v_{n}^{T}C^{+}v_{n}}$
(2.8)

where $C^{+}$ denotes the pseudo-inverse of covariance matrix $C$ and
$\det^{*}(C)$ denotes its pseudo-determinant (the product of its non-zero
singular values). We will verify in Section 2.4 that (2.8) is equivalent to

$\displaystyle p_{\text{off}}^{\perp}(v_{n})=\frac{\det(\nabla q(x)^{T}\nabla
q(x)^{T})^{1/2}}{(2\pi)^{m/2}\sigma_{\text{prp}}^{m}}\exp\Big{(}-v_{n}^{T}\frac{\nabla
q(x)\nabla q(x)^{T}}{2\sigma_{\text{prp}}^{2}}v_{n}\Big{)}.$ (2.9)

We will also verify in Section 2.4 that (2.9) can be written in a more compact
form as

$\displaystyle
p_{\text{off}}^{\perp}(r_{n})=\frac{\Pi_{i=1}^{m}\sigma_{i}}{(2\pi)^{m/2}\sigma_{\text{prp}}^{m}}\exp\Big{(}-\frac{|r_{n}|^{2}}{2\sigma_{\text{prp}}^{2}}\Big{)},$
(2.10)

where $v_{n}=N_{x}r_{n}$. By similar arguments, we can write the density for
$V_{t}=T_{x}R_{t}$, the step in the tangent direction, as

$\displaystyle
p_{\text{off}}^{\|}(v_{t})=\frac{1}{(2\pi)^{d/2}\sigma_{\text{tan}}^{d}}\exp\Big{(}-\frac{|v_{t}|^{2}}{2\sigma_{\text{tan}}^{2}}\Big{)},$
(2.11)

or equivalently

$\displaystyle
p_{\text{off}}^{\|}(r_{t})=\frac{1}{(2\pi)^{d/2}\sigma_{\text{tan}}^{d}}\exp\Big{(}-\frac{|r_{t}|^{2}}{2\sigma_{\text{tan}}^{2}}\Big{)}.$
(2.12)

Then, the overall point proposal density for an Off move is given by

$\displaystyle\begin{split}h_{21}(x,y)&=p_{\text{off}}^{\perp}(r_{n})p_{\text{off}}^{\|}(r_{t})\\\
\\\
&=\frac{\Pi_{i=1}^{m}\sigma_{i}}{(2\pi)^{(m+d)/2}\sigma_{\text{prp}}^{m}\sigma_{\text{tan}}^{d}}\exp\Big{(}-\frac{|r_{n}|^{2}}{2\sigma_{\text{prp}}^{2}}-\frac{|r_{t}|^{2}}{2\sigma_{\text{tan}}^{2}}\Big{)}.\end{split}$
(2.13)

Before accepting $y$, we must check that the reverse move is possible. That
is, we must check whether it is possible to propose $x$ starting from $y$.
This requires knowledge of the On move proposal process, therefore we
momentarily postpone this discussion to Section 2.3. In any case, if the
reverse check is successful, we set the proposal to $(y,1)$, otherwise we stop
and set $(X_{k+1},i_{k+1})=(x,2)$.

#### On move: $x\in\Omega_{1},\ y\in\Omega_{2}.$

Given $x\in S^{c}$, an Off move generates a point proposal $y\in S$. We start
by projecting $x$ onto the hard constraint $S$. The projection step finds a
vector $a_{1}\in\mathbb{R}^{m}$ such that

$\displaystyle q\big{(}x+\nabla q(x)a_{1}\big{)}=0.$ (2.14)

We solve equation (2.14) by standard Newton’s method. If the projection fails,
we must immediately stop and set $(X_{k+1},i_{k+1})=(x,1)$. If the projection
is successful, we take the step $x\rightarrow x_{s}=x+\nabla q(x)a_{1}\in S$
and move on to the rest of the proposal generation process. Note that in the
soft constraints limit we can assume that $\text{dist}(x,S)=O(\varepsilon)$,
so that the column space of $\nabla q(x)$ will give a good approximation to
the normal space of $S$ at $x_{s}$, $\mathcal{N}_{x_{s}}$. If we were to stop
here and simply propose $x_{s}$, we would have a transition kernel that
involves only delta-function components. This is because, conditional on $x$,
$x_{s}$ has been generated in a completely deterministic way, which makes
detailed balance impossible. Hence, we need to introduce some non-
deterministic component to the proposal generation process. Given $x_{s}$, we
generate the proposal by taking a surface move to a point $y=x_{s}+v+w\in S$,
where $v\in\mathcal{T}_{x_{s}}$ and $w\in\mathcal{N}_{x_{s}}$. As with the
hard move, we first take a step $x_{s}\rightarrow x_{s}+v$ with
$v\in\mathcal{T}_{x_{s}}$. We generate $v$ by sampling an isotropic Gaussian
in $\mathcal{T}_{x_{s}}$ with variance $\sigma_{\text{on}}^{2}$. This tangent
space move has a density

$\displaystyle
p_{\text{on}}(v)=\frac{1}{(2\pi)^{d/2}\sigma_{\text{on}}^{d}}e^{-\frac{|v|^{2}}{2\sigma_{\text{on}}^{2}}}$

We then project to $x+v$ back onto $S$ by finding a $w\in\mathcal{N}_{x_{s}}$
such that $x_{s}+v+w\in S$. This projection step is carried out by using
Newton’s method to find an $a_{2}\in\mathbb{R}^{m}$ such that
$q\big{(}x_{s}+v+\nabla q(x_{s})a_{2}\big{)}=0$, and then setting $w=\nabla
q(x_{s})a_{2}\in\mathcal{N}_{x_{s}}$. Note that here, as for the Hard move, we
need to reject the proposal in case the Newton projection step is
unsuccessful. The only difference with the Hard move is that in case of an On
move there is no need to do the reverse check. In fact, given $x$ and $y$, an
Off move from $y$ to $x$ is always possible. We can construct the reverse Off
move as $x=y+V^{\prime}_{n}+V^{\prime}_{t}$ (where
$V^{\prime}_{n}\in\mathcal{N}_{y}$ and $V^{\prime}_{t}\in\mathcal{T}_{y}$) by
projecting $x-y$ onto $\mathcal{T}_{y}$ and
$\mathcal{N}_{y}=\mathcal{T}_{y}^{\perp}$.

The point proposal density for the On move is given by

$\displaystyle\begin{split}h_{12}(x,y)&=p_{\text{on}}(v)J(x_{s},y)\\\ \\\
&=\frac{1}{(2\pi)^{d/2}\sigma_{\text{on}}^{d}}e^{-\frac{|v|^{2}}{2\sigma_{\text{on}}^{2}}}\det(T_{x_{s}}^{T}T_{y}),\end{split}$
(2.15)

where, as before, $J(x_{s},y)=\det(T_{x_{s}}^{T}T_{y})$ is the Jacobian factor
that accounts for the area distortion induced by the projection map
$v\rightarrow y$ (see Section 2.5 for details on how to compute the Jacobian).

#### Soft move: $x,\ y\in\Omega_{1}.$

Given a point $x\in S^{c}$, we propose $y\in S^{c}$. For this, we sample an
isotropic Gaussian in $\mathbb{R}^{d_{a}}$ with mean zero and variance
$\sigma_{\text{sft}}^{2}$. We denote this sample by $v$ and set $y=x+v$. The
point proposal density for the Soft move is just

$\displaystyle
h_{11}(x,y)=\frac{1}{(2\pi)^{d_{a}/2}\sigma_{\text{sft}}^{d_{a}}}e^{-\frac{|v|^{2}}{2\sigma_{\text{sft}}^{2}}}.$
(2.16)

Note that the Soft move is just a special case of a Hard move in
$\mathbb{R}^{d_{a}}$.

### 2.3 Reverse check for the Off move

Here we discuss the reverse check process for the Off move. Given $x\in S$,
assume that the Off proposal has generated a point $y=x+V_{n}+V_{t}\in S^{c}$,
where $V_{n}\in\mathcal{N}_{x}$ and $V_{t}\in\mathcal{T}_{x}$. We must check
whether the On proposal can generate $x$, starting from $y$. The first step is
to project $y$ back onto $S$. The reverse projection step attempts to find a
vector $a\in\mathbb{R}^{m}$ such that $q\big{(}y+\nabla q(y)a\big{)}=0$. If
Newton’s method fails to converge, we must stop immediately and reject the Off
move point proposal $y$. On the other hand, if Newton’s method converges to
some $y_{s}=x+\nabla q(y)a\in S$, we need to check whether we can get from
$y_{s}$ to $x$ with a reverse surface move. The reverse surface move is
constructed as $x=y_{s}+v^{\prime}+w^{\prime}$ where
$v^{\prime}\in\mathcal{T}_{y_{s}}$ and $w^{\prime}\in\mathcal{N}_{y_{s}}$.
Since we already know both $x$ and $y$, we can always find $v^{\prime}$ and
$w^{\prime}$ by projecting $x-y_{s}$ onto $\mathcal{T}_{y_{s}}$ and
$\mathcal{N}_{y_{s}}=\mathcal{T}_{y_{s}}^{\perp}$. Thus, we only must check
whether the Newton solver would find $x$ starting from $y_{s}+v^{\prime}$. If
not, we reject $(y,1)$ and set $(X_{k+1},i_{k+1})=(x,2)$.

### 2.4 Density for the Off move

In this section we will verify that (2.8), (2.9) and (2.10) are all
equivalent. This will follow from the following four facts. First,

$\displaystyle\begin{split}C^{+}&=\frac{1}{\sigma_{\text{prp}}^{2}}\nabla
q(x)\nabla q(x)^{T},\end{split}$ (2.17)

second,

$\displaystyle\det(\nabla q(x)^{T}\nabla
q(x)^{T})^{1/2}=\Pi_{i=1}^{m}\sigma_{i},$ (2.18)

where the $\sigma_{i}$’s are the singular values in $\Sigma$. Third,

$\displaystyle\text{det}^{*}(C)=\frac{\sigma_{\text{prp}}^{2m}}{\det(\nabla
q(x)^{T}\nabla q(x))},$ (2.19)

and fourth,

$\displaystyle|\nabla q(x)^{T}v_{n}|=|r_{n}|,$ (2.20)

where $v_{n}=N_{x}r_{n}$. Recall that $\nabla q(x)=U\Sigma V^{T}$,
$N_{x}=U\Sigma^{-1}V^{T}$ and $C=\sigma_{\text{prp}}^{2}N_{x}N_{x}^{T}$. We
can verify (2.17) by

$\displaystyle\begin{split}C^{+}&=\frac{1}{\sigma_{\text{prp}}^{2}}U\Sigma\Sigma^{T}U^{T}\\\
&=\frac{1}{\sigma_{\text{prp}}^{2}}\nabla q(x)\nabla q(x)^{T}.\end{split}$

Let the $\sigma_{i}$’s be the singular values in $\Sigma$. (2.18) holds since

$\displaystyle\det(\nabla q(x)^{T}\nabla q(x))=\det(V^{T}\Sigma^{T}\Sigma
V)=\Pi_{i=1}^{m}\sigma_{i}^{2},$

and for (2.19) we have

$\displaystyle\begin{split}\text{det}^{*}(C)&=\text{det}^{*}(\sigma_{\text{prp}}^{2}U\Sigma^{-1}\Sigma^{-T}U^{T})\\\
&=\sigma_{\text{prp}}^{2m}\text{det}^{*}(\Sigma^{-1}\Sigma^{-T})=\frac{\sigma_{\text{prp}}^{2m}}{\Pi_{i=1}^{m}\sigma_{i}^{2}}=\frac{\sigma_{\text{prp}}^{2m}}{\det(\nabla
q(x)^{T}\nabla q(x))}.\end{split}$

We will now verify (2.20). We have

$\displaystyle-v_{n}^{T}$ $\displaystyle\nabla q(x)\nabla
q(x)^{T}v_{n}=-r_{n}^{T}N_{x}^{T}U\Sigma\Sigma^{T}U^{T}N_{x}r_{n}=r_{n}^{T}r_{n}.$

### 2.5 Calculating the Jacobian

In this section we will show how calculate the Jacobian factor appearing in
the point proposal densities for the Hard and On moves. Similar calculations
based on linearization of the projection map are presented in [1] and [4]. Let
$y=x+v+w$ be the surface sampler proposal point, and suppose $v\rightarrow
v+\Delta v$ and $w\rightarrow w+\Delta w$, where $\Delta v\in\mathcal{T}_{x}$,
$\Delta w\in\mathcal{N}_{x}$ and $\Delta y\in\mathcal{T}_{y}$. Then, there
exists vectors $a,b,c$ such that $\Delta v=T_{x}a$, $\Delta w=N_{x}b$ and
$\Delta y=T_{y}c$. We have that $\Delta y=\Delta v+\Delta w$, which can be
rewritten as $T_{y}c=T_{x}a+N_{x}b$. If we multiply the last equation by
$T_{x}^{T}$ we get

$\displaystyle T_{x}^{T}T_{y}c=a,$ (2.21)

where we used that $T_{x}^{T}T_{x}=I$ (by orthonormality of $T_{x}$) and that
$T_{x}^{T}N_{x}=0$ (by $\mathcal{N}_{x}=\mathcal{T}_{x}^{\perp}$). The
Jacobian of the map $v\rightarrow y$ is thus given by $(T_{x}^{T}T_{y})^{-1}$.
From this, we conclude that
$J(x,y)=(\det(T_{x}^{T}T_{y})^{-1})^{-1}=\det(T_{x}^{T}T_{y})$.

## 3 Choosing the Correct Scaling: the Flat Surface Case

In Section 2.2 we found the point-proposal densities for each move type.
However, we did not provide any details on how to choose the scale parameters
$\sigma_{\text{prp}},\sigma_{\text{tan}},\sigma_{\text{on}},\sigma_{\text{hrd}},\sigma_{\text{sft}}$.
In this section we will specialize to the setting where $S$ is a flat manifold
defined by linear constraints (such as a line or a plane). In this case, we
will see that it is possible to choose a scaling for the On and Off moves that
produces a $100\%$ acceptance probability. In other words, we will determine
$\sigma_{\text{prp}},\sigma_{\text{tan}}$, and $\sigma_{\text{on}}$ such that
the detailed balance condition

$\displaystyle\begin{split}f_{2}(x)&\lambda_{21}(x)h_{21}(x,y)=f_{1}(y)\lambda_{12}(y)h_{12}(y,x),\quad\forall\
x\in\Omega_{2},y\in\Omega_{1}\end{split}$ (3.1)

holds. This will imply that $a_{21}(x,y)=a_{12}(y,x)\equiv 1$. When $S$ is not
flat, the smoothness of the constraint function $q$ ensures that the surface
is locally well-approximated by a plane. Therefore, we can expect that the
same choice of $\sigma_{\text{prp}},\sigma_{\text{tan}}$, and
$\sigma_{\text{on}}$ will lead to $a_{12}(x,y),a_{21}(y,x)\nearrow 1$ as
$\varepsilon\downarrow 0$.

We mention here that the scale parameters for the Hard and Soft moves
($\sigma_{\text{hrd}}$ and $\sigma_{\text{sft}}$) should be $O(1)$ and
$O(\varepsilon)$, respectively. The reason for having a Hard move is that its
scaling is not affected by the size of $\varepsilon$, that is, we can have
$\sigma_{\text{hrd}}=O(1)$. This move should substantially de-correlate the
samples and produce a method that is effective in the sense of Section 1.
Nothing much can be done about the Soft move. As we mentioned before, a
typical sample $X$ from (1.2) will have $\text{dist}(X,S)=O(\varepsilon)$.
Therefore, the size for the isotropic Gaussian proposal (2.16) will have to be
to be small, that is, we need to have $\sigma_{\text{sft}}=O(\varepsilon)$.

### 3.1 Detailed balance with a flat surface

In this section we will specialize to the flat surface case and find a scaling
for the On and Off moves for which condition (3.1) holds. Suppose we have a
linear constraint function and a surface of the form

$\displaystyle q(x)=A^{T}x,\quad S=\big{\\{}x\in\mathbb{R}^{d_{a}}\ :\
q(x)=0\big{\\}}.$ (3.2)

where $A\in\mathbb{R}^{d_{a}\times m}$ is assumed to have full column rank
$m$. Let $A=U\Sigma V^{T}$ be the singular value decomposition for $\nabla
q(x)=A$, and define the matrix $N=U\Sigma^{-1}V^{T}$. Note that in this case
$\mathcal{N}_{x}=\text{Col}(A)=\text{Col}(N)$ and
$\mathcal{T}_{x}=\text{Null}(A^{T})=S$, for any $x\in S$. We can rewrite the
distribution in (1.2) by,

$\displaystyle\begin{split}\pi_{\varepsilon}(dx)&=\frac{1}{Z_{\varepsilon}}\exp\Big{(}-x^{T}\frac{AA^{T}}{2\varepsilon^{2}}x\Big{)}dx.\end{split}$
(3.3)

We see that $\pi_{\varepsilon}$ is just a degenerate (since $AA^{T}$ has rank
$m<d_{a}$) Gaussian distribution with covariance matrix
$\varepsilon^{2}NN^{T}$ and normalizing constant

$\displaystyle\begin{split}Z_{\varepsilon}&=\sqrt{(2\pi)^{m}\varepsilon^{2m}\text{det}^{*}(NN^{T})}\\\
&=\sqrt{(2\pi)^{m}\varepsilon^{2m}\text{det}^{*}(\Sigma^{-1}\Sigma^{-T})}\\\
&=\sqrt{(2\pi)^{m}\varepsilon^{2m}\Pi_{i=1}^{m}\sigma_{i}^{-2}}\\\
&=\sqrt{\frac{(2\pi)^{m}\varepsilon^{2m}}{\det(A^{T}A)}}.\end{split}$ (3.4)

Therefore, we can rewrite (3.3) as

$\displaystyle\begin{split}\pi_{\varepsilon}(dx)&=\frac{\det(A^{T}A)^{1/2}}{(2\pi)^{m/2}\varepsilon^{m}}\exp\Big{(}-x^{T}\frac{AA^{T}}{2\varepsilon^{2}}x\Big{)}dx.\end{split}$
(3.5)

Now consider an Off move from $x\in\Omega_{2}=S$ to point
$y\in\Omega_{1}=S^{c}$. The reverse On move has the form
$x=y_{s}+v_{\text{rev}}+w_{\text{rev}}$, where $y_{s}\in S$,
$v_{\text{rev}}\in\mathcal{T}_{y_{s}}$ and
$w_{\text{rev}}\in\mathcal{N}_{y_{s}}$. Recall that $y_{s}$ is determined by
solving the linear system $q\big{(}y+\nabla q(y)a\big{)}=0$ for $a$ via
Newton’s method. Newton’s method finds the solution in one step222The solution
is $a=-(A^{T}A)^{-1}A^{T}y$. and returns
$y_{s}=\text{proj}_{S}(y)=y+z_{\text{rev}}$, where $z_{\text{rev}}=\nabla
q(y)a\in\mathcal{N}_{y_{s}}$. Since $S$ is flat, we know that the tangent step
in the surface move, $v_{\text{rev}}$, will not bring us off the surface. We
conclude that the surface move projection step is unnecessary and set
$w_{\text{rev}}=0$. To sum up, we have $x=y+z_{\text{rev}}+v_{\text{rev}}$
with $z_{\text{rev}}\in\mathcal{N}_{y_{s}}$ and
$v_{\text{rev}}\in\mathcal{T}_{y_{s}}$. The point proposal density for the
reverse On move is given by $h_{12}(y,x)=p_{\text{on}}(v_{\text{rev}})$, where
there is no Jacobian factor since $\det(T_{y_{s}}^{T}T_{x})\equiv 1$333Since
$S$ is flat, we have $T_{x}=T$ and by orthonormality
$T_{x}^{T}T_{y}=T_{y_{s}}^{T}T_{x}=T^{T}T=I$.. Detailed balance condition
(3.1) can be written as

$\displaystyle\begin{split}h_{21}(x,y)&=\frac{\lambda_{12}(y)f_{1}(y)}{\lambda_{21}(x)f_{2}(x)}\
h_{12}(y,x)\\\ \\\ &=\frac{\lambda_{12}}{\lambda_{21}}\
\frac{f_{1}(y)}{f_{2}(x)}\ p_{\text{on}}(v_{\text{rev}})\\\ \\\
&=\frac{\exp\big{(}-y^{T}\frac{AA^{T}}{2\varepsilon^{2}}y\big{)}}{\frac{\lambda_{21}}{\lambda_{12}k_{1}}f_{2}(x)}\
\frac{\exp\big{(}-\frac{|v_{\text{rev}}|^{2}}{2\sigma_{\text{on}}^{2}}\big{)}}{(2\pi)^{d/2}\sigma_{\text{on}}^{d}},\end{split}$
(3.6)

where we used that
$f_{1}(y)=k_{1}\exp\big{(}-y^{T}\frac{AA^{T}}{2\varepsilon^{2}}y\big{)}$. Note
that (3.6) suggests how to construct the Off move point proposal so that the
detailed balance condition (3.1) holds. To balance the two Gaussian terms on
the right hand side of (3.6), the Off move should consist of a combination of
$(i)$ a Gaussian step in the normal space with basis $N_{x}=N$ and $(ii)$ an
isotropic Gaussian step in tangent space (with any choice of orthonormal basis
$T_{x}=T$). This corresponds to the Off proposal mechanism we presented in
Section 2.2, where the proposal is set to $y=x+v$ with $v=v_{n}+v_{t}$ and
$v_{n}=Nr_{n}\in\mathcal{N}_{x}$ and $v_{t}=Tr_{t}\in\mathcal{T}_{x}$. We can
use (2.9), (2.11) to write the Off proposal density as

$\displaystyle\begin{split}h_{21}(x,y)&=p_{\text{off}}^{\perp}(v_{n})p_{\text{off}}^{\|}(v_{t})\\\
\\\
&=\frac{\exp\big{(}-v_{n}^{T}\frac{AA^{T}}{2\sigma_{\text{prp}}^{2}}v_{n}\big{)}}{\frac{(2\pi)^{m/2}\sigma_{\text{prp}}^{m}}{\det(A^{T}A)^{1/2}}}\
\frac{\exp\big{(}-\frac{|v_{t}|^{2}}{2\sigma_{\text{tan}}^{2}}\big{)}}{(2\pi)^{d/2}\sigma_{\text{tan}}^{d}},\end{split}$
(3.7)

where we omitted the Jacobian factor since $\det(T_{x}^{T}T_{y})\equiv 1$. If
we substitute (3.7) into (3.6) we obtain

$\displaystyle\begin{split}\frac{\exp\big{(}-v_{n}^{T}\frac{AA^{T}}{2\sigma_{\text{prp}}^{2}}v_{n}\big{)}}{\frac{(2\pi)^{m/2}\sigma_{\text{prp}}^{m}}{\det(A^{T}A)^{1/2}}}\
\frac{\exp\big{(}-\frac{|v_{t}|^{2}}{2\sigma_{\text{tan}}^{2}}\big{)}}{(2\pi)^{d/2}\sigma_{\text{off}}^{d}}=\frac{\exp\big{(}-y^{T}\frac{AA^{T}}{2\varepsilon^{2}}y\big{)}}{\frac{\lambda_{21}}{\lambda_{12}k_{1}}f_{2}(x)}\
\frac{\exp\big{(}-\frac{|v_{\text{rev}}|^{2}}{2\sigma_{\text{on}}^{2}}\big{)}}{(2\pi)^{d/2}\sigma_{\text{on}}^{d}}.\end{split}$
(3.8)

Now observe that we must have $v_{\text{rev}}=-v_{t}$, so that if we choose
$\sigma_{\text{tan}}=\sigma_{\text{on}}$, then the Gaussian densities for the
step in tangent space simplify in (3.8). Moreover, since $y=x+v_{n}+v_{t}$
with $x,v_{t}\in S=\text{Null}(A^{T})$, we have $A^{T}y=A^{T}v_{n}$.
Therefore, (3.8) becomes

$\displaystyle\begin{split}\frac{\exp\big{(}-v_{n}^{T}\frac{AA^{T}}{2\sigma_{\text{prp}}^{2}}v_{n}\big{)}}{\frac{(2\pi)^{m/2}\sigma_{\text{prp}}^{m}}{\det(A^{T}A)^{1/2}}}=\frac{\exp\big{(}-v_{n}^{T}\frac{AA^{T}}{2\varepsilon^{2}}v_{n}\big{)}}{\frac{\lambda_{21}}{\lambda_{12}k_{1}}f_{2}(x)}.\end{split}$
(3.9)

This last equation tells us the scaling for the normal step component of an
Off move. Indeed, to simplify the exponential factors in (3.9) we have to take
$\sigma_{\text{prp}}=\varepsilon$. Equation (3.9) also suggests what the
density on the surface, $f_{2}$, should be. In fact, we can rewrite (3.9) as

$\displaystyle\begin{split}\frac{f_{2}(x)}{k_{1}}&=\frac{\lambda_{12}(2\pi)^{m/2}\varepsilon^{m}}{\lambda_{21}}\det(A^{T}A)^{-1/2}.\end{split}$
(3.10)

Thanks to our choice of surface density $f_{2}(x)=k_{2}\det(\nabla
q(x)^{T}\nabla q(x))^{-1/2}=k_{2}\det(A^{T}A)^{-1/2}$, the gradient factors
simplify. Therefore, (3.10) becomes

$\displaystyle\begin{split}\frac{k_{2}}{k_{1}}=\frac{\lambda_{12}(2\pi)^{m/2}\varepsilon^{m}}{\lambda_{21}},\end{split}$
(3.11)

which can obviously be reduced to an identity since we are free to choose
constants $k_{1}$ and $k_{2}$.

To summarize, we showed that when $S$ is a flat surface, if we set
$\sigma_{\text{prp}}=\varepsilon$, $\sigma_{\text{tan}}=\sigma_{\text{on}}$
then detailed balance condition (3.1) holds, provided that our choice of
$k_{1},k_{2}$ satisfies (3.11).

## 4 Experiments

In this section we illustrate the ith two examples of hard constraint surface:
the intersection of two spheres and the intersection of an ellipsoid with a
sphere.

For the rest of this section, we will use $X=(X^{1},X^{2},X^{3})$ to denote a
random variable in $\mathbb{R}^{3}$ with distribution $\pi_{\varepsilon}$ as
in (1.2). We will also use $x=(x_{1},x_{2},x_{3})$ to denote a generic point
in $\mathbb{R}^{3}$.

In what follows we will estimate the auto-covariance function and the auto-
correlation time of MCMC samples $X_{1},X_{2},...,X_{N}$ drawn from
$\pi_{\varepsilon}$. Consider a general estimator

$\displaystyle\hat{S}_{N}=\frac{1}{N}\sum_{j=1}^{N}F(X_{j}),$ (4.1)

for some given function $F(\cdot)$. Assuming that $\\{X_{j}\\}_{j=1}^{N}$ are
sampled from (1.2), the equilibrium lag $t$ auto-covariance is

$\displaystyle C_{t}=\text{Cov}[F(X_{j}),F(X_{j+t})].$ (4.2)

Note that this is an even function of $t$ since $C_{t}=C_{-t}$. We estimate
the auto-covariance function from the MCMC samples as

$\displaystyle
C_{t}\approx\frac{1}{N-t}\sum_{j=1}^{N-t}\big{(}F(X_{j})-\hat{S}_{N}\big{)}\big{(}F(X_{j+t})-\hat{S}_{N}\big{)}$
(4.3)

The integrated auto-correlation time, denoted by $\tau$, is defined as

$\displaystyle\tau=\sum_{t=-\infty}^{\infty}\frac{C_{t}}{C_{0}}=1+\frac{2}{\text{Var}(F(X_{j}))}\sum_{t=1}^{\infty}C_{t}$
(4.4)

We use the self-consistent window approach (described in [3]) to estimate
$\tau$. The self-consistent window approach estimates $\tau$ by using a number
of samples $M\ll N$. This produces an estimator with lower variance at the
cost of increasing the bias. The point of having an accurate estimate for the
integrated auto-correlation time is that

$\displaystyle\text{Var}(\hat{S}_{N})\approx\frac{\tau}{N}C_{0},$ (4.5)

so that the error in the estimate decreases as $1/\sqrt{N_{\text{eff}}}$,
where $N_{\text{eff}}=N/\tau$ is the effective sample size. The integrated
auto-correlation time is a measure of how many MCMC steps one has to take to
generate approximately independent samples.

Note that instead of estimating the auto-covariance function from the sum in
(4.3) (this takes $O(N^{2})$ operations), it is computationally more efficient
to take the fast Fourier transform of the entire sample (this takes $O(N\log
N)$ operations). This uses the fact that the auto-correlation of a continuous-
time deterministic signal can be written as a convolution, which in Fourier
space corresponds to multiplication.

Note that for most of our experiments we generated approximately $5\times
10^{6}$ samples from the augmented target distribution (2.1). We used
parameters $\lambda_{11}=0.2$ (Soft move), $\lambda_{12}=0.8$ (On move),
$\lambda_{22}=0.8$ (Hard move), $\lambda_{21}=0.2$ (Off move) and
$\sigma_{\text{prp}}=\sigma_{\text{tan}}=\sigma_{\text{on}}=\varepsilon$,
$\sigma_{\text{hrd}}=1$ and $\sigma_{\text{sft}}=0.7\times\varepsilon$. About
$20\%$ of the samples were off the hard constraint surface $S$. That is, with
these parameters we obtained $N\approx 10^{6}$ samples from
$\pi_{\varepsilon}$, which we then used for analysis (samples on the surface
$S$ were discarded since these are drawn from $\pi_{s}$).

### Two Spheres model

Consider two soft constraints $q_{i}:\mathbb{R}^{3}\rightarrow\mathbb{R}$
given by

$\displaystyle q_{i}(x)=|x-c_{i}|^{2}-r_{i}^{2},$

where $c_{1},c_{2}\in\mathbb{R}^{3}$ and $r_{1},r_{2}$ are positive real
numbers. Geometrically, the hard constraint surface
$S=\big{\\{}x\in\mathbb{R}^{3}\ :\ q(x)=0\big{\\}}$ is the intersection of two
spheres with centers $c_{1}$, $c_{2}$ and radii $r_{1}$, $r_{2}$,
respectively. In our experiments, we use centers $c_{1}=(0,0,1)$,
$c_{2}=(0,-1,0)$ and radii $r_{1}=r_{2}=\sqrt{2}$. In this way, the point
$s=(1,0,0)\in S$. The surface $S$ consists of a circle with center
$c=(0,-1/2,1/2)$ (the midpoint between $c_{1}$ and $c_{2}$) and radius
$r=|s-c|=\sqrt{3/2}$. The circle lies on the plane perpendicular to the line
segment between $c_{1}$ and $c_{2}$, that is, it lies on the plane
$P=\big{\\{}x\in\mathbb{R}^{3}\ :\ (x-p)\cdot\vec{n}=0\big{\\}}$, where
$\vec{n}=(0,1,1)$.

In general, this problem is too simple since the gradient factor $\det(\nabla
q(x)^{T}\nabla q(x))$ is constant and the Hard move does not have any
Metropolis or projection failures. However, this setting allows for a simple
correctness check which involves the angle $\theta$ between the projection of
$X\sim\pi_{\varepsilon}$ onto $P$ and any reference vector perpendicular to
$\vec{n}$ (say $w=s-c$). Since in this case the gradient factor in $\pi_{s}$
is constant, the angle $\theta$ is uniformly distributed.

We performed the $\theta$-check with $\varepsilon=0.022$ (inverse temperature
$\beta=1/2\varepsilon^{2}=10^{3}$). We computed the angle between the
projection of each sample onto $P$ and the reference vector $w$. Figure 1
displays an histogram for the sampled angles. The empirical distribution of
the angle $\theta$ looks correct within statistical sampling error.

We also computed the sample auto-covariance function and auto-correlation time
for decreasing values of $\varepsilon$. In particular, we took inverse
temperatures $\beta=1/2\varepsilon^{2}=10^{2},10^{3},10^{4},2\times 10^{4}$
and used $F(x)=x_{1}$, the projection map onto the first coordinate. Figure 1
contains plots of the auto-covariance structure and lists the integrated auto-
correlation time for various values of $\varepsilon$. We see that the auto-
correlation time stays approximately constant as $\varepsilon$ decreases.

(a) Two-Sphere Model

(b) Two-Sphere Model

Figure 1. (A): Sample distribution and theoretical density for $\theta$ with
$\varepsilon=0.022$. (B): Sample auto-covariance function and auto-correlation
time with $F(x)=x_{1}$ for decreasing values of $\varepsilon$.

As another check, we also estimate the distribution of $X^{1}$, the first
coordinate of $X\sim\pi_{\varepsilon}$. To do this, we divide the $x_{1}$-axis
into $L$ bins. Note that we can compute the un-normalized density of $X^{1}$,

$\displaystyle p_{X^{1}}(x_{1})=\iint p_{X}(x_{1},x_{2},x_{3})dx_{2}dx_{3},$

via numerical quadrature444We use the midpoint rule.. We use estimators

$\displaystyle\hat{R}_{N}^{i}=\frac{\frac{1}{N}\sum_{j=1}^{N}\mathds{1}_{B_{i}}(X^{1}_{j})}{p_{X^{1}}(b_{i})|B_{i}|},\quad\quad
i=1,...,L$ (4.6)

where $B_{i}$ is the $i$-th bin with center $b_{i}$ and $X^{1}_{j}$ is the
first coordinate of the $j$-th sample from the MCMC sampler. If the sampler is
correct, we should find that the random variables $\hat{R}_{N}^{i}$ are
roughly constant over the bin index $i$ since

$\displaystyle\begin{split}\hat{R}_{N}^{i}&=\frac{\frac{1}{N}\sum_{j=1}^{N}\mathds{1}_{B_{i}}(X^{1}_{j})}{p_{X^{1}}(b_{i})|B_{i}|}\\\
\\\ &\approx\frac{\text{Pr}\big{(}X^{1}\in
B_{i}\big{)}}{\int_{B_{i}}p_{X^{1}}(x_{1})dx_{1}}=\frac{1}{Z},\end{split}$

where $Z=\int p_{X^{1}}(x_{1})dx_{1}$. This is what we see in Table (A) of
Figure 2.

(a) Two-Sphere Model

(b) Ellipsoid-Sphere Model

Figure 2. Tables (A) and (B) display bin index (”bin”), center of the bin
(”center”), number of bin hits (”count”), exact marginal density at bin center
point (”pdf”), and the value of the random variable $\hat{R}^{i}_{N}$ (”1/Z”).
For both tables, we generated $N\approx 2\times 10^{6}$ samples from
$\pi_{\varepsilon}$ with $\varepsilon=0.022$ (i.e., inverse temperature
$\beta=1/2\varepsilon^{2}=1000$).

### Ellipsoid-Sphere model

Consider two soft constraints $q_{i}:\mathbb{R}^{3}\rightarrow\mathbb{R}$
given by

$\displaystyle\begin{split}q_{1}(x)&=|x-c_{1}|^{2}-r_{1}^{2},\\\ \\\
q_{2}(x)&=\frac{(x_{1}-c_{2,1})^{2}}{s_{1}^{2}}+\frac{(x_{2}-c_{2,2})^{2}}{s_{2}^{2}}+\frac{(x_{3}-c_{2,3})^{2}}{s_{3}^{2}}-1,\end{split}$

where $c_{1},c_{2}\in\mathbb{R}^{3}$ and $r_{1}$, $s_{1},s_{2},s_{3}$ are
positive real numbers. Geometrically, the hard constraint surface
$S=\big{\\{}x\in\mathbb{R}^{3}\ :\ q(x)=0\big{\\}}$ is the intersection of a
sphere centered at $c_{1}$ with radius $r_{1}$ and an ellipsoid centered at
$c_{2}$ (where $s_{1},s_{2},s_{3}$ are half the length of its principal axes).
In our experiments, we use centers $c_{1}=(0,0,1)$, $c_{2}=(0,-1,0)$, radius
$r_{1}=\sqrt{2}$ and $s_{1}=2,s_{2}=3,s_{3}=5$.

This problem is more complex than the two-spheres model discussed above
because in this case $\det(\nabla q(x)^{T}\nabla q(x))$ is non-constant. This
means that the Hard move will have both Metropolis and reverse check failures.

We estimated the $X^{1}$ marginal density in the same way as for the two-
sphere model, see Table (B) in Figure 2.

We also computed the sample auto-covariance function and auto-correlation time
for decreasing values of $\varepsilon$. As before, we took inverse
temperatures $\beta=1/2\varepsilon^{2}=10^{2},10^{3},10^{4},2\times 10^{4}$
and used $F(x)=x_{1}$, the projection map onto the first coordinate. See
Figure 3 part (A) for plots of the auto-covariance function and for the auto-
correlation time estimates. As expected, we found that the auto-correlation
time stays approximately constant as $\varepsilon$ decreases.

As a motivation for the n Figure 3 part (B) we display the auto-correlation
structure for the isotropic Gaussian proposal Metropolis sampler with an
$O(\varepsilon)$ step size. We tested this sampler for decreasing values of
$\varepsilon$. In particular, we took inverse temperatures
$\beta=1/2\varepsilon^{2}=5,10,15,20$ and used $F(x)=x_{1}$. We set a proposal
size to obtain an average acceptance probability of about $40\%$. Not
surprisingly, we observe a blow-up of the auto-correlation time even for
relatively small inverse temperatures (i.e., relatively large $\varepsilon$).

As a last display, in Table 1 we show the average acceptance probability for
Off and On moves for a decreasing sequence of $\varepsilon$-values. The data
are generated from the sing the ellipsoid-sphere model. We used step-size
parameters
$\sigma_{\text{prp}}=\sigma_{\text{tan}}=\sigma_{\text{on}}=\varepsilon$. As
expected, the average acceptance probabilities converge to $1$ as
$\varepsilon\downarrow 0$, confirming the analysis presented in Section 3.

(a) Ellipsoid-Sphere Model

(b) Ellipsoid-Sphere Model

Figure 3. (A): Auto-correlation structure for or decreasing values of $\varepsilon$. (B): Auto-correlation structure for Isotropic Gaussian Metropolis sampler with $O(\varepsilon)$ step size for decreasing values of $\varepsilon$. Value of $\varepsilon$ |  Av. Acc. Pr.
Off Move |  Av. Acc. Pr.
On Move
---|---|---
0.223 |  0.768781 |  0.763992
0.070 |  0.919511 |  0.919816
0.022 |  0.974695 |  0.974659
0.007 |  0.992233 |  0.99188
0.002 |  0.997509 |  0.997627
Table 1. Average acceptance probability for Off and On moves as
$\varepsilon\downarrow 0$. The chains are generated by the sing the ellipsoid-
sphere model. We used parameters
$\sigma_{\text{prp}}=\sigma_{\text{tan}}=\sigma_{\text{on}}=\varepsilon$.

## 5 Verifying the Invariant Measure

In this section we verify that the Markov chain generated by the reserves the
target distribution $\rho$ as in (2.1), which we rewrite here for convenience
of the reader

$\displaystyle\rho(dx)$
$\displaystyle=\frac{1}{Z}\Big{[}f_{1}(x)\mu_{1}(dx)+f_{2}(x)\mu_{2}(dx)\Big{]}.$
(5.1)

Recall that $\mu_{1}$ is the reference measure for the space
$(\Omega_{1},\mathcal{F}_{1})$ and that $\mu_{2}$ is the reference measure for
the space $(\Omega_{2},\mathcal{F}_{2})$, where $\mathcal{F}_{1}$,
$\mathcal{F}_{2}$ are the Borel sigma algebras for $\Omega_{1}$, $\Omega_{2}$,
respectively. The target $\rho$ defines a probability distribution in ambient
space $(\Omega,\mathcal{F})$, where $\Omega=\Omega_{1}\cup\Omega_{2}$ and
$\mathcal{F}$ is its Borel sigma algebra.

Note that we can extend the measures $\mu_{1}$ and $\mu_{2}$ to the entire
ambient space $(\Omega,\mathcal{F})$ so that these are mutually singular
measures on $(\Omega,\mathcal{F})$ which ”live” on $\Omega_{1}$ and
$\Omega_{2}$, respectively. We will assume this throughout the rest of the
section.

We will prove that the Markov transition kernel for the s self-adjoint in
$L_{\rho}^{2}(\Omega)$. As we will see, this implies that $\rho$ is the
invariant measure for the For each pair of labels $i,j\in\\{1,2\\}$, we use
$R_{ij}(x,dy)$ to denote the Markov transition kernel for the point-update
mechanism which generates a new sample $y\in\Omega_{j}$, given that the
current state is $x\in\Omega_{i}$ and that we’ve chosen new label $j$. This
kernel can be represented as

$\displaystyle
R_{ij}(x,dy)=r_{ij}(x)\delta_{x}(dy)+a_{ij}(x,y)h_{ij}(x,y)\mu_{j}(dy),$ (5.2)

where $a_{ij}(x,y)$ is the Metropolis-Hastings acceptance probability (2.2)
and $h_{ij}(x,y)$ is the point proposal density as defined in Section 2.2. We
also have

$\displaystyle\begin{split}r_{ij}(x)&=1-\int_{\Omega}a_{ij}(x,y)h_{ij}(x,y)\mu_{j}(dy)\\\
&=1-\int_{\Omega_{j}}a_{ij}(x,y)h_{ij}(x,y)\mu_{j}(dy),\end{split}$

the overall probability of rejecting a proposal starting from $x$. The kernel
$R_{ij}(x,dy)$ has the following two properties: $(i)$ for each
$x\in\Omega_{i}$, $R_{ij}(x,dy)$ defines a probability distribution on
$(\Omega,\mathcal{F})$ which ”lives” in $\Omega_{j}$ and $(ii)$ for each set
$A\in\mathcal{F}$, the map $x\mapsto R_{ij}(x,A)$ is measurable on
$(\Omega_{i},\mathcal{F}_{i})$. With this, we can write the overall transition
kernel for the s

$\displaystyle\begin{split}R(x,dy)=\lambda_{11}(x)&R_{11}(x,dy)+\lambda_{12}(x)R_{12}(x,dy)\\\
\\\ &+\lambda_{21}(x)R_{21}(x,dy)+\lambda_{22}(x)R_{22}(x,dy),\end{split}$
(5.3)

where $\lambda_{ij}(x)=\lambda_{ij}\mathds{1}_{\Omega_{i}}(x)$555Formally,
$R_{ij}(x,dy)$ is undefined when $x\notin\Omega_{i}$. However, we can extend
the map $x\mapsto R_{ij}(x,dy)$ to be measurable on the entire space
$(\Omega,\mathcal{F})$ by setting (say) $a_{ij}(x,y)=h_{ij}(x,y)\equiv 0$
whenever $x\notin\Omega_{i}$. In any case, if $x\notin\Omega_{i}$, then
transition kernel $R_{ij}(x,dy)$ doesn’t contribute to the overall kernel
$R(x,dy)$ in (5.3) since $R_{ij}(x,dy)$ is multiplied by
$\lambda_{ij}(x)=\lambda_{ij}\mathds{1}_{\Omega_{i}}(x)=0$. . Note that since
$\sum_{i,j}\lambda_{ij}(x)\equiv 1$, we can conclude that $R(x,dy)$ also
defines a Markov transition kernel. Equivalently, we can write

$R(x,dy)=\begin{cases}\lambda_{11}R_{11}(x,dy)+\lambda_{12}R_{12}(x,dy),&\text{if
}x\in\Omega_{1}\\\ \\\
\lambda_{21}R_{21}(x,dy)+\lambda_{22}R_{22}(x,dy),&\text{if
}x\in\Omega_{2}.\end{cases}$

We are now ready to show that $R$ is self-adjoint in $L_{\rho}^{2}(\Omega)$.
First, note that the Metropolis-Hastings acceptance rule (2.2) implies the
detailed balance condition

$\displaystyle
f_{i}(x)\lambda_{ij}(x)h_{ij}(x,y)a_{ij}(x,y)=f_{j}(y)\lambda_{ji}(y)h_{ji}(y,x)a_{ji}(y,x)$
(5.4)

for all pairs of labels $i,j\in\\{1,2\\}$ and for any $x\in\Omega_{i}$ and
$y\in\Omega_{j}$. We will use this to show that for any $f,g\in C_{b}(\Omega)$
we have

$\displaystyle(f,Rg)_{\rho}$
$\displaystyle=\sum_{i,j}(f,R_{ij}^{\lambda}g)_{\rho}$
$\displaystyle=\sum_{i,j}(R_{ij}^{\lambda}f,g)_{\rho}=(Rf,g)_{\rho},$

where $R_{ij}^{\lambda}$ is notation for the kernel

$\displaystyle R_{ij}^{\lambda}(x,dy)=\lambda_{ij}(x)R_{ij}(x,dy),$

and

$\displaystyle
R_{ij}^{\lambda}f(x)=\lambda_{ij}(x)\int_{\Omega}R_{ij}(x,dy)f(y).$

Let $f,g\in C_{b}(\Omega)$. In what follows, all the integrals are intended to
be over $\Omega$, so we will omit the domain subscripts for ease of notation.
For any $i,j\in\\{1,2\\}$, we have

$\displaystyle\begin{split}(f,R_{ij}^{\lambda}g)_{\rho}&=\iint
f(x)g(y)\lambda_{ij}(x)R_{ij}(x,dy)\rho(dx)\\\ \\\ &=\iint
f(x)g(y)\lambda_{ij}(x)\big{[}r_{ij}(x)\delta_{x}(dy)+h_{ij}(x,y)a_{ij}(x,y)\mu_{j}(dy)\big{]}\rho(dx)\\\
\\\ &=\int f(x)g(x)\lambda_{ij}(x)r_{ij}(x)\rho(dx)\ +\\\
&\quad\quad\quad\quad\quad\frac{1}{Z}\iint
f(x)g(y)f_{i}(x)\lambda_{ij}(x)h_{ij}(x,y)a_{ij}(x,y)\mu_{j}(dy)\mu_{i}(dx)\\\
\\\ &=\int f(x)g(x)\lambda_{ij}(x)r_{ij}(x)\rho(dx)\ +\\\
&\quad\quad\quad\quad\quad\frac{1}{Z}\iint
f(x)g(y)f_{j}(y)\lambda_{ji}(y)h_{ji}(y,x)a_{ji}(y,x)\mu_{j}(dy)\mu_{i}(dx)\\\
\\\ &=\int f(x)g(x)\lambda_{ij}(x)r_{ij}(x)\rho(dx)\ +\\\
&\quad\quad\quad\quad\quad\frac{1}{Z}\iint
f(y)g(x)f_{j}(x)\lambda_{ji}(x)h_{ji}(x,y)a_{ji}(x,y)\mu_{j}(dx)\mu_{i}(dy),\end{split}$
(5.5)

where for the third equality we used that $\lambda_{ij}(x)=0$ when
$x\notin\Omega_{i}$ (so that only the $f_{i}(x)\mu_{i}(dx)$ part of the target
$\rho(dx)$ matters in the integral), for the fourth equality we used the
detailed balance condition (5.4), and for the fifth equality we interchanged
the variables $x$ and $y$ in the double integral. When $i=j=1$, we see that

$\displaystyle\begin{split}(f,R_{11}^{\lambda}g)_{\rho}&=\iint
f(x)g(y)\lambda_{11}(x)R_{11}(x,dy)\rho(dx)\\\ \\\ &=\int
f(x)g(x)\lambda_{11}(x)r_{11}(x)\rho(dx)\ +\\\
&\quad\quad\quad\quad\quad\frac{1}{Z}\iint
f(y)g(x)f_{i}(x)\lambda_{11}(x)h_{11}(x,y)a_{11}(x,y)\mu_{i}(dx)\mu_{i}(dy)\\\
\\\ &=\iint
f(y)g(x)\lambda_{11}(x)R_{11}(x,dy)\rho(dx)=(R_{11}^{\lambda}f,g)_{\rho},\end{split}$
(5.6)

where for the second equality we used (5.5). Therefore,
$R_{11}^{\lambda}(x,dy)$ is self-adjoint. The case when $i=j=2$ is similar.

We are left to verify that

$\displaystyle(f,R_{12}^{\lambda}g)_{\rho}+(f,R_{21}^{\lambda}g)_{\rho}=(R_{12}^{\lambda}f,g)_{\rho}+(R_{21}^{\lambda}f,g)_{\rho}.$
(5.7)

Note that by (5.5), we can write

$\displaystyle\begin{split}(f,R_{12}^{\lambda}g)_{\rho}&=\int
f(x)g(x)\lambda_{12}(x)r_{12}(x)\rho(dx)\ +\\\ &\ \ \ \ \ \ \ \ \ \ \ \
\frac{1}{Z}\iint
f(y)g(x)f_{2}(x)\lambda_{21}(x)h_{21}(x,y)a_{21}(x,y)\mu_{2}(dx)\mu_{1}(dy),\end{split}$
(5.8)

and

$\displaystyle\begin{split}(f,R_{21}^{\lambda}g)_{\rho}&=\int
f(x)g(x)\lambda_{21}(x)r_{21}(x)\rho(dx)\ +\\\ &\ \ \ \ \ \ \ \ \ \ \ \
\frac{1}{Z}\iint
f(y)g(x)f_{1}(x)\lambda_{12}(x)h_{12}(x,y)a_{12}(x,y)\mu_{1}(dx)\mu_{2}(dy).\end{split}$
(5.9)

Finally, by simply applying the definitions we get

$\displaystyle\begin{split}(R_{12}^{\lambda}f,g)_{\rho}&=\iint
g(x)f(y)\lambda_{12}(x)R_{12}(x,dy)\rho(dx)\\\ \\\ &=\int
g(x)f(x)\lambda_{12}(x)r_{12}(x)\rho(dx)\ +\\\ &\ \ \ \ \ \ \ \ \ \ \ \
\frac{1}{Z}\iint
f(y)g(x)f_{1}(x)\lambda_{12}(x)h_{12}(x,y)a_{12}(x,y)\mu_{2}(dy)\mu_{1}(dx),\end{split}$
(5.10)

and

$\displaystyle\begin{split}(R_{21}^{\lambda}f,g)_{\rho}&=\iint
g(x)f(y)\lambda_{21}(x)R_{21}(x,dy)\rho(dx)\ +\\\ \\\ &=\int
g(x)f(x)\lambda_{21}(x)r_{21}(x)\rho(dx)\\\ &\ \ \ \ \ \ \ \ \ \ \ \
\frac{1}{Z}\iint
f(y)g(x)f_{2}(x)\lambda_{21}(x)h_{21}(x,y)a_{21}(x,y)\mu_{1}(dy)\mu_{2}(dx),\end{split}$
(5.11)

By combining (5.8), (5.9), (5.10) and (5.11), we see that (5.7) holds. We thus
proved that $(f,Rg)_{\rho}=(Rf,g)_{\rho}$ for any $f,g\in C_{b}(\Omega)$.
Since $C_{b}(\Omega)$ is dense in $L^{2}_{\rho}(\Omega)$, we conclude that $R$
is self-adjoint in $L^{2}_{\rho}(\Omega)$.

As a last step, we show that self-adjointness of $R$ in $L^{2}_{\rho}(\Omega)$
implies that $R$ preserves the target distribution $\rho$ as in (5.1). Indeed,
by self-adjointness of $R$ we have

$\displaystyle(f,1)_{\rho}=(f,R1)_{\rho}=(Rf,1)_{\rho},$ (5.12)

for any $f\in L^{2}_{\rho}(\Omega)$. In particular, if we take any measurable
set $A$ and set $f=\mathds{1}_{A}$ in (5.12), we get

$\displaystyle\rho(A)=\int R(x,A)\rho(dx).$ (5.13)

##  Appendix A Target Distribution in the Soft Constraints Limit

We will give a semi-rigorous argument for (1.4), (1.5) that uses the co-area
formula. Take any function $\varphi\in C_{c}^{\infty}(\mathbb{R}^{d_{a}})$.
Use $S_{a}$ to denote the set
$S_{a}=\big{\\{}x\in\mathbb{R}^{d_{a}}:q(x)=a\big{\\}}$ and $\sigma_{a}(dx)$
to denote the $d$-dimensional Hausdorff measure on $S_{a}$. Let $R(x)=\nabla
q(x)^{T}\nabla q(x)$, and define the functions
$g(a)=\int_{S_{a}}\varphi(x)|R(x)|^{-1/2}\sigma_{a}(dx)$ and
$h(a)=\int_{S_{a}}|R(x)|^{-1/2}\sigma_{a}(dx)$. First, we observe that

$\displaystyle\begin{split}Z_{\varepsilon}=\int_{\mathbb{R}^{d_{a}}}&e^{\frac{U(x)}{2\varepsilon^{2}}}dx=\int_{\mathbb{R}^{m}}\int_{S_{a}}e^{\frac{U(x)}{2\varepsilon^{2}}}|R(x)|^{-1/2}\sigma_{a}(dx)da\\\
\\\ &=\int_{\mathbb{R}^{m}}e^{-\frac{|a|^{2}}{2\varepsilon^{2}}}h(a)da\ \sim\
h(0)(2\pi)^{m/2}\varepsilon^{m},\quad\text{as}\ \varepsilon\downarrow
0,\end{split}$ (A.1)

where the second equality is just the co-area formula, in the third equality
we used that $U(x)=|a|^{2}$ when $x\in S_{a}$, and the last step follows from
Laplace’s method. Then,

$\displaystyle\begin{split}\int_{\mathbb{R}^{d_{a}}}\varphi(x)\pi_{\varepsilon}(dx)&=\int_{\mathbb{R}^{m}}\int_{S_{a}}\varphi(x)\frac{1}{Z_{\varepsilon}}e^{-\frac{U(x)}{2\varepsilon^{2}}}|R(x)|^{-1/2}\sigma_{a}(dx)da\\\
\\\
&=\frac{1}{Z_{\varepsilon}}\int_{\mathbb{R}^{m}}e^{-\frac{|a|^{2}}{2\varepsilon^{2}}}g(a)da\\\
\\\
&\approx\frac{1}{h(0)(2\pi)^{m/2}\varepsilon^{m}}\int_{\mathbb{R}^{m}}e^{-\frac{|a|^{2}}{2\varepsilon^{2}}}g(a)da\\\
\\\ &\xrightarrow[]{\varepsilon\downarrow 0}\ \frac{1}{h(0)}g(0)\\\ \\\
&=\frac{1}{h(0)}\int_{S}\varphi(x)|R(x)|^{-1/2}\sigma(dx)\\\ \\\
&=\int_{S}\varphi(x)\pi_{s}(dx),\end{split}$ (A.2)

where $\pi_{s}(x)=|R(x)|^{-1/2}\sigma(dx)/h(0)$. The first equality in (A.2)
is (once again) the co-area formula. In the third line of (A.2) we used the
asymptotic approximation (A.1). The limit in the fourth line holds by standard
approximate identity arguments.

## References

  * [1] Miranda Holmes-Cerfon  “Simulating sticky particles: A Monte Carlo method to sample a stratification”  In _The Journal of Chemical Physics_ 153.16 AIP Publishing, 2020, pp. 164112  DOI: 10.1063/5.0019550
  * [2] Miranda Holmes-Cerfon, Louis Theran and Steven J. Gortler  “Almost-rigidity of frameworks”  arXiv, 2019  DOI: 10.48550/ARXIV.1908.03802
  * [3] Alan Sokal  “Monte Carlo methods in statistical mechanics: Foundations and new algorithms”  In _Functional integration: basics and applications_ Springer, Boston, 1997, pp. 131–192
  * [4] Emilio Zappa, Miranda Holmes-Cerfon and Jonathan Goodman  “Monte Carlo on manifolds: sampling densities and integrating functions”  arXiv, 2017  DOI: 10.48550/ARXIV.1702.08446