Monte Carlo on manifolds: sampling densities and integrating functions

TL;DR

This paper introduces Monte Carlo methods for sampling and integrating over manifolds defined by constraints, using tangent space projections to efficiently handle probability distributions and compute integrals without curvature information.

Contribution

It presents a tangent-space-based MCMC sampler and a multi-stage integral computation algorithm for constrained manifolds, with practical error estimates and applications to physical systems.

Findings

Effective sampling and integration algorithms demonstrated on test problems.

Error estimates enable single-run accuracy assessment.

Application to sticky hard sphere systems predicts phase transition points.

Abstract

We describe and analyze some Monte Carlo methods for manifolds in Euclidean space defined by equality and inequality constraints. First, we give an MCMC sampler for probability distributions defined by un-normalized densities on such manifolds. The sampler uses a specific orthogonal projection to the surface that requires only information about the tangent space to the manifold, obtainable from first derivatives of the constraint functions, hence avoiding the need for curvature information or second derivatives. Second, we use the sampler to develop a multi-stage algorithm to compute integrals over such manifolds. We provide single-run error estimates that avoid the need for multiple independent runs. Computational experiments on various test problems show that the algorithms and error estimates work in practice. The method is applied to compute the entropies of different sticky hard sphere systems. These predict the temperature or interaction energy at which loops of hard sticky spheres become preferable to chains.

Figures

Tables

Table 1. Table 1: Numerical vs analytical results for the volume of the special orthogonal group S ​ O ​ ( n ) 𝑆 𝑂 𝑛 SO(n) , for 2 ≤ n ≤ 7 2 𝑛 7 2\leq n\leq 7 . Here d 𝑑 d is the dimension of S ​ O ​ ( n ) 𝑆 𝑂 𝑛 SO(n) , V ¯ ¯ 𝑉 \overline{V} is the sample mean obtained from a sample of size 50 50 50 and the relative standard deviation s V ~ ~ subscript 𝑠 𝑉 \widetilde{s_{V}} and relative error ϵ r subscript italic-ϵ 𝑟 \epsilon_{r} are given in ( 47 ).

$n$	$d$	$V_t$	$\overline{V}$	$\widetilde{s_V}$		${ϵ}_r$
2	1	8.89	8.9	0.	0038	0.	0047
3	3	223.3	223	0.	0078	0.	0024
4	6	$1.24\cdot {10}^4$	$1.2\cdot {10}^4$	0.	016	0.	012
5	10	$1.31\cdot {10}^6$	$1.3\cdot {10}^6$	0.	022	0.	024
6	15	$2.30\cdot {10}^8$	$2.2\cdot {10}^8$	0.	0354	0.	0519
7	21	$6.09\cdot {10}^{10}$	$6.0\cdot {10}^{10}$	0.	0547	0.	0094

Table 2. Table 2: Geometrical partition functions z 𝑧 z of the configuration spaces of clusters of N = 6 𝑁 6 N=6 unit spheres with m = 10 𝑚 10 m=10 contacts (the numbering is as in [ 1 ] , where the manifolds are called “modes.”) Each computation is done using a single MCMC run of n t = 10 8 subscript 𝑛 𝑡 superscript 10 8 n_{t}=10^{8} points. The error bars (such as ± 0.12 plus-or-minus 0.12 \pm 0.12 for manifold 8) are one standard deviation, estimated using ( 43 ). We compare the calculations with the value z t subscript 𝑧 𝑡 z_{t} obtained in [ 1 ] .

Manifold	$z$				$z_t$
$8$	6.	99	$\pm$	0.12	7.	014
$9$	44.	4	$\pm$	0.33	44.	402
$10$	57.	1	$\pm$	0.26	57.	07
$11$	18.	9	$\pm$	0.20	19.	18
$12$	3.	30	$\pm$	0.020	3.	299
$13$	62.	4	$\pm$	0.23	62.	34
$14$	3.	35	$\pm$	0.039	3.	299
$15$	27.	7	$\pm$	0.14	27.	78
$16$	11.	7	$\pm$	0.058	11.	67
$17$	28.	5	$\pm$	0.31	28.	41
$18$	8.	44	$\pm$	0.088	8.	238
$19$	69.	2	$\pm$	0.33	69.	43
$20$	50.	9	$\pm$	0.27	50.	86

Table 3. Table 3: Volume V 𝑉 V , geometrical partition function z 𝑧 z and mean vibrational contribution h ¯ = z / V ¯ ℎ 𝑧 𝑉 \overline{h}=z/V for the manifolds M i , N subscript 𝑀 𝑖 𝑁 M_{i,N} corresponding to the loop/chains of N 𝑁 N particles. The computations were done using n t = 10 8 subscript 𝑛 𝑡 superscript 10 8 n_{t}=10^{8} points. Error bars are calculated using ( 43 ). The results are compared with the average vibrational contribution h ~ ~ ℎ \widetilde{h} estimated from our sampling algorithm. Specifically, we generate n p subscript 𝑛 𝑝 n_{p} points on M i , N subscript 𝑀 𝑖 𝑁 M_{i,N} (and repeat this n m = 10 subscript 𝑛 𝑚 10 n_{m}=10 times) and estimate the average as h ~ = 1 n m ​ ∑ j = 1 n m h ~ j ~ ℎ 1 subscript 𝑛 𝑚 superscript subscript 𝑗 1 subscript 𝑛 𝑚 subscript ~ ℎ 𝑗 \widetilde{h}=\frac{1}{n_{m}}\sum_{j=1}^{n_{m}}\widetilde{h}_{j} where h ~ j = 1 n p ​ ∑ i = 1 n p f ​ ( X i ) subscript ~ ℎ 𝑗 1 subscript 𝑛 𝑝 superscript subscript 𝑖 1 subscript 𝑛 𝑝 𝑓 subscript 𝑋 𝑖 \widetilde{h}_{j}=\frac{1}{n_{p}}\sum_{i=1}^{n_{p}}f(X_{i}) . The error bars of the latter are computed using the sample standard deviation, using a sample of size n m = 10 subscript 𝑛 𝑚 10 n_{m}=10 .

	$V$			$z$			$\overline{h}$					$\tilde{h}$
4 Chain	$3.51\cdot {10}^2$	$\pm$	1.03	$1.32\cdot {10}^2$	$\pm$	0.38	0.	374	$\pm$	0.	0022	0.	379	$\pm$	$1.5\cdot {10}^{-5}$
4 Loop	$8.27\cdot {10}^1$	$\pm$	0.20	$2.18\cdot {10}^1$	$\pm$	0.054	0.	264	$\pm$	0.	0013	0.	261	$\pm$	$7.6\cdot {10}^{-5}$
5 Chain	$3.03\cdot {10}^3$	$\pm$	16	$8.31\cdot {10}^2$	$\pm$	4.3	0.	274	$\pm$	0.	0029	0.	279	$\pm$	$1.1\cdot {10}^{-4}$
5 Loop	$3.19\cdot {10}^2$	$\pm$	4.3	$6.19\cdot {10}^1$	$\pm$	0.85	0.	194	$\pm$	0.	0053	0.	194	$\pm$	$2.8\cdot {10}^{-4}$
6 Chain	$2.70\cdot {10}^4$	$\pm$	$2.7\cdot {10}^2$	$5.33\cdot {10}^3$	$\pm$	52	0.	197	$\pm$	0.	0039	0.	206	$\pm$	$1.9\cdot {10}^{-4}$
6 Loop	$1.66\cdot {10}^3$	$\pm$	13	$2.32\cdot {10}^2$	$\pm$	1.9	0.	143	$\pm$	0.	0024	0.	147	$\pm$	$3.1\cdot {10}^{-5}$
7 Chain	$2.42\cdot {10}^5$	$\pm$	$1.3\cdot {10}^4$	$3.84\cdot {10}^4$	$\pm$	$1.9\cdot {10}^3$	0.	158	$\pm$	0.	0017	0.	151	$\pm$	$1.7\cdot {10}^{-4}$
7 Loop	$1.01\cdot {10}^4$	$\pm$	$1.2\cdot {10}^2$	$1.09\cdot {10}^3$	$\pm$	$13$	0.	109	$\pm$	0.	0026	0.	110	$\pm$	$2.6\cdot {10}^{-5}$
8 Chain	$2.19\cdot {10}^6$	$\pm$	$6.4\cdot {10}^4$	$2.77\cdot {10}^5$	$\pm$	$8.3\cdot {10}^3$	0.	127	$\pm$	0.	0075	0.	111	$\pm$	$9.6\cdot {10}^{-4}$
8 Loop	$6.76\cdot {10}^4$	$\pm$	$1.6\cdot {10}^3$	$5.65\cdot {10}^3$	$\pm$	$1.3\cdot {10}^2$	0.	0836	$\pm$	0.	0039	0.	0824	$\pm$	$2.5\cdot {10}^{-5}$
9 Chain	$2.12\cdot {10}^7$	$\pm$	$1.1\cdot {10}^6$	$1.73\cdot {10}^6$	$\pm$	$8.9\cdot {10}^4$	0.	0816	$\pm$	0.	0084	0.	0820	$\pm$	$1.1\cdot {10}^{-5}$
9 Loop	$4.94\cdot {10}^5$	$\pm$	$1.7\cdot {10}^4$	$2.87\cdot {10}^4$	$\pm$	$9.2\cdot {10}^2$	0.	0582	$\pm$	0.	0038	0.	0612	$\pm$	$1.2\cdot {10}^{-5}$
10 Chain	$2.07\cdot {10}^8$	$\pm$	$1.2\cdot {10}^7$	$1.43\cdot {10}^7$	$\pm$	$8.1\cdot {10}^5$	0.	0689	$\pm$	0.	0080	0.	0604	$\pm$	$3.5\cdot {10}^{-4}$
10 Loop	$3.81\cdot {10}^6$	$\pm$	$2.2\cdot {10}^5$	$1.71\cdot {10}^5$	$\pm$	$9.7\cdot {10}^3$	0.	0448	$\pm$	0.	0051	0.	0453	$\pm$	$4.7\cdot {10}^{-5}$

Table 4. Table 4: Ratios between the geometrical partition functions of a chain and a loop of N 𝑁 N particles, for 4 ≤ N ≤ 10 4 𝑁 10 4\leq N\leq 10 . We report both the ratio z C / z L subscript 𝑧 𝐶 subscript 𝑧 𝐿 z_{C}/z_{L} of a single chain vs. loop, and the ratio ( n C ​ z C ) / ( n L ​ z L ) subscript 𝑛 𝐶 subscript 𝑧 𝐶 subscript 𝑛 𝐿 subscript 𝑧 𝐿 (n_{C}z_{C})/(n_{L}z_{L}) when spheres are treated as indistinguishable. The latter may be used to estimate the sticky parameter κ 𝜅 \kappa as in ( 54 ).

$N$	$z_C/z_L$				$\left(n_L{z}_L\right)/\left(n_C{z}_C\right)$
$4$	6.	02	$\pm$	0.032	24.	1	$\pm$	0.	13
$5$	13.	5	$\pm$	0.13	67.	5	$\pm$	0.	65
$6$	22.	9	$\pm$	0.41	138		$\pm$	2.	5
$7$	35.	1	$\pm$	1.0	246		$\pm$	7.	0
$8$	49.	0	$\pm$	2.6	392		$\pm$	21
$9$	60.	3	$\pm$	5.0	543		$\pm$	45
$10$	83.	6	$\pm$	9.5	836		$\pm$	95

Equations

$$Z={\int }_{M}f(x)\sigma (dx),$$

$$Z_B={\int }_{M\cap B}f(x)\sigma (dx).$$

$${\rho }_B(dx)=\frac{1}{Z_B}f(x)\sigma (dx).$$

$$R_i=\frac{Z_i}{Z_{i+1}}.$$

$$\hat{Z}={\hat{Z}}_k\prod _{i=0}^{k-1}{\hat{R}}_i,$$

$${\sigma }_{\hat{Z}}^2=\text{mbox}var\left(\hat{Z}\right),$$

$$M=\left\{x\in {\mathbb{R}}^{d_a}:q_i(x)=0,i=1,\dots ,m,h_j(x)>0,j=1,\dots ,l\right\}.$$

$${\left(Q_x\right)}_{ij}=\frac{\partial q_i(x)}{\partial x_j}.$$

$$\rho (dx)=\frac{1}{Z}f(x)\sigma (dx).$$

$$\text{mbox}if{X}_k\sim \rho \text{mbox}then{X}_{k+1}\sim \rho .$$

$$p(v|x)=\frac{1}{(2\pi {)}^{d/2}{s}^d}{e}^{-\frac{|v{|}^2}{2s^2}}.$$

$$w=\sum _{j=1}^m{a}_j\nabla q_j(x)=Q_{x}a.$$

$$q_i\left(x+v+Q_{x}a\right)=0,i=1,\dots m.$$

$$J_{ij}={\partial }_{a_j}{q}_i(x+v+Q_{x}a)={\left(\nabla q_i(x+v+Q_{x}a)\right)}^t\nabla q_j(x).$$

$$|q(x+v+Q_{x}a)|\le \text{mbox}tol.$$

$$\underset{y\in M}{\mathrm{argmin}}{\left|x+v-y\right|}_2.$$

$$a(y|x)=\min \left(1,\frac{f(y)p(v^{\prime }|y)}{f(x)p(v|x)}\right).$$

$$a(y|x)=\min \left(1,\frac{f(y)p(v^{\prime }|y)J(x|y)}{f(x)p(v|x)J(y|x)}\right),$$

$$J(y|x)=J(x|y).$$

$$\sum _j{u}_{x,i}^t{u}_{y,j}{b}_j=a_i.$$

$$J(y|x{)}^{-1}=\text{mbox}det\left(P^{-1}\right)=\text{mbox}det\left((U_x^t{U}_y{)}^{-1}\right).$$

$$J(x|y{)}^{-1}=\text{mbox}det\left((U_y^t{U}_x{)}^{-1}\right).$$

$$P_x(dy)=h(y|x)\sigma (dy)+\zeta (x){\delta }_x(dy).$$

$$h(y|x)=(1-1_{F_x}(y))\cdot p(v|x)\cdot J(y|x),$$

$$R_x(dy)=a(y|x)h(y|x)\sigma (dy)+\xi (x){\delta }_x(dy).$$

$$\xi (x)=\zeta (x)+{\int }_M\left(1-a(y|x)\right)h(y|x)\sigma (dy).$$

$$f(x)a(y|x)h(y|x)=f(y)a(x|y)h(x|y).$$

$$X_k\sim {\rho }_k(dx)=g_k(x)\sigma (dx).$$

$$g_{k+1}(x)={\int }_M{g}_k(y)a(x|y)h(x|y)d\sigma (y)+\xi (x)g_k(x).$$

$$g_{k+1}(x)={\int }_{M}f(x)a(y|x)h(y|x)d\sigma (y)+\xi (x)f(x).$$

Full text

1 Introduction

Many physical and engineering systems involve motion with equality constraints. For example, in physics and chemistry, bonds between atoms or colloidal particles may be modeled as having a fixed length [1], while in robotics, joints or hinges connecting moving components may be written as constraints on distances or angles [2]. Constraints also arise in statistics on the set of parameters in parametric models [3]. In such problems the space of accessible configurations is lower-dimensional than the space of variables which describe the system, often forming a manifold embedded in the full configuration space. One may then be interested in sampling a probability distribution defined on the manifold, or calculating an integral over the manifold such as its volume.

This paper presents Monte Carlo methods for computing a $d$-dimensional integral over a $d$-dimensional connected manifold $M$ embedded in a $d_{a}$-dimensional Euclidean space, with $d_{a}\geq d$. The manifold is defined by sets of constraints. We first describe an MCMC sampler for manifolds defined by constraints in this way. Then we describe a multi-phase procedure that uses the sampler to estimate integrals. It should go without saying (but rarely does) that deterministic algorithms to compute integrals are impractical except in very low dimensions or for special problems (e.g.[4]).

Our MCMC sampler has much in common with other algorithms for constraint manifolds, including other samplers, optimizers, differential algebraic equation (DAE) solvers (see, e.g., [5]), etc. A move from a point $x\in M$ starts with a move in the tangent space at $x$, which is $v\in T_{x}$. This is followed by a projection back to $y\in M$, which we write as $y=x+v+w$. The sampler simplifies if we require $w\perp T_{x}$, so the projection is perpendicular to the tangent space at the original point. This idea was first suggested in [6]111We learned of this work after we submitted the original version of this paper for publication.. This is different from other natural choices, such as choosing $y\in M$ to minimize the projection distance $\left\|y-(x+v)\right\|$. The choice of $w\perp T_{x}$ makes the sampler easy to implement as it requires only first derivatives of the constraint surfaces. Neither solving for $y$ using Newton’s method nor the Metropolis Hastings detailed balance condition require any higher derivatives. Other surface sampling methods that we know of require second derivatives, such as the explicit parameterization method of [3], the Hamiltonian method such as of [7], geodesic methods such as [8, 9], and discretized SDE methods such as [10]. It is common in our computational experiments that there is no $w\perp T_{x}$ so that $x+v+w\in M$, or that one or more $w$ values exist but the Newton solver fails to find one. Section 2 explains how we preserve detailed balance even in these situations.

Our goal is to compute integrals of the form

[TABLE]

where $\sigma(dx)$ is the $d-$dimensional surface area measure (Hausdorff measure) and $f$ is a positive smooth function. Writing $B\subset\mathbb{R}^{d_{a}}$ for the ball centered at a point $x_{0}\in M$, we also estimate certain integrals of the form

[TABLE]

For each $B$, we define a probability distribution on $M\cap B$:

[TABLE]

The sampling algorithm outlined above, and described in more detail in Section 2, allows us to draw samples from $\rho_{B}$.

To calculate integrals as in (1) we use a multi-phase strategy related to nested sampling and thermodynamic integration (e.g.[11]). Similar strategies were applied to volume estimation problems by the Lovasz school, see, e.g., [12]. Let $B_{0}\supset B_{1}\supset\cdots\supset B_{k}$ be a nested contracting family of balls with a common center. With an abuse of notation, we write $Z_{i}$ for $Z_{B_{i}}$. We use MCMC sampling in $M\cap B_{i}$ to estimate the ratio

[TABLE]

We choose the last ball $B_{k}$ small enough that we may estimate $Z_{k}$ by direct Monte-Carlo integration. We choose $B_{0}$ large enough that $M\subset B_{0}$ so $Z=Z_{0}$. Our estimate of $Z$ is

[TABLE]

where hatted quantities such as $\widehat{Z}$ are Monte Carlo estimates. Section 3 describes the procedure in more detail.

Error estimates should be a part of any Monte Carlo computation. Section 3.2 describes our procedure for estimating

[TABLE]

a procedure adapted from [13]. If we can estimate $\sigma^{2}_{\widehat{R}_{i}}$, and if these are small enough, then we can combine them to estimate $\sigma^{2}_{\widehat{Z}}$. We explore how $\sigma_{\widehat{Z}}^{2}$ depends on parameters of the algorithm via several simplified model problems, in section 5. We specifically analyze how the error depends on the amount by which we shrink the balls for each ratio, and make predictions that could lead to good systematic choices of this parameter.

We illustrate our methods on a number of examples. We apply the sampler to the surface of a torus and cone in 3D, as well as the special orthogonal group $SO(n)$ for $n=11$, which is a manifold of dimension $55$ in an ambient space of dimension 121 (section 2.3) . All of these examples have certain marginal distributions which are known analytically, so we verify the sampler is working correctly. We verify the integration algorithm is correct by calculating the surface areas of a torus and $SO(n)$ for $2\leq n\leq 7$ (sections 4.1, 4.2.)

Finally, we apply our algorithm to study clusters of sticky spheres (section 4.3), a model that is studied partly to understand how colloidal particles self-assemble into various configurations [14, 1, 15]. A simple question in that direction concerns a linear chain of $n$ spheres – whether it is more likely to be closed in a loop or open in a chain. According to equilibrium statistical mechanics, the ratio of the probabilities of being open versus closed depends on a parameter ( the “sticky parameter”) characterizing the binding energy and temperature, and the entropies, which are surface areas of the kind we are computing. We compute the entropies for open and closed chains of lengths 4–10. This data can be used to predict the relative probabilities of being in a loop or a chain, if the sticky parameter is known, but if it is not known our calculations can be used to measure it by comparing with data.

Notation and assumptions.

Throughout the paper, we will denote by $M$ a $d$-dimensional connected manifold embedded in an ambient space $\mathbb{R}^{d_{a}}$, which is defined by equality and inequality constraints. There are $m$ equality constraints, $q_{i}(x)=0$, $i=1,\ldots,m$, where the $q_{i}$ are smooth functions $q_{i}:\mathbb{R}^{d_{a}}\rightarrow\mathbb{R}$. There are $l$ inequality constraints of the form $h_{j}(x)>0$, for $j=1,\ldots,l$. The manifold $M$ is

[TABLE]

If $M$ has several connected components, then the sampler may sample only one component, or it may hop between several components. It is hard to distinguish these possibilities computationally.

We let $Q_{x}$ be the matrix whose columns are the gradients $\{\nabla q_{i}(x)\}_{i=1}^{m}$. This makes $Q_{x}$ the transpose of the Jacobian of the overall constraint function $q:\mathbb{R}^{d_{a}}\to R^{m}$. The entries of $Q_{x}$ are

[TABLE]

We assume that $Q_{x}$ has full rank $m$ everywhere on $M$. By the implicit function theorem, this implies that the dimension of $M$ is $d=d_{a}-m$ and that the tangent space $T_{x}\equiv T_{x}M$ at a point $x\in M$ is well defined. In this case $\{\nabla q_{i}(x)\}_{i=1}^{m}$ form a basis of the orthogonal space $T^{\perp}_{x}\equiv T_{x}M^{\perp}$. Note that $M$ inherits the metric from the ambient space $\mathbb{R}^{d_{a}}$ by restriction. The corresponding $d-$dimensional volume element is $d-$dimensional Hausdorff measure, which we denote by $\sigma(dx)$.

2 Sampling on a constraint manifold

2.1 The algorithm

We describe the MCMC algorithm we use to sample a probability distribution of the form

[TABLE]

Here $\sigma$ is the $d-$dimensional surface measure on $M$ given by (6). As usual, the algorithm requires repeated evaluation of $f$ but does not use the unknown normalization constant $Z$. The algorithm is described in pseudocode in Figure 1.

Our Metropolis style MCMC algorithm generates a sequence $X_{k}\in M$ with the property that

[TABLE]

Once $X_{k}=x$ is known, the algorithm makes a proposal $y\in M$. If the proposal fails (see below) or if it is rejected, then $X_{k+1}=X_{k}$. If the algorithm succeeds in generating a $y$, and if $y$ is accepted, then $X_{k+1}=y$. If $M$ is connected, compact, and smooth, then the algorithm is geometrically ergodic. See e.g. [16] for an efficient summary of the relevant MCMC theory.

The proposal process begins with a tangential move $x\to x+v$ with $v\in T_{x}$. We generate $v$ by sampling a proposal density $p(v|x)$ defined on $T_{x}$. All the numerical experiments reported here use an isotropic $d-$dimensional Gaussian with width $s$ centered at $x$:

[TABLE]

We generate $v$ using an orthonormal basis for $T_{x}$, which is the orthogonal complement of the columns of the constraint gradient matrix $Q_{x}$, see (7). This orthonormal basis is found as the last $d$ columns of the $d_{a}\times d_{a}$ matrix $Q$ in the QR decomposition222Please be aware of the conflict of notation here. The $d_{a}\times d_{a}$ orthogonal $Q$ of the QR decomposition is not the $d_{a}\times d$ gradient $Q_{x}$. of $Q_{x}$, which we find using dense linear algebra.

Given $x$ and $v$, the projection step looks for $w\in T_{x}^{\perp}$ with $y=x+v+w\in M$. See Figure 2(a) for an illustration. It does this using an $m-$component column vector $a=(a_{1},\ldots,a_{m})^{t}$ and

[TABLE]

The unknown coefficients $a_{j}$ are found by solving the nonlinear equations

[TABLE]

This can be done using any nonlinear equation solver. In our code, we solve (12) using simple Newton’s method (no line search, regularization, etc.) and initial guess $a=0$. The necessary Jacobian entries are

[TABLE]

We iterate until the first time a convergence criterion is satisfied

[TABLE]

where $|q|$ is some norm. In our implementation, we used the $l_{2}$ norm. When this criterion is satisfied, the projection is considered a success and $y=x+v+w$ is the proposal. If the number of iterations reaches nmax before the convergence criterion is satisfied, then the projection phase is considered a failure. It is possible to change these details, but if you do you must also make the corresponding changes in the detailed balance check below.

There are other possible projections from $x+v$ to $M$, each with advantages and disadvantages. It would seem natural, for example, to find $y$ by computing

[TABLE]

This produces a projection $w\in T_{y}^{\perp}$, rather than our $w\in T_{x}^{\perp}$. It also has the advantage that such a $y$ always exists and is almost surely unique. The optimization algorithm that computes (15) may be more robust than our Newton-based nonlinear equation solver. The disadvantage of (15) is that there are curvature effects in the relative surface area computation that are part of detailed balance. Computing these curvatures requires second derivatives of the constraint functions $q_{i}$. This is not be an intrinsic limitation, but requires more computation and is not as straightforward to implement as our algorithm.

The proposal algorithm just described produces a random point $y\in M$. There are two cases where we must reject it immediately so the detailed balance relation may hold (which we verify in section 2.2.) First, we must check whether an inequality constraint is violated, so $h_{i}(y)\leq 0$ for some $i$. If so, then $y$ is rejected. Second, we must check whether the reverse step is possible, i.e. whether it is possible to propose $x$ starting from $y$. To make the reverse proposal, the algorithm would have to choose $v^{\prime}\in T_{y}$ so that $x=y+v^{\prime}+w^{\prime}$ with $w^{\prime}\perp T_{y}$. Figure 2(b) illustrates this. Since $x$ and $y$ are known, we may find $v^{\prime}$ and $w^{\prime}$ by the requirement that $x-y=v^{\prime}+w^{\prime}$, $v^{\prime}\in T_{y}$, $w^{\prime}\in T_{y}^{\perp}$. This is always possible by projecting $x-y$ onto $T_{y}$ and $T_{y}^{\perp}$, which are found using the $QR$ decomposition of $Q_{y}$. However, we must additionally determine whether the Newton solver would find $x$ starting from $y+v^{\prime}$, a step that was neglected in earlier studies, e.g. [6]. Even though we know that $x$ is a solution to the equations, we must still verify that the Newton solver produces $x$. With that aim, we run the Newton algorithm from initial guess $y+v^{\prime}$ and see whether the the converge criterion (14) within nmax iterations. If not the proposal $y$ is rejected. Whenever $y$ is rejected we set $X_{k+1}=x$.

If we come to this point, we compute an acceptance probability $a(y|x)$ using the Metropolis Hastings formula

[TABLE]

The reader might expect that we should use the more complicated formula

[TABLE]

where $J(y|x)$ is the inverse of the determinant of the linearization of the projection $v\rightarrow y$, formally, $\left|\frac{\partial v}{\partial y}\right|$. This would account for how the area of a small patch of surface near $x$, is distorted upon mapping it to $y$. However, it turns out that the Jacobian factors cancel:

[TABLE]

This is evident in part b of Figure 2, where the angle between $w$ and $T_{y}$ is the same as the angle between $w^{\prime}$ and $T_{x}$. The proof in two dimensions is an exercise in Euclidean geometry. In more dimensions, $J(y|x)$ and $J(x|y)$ are products of cosines of principle angles between $T_{x}$ and $T_{y}$, see [17].

To prove (17), let $U_{x}$ and $U_{y}$ be matrices whose columns are orthonormal bases for $T_{x}$ and $T_{y}$ respectively. Column $i$ of $U_{x}$ will be denoted $u_{x,i}$, and similarly for $U_{y}$. Consider a vector $\xi\in T_{x}$, which may be written $\xi=\sum a_{i}u_{x,i}$. Let $\eta\in T_{y}$ be the projection of $\xi$ to $T_{y}$ normal to $T_{x}$, and write $\eta=\sum b_{j}u_{y,j}$. The orthogonality condition is $u_{x,i}^{t}\eta=u_{x,i}^{t}\xi=a_{i}$. This leads to the equations

[TABLE]

These take the form $Pb=a$, or $b=P^{-1}a$, where $P=U_{x}^{t}U_{y}$. Since $U_{x}$ and $U_{y}$ are orthonormal bases, the volume element going from $x\to y$ is expanded by a factor

[TABLE]

Similar reasoning shows that the projection from $T_{y}$ to $T_{x}$ perpendicular to $T_{y}$ expands volumes by a factor of

[TABLE]

These factors are equal because the determinants of a matrix and its transpose are the same, proving (17).

We close this section with a general remark concerning the choice of parameters. A key part of the algorithm is the projection step, which requires using Newton’s method or some other algorithm to solve a nonlinear system of equations. This algorithm called twice at every step of the sampler and so making it efficient is critical to the overall efficiency of the sampler. Profiling reveals that most of the computing time is spent in the ** project** part (see Figure 1). We found that setting the parameter nmax low, or otherwise terminating the solver rapidly, made a significant difference to the overall efficiency. Unlike most problems that require solving equations, we are not interested in guaranteeing that we can find a solution if one exists – we only need to find certain solutions rapidly; our reverse projection check and subsequent rejection ensures that we correct for any asymmetries in the solver and still satisfy detailed balance.

2.2 Verifying the balance relation

In this section we verify what may already be clear, that the algorithm described above satisfies the desired balance relation (9). In the calculations below, we assume that the constraints hold exactly, an assumption that will not be true when they are solved for numerically, as we discuss in the remark at the end of the section.

To start, let us express the probability distribution of the proposal point $y$ as the product of a surface density and $d-$dimensional surface measure,

[TABLE]

Here $\delta_{x}(dy)$ is the unit point mass probability distribution at the point $x$, $\zeta(x)$ is the probability that the projection step failed, i.e. Newton’s method failed to produce a $y$ or it lay outside the boundary, and $h(y|x)$ is the proposal density, i.e. the probability density of successfully generating $y\in M$ from $x$. The proposal density is the product of three factors:

[TABLE]

where the set $F_{x}\subset M$ consists of those points $y\in M$ that cannot be reached using our Newton solver starting from $x$, and ${\bf 1}_{F_{x}}(y)$ denotes the characteristic function of the set $F_{x}$, which equals 1 if $y\in F_{x}$ and zero if $y\notin F_{x}$. Note that we must have $\int_{M}h(y|x)\,d\sigma(y)+\zeta(x)=1$.

If the acceptance probability is $a(y|x)$, then the overall probability distribution of $y$ takes the form

[TABLE]

Here, $\xi(x)$ is the overall probability of rejecting a move when starting at $x$, equal to the sum of the probability of not generating a proposal move, and the probability of successfully generating a proposal move that is subsequently rejected. Since $1-a(y|x)$ is the probability of rejecting a proposal at $y$, the overall rejection probability is

[TABLE]

We quickly verify that the necessary balance relation (9) is satisfied if $a(y|x)$ is chosen to satisfy the detailed balance formula

[TABLE]

For this, suppose $X_{k}$ has probability distribution on $M$ with density $g_{k}(x)$ with respect to surface area measure:

[TABLE]

We check that if (22) is satisfied, and $g_{k}=f$, then $g_{k+1}=f$. We do this in the usual way, which is to write the integral expression for $g_{k+1}$ and then simplify the result using (22). To start,

[TABLE]

The integral on the right represents jumps to $x$ from $y\neq x$ on $M$. The second term on the right represents proposals from $x$ that were unsuccessful or rejected. If $g_{k}=f$, using the detailed balance relation (22) this becomes

[TABLE]

The integral on the right now represents the probability of accepting a proposal from $x$:

[TABLE]

When combined with the second term, the right side of (23) is equal to $f(x)$, as we wished to show.

We use the standard Metropolis Hastings formula to enforce (22):

[TABLE]

If we succeed in proposing $y$ from $x$, then we know $y\notin F_{x}$ so we are able to evaluate $h(y|x)$. However, it is still possible that $x\notin F_{y}$, in which case $h(x|y)=0$. This is why we must apply the Newton solver to the reverse move, to see whether it would succeed in proposing $x$ from $y$. If not, we know $h(x|y)=0$ and therefore $a(y|x)=0$ and we must reject $y$.

Remark. We have assumed so far the constraints hold exactly, which is not true when we solve for $y$ numerically. Rather, the proposal move will lie within some tolerance tol, so we are sampling a “fattened” region $M_{\texttt{tol}}=\left\{y:\;|q_{i}(y)|\leq\texttt{tol}\;i=1,\ldots,m,\;h_{j}(x)>0,\;j=1,\ldots,l\right\}$. It is well-known (e.g. [18, 19]) that if the fattened region is uniformly sampled, the distribution near $M$ will not be the uniform distribution, but rather will differ by an $O(1)$ amount which does not vanish as $\texttt{tol}\to 0$.

We do not know the distribution by which our numerical solver produces a point in $M_{\texttt{tol}}$. However, since we observe the correct distributions on $M$ in the examples below (for small enough tol), we infer that we are not sampling the fattened region uniformly, but rather in a way that is consistent with the surface measure on the manifold, in the limit when $\texttt{tol}\to 0$. Moreover, it is practical to take tol to be extremely small, on the order of machine precision, because of the fast local convergence of Newton’s method.

2.3 Examples

In this section we illustrate the MCMC sampler with three examples: a torus, a cone and the special orthogonal group $SO(n)$.

2.3.1 Torus

Consider a torus $\mathbb{T}^{2}$ embedded in $\mathbb{R}^{3}$, implicitly defined by

[TABLE]

where $R$ and $r$ are real positive numbers, with $R>r$. Geometrically, $\mathbb{T}^{2}$ is the set of points at distance $r$ from the circle of radius $R$ in the $(x,y)$ plane centered at $(0,0,0)$. An explicit parameterization of $\mathbb{T}^{2}$ is given by

[TABLE]

We ran $N=10^{6}$ MCMC steps to sample uniform measure ($f(x)=1$) on $\mathbb{T}^{2}$ with toroidal radius $R=1$, poloidal radius $r=.5$, and step size scale $s=.5$. Figure 3 shows our empirical marginal distributions and the exact theoretical distributions. The empirical distributions of the toroidal angle $\theta$ and the poloidal angle $\phi$ are correct to within statistical sampling error. Around $6\%$ of proposed moves were rejected because of failure in the reverse projection step. When we didn’t implement this step, the marginals were not correct.

2.3.2 Cone

We consider the circular right-angle cone $\mathcal{C}\subset\mathbb{R}^{3}$ with vertex $(0,0,0)$ given by

[TABLE]

We study this example in order to determine how our algorithm might behave near singular points of a hypersurface. We have written strict inequalities to be consistent with the general formalism (6). In principle, the algorithm should do the same thing with non-strict inequalities, since the set of points with equality (e.g. $z=0$) has probability zero. If we allowed $z=0$, then $\mathcal{C}$ would not be a differentiable manifold because of the singularity at the vertex.

We suspected that the projection steps (see Figure 1) would be failure prone near the vertex because of high curvature, and that the auto-correlation times would be large or even unbounded. Nevertheless, Figure 4 shows satisfactory results, at least for this example. We ran for $N=10^{6}$ MCMC steps with step size parameter $s=.9$. The theoretical marginal densities for $X$, $Y$, and $Z$ are easily found to be

[TABLE]

2.3.3 Special orthogonal group

We apply our surface sampler to the manifold $SO(n)$, which is the the group of $n\times n$ orthogonal matrices with determinant equal to 1. We chose $SO(n)$ because these manifolds have high dimension, $\frac{1}{2}n(n-1)$, and high co-dimension, $\frac{1}{2}n(n+1)$ [20]. We view $SO(n)$ as the set of $n\times n$ matrices, $x\in\mathbb{R}^{n\times n}$ that satisfy the row ortho-normality constraints for $k=1,\ldots,n$ and $l>k$:

[TABLE]

This is a set of $\frac{1}{2}n(n+1)$ equality constraints. The gradient matrix $Q_{x}$ defined by (7) may be seen to have full rank for all $x\in SO(n)$. Any $x$ satisfying these constraints has $\mbox{det}(x)=\pm 1$. The set with $\mbox{det}(x)=1$ is connected. It is possible that our sampler would propose an $x$ with $\mbox{det}(x)=-1$, but we reject such proposals. We again chose density $f(x)=1$, which gives $SO(n)$ a measure that is, up to a constant, the same as the natural Haar measure on $SO(n)$.

One check of the correctness of the samples is the known large-$n$ distribution of $T=\mbox{Tr}(x)$, which converges to a standard normal as $n\to\infty$ [21]. This is surprising in view of the fact that $\mbox{Tr}(x)$ is the sum of $n$ random numbers each of which has $\mbox{var}(x_{kk})=1$. If the diagonal entries were independent with this variance, then $\mbox{var}(T)=n^{1/2}$ instead of the correct $\mbox{var}(T)=1$. This makes the correctness of the $T$ distribution an interesting test of the correctness of the sampler.

Figure 5 presents results for $n=11$ and $N=10^{6}$ MCMC steps, and proposal length scale $s=.28$. The distribution of $T$ seems correct but there are non-trivial auto-correlations. The acceptance probability was about $35\%$.

Note that there are better ways to sample $SO(n)$ that to use our sampler. For example, one can produce independent samples by choosing $y\in\mathbb{R}^{n\times n}$ with i.i.d. entries from ${\cal N}(0,1)$ and use the QR decomposition $y=Rx$, where $R$ is upper triangular and $x$ is orthogonal.

3 Integrating over a manifold

3.1 Algorithm

In this section we describe an algorithm to calculate integrals over a manifold $M$ of the form (1). We assume that $M$ given in (6) is bounded. The strategy as outlined in the Introduction is to consider a sequence of $d_{a}$-dimensional balls $B_{0}\supset B_{1}\supset\ldots\supset B_{k}$, centered at $x_{0}\in M$ with radii $r_{i}$, for $i=1,\ldots,k$, where $B_{0}\supset M$ entirely contains the manifold, and $B_{k}$ is “small enough” in a way we describe momentarily. We assume that each $M\cap B_{i}$ is connected, an assumption that is critical for our algorithm to work. We define the collection of integrals (see (2))

[TABLE]

The final estimate is $Z\approx\widehat{Z}$ according to (5).

To estimate the ratio $R_{i}=Z_{i}/Z_{i+1}$ we generate $n_{i}$ points in $B_{i}\cap M$ with the probability distribution $\rho_{i}$ as in (18). This is done using the MCMC algorithm given in Section 2 by including the extra inequality $|x-x_{0}|^{2}<r_{i}^{2}$. An estimator for $R_{i}$ is then given by

[TABLE]

where, for $i\leq j$, we define

[TABLE]

We could use additional elements $N_{i,j}$ with $j>i+1$ to enhance the estimate of the ratios, but we leave exactly how to do this a question for future research. Note that $N_{i,i}=n_{i}$.

The integral over the smallest set $M\cap B_{k}$ is computed in a different way. If $M$ had no equality constraints, so that the dimension of $M$ were equal to the dimension of $B_{k}$, then we could choose $B_{k}\subset M$ to lie entirely inside $M$ and so the integral (for constant $f$) would be known analytically. Indeed, this idea was used in several studies that developed efficient algorithms to calculate the volume of a convex body [4, 12]. With the equality constraints, in general there is no $B_{k}$ with radius $r_{k}>0$ for which the integral $Z_{k}$ is known in closed form exactly.

Instead, we choose $r_{k}$ small enough that that there is an easily computed one-to-one correspondence between points in $M\cap B_{k}$ and points in the tangent plane $D_{k}=T_{x_{0}}\cap B_{k}$. This happens, for example, when the manifold $M\cap B_{k}$ is a graph above $D_{k}$. We use the coordinates in the flat disk $D_{k}$ to integrate over $M\cap B_{k}$, which is simply the surface integral as illustrated in Figure 6:

[TABLE]

This uses the following notation. The integral on the right is over the $d$-dimensional disk. For any $x\in D_{k}$, the projection to $M$ perpendicular to $T_{x_{0}}$ is called $y(x)\in M$. For small enough $r_{k}$, the projection is easily calculated using the projection algorithm of Figure 1. The Jacobian factor $J(x)=\text{det}(U_{x_{0}}^{t}U_{y(x)})$ is known from (18), where as before $U_{x}$ is a matrix whose columns are an orthonormal basis of $T_{x}$.

The integral over $D_{k}$ may be computed deterministically if $d$ is small enough, but otherwise we estimate it by direct Monte Carlo integration. We choose $n_{k}$ independent points $x_{i}\in D_{k}$, which is easy for a $d-$dimensional disk. We compute the corresponding projections and Jacobian weight factors, $y_{i}=y(x_{i})$ and $J_{i}=J(x_{i})$. Some of the $y_{i}$ will be outside of $M\cap B_{k}$, so we need an extra indicator function factor, ${\bf 1}_{B_{k}}(y_{i})$, that is equal to 1 if $y_{i}\in M\cap B_{k}$. The $d$ dimensional volume of $D_{k}$ is known, so we have the direct Monte Carlo estimate

[TABLE]

The estimators (28) and (31) are substituted in (5) to give our overall estimate of the integral $Z$.

There is one additional consideration which is important to make this work in practice. That is, the projection from $D_{k}$ to $M\cap B_{k}$ must never fail. If it does, the estimate (31) will not be convergent. Empirically we found the projection always worked when $r_{k}$ was small enough, which proved to be a more limiting criterion than the criterion that $M\cap B_{k}$ be a graph above $D_{k}$.

3.2 Variance estimate

This section serves two purposes. One is to present the method we use to estimate the variance, $\sigma^{2}_{\widehat{Z}}$, of the estimator (5). This is similar to the strategy of [13], but adapted to this setting. This estimate is accurate only if the MCMC runs are long enough so that $Z_{k}$ and the $R_{i}$ are estimated reasonably accurately. Our computational tests below show that the estimator can be useful in practice. The other purpose is to set up notation for some heuristics for choosing parameters, a question dealt with in more depth in Section 5.

We start with the simple estimator $\widehat{Z}_{k}$ in (31) for the integral over the smallest ball. It is convenient to write the estimator as the sum of the exact quantity and a statistical error, and then to work with relative error rather than absolute error. We neglect bias, as we believe that error from bias is much smaller than statistical error in all our computations. In the present case, we write

[TABLE]

where $\zeta_{k}$ is a zero mean random variable with $\mbox{var}(\zeta_{k})=1$, $\sigma_{k}$ is the absolute standard deviation of $\widehat{Z}_{k}$, and $\rho_{k}=\sigma_{k}/Z_{k}$ is the relative standard deviation. Since the samples in (31) are independent, the standard deviation satisfies

[TABLE]

We estimate $\text{Var}(G(y))$ with the sample variance and obtain, for the absolute and relative variances,

[TABLE]

Here $\overline{G}=\sum_{i=1}^{n_{k}}\frac{1}{n_{k}}G(y_{i})$ is the sample mean.

Now consider the estimator $\widehat{R}_{i}$ in (28) for the $i$th volume ratio. This involves the number count $N_{i,i+1}$. We start by approximating $\mbox{var}(N_{i,i+1})$ and then incorporate this in the approximation for $\mbox{var}(\widehat{R}_{i})$. The number count may be written

[TABLE]

where $X_{j}\in M\cap B_{i}$ is our sequence of samples. In the invariant probability distribution of the chain, the probability of $X_{j}\in B_{i+1}$ is

[TABLE]

Since ${\bf 1}_{B_{i+1}}(X_{j})$ is a Bernoulli random variable, we have (in the invariant distribution)

[TABLE]

We will call this the static variance. It alone is not enough to determine $\mbox{var}(N_{i,i+1})$, since the steps of the Markov chain are correlated. To account for this correlation, we use the standard theory (see e.g. [22]) of error bars for MCMC estimators of the form (34), which we now briefly review.

Consider a general sum over a function of a general MCMC process

[TABLE]

The equilibrium lag $t$ auto-covariance is (assuming $X_{j}$ is in the invariant distribution and $t\geq 0$)

[TABLE]

This is extended to be a symmetric function of $t$ using $C_{t}=C_{\left|t\right|}$. It is estimated from the MCMC run using

[TABLE]

The Einstein Kubo sum is

[TABLE]

The auto-correlation time is

[TABLE]

The sum may be estimated from data using the self-consistent window approximation described in [22]. The static variance $C_{0}$ in (38) may be estimated using (36) rather than (37). The point of all this is the large $n$ variance approximation

[TABLE]

For our particular application, this specializes to

[TABLE]

where $\tau_{i}$ is the correlation time for the indicator function ${\bf 1}_{B_{i+1}}(X_{j})$, assuming points are generated in $M\cap B_{i}$. It may be calculated using (38).

We use this to write $N_{i,i+1}$ in terms of a random variable $\xi_{i}$ with unit variance and mean zero as

[TABLE]

Substituting this approximation into (28) gives the approximation

[TABLE]

where we expanded the term $1/N_{i,i+1}$ assuming that $n_{i}$ was large. Continuing, we substitute into (5) to get

[TABLE]

If the relative errors are all small, then products of them should be smaller still. This suggests that we get a good approximation of the overall relative error keeping only terms linear in $\zeta_{k}$ and the $\xi_{i}$. This gives

[TABLE]

Finally, we assume that the random variables $\zeta_{k}$ and $\xi_{i}$ are all independent. As with the earlier approximations, this becomes valid in the limit of large $n_{i}$. In that case, we get

[TABLE]

The one standard deviation error bar for $\widehat{Z}$ is $Z\sigma_{r}$, where the relative standard deviation is

[TABLE]

All the quantities on the right in (43) are estimated from a single run of the sampler in each $M\cap B_{i}$, with $p_{i}$ estimated as $1/\widehat{R}_{i}$ (see (28)), $\tau_{i}$ estimated from the sequence ${\bf 1}_{B_{i+1}}(X_{j})$ using a self-consistent window approximation to (38), and $\rho_{k}$ estimated as in (33).

3.3 Algorithm parameters

This section describes how we choose the base point, $x_{0}$, and the sequence of ball radii $r_{0},r_{1},\ldots,r_{k}$. A guiding motivation is the desire to minimize the variance of $\widehat{Z}$, estimated by (42), for a given amount of work. Our parameter choices are heuristic and certainly may be improved. But they lead to an overall algorithm that is “correct” in that it gives convergent estimates of the integral $Z$.

Our aim is to choose $r_{0}$ as small as possible, and $r_{k}$ as large as possible, to have a smaller overall ratio to estimate. Additionally, we expect the estimate of the final integral $\widehat{Z}_{k}$ to be significantly faster and more accurate than each ratio estimate $\widehat{R}_{i}$, since $\widehat{Z}_{k}$ is estimated by independently chosen points, so we wish $r_{k}$ to be large. With this aim in mind, the first step is to learn about $M$ through sampling. We generate $n$ points $x_{i}\in M$ by sampling the surface measure $\rho(x)=\frac{1}{Z}f(x)d\sigma(x)$, using $n$ MCMC steps of the surface sampling algorithm of Section 2. Then we choose $x_{0}$ from among these samples. Usually we would like it to be near the center of $M$, to reduce $r_{0}$ or increase $r_{k}$. Unless otherwise noted, for our implementation we choose it far from an inequality constraint boundary using

[TABLE]

The $h_{j}$ are the functions that determine the boundary of $M$ (cf. (6)). This choice may allow the smallest ball, $B_{k}$ to be large, or make it so that most points of $B_{k}$ project to points that satisfy the inequality constraints. This is heuristic in that there is no quantitative relation between $\min_{j}h_{j}(x_{i})$ and the distance (Euclidean or geodesic or whatever) from $x_{i}$ to the boundary of $M$, if $M$ has a boundary.

This is not the only possibility for choosing $x_{0}$. One could also choose it randomly from the set of sampled points. Another possibility would be to center $x_{0}$ by minimizing the “radius” of $M$ about $x_{0}$, which is the maximum distance from $x_{0}$ to another $x\in M$. We could do this by setting $x_{0}=\operatorname*{arg\,min}_{x_{i}}\max_{j}\left|x_{i}-x_{j}\right|\;.$ Having a smaller radius may mean that we need fewer balls $B_{j}$ in the integration algorithm.

Once $x_{0}$ is fixed, the radius $r_{0}$ of the biggest ball $B_{0}$ can be found using the sample points:

[TABLE]

Next consider the minimum radius $r_{k}$. We need this to be small enough that (a) $B_{k}\cap M$ is single-valued over $T_{x_{0}}$, and (b) the projection from $D_{k}\to M\cap B_{k}$ never fails. We test (b) first because it is the limiting factor. We choose a candidate radius $\tilde{r}$, sample $n\approx 10^{5}$ points uniformly from $D_{k}$, and project to $M\cap B_{k}$. If any single point fails, we shrink $\tilde{r}$ by a factor of 2 and test again. When we have a radius $\tilde{r}$ where all of the $10^{5}$ projections succeed, our error in estimating (31) should be less than about 0.001%.

Then, we proceed to test (a) approximately as follows. We sample $n$ points $x_{1},\ldots,x_{n}$ uniformly from $B_{\tilde{r}}(x_{0})\cap M$, where $B_{\tilde{r}}(x_{0})$ is the ball of radius $\tilde{r}$ centered at $x_{0}$. We then consider every pair $(x_{i},x_{j})$ and check the angle of the vector $v_{ij}=x_{i}-x_{j}$ with the orthogonal complement $T_{x_{0}}^{\perp}$. Specifically, we construct an orthogonal projection matrix onto the normal space $T_{x_{0}}^{\perp}$, for example as $P=U_{x_{0}}U_{x_{0}}^{T}$. Then, if there exists a pair for which $|v-Pv|<\texttt{tol}$, where tol is a tolerance parameter, we shrink $\tilde{r}$ and repeat the test. In practice, we shrink $\tilde{r}$ by a constant factor each time (we used a factor $C=2$). We expect this test to detect regions where $B_{\tilde{r}}(x_{0})\cap M$ is multiple-valued as long as it is sampled densely enough. Of course, if the overlapping regions are small, we may miss them, but then these small regions will not contribute much to the integral $Z_{k}$.

An alternative test (which we have not implemented) would be to sample $B_{\tilde{r}}(x_{0})\cap M$ and check the sign of the determinant of the projection to $T_{x_{0}}$, namely $\mbox{sgn}\left(\mbox{det}(U_{x_{0}}^{t}U_{y})\right)$ where $y\in B_{\tilde{r}}(x_{0})\cap M$. If $B_{\tilde{r}}(x_{0})\cap M$ is multiple-valued, we expect there to be open sets where this determinant is negative, and therefore we would find these sets by sampling densely enough.

Having chosen $x_{0}$, $r_{0}$, and $r_{k}$, it remains to choose the intermediate ball radii $r_{i}$. We use a fixed $\nu>1$ and take the radii so the balls have fixed $d$-dimensional volume ratios (note $d$ is typically different from the ball’s actual dimension $d_{a}$), as

[TABLE]

Such fixed $\nu$ strategies are often used for volume computation [12]. Given (46), the radii are $r_{i}=r_{0}\nu^{-i/d}.$ Since we have already chosen the smallest radius, this gives a relationship between $\nu$ and $k$ as $\nu=\left(r_{0}/r_{k}\right)^{d/k}$. Note that we choose $r_{k}$, the smallest radius, before we choose $k$, the number of stages.

It is not clear what would be the best choice for $\nu$. If $\nu$ is too large then the ratios $R_{i}$ (see (28)) will not be estimated accurately – indeed they will be large so $p_{i}$ in (43) will be small, increasing the relative standard deviation of the estimator. If $\nu$ is too small, then $k$ will be large and the errors at each stage will accumulate. In addition, the error at each single stage will be large if we fix the total number of MCMC points, since there will be fewer points per stage. Clearly, these issues are important for the quality of the algorithm in practice, but they are issues that are seriously under-examined in the literature. We return to them in section 5.

4 Integration: examples

This section presents some experiments with the integration algorithm of Section 3. There are three sets of computations. We first compute the surface area of the torus $\mathbb{T}^{2}$ (section 4.1.) This computation is cheap enough that we can repeat it many times and verify that the error bar formula (43) is reasonably accurate. We then calculate the volume of $SO(n)$ (section 4.2), to verify the algorithm can work on surfaces with high dimension and high co-dimension. Finally, we apply the algorithm to study clusters of sticky spheres (section 4.3), and make new predictions that could be tested experimentally.

4.1 Surface of a torus

Consider the torus of Subsection 2.3.1 with $R=1$ and $r=.5$. The exact area is $Z=4\pi^{2}rR\approx 19.7392$. The geometry of the manifold is simple enough that we can choose the parameters of the algorithm analytically. We take the initial point to be $x_{0}=(R+r,0,0)$. The smallest ball $B_{0}$ centered at $x_{0}$ such that $\mathbb{T}^{2}\subset\mathbb{T}^{2}\cap B_{0}$ has radius $r_{0}=2(R+r)$. We take the radius of the smallest ball to be $r_{k}=r$, for which the projection $T_{x_{0}}\cap B_{k}\to\mathbb{T}^{2}\cap B_{k}$ is clearly single valued. For each computation we generate $n$ points at each stage, with the total number of points fixed to be $n_{t}=kn=10^{5}$.

We calculated the volume for each number of stages $k=1,\ldots,20$. We repeated each calculation $50$ times independently. This allowed us to get a direct estimate of the variance of the volume estimator $\widehat{Z}$. Figure 7 shows that the standard deviation estimated from a single run using (43) is in reasonable agreement with the direct estimate from $50$ independent runs.333We do not strive for very precise agreement. An early mentor of one of us (JG), Malvin Kalos, has the wise advice: “Don’t put error bars on error bars.” It also shows that for this problem $k=2$ seems to be optimal. It does not help to use a lot of stages, but it is better to use more than one.

4.2 The special orthogonal group

The volume of $SO(n)$ as described in Subsection 2.3.3 is (see, e.g., [23])

[TABLE]

is the volume of the $i$th dimensional sphere,. This formula may be understood as follows. An element of $SO(n)$ is a positively oriented orthonormal basis of $\mathbb{R}^{n}$. The first basis element may be any element of the unit sphere in $\mathbb{R}^{n}$, which is $S^{n-1}$. The next basis element may be any element orthogonal to the first basis element, so it is in $S^{n-2}$, and so on. The prefactor $2^{\frac{n(n-1)}{4}}$ arises from the embedding of $SO(n)$ into $\mathbb{R}^{n\times n}$. For example, $SO(2)$ is a circle that would seem to have length $2\pi$, but in $\mathbb{R}^{4}$ it is the set $(\cos(\theta),\sin(\theta),-\sin(\theta),\cos(\theta))$, $0\leq\theta<2\pi$. This curve has length $2\sqrt{2}\pi$. The difference is the ratio $2^{\frac{2\times 1}{4}}=\sqrt{2}$.

We applied the integration algorithm to compute the volume of $SO(n)$, for $2\leq n\leq 7$. For each $n$, we ran the algorithm $N=50$ independent times and computed the mean of the estimates $\overline{V}$, and the sample standard deviation $s_{V}$. We used a total of $n_{t}=10^{5}$ points and $k=4$ steps for each $n$. Since the true value $V_{t}$ is known, we can compute the relative standard deviation $\widetilde{s_{v}}$ and the relative error $\epsilon_{r}$ via

[TABLE]

We give the results in Table 1, which shows our relative errors are small, and close the the estimated standard deviations. We notice that the errors increase as the dimension gets higher.

4.3 Sticky-sphere clusters

Next we consider a system of particles interacting with a so-called “sticky” potential, which can be thought of as a delta-function when the surfaces of two particles are in contact. This is a model for particles whose range of interaction is very short compared to the particles’ diameters, as is the case for many kinds of colloidal particles. The probabilities of finding the system in different states in equilibrium can be calculated using statistical mechanics [15] and we briefly summarize the relevant ideas.

Consider a system of $N$ unit spheres in $\mathbb{R}^{3}$, represented as a vector $x=(x_{1},\ldots,x_{N})^{t}\in\mathbb{R}^{3N}$, where $x_{i}\in\mathbb{R}^{3}$ is the center of the $i$-th sphere. We suppose there are $m$ pairs of spheres in contact $E=\{(i_{1},j_{1}),\ldots,(i_{m},j_{m})\}$. For each pair in contact there is a constraint

[TABLE]

and for each pair not in contact there is an inequality which says the spheres cannot overlap:

[TABLE]

To remove the translational degrees of freedom of the system, we fix the center of mass of the cluster at the origin. We do this by imposing three extra constraints, which together are written in vector form as

[TABLE]

We define $M$ to be the set of points $x$ satisfying the constraints (48), (50), and the inequalities (49). It is possible that $M$ has singular points where the equality constraints (48) are degenerate (have gradients that are not linearly independent), but we ignore this possibility (note that the cone in section 2.3.2 shows our algorithm may work even near singular points.) Therefore $M$ is a manifold of dimension $d=3N-m-3$ embedded in an ambient space of dimension $d_{a}=3N$, with surface area element $d\sigma$. If $M$ is not connected, our algorithm samples a single connected component.

The equilibrium probability to find the spheres in a particular cluster defined by contacts $E$ is proportional to the partition function $Z$. In the sticky-sphere limit, and where the cluster is in a container much larger than the cluster itself, the partition function is calculated (up to constants that are the same for all clusters with the same number of spheres) as [1, 15]

[TABLE]

Here $\kappa$ is the “sticky parameter”, a single number depending on both the strength of the interaction between the spheres and the temperature. This parameter is something that must be known or measured for a given system. The factor $z$ is called the “geometrical partition function” because it does not depend on the nature of the interaction between the spheres, nor the temperature. The function $f(x)$ is given by

[TABLE]

which is the product of the non-zero eigenvalues $\lambda_{i}$ of $R^{t}R$, where $R$ is half the Jacobian of the set of contact constraints (48). In the mathematical theory of rigidity, the matrix $R$ is called the rigidity matrix [24].

As a first check that our algorithm is working we calculate geometrical partition functions $z$ that are already known, by considering the 13 different clusters of $6$ spheres with $m=10$ contacts. All manifolds formed in this way are five-dimensional. For each manifold we compute $z$ as in (51) using $n_{t}=10^{8}$ points and $k=4$ steps. For each estimate we compute its error bar using formula (43). We compare these values with the values $z_{t}$ calculated in [1], which were obtained by parameterizing and triangulating each manifold and calculating the integrals using finite elements.444 In [1], the calculations were done on the quotient manifold obtained by modding out the rotational degrees of freedom. We originally tried working in a quotient space here, by fixing six coordinates; this worked for these low-dimensional manifolds but we had problems with higher-dimensional ones, particularly chains. We presume this was because of the singularity in the quotient space. To be consistent, all calculations in this paper are performed on the original manifold $M$, which contains rotational degrees of freedom. Therefore, we compare our results to the value $z_{t}=8\pi^{2}\widetilde{z}$, where $\widetilde{z}$ is the value reported in [1]. The constant $8\pi^{2}$ is the volume of the rotational space $SO(3)$.

Table 2 shows the two calculations agree.

Next, we consider a cluster of $N$ identical spheres which may form either a chain (C) or a loop (L), as in Figure 8. We are interested in the conditions under which chains are preferable to loops. Let $M_{C,N},M_{L,N}$ denote the manifolds corresponding to a chain, loop respectively; these have dimensions $2N-2$, $2N-3$. For example, for $N=10$ spheres the dimensions are 18 and 17. Integrals over these manifolds are not easily calculated by deterministic parameterizations.

For each $N=4-10$, and for each $i\in\{L,C\}$, we compute the volume $V$ and the geometrical partition function $z$ of the manifold $M_{i,N}$, using $n_{t}=10^{8}$ points and $k=4$ steps each. Because we don’t have other results to compare to, we check that our calculations are correct by comparing our calculated value of the average of $f$, namely $\overline{h}=z/V$, to that estimated by our sampling algorithm, $\tilde{h}$. The values of $\tilde{h}$, $\bar{h}$ agree within statistical error as shown in Table 3.

With the values of $z_{C},z_{L}$ we can make experimental predictions. Indeed, the ratio of the probabilities of finding a system of $N$ sticky sphere in a chain versus a loop in equilibrium is

[TABLE]

where $n_{C},n_{L}$ are the number of distinct copies of each manifold that one obtains by considering all permutations of spheres. If the spheres are indistinguishable then one can calculate555 The number of distinct copies is $N!/o$, where $o$ is the symmetry number of the clusters, i.e. the number of permutations that are equivalent to an overall rotation. For the chain, the symmetry group is $C_{2}$, the cyclic group of order 2, and therefore the symmetry number is $2$. For the loop, the symmetry group is $C_{N}$, the rotational group of a regular N-gon of order $N$, so the symmetry number is $2N$. It is also possible that the spheres are distinguishable, for example if the bonds forming the chain are unbreakable and only one bond can break and reform, in which case we would have $n_{C}{=}n_{L}{=}1$. Note that for three-dimensional clusters, the number of distinct clusters of $N$ spheres should usually be computed as $2N!/o$, where $o$ is the number of permutations that are equivalent to a rotation of the cluster (e.g. [25].) In our case, however, the clusters can be embedded into a two-dimensional plane, so a reflection of the cluster is equivalent to a rotation, hence we do not need to include the factor 2.

that $n_{C}=N!/2$ and $n_{L}=(N{-}1)!/2$. From the calculated ratios $(n_{C}z_{C})/(n_{L}z_{L})$, shown in Table 4, and the value of the sticky parameter, which must be measured or inferred in some other way, we obtain a prediction for the fraction of chains versus loops, which could be compared to experimental data. Of course, the spheres could assemble into other configurations, whose entropies we could in principle calculate, but (53) still holds without calculating these additional entropies because the normalization factors cancel.

If we don’t know the sticky parameter $\kappa$ for the system, then (53) can actually be used to measure it for a given system. Indeed, an estimate for $\kappa$ using a fixed value of $N$ would be

[TABLE]

Such a measurement could be useful, for example, in systems of emulsion droplets whose surfaces have been coated by strands of sticky DNA [26]. The DNA attaches the droplets together, like velcro, but the sticky patch can move around on the surface of a droplet, so a chain of attached droplets can sample many configurations. Understanding the strength of the droplet interaction from first principles, or even measuring it directly, is a challenge because it arises from the microscopic interactions of many weak DNA bonds as well as small elastic deformation of the droplet. But, (53) could be used to measure the stickiness simply by counting clusters in a microscope.

Finally, we note that the ratio $(n_{C}z_{C})/(n_{L}z_{L})$ is related to the loss of entropy as a chain forms a bond to become a loop. Even though a loop has more bonds, hence lower energy, it still might be less likely to be observed in equilibrium because it has lower entropy. The ratio $(n_{C}z_{C})/(n_{L}z_{L})$ is exactly the value of stickiness $\kappa$ above which a loop becomes more likely than a chain, which corresponds to a threshold temperature or binding energy. Note that as $N$ increases the ratio does too, so we expect fewer loops for longer polymers than for shorter ones.

5 Minimizing the error bars: some toy models to understand parameter dependence

Equation (43) estimates the relative standard deviation of the estimator $\widehat{Z}$ for $Z$. We would like to choose parameters in the algorithm to minimize this relative error. Previous groups studying similar algorithms in Euclidean spaces have investigated how the overall number of operations scales with the dimension [[, e.g.]]vempala, but to our knowledge there is little analysis of how such an algorithm behaves for fixed dimension, and how to optimize its performance with fixed computing power.

To this aim we construct a hierarchy of toy models that suggest how the error formula may depend on parameters. We focus exclusively on the sum $\sum_{i}(1-p_{i})\tau_{i}/(n_{i}p_{i})$ at the right of (43) which estimates the relative variance of a product of ratio estimates. This is because the estimator $\widehat{Z}_{k}$ for $\rho_{k}$ usually has a much smaller error, since it is estimated by straightforward Monte Carlo integration with independently chosen points. Our models may apply to other multi-phase MCMC algorithms.

We consider an algorithm which chooses the same number of points $n$ at each stage and has a fixed number of total points $n_{t}=kn$, where $k$ is the number of stages. Additionally, we fix the largest and smallest radii $r_{0}$, $r_{k}$, and take the $d$-dimensional volume ratios between balls to be the same constant $\nu_{i}=\nu$ at each stage. The ratio $\nu$ and number of stages $k$ are related by

[TABLE]

where we will write $C$ to denote a constant (not always the same) which does not depend on $\nu$ or $k$.

We now make the approximation that $R_{i}\propto\nu$. This would be exactly true if $M$ were a hyperplane of dimension $d$. The relative variance we wish to minimize is

[TABLE]

Recall that $\tau_{i}$ is the correlation time for the indicator function ${\bf 1}_{B_{i+1}}(X_{j})$, assuming points are generated in $M\cap B_{i}$ (see (38), (39).)

We now consider several models for the correlation times $\tau_{i}$ that will allow us to investigate the behavior of $g(\nu)$. Ultimately, we would like to find the value $\nu^{*}$ which minimizes this function.

5.1 Constant $\tau$.

Suppose the $\tau_{i}=\tau$ are equal and independent of $k$. This would happen, for example, if we were able to make independent samples and achieve $\tau_{i}=1$. It also could happen if there were a scale invariant MCMC algorithm such as the hit-and-run method or affine invariant samplers. It seems unlikely that such samplers exists for sampling embedded manifolds. In that case, we have

[TABLE]

Optimizing this over $\nu$ gives $\nu^{*}\approx 4.9$, independent of dimension. This result is interesting in that it suggests a more aggressive strategy than the previously-suggested choice $\nu=2$ [12], i.e. one should shrink the balls by a larger amount and hence have a smaller number of stages.

5.2 Diffusive scaling

As a second scenario, we model $\tau_{i}$ as

[TABLE]

This is a “diffusive scaling” because it represents the time scale of a simple diffusion or random walk in a ball of radius $r_{i}$ (see (72) in section 5.3 below.) From (46), we have $\tau_{i}\propto\nu^{-2i/d}$. Hence (56) becomes, after some algebra and using (55):

[TABLE]

Plots of $g_{d}(\nu)$ are given in Figure 9. Examples of minimizers $\nu_{d}^{*}$ are $\{2.6,2.7,3.1,3.4,3.6,3.7,4.1,4.5,4.7\}$ for $d=\{1,2,3,4,5,6,10,20,50\}$ respectively. Some analysis shows that $\nu^{*}_{d}\to\nu^{*}$ as $d\to\infty$, where $\nu^{*}$ is the optimal ratio assuming constant $\tau$ above. The same analysis shows that $g_{d}(\nu^{*}_{d})\approx\frac{d(\nu^{*}-1)}{2\log^{2}(\nu^{*})}$ for large $d$, so the relative variance in this model is proportional to the manifold’s intrinsic dimension $d$. The figure shows that $g_{d}(\nu)$ increases more slowly for $\nu>\nu^{*}$ than for $\nu<\nu^{*}$. This also suggests taking larger values of $\nu$ in computations.

We found it surprising that the model with diffusive scaling predicted smaller ratios than the one with constant scaling, since it has a correlation time which decreases as the balls shrink so one might prefer to shrink them rapidly. Some insight into why comes from thinking about the marginal cost of adding an additional ball, given we already have $k^{\prime}$ balls. For the first model, this marginal cost is constant no matter what $k^{\prime}$ is, and equal to the cost of just the first ball. For the second model, the marginal cost can decrease as $k^{\prime}$ increases, since the contribution to the cost from the $j$th ball is a decreasing function of $j$. So, with diffusive scaling, each additional ball does not increase the error by as much. Another way of seeing this is to plot the correlation times as a function of ball number, and note that the first and last balls have fixed correlation times. The total error is roughly the area under this plot. But, the area under an exponential function connecting $(0,Cr_{0}^{d})$ to $(k^{\prime},Cr_{k}^{d})$ increases less, when $k^{\prime}\to k^{\prime}+1$, than the area under a constant function connecting $(0,C^{\prime})$ to $(k^{\prime},C^{\prime})$ (here $C$, $C^{\prime}$ are constants.)

5.3 A Brownian motion model in balls.

As a third model, we approximate the MCMC sequence $X_{0},X_{1},\ldots$ with a continuous-time Brownian motion, and simplify the geometry to two concentric balls in a $d-$dimensional plane. This allows us to find an analytic formula for $\tau_{i}$. While this is a natural model to consider, we will describe later why it fails to capture the behavior of almost any system of interest. Nevertheless, the reasons behind its failure help build intuition into high-dimensional geometry, so we feel it useful to describe this mathematical exercise.

What follows is an exact derivation of the correlation times in this model. Readers interested in only the result may skip to equation (75) and the subsequent discussion.

Let $X_{t}$ be a $d$-dimensional Brownian motion on a $d$-dimensional ball $B_{0}$ of radius $R_{0}$ centered at the origin, with reflecting boundary conditions at $\left|x\right|=R_{0}$. We consider the stochastic process

[TABLE]

where $B_{1}$ is a ball with radius $R_{1}<R_{0}$, also centered at the origin. Our goal is to derive an analytic formula for the correlation time of the process, defined for a continuous-time process to be (compare with (38)) [22]

[TABLE]

This is expressed in terms of the stationary covariance function of $F(X_{t})$:

[TABLE]

The expectation is calculated assuming that $X_{0}$ is drawn from the stationary distribution, which for a Brownian motion is the uniform distribution.

We know that

[TABLE]

(as for a Bernoulli random variable), so we focus on computing the integral in (61). Notice that $\overline{F}=(R_{1}/R_{0})^{d}$.

We start by writing (62) in terms of the probability density $p(t,x,y)$ for $X_{t}=x$ given $X_{0}=y$. If $X_{0}$ is in the stationary distribution, then the joint density of $X_{0}=y$ and $X_{t}=x$ is $Z_{0}^{-1}p(t,x,y)$, where $Z_{0}=\mbox{vol}(B_{0})$. Using this expression to write (62) as an integral over $dx,dy$, integrating over time, and rearranging the integrals gives

[TABLE]

To calculate the time integral we must evaluate $\overline{p}(x,y)=\int_{0}^{+\infty}p(t,x,y)\,dt.$ Notice that $p(t,x,y)$ satisfies the forward Kolmogorov equation with boundary and initial conditions (see, e.g., [27]):

[TABLE]

where $\hat{\boldsymbol{n}}$ is the vector normal to the boundary. Therefore $\overline{p}(x,y)$ satisfies the equation

[TABLE]

Therefore $\overline{p}(x,y)$ is proportional to the Green’s function for the Laplace equation in a ball. We now wish to evaluate the integral over the $y$ coordinate

[TABLE]

where we have defined $f(x)={\bf 1}_{B_{1}}(x)-\overline{F}$ to be the centered function whose covariance we wish to evaluate. One can check that $u(x)$ solves the Poisson equation

[TABLE]

We now solve this equation analytically for $u(x)$. For $d=1$, (66) reduces to an ODE, whose solution is

[TABLE]

For $d\geq 2$, we solve (66) by looking at radial solutions $u=u(|x|)=u(r)$. In polar coordinates (66) is

[TABLE]

We solve this equation using the elementary radial solutions $\Phi(r)$ of the Laplace equation ($\Delta u=0$) in the disk,666These are given by

$\Phi(r)=\left\{\begin{aligned} &C\;\text{log}(r)+K&d=2\\ &\frac{C}{2-d}r^{2-d}+K&d\geq 3\end{aligned}\right.$

where $C$ and $K$ are constants [28].

and requiring the solution to satisfy $\lim_{r\to 0}u(r)=0$ and to be continuous at $r=R_{1}$. For $d=2$ we have

[TABLE]

and for $d\geq 3$ we have

[TABLE]

We can now evaluate the integral over the $x$ coordinate to obtain

[TABLE]

where in the last step we have written the integral using polar coordinates. Here $\mbox{vol}(S_{d-1})$ is the surface area of the boundary of the unit $d$-dimensional sphere. Substituting our analytic expression for $u(x)$, and using (63) in (61), we obtain

[TABLE]

Here $\nu=(R_{o}/R_{1})^{d}$ is the volume ratio as in (46), and the functions $h_{d}(\nu)$ are given by

[TABLE]

The functions $h_{d}(\nu)$ are proportional to the correlation time and are plotted in Figure 10 for $d=1,2,3,4,5$. For fixed $\nu$ and $d\geq 3$ the correlation time decreases with the dimension, and for all $d$ we have

[TABLE]

We found all of these observations to be surprising and will discuss them further in section 5.4.

We use (72) to approximate the correlation time $\tau_{i}$ in (56) as

[TABLE]

After some algebra, we obtain

[TABLE]

One can verify from (59), (75) that $l_{d}(\nu)=g_{d}(\nu)h_{d}(\nu)$, so this model’s only difference with that of section 5.2 is that it additionally models how $\tau_{i}$ depends on $\nu$.

Plots of the function $l_{d}(\nu)$ for dimension up to $d=5$ are given in Figure 10. The function $l_{1}(\nu)$ is strictly decreasing, while $l_{d}(\nu)$ is strictly increasing for $d\geq 2$. It follows that the ratio $\nu^{*}$ that minimizes (75) is $\nu^{*}=1$ for $d\geq 2$.

5.4 Comparing the toy models with data

Our first two models for the correlation time predicted an error that was minimized at an intermediate value $\nu^{*}$, equal to about $\nu^{*}\approx 5$ for the model with constant $\tau$, and slightly smaller for the model with diffusive scaling though approaching it as $d\to\infty$. The third was totally different – it predicted taking as small $\nu$ as possible. What does the empirical evidence suggest?

Figure 11 plots the relative standard deviation of $\widehat{Z}$ obtained empirically from averaging over several simulations, as a function of $\log(\nu)$, for several different manifolds of varying dimensions and co-dimensions. For this figure we calculated the volume using $k=1,\ldots,40$ and computed $\nu$ using (55), assuming fixed $r_{0},r_{k}$. All curves have the same basic shape, which is qualitatively the same as that predicted by the first two models but not the third.

What went wrong with the third model? The only difference it has with the diffusive scaling model is an additional assumption for how the correlation time depends on $\nu$. This correlation time was small for large $\nu$, something to be expected since a Brownian motion should leave $B_{1}$ very quickly when $B_{1}$ is much smaller than $B_{0}$. But, importantly, the correlation time was also small for small $\nu$; small enough to make the overall relative variance decrease monotonically as $\nu$ decreases. The fact that the correlation time approaches 0 as $\nu\to 1$ implies that a Brownian motion leaves $B_{1}$ very quickly even when $B_{1}$ is nearly as big as $B_{0}$.

We believe this is an artifact of our highly restrictive geometry, and almost any problem of interest will not have this property. In high dimensions, virtually all of the volume of an object will be concentrated near the extremal parts of its boundary, a fact that can be checked for example by calculating the volume of a spherical shell of width $\epsilon\ll 1$. A sphere’s boundary therefore is concentrated near the outermost spherical shell, so when a Brownian motion enters $B_{1}$ it will with overwhelming probability stay near the boundary, where it is easy for it to leave $B_{1}$. Therefore, we expect the Brownian motion to quickly enter and leave $B_{1}$ repeatedly, like a Brownian motion jumping back and forth between the boundaries of a thin strip. The correlation time therefore should be short.

When $B_{1}\cap M$, $B_{0}\cap M$ are not spherical, then some parts of their boundaries will lie further from the origin than others, and it is in these “corners” that we expect a Brownian motion to get stuck for longer times. In convex geometry such corners are known to control the computational complexity of Monte-Carlo algorithms [4], an effect that is simply not possible to capture with Brownian motion in spheres. We expect that inhomogeneous boundary densities are the norm rather than the exception in any problem of scientific interest, and therefore do not believe the predictions of this third model.

6 Conclusions

We introduced algorithms to sample and integrate over a connected component of a manifold $M$, defined by equality and inequality constraints in Euclidean space. Our sampler is a random walk on $M$ that does not require an explicit parameterization of $M$ and uses only the first derivatives of the equality constraint functions. This is simpler than other samplers for manifolds that use second derivatives. Our integration algorithm was adapted from the multi-phase Monte Carlo methods for computing high dimensional volumes, based on intersecting the manifold with a decreasing sequence of balls.

We discussed several toy models aimed at understanding how the integration errors depend on the factor by which the balls are shrunk at each step. Our models and data suggested taking larger factors than those suggested in the literature. We also estimated these errors by computing results analytically for a Brownian motion in balls, but the results of these calculations disagree with computational experiments with the sampler. We discussed the geometric intuition behind this failure.

We tested the algorithms on several manifolds of different dimensions whose distributions or volumes are known analytically, and then we used our methods to compute the entropies of various clusters of hard sticky spheres. These calculations make specific predictions that could be tested experimentally, for example using colloids interacting via sticky DNA strands.

We expect our methods to apply to other systems of objects that can be effectively modeled as “sticky,” i.e. interacting with a short-range potential [29]. For example, one could calculate the entropies of self-assembling polyhedra, a system that has been realized experimentally in the hope that it will lead to efficient drug-delivery methods [30]; one could calculate the entire free energy landscape of small clusters (building on results in [1]); or one could more accurately calculate the phase diagrams of colloidal crystals [31, 32] by computing entropic contributions to the free energy.

While our methods have worked well on the systems considered in this paper, we expect that applying them to certain larger systems will require additional techniques and ideas. For example, the manifold corresponding to the floppy degrees of freedom of a body-centered cubic crystal is thought to have a shape like a bicycle wheel, with long thin spokes leading out of a small round hub near the center of the crystal [33]. Sampling efficiently from a manifold with such different length scales will be a challenge, though one may draw on techniques developed to sample the basins of attraction of jammed sphere packings, for which the same issues arise [34]. A different, though possibly related problem, is to sample from manifolds whose constraint functions are not linearly independent. Such manifolds arise for example in clusters of identical spheres [35] or origami-based models [36], and could be important for understanding how crystals nucleate.

There are many possible improvements to the Monte Carlo methods described here. For the sampling algorithm, one improvement would be to use gradient information of the density function $f$ to create a more effective proposal distribution. It may be possible to find a manifold version of exact discrete Langevin methods described, for example, in [22] and [16]. An interesting extension to the sampler would be to allow it to jump between manifolds of different dimensions. With such a sampler, one could calculate the ratio of entropies in (53) directly, without first calculating the volumes.

Our integration algorithm could also be enhanced in several ways. It could be extended to the case when the intermediate manifolds $M\cap B_{i}$ are not connected, for example by using umbrella sampling techniques. It could be more adaptive: the radii $R_{i}$ could be chosen more systematically using variance estimators, as is done in thermodynamic integration [13]. We could use additional ratios such as $N_{i,i+2}$ (see (29)) to improve the ratio estimates. We could also let $n_{i}$ vary with ball radius so the errors at each ratio more equally distributed; according to the diffusive scaling model of section 5.2, we should use more points to estimate the larger ratios. Finally, it might be possible to change the probability distribution using smoother “steps” when breaking up the volume into a product of ratios, as in [13] and thermodynamic integration [11], though how to adapt this strategy to manifolds is not immediately clear.

Acknowledgements

E. Z. and M. H.-C. acknowledge support from DOE grant DE-SC0012296. The authors would additionally like to thank Gabriel Stoltz and Tony Lelievre for interesting discussions that helped to improve this article.