Monte Carlo on manifolds: sampling densities and integrating functions
Emilio Zappa Miranda Holmes-Cerfon Jonathan Goodman






1]Emilio Zappa  
1]Miranda Holmes-Cerfon
1]Jonathan Goodman

[1]Courant Institute of Mathematical Sciences, New York University, NY.


Monte Carlo on manifolds: sampling densities and integrating functions
    [
    
======================================================================




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.





§ 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 <cit.>, while in robotics, joints or hinges connecting moving components may be written as constraints on distances or angles <cit.>. Constraints also arise in statistics on the set of parameters in parametric models <cit.>. 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≥ 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.<cit.>).

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., <cit.>), etc.  
A move from a point x ∈ M starts with a move in the tangent space at x, which is v ∈ T_x.
This is followed by a projection back to y ∈ M, which we write as y = x + v + w.
The sampler simplifies if 
we require w ⊥ T_x, so the projection is 
perpendicular to the tangent space at the original point. 
This idea was first suggested in <cit.>[We 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 ∈ M
to minimize the projection distance y-(x+v). 




The choice of w⊥ 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 <cit.>, the Hamiltonian method
such as of <cit.>, 
geodesic methods such as <cit.>, 
and discretized SDE methods such as <cit.>.
It is common in our computational experiments that there is no w ⊥ T_x so that x+v+w ∈ M, or that one or more w values exist 
but the Newton solver fails to find one.
Section <ref> explains how we preserve detailed balance even in these situations.

Our goal is to compute integrals of the form 

    Z =   ∫_M f(x)  σ(dx)   ,

where σ(dx) is the d-dimensional surface area measure (Hausdorff measure) and 
f is a positive smooth function.
Writing B ⊂^d_a for the ball centered at a point x_0∈ M, we also estimate
certain integrals of the form 

    Z_B =   ∫_M∩ B f(x)  σ(dx)   .

For each B, we define a probability distribution on M∩ B:

    ρ_B(dx) = 1/Z_B f(x) σ(dx)   .

The sampling algorithm outlined above, and described in more detail in Section <ref>, 
allows us to draw samples from ρ_B.

To calculate integrals as in (<ref>) we use a multi-phase strategy related to nested sampling and thermodynamic integration (e.g.<cit.>).
Similar strategies were applied to volume estimation problems by the Lovasz school, 
see, e.g., <cit.>.
Let B_0 ⊃ B_1 ⊃⋯⊃ 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∩ B_i to estimate the ratio

    R_i = Z_i/Z_i+1  .

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 ⊂ B_0 
so Z = Z_0.
Our estimate of Z is 

    Z = Z_k  ∏_i=0^k-1R_i   ,

where hatted quantities such as Z are Monte Carlo estimates.
Section <ref> describes the procedure in more detail.

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

    σ^2_Z = ( Z)   ,

a procedure adapted from  <cit.>.
If we can estimate σ^2_R_i, and if these are small enough, then we can combine
them to estimate σ^2_Z.
We explore how σ_Z^2 depends on parameters of the algorithm via several simplified model problems, in section <ref>. 
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 <ref>) . 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≤ n ≤ 7 (sections <ref>, <ref>.) 

Finally, we apply our algorithm to study clusters of sticky spheres (section <ref>), a model that is studied partly to understand how colloidal particles self-assemble into various configurations <cit.>.
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 ^d_a, which is defined by equality and 
inequality constraints.
There are m equality constraints, q_i(x) = 0, i = 1, …, m,
where the q_i are smooth functions q_i : ^d_a→.
There are l inequality constraints of the form h_j(x) > 0, for j =1, …, l.
The manifold M is 

    M = { x ∈^d_a : q_i(x) = 0,   i=1, …, m,   h_j(x) > 0,   j = 1, …, l }.

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 {∇ q_i (x) }_i=1^m.
This makes Q_x the transpose of the Jacobian of the overall constraint function 
q: ^d_a→ R^m.
The entries of Q_x are 

    ( Q_x )_ij  = ∂ q_i(x)/∂ x_j  .

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≡ T_xM at a point x ∈ M is well defined. In this case {∇ q_i (x) }_i=1^m form a basis of the orthogonal space T^⊥_x ≡ T_xM^⊥. 
Note that M inherits the metric from the ambient space ^d_a by restriction. 
The corresponding d-dimensional volume element is d-dimensional Hausdorff measure, which we
denote by σ(dx).








§ SAMPLING ON A CONSTRAINT MANIFOLD
  













 §.§ The algorithm


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

    ρ(dx) = 1/Z f(x) σ(dx)   .

Here σ is the d-dimensional surface measure on M given by (<ref>).
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 <ref>.

Our Metropolis style MCMC algorithm generates a sequence X_k ∈ M with the property that 

    X_k ∼ρ                 X_k+1∼ρ  .

Once X_k = x is known, the algorithm makes a proposal y ∈ 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. <cit.> for an efficient summary of the relevant MCMC theory.

The proposal process begins with a tangential move x → x+v with v∈ 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:

    p(v|x) = 1/(2π)^d/2 s^d e^-|v|^2/2s^2  .


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 (<ref>).
This orthonormal basis is found as the last d columns of the d_a× d_a matrix Q in the QR decomposition[Please be aware of the
   conflict of notation here.
   The d_a× d_a orthogonal Q of the QR decomposition is not the d_a× d gradient Q_x.] 
of Q_x, which we find using dense linear algebra.

Given x and v, the projection step looks for w ∈ T_x^⊥ with y = x + v + w ∈ M. 
See Figure <ref>(a) for an illustration. 
It does this using an m-component column vector a = (a_1, …, a_m)^t and

    w = ∑_j=1^m a_j ∇ q_j(x) = Q_x a   .

The unknown coefficients a_j are found by solving the nonlinear equations

    q_i(x+v+Q_xa) =0 ,      i=1,… m.

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

    J_ij = ∂_a_j q_i(x + v + Q_xa) 
                  = ( ∇ q_i(x + v + Q_xa))^t ∇ q_j(x)   .

We iterate until the first time a convergence criterion is satisfied

    |q(x+v+Q_xa)| ≤  .

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  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

    _y ∈ M| x + v - y |_2   .

This produces a projection w ∈ T_y^⊥, rather than our w ∈ T_x^⊥.
It also has the advantage that such a y always exists and is almost surely unique.
The optimization algorithm that computes (<ref>) may be more robust than our 
Newton-based nonlinear equation solver.
The disadvantage of (<ref>) 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 ∈ M.
There are two cases where we must reject it immediately so the detailed balance relation may hold (which we verify in section <ref>.) 
First, we must check whether an inequality  constraint is violated, so h_i(y) ≤ 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^'∈ T_y so that 
x = y + v^' + w^' with w^'⊥ T_y.
Figure <ref>(b) illustrates this.
Since x and y are known, we may find v^' and w^' by the requirement that
x-y = v^' + w^', v^'∈ T_y, w^'∈ T_y^⊥.
This is always possible by projecting x-y onto T_y and T_y^⊥, 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^', a step that was neglected in earlier studies, e.g. <cit.>. 
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^' and see whether the 
the converge criterion (<ref>) within 
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

    a(y|x) = min( 1, f(y)p(v^'|y)/f(x)p(v|x))   .

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

    a(y|x) = min( 1, f(y)p(v^'|y) J(x|y)/f(x)p(v|x) J(y|x))   ,

where J(y|x) is the inverse of the determinant of the linearization of the projection 
v → y,  formally, |∂ v/∂ y|. 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:

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

This is evident in part b of Figure <ref>, where the angle between w and T_y is the 
same as the angle between w^' 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 <cit.>.

To prove (<ref>), 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 ξ∈ T_x, which may be written ξ = ∑ a_i u_x,i.
Let η∈ T_y be the projection of ξ to T_y normal to T_x, and write 
η = ∑ b_j u_y,j.
The orthogonality condition is u_x,i^t η = u_x,i^tξ=a_i.
This leads to the equations

    ∑_j u_x,i^t u_y,jb_j = a_i   .

These take the form Pb = a, or b = P^-1a, where P = U_x^t U_y.
Since U_x and U_y are orthonormal bases, the volume element going from x→ y is expanded by a factor

    J(y|x)^-1 = (P^-1)= ((U_x^tU_y)^-1)  .

Similar reasoning shows that the projection from T_y to T_x perpendicular to T_y 
expands volumes by a factor of

    J(x|y)^-1 =  ((U_y^tU_x)^-1)  .

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

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  
part (see Figure <ref>).
We found that setting the parameter   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. 






 §.§ Verifying the balance relation


In this section we verify what may already be clear, that the algorithm described above 
satisfies the desired balance relation (<ref>). 
 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, 

    P_x(dy) = h(y|x) σ(dy) + ζ(x) δ_x(dy).

Here δ_x(dy) is the unit point mass probability distribution at the point x, 
ζ(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∈ M from x. 
The proposal density is the product of three factors:

    h(y|x) = (1- 1_F_x(y))  ·  p(v|x)  ·  J(y|x)   ,

where the set F_x⊂ M consists of those points y∈ M that cannot be reached using our Newton solver starting from x,
and 1_F_x(y) denotes the characteristic function of the set F_x, which 
equals 1 if y ∈ F_x and zero if y ∉ F_x. 
Note that we must have ∫_M h(y|x)  dσ(y)  + ζ(x) = 1.


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

    R_x(dy) = a(y|x) h(y|x) σ(dy) + ξ(x)δ_x(dy)   .

Here, ξ(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

    ξ(x) = ζ(x) + ∫_M ( 1 - a(y|x) )h(y|x)σ(dy)   .


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

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

For this, suppose X_k has probability distribution on M with density g_k(x) with respect 
to surface area measure:

    X_k ∼ρ_k(dx) = g_k(x) σ(dx)   .

We check that if (<ref>) 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 (<ref>).
To start, 

    g_k+1(x) = ∫_M g_k(y) a(x|y) h(x|y) dσ(y)  + ξ(x)g_k(x)   .

The integral on the right represents jumps to x from y ≠ 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 (<ref>) this becomes

    g_k+1(x) = ∫_M f(x) a(y|x) h(y|x) dσ(y)  + ξ(x)f(x)   .

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

    ( ) = ∫_M a(y|x)h(y|x) dσ(y)   .

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



We use the standard Metropolis Hastings formula to enforce (<ref>):

    a(y|x) = min(1, f(y)h(x|y)/f(x)h(y|x))   .

If we succeed in proposing y from x, then we know y ∉ F_x so we are able to evaluate h(y|x). 
However, it is still possible that x ∉ 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 , so we are sampling a “fattened” region M_ = {y: |q_i(y)| ≤  i=1,…,m,  h_j(x) > 0,   j = 1, …, l }.  It is well-known (e.g. <cit.>) 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 → 0. 

We do not know the distribution by which our numerical solver produces a point in M_. However, since we observe the correct distributions on M in the examples below (for small enough ), 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 → 0.
Moreover, it is practical to take  to be extremely small, on the order of machine
precision, because of the fast local convergence of Newton's method.





 §.§ Examples
 

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





  §.§.§ Torus
  
Consider a torus 𝕋^2 embedded in ^3, implicitly defined by

    𝕋^2 = { (x,y,z) ∈^3 : ( R- √(x^2+y^2))^2 + z^2-r^2 = 0 },

where R and r are real positive numbers, with R > r. 
Geometrically, 𝕋^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 𝕋^2 is given by

    𝕋^2 = {(R + r cos(ϕ) cos(θ), (R + r cos(ϕ) sin(θ), r sin(ϕ) ) : θ, ϕ∈ [0,2π] }.


We ran N = 10^6 MCMC steps to sample uniform measure (f(x) = 1) on 𝕋^2
with toroidal radius R = 1, poloidal radius r = .5, and step size scale s = .5. 
Figure <ref> shows our empirical marginal distributions and the exact theoretical distributions. 
The empirical distributions of the toroidal angle θ and the poloidal angle ϕ
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. 







  §.§.§ Cone

 We consider the circular right-angle cone 𝒞⊂^3
with vertex (0,0,0) given by

    𝒞 = { (x,y,z) ∈^3 : z-√(x^2+y^2) = 0 ,   x^2+y^2 < 1,   z > 0}.

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 (<ref>).
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 𝒞 would not be a differentiable manifold because of 
the singularity at the vertex.

We suspected that the projection steps (see Figure <ref>) 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 <ref> 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

    g_X(x) = 2/π√(1-x^2) ,  
          g_Y(y) = 2/π√(1-y^2) ,  
          g_Z(z) =   2z   .


 






  §.§.§ Special orthogonal group
 
We apply our surface sampler to the manifold SO(n), which is the
the group of n × n orthogonal matrices with determinant equal to 1. 
We chose SO(n) because these manifolds have high dimension, 1/2n(n-1), and high co-dimension, 1/2n(n+1) <cit.>. 
We view SO(n) as the set of n × n matrices, x ∈^n× n that satisfy the row 
ortho-normality constraints for k = 1, …, n and l > k:

    g_kk(x) = ∑_m=1^n x_km^2 = 1    ,      g_kl(x) = ∑_m=1^n x_kmx_lm = 0   .

This is a set of 1/2n(n+1) equality constraints.
The gradient matrix Q_x defined by (<ref>) may be seen to have full rank for all x ∈ SO(n).
Any x satisfying these constraints has (x) = ± 1.
The set with (x) = 1 is connected.
It is possible that our sampler would propose an x with (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 = (x),
which converges to a standard normal as n →∞ <cit.>.
This is surprising in view of the fact that (x) is the sum of n random numbers each of
which has (x_kk) = 1.
If the diagonal entries were independent with this variance, then (T) = n^1/2
instead of the correct (T) = 1.
This makes the correctness of the T distribution an interesting test of the correctness 
of the sampler.

Figure <ref> 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 ∈^n× n with i.i.d. entries from N(0,1) and use the 
QR decomposition y = Rx, where R is upper triangular and x is orthogonal.















§ INTEGRATING OVER A MANIFOLD
  



 §.§ Algorithm


In this section we describe an algorithm to calculate integrals over a manifold M of the form (<ref>).
We assume that M given in (<ref>) is bounded. 
The strategy as outlined in the Introduction is to 
 consider a sequence of d_a-dimensional balls B_0⊃ B_1⊃…⊃ B_k, centered at 
x_0 ∈ M with radii r_i, for i = 1,…, k, where B_0⊃ M entirely contains the manifold, 
and B_k is “small enough" in a way we describe momentarily. 
We assume that each M∩ B_i is connected, an assumption that is critical for our algorithm to work. 
We define the collection of integrals (see (<ref>))

    Z_i = ∫_B_i ∩ M f(x) dσ(x)   .

The final estimate is Z ≈Z according to (<ref>).

To estimate the ratio R_i = Z_i/Z_i+1 we generate n_i points in B_i∩ M with the probability distribution ρ_i as in (<ref>). 
This is done using the MCMC algorithm given in Section <ref> by including the extra inequality 
|x-x_0|^2 < r_i^2. 
An estimator for R_i is then given by

    R_i  = n_i/N_i,i+1,

where, for i≤ j, we define

    N_i,j = # points generated in B_i ∩ M that also lie in B_j ∩ M.

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∩ 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⊂ 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 <cit.>.  
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 ∩ B_k and points in the tangent plane D_k = T_x_0∩ B_k.
This happens, for example, when the manifold M∩ B_k is a graph above D_k. 
We use the coordinates in the flat disk D_k to integrate over M∩ B_k, which is simply the surface integral as illustrated 
in Figure <ref>:

    Z_k =  ∫_M∩ B_k f(y)  dσ(y) = ∫_D_k f(y(x)) J^-1(x)  dx   .

This uses the following notation.
The integral on the right is over the d-dimensional disk.
For any x ∈ D_k, the projection to M perpendicular to T_x_0 is called y(x)∈ M. 
For small enough r_k, the projection is easily calculated using the projection
algorithm of Figure <ref>.
The Jacobian factor J(x)  = det(U_x_0^tU_y(x)) is known from (<ref>), 



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 ∈ 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 ∩ B_k, so we need an extra indicator function
factor, 1_B_k(y_i), that is equal to 1 if y_i ∈ M∩ B_k.
The d dimensional volume of D_k is known, so 



 we have the direct Monte Carlo estimate

    Z_k = 
               _d(D_k)/n_k ∑_i=1^n_k 1_B_k(y_i)f(y_i) J_i^-1  ,
                   _d(D_k) = π^d/2/Γ(d/2+1)  r_k^d   .

The estimators (<ref>) and (<ref>) are substituted in (<ref>) 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∩ B_k must never fail. If it does, the estimate (<ref>) 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∩ B_k be a graph above D_k. 

   



 §.§ Variance estimate


This section serves two purposes.
One is to present the method we use to estimate the variance, σ^2_Z, 
of the estimator (<ref>).
This is similar to the strategy of <cit.>, 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 <ref>.

We start with the simple estimator Z_k in (<ref>) 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 

    Z_k = Z_k + σ_k ζ_k 
                       = Z_k ( 1 + ρ_kζ_k)   ,

where ζ_k is a zero mean random variable with (ζ_k) = 1,
σ_k is the absolute standard deviation of Z_k, and  
ρ_k=σ_k/Z_k is the relative standard deviation.
Since the samples in (<ref>) are independent, the standard deviation satisfies

    σ_k^2 =  ^2_d(D_k)/n_kVar(G(y)),      G(y) =   1_B_k(y)f(y) J^-1(y).

We estimate Var(G(y)) with the sample variance and 



 obtain, for the absolute and relative variances, 

    σ_k^2 ≈^2_d(D_k)/n_k^2∑_i=1^n_k (G(y_i)-G)^2,     ρ_k^2 ≈σ_k^2/Z_k^2.

Here G = ∑_i=1^n_k1/n_k G(y_i) is the sample mean. 

Now consider the estimator R_i in (<ref>) for the ith volume ratio. This involves the number count N_i,i+1. We start by approximating  (N_i,i+1) and then incorporate this in the approximation for (R_i). 



The number count may be written

    N_i,i+1 = ∑_j=1^n_i 1_B_i+1(X_j)   ,

where X_j ∈ M ∩ B_i is our sequence of samples. 
In the invariant probability distribution of the chain, the probability of X_j ∈ B_i+1 is

    p_i = 𝔼[  1_B_i+1(X_j)] = 1/R_i = Z_i+1/Z_i  .

Since 1_B_i+1(X_j) is a Bernoulli random variable, we have (in the invariant distribution)

    (  1_B_i+1(X_j)) = p_i (1-p_i)   .

We will call this the static variance. It alone is not enough to determine (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. <cit.>) of error bars for MCMC estimators of 
the form (<ref>), which we now briefly review. 

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

    S = ∑_j=1^n F(X_j)   .

The equilibrium lag t auto-covariance is (assuming X_j is in the invariant distribution
and t ≥ 0)

    C_t = [ F(X_j),F(X_j+t)]   .

This is extended to be a symmetric function of t using C_t = C_|t|.
It is estimated from the MCMC run using

    F≈1/n_i∑_j=1^n_i F(X_j)   ,     
        C_t ≈1/n_i-t∑_j=1^n_i-t( F(X_j)- F)( F(X_j+t)- F)   .

The Einstein Kubo sum is

    D = ∑_t = - ∞^∞ C_t   .

The auto-correlation time is

    τ = D/C_0 = 1 + 2/(F(X_j))∑_t=1^∞ C_t   .

The sum may be estimated from data using the self-consistent window approximation described in <cit.>.
The static variance C_0 in (<ref>) may be estimated using 
(<ref>) rather than (<ref>).
The point of all this is the large n variance approximation

    (S) ≈ nD = nC_0τ  .


For our particular application, this specializes to

    (N_i,i+1) ≈  np_i(1-p_i)τ_i   ,

where τ_i is the correlation time for the indicator function 1_B_i+1(X_j), assuming points are generated in M∩ B_i. It may be calculated using (<ref>). 

We use this to write N_i,i+1 in terms of a random variable ξ_i 
with unit variance and mean zero as 

    N_i,i+1
                    ≈ n_ip_i ( 1 + √((1-p_i)τ_i/n_ip_i) ξ_i )   .

Substituting this approximation into (<ref>) gives the approximation




    R_i≈1/p_i( 1 -   √((1-p_i)τ_i/n_ip_i) ξ_i )   .

where we expanded the term 1/N_i,i+1 assuming that n_i was large. 

Continuing, we substitute into (<ref>) to get

    Z
                      
                     ≈ Z [( 1 + ρ_k ζ_k)
                         ∏_i=0^k-1( 1 -   √((1-p_i)τ_i/n_ip_i) ξ_i )] .

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 ζ_k and the ξ_i.
This gives

    Z≈ Z 
             [ 1 + ρ_k ζ_k + ∑_i=0^k-1√((1-p_i)τ_i/n_ip_i) ξ_i ]   .

Finally, we assume that the random variables ζ_k and ξ_i are all independent.
As with the earlier approximations, this becomes valid in the limit of large n_i.
In that case, we get

    (Z) 
           = σ_Z^2      ≈     Z^2 ( ρ_k^2 + ∑_i (1-p_i)τ_i/n_ip_i)    .

The one standard deviation error bar for Z is Zσ_r, where the relative 
standard deviation is

    σ_r ≈√(ρ_k^2 + ∑_i=0^k-1(1-p_i)τ_i/n_ip_i)  .


All the quantities on the right in (<ref>) are estimated from a single run of the sampler in each M∩ B_i, with p_i estimated as 1/R_i  (see (<ref>)), τ_i estimated from 
the sequence  1_B_i+1(X_j) using a self-consistent window approximation to (<ref>), 
and ρ_k estimated as in (<ref>).





 §.§ Algorithm parameters


This section describes how we choose   
the base point, x_0, and the sequence of ball radii r_0, r_1,…,r_k. 
A guiding motivation is the desire to minimize the variance of Z, 
estimated by (<ref>), 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 Z_k to be significantly faster and more accurate than each ratio estimate R_i, since 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 ∈ M by sampling the surface measure 
ρ(x) = 1/Z f(x) dσ(x), using n MCMC steps of the surface sampling algorithm of Section <ref>.
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 

    x_0 = _x_i = x_1, …, x_n  min_j   h_j(x_i)   .

The h_j are the functions that determine the boundary of M (cf. (<ref>)).
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 ∈ M. 
We could do this by setting x_0 = _x_i
   max_j | x_i - x_j|   .
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: 

    r_0 = max{ |x_0-x_i| : i = 1, …, n }.









Next consider the minimum radius r_k. We need this to be small enough that (a) B_k∩ M is 
single-valued over T_x_0,  and (b) the projection from D_k→ M∩ B_k never fails. 
We test (b) first because it is the limiting factor. 
We choose a candidate radius r̃, sample n≈ 10^5 points uniformly from D_k, and project to M∩ B_k. If any single point fails, we shrink r̃ by a factor of 2 and test again. 
When we have a radius r̃ where all of the 10^5 projections succeed, our error in estimating (<ref>) should be less than about 0.001%. 

Then, we proceed to test (a) approximately as follows. 
We sample n points x_1, …, x_n uniformly from B_r̃(x_0) ∩ M, where B_r̃(x_0) is the ball of radius 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^⊥. 
Specifically, we construct an orthogonal projection matrix onto the normal space T_x_0^⊥, for example as P = U_x_0U_x_0^T. 
Then, if there exists a pair for which |v-Pv| <, where  is a tolerance parameter, we shrink r̃ and repeat the test. In practice, we shrink r̃ by a constant factor each time (we used a factor C = 2). 
We expect this test to detect regions where B_r̃(x_0) ∩ 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_r̃(x_0) ∩ M and check the sign of the determinant of the projection to T_x_0, namely ( (U_x_0^tU_y)) where y∈ B_r̃(x_0) ∩ M. If B_r̃(x_0) ∩ 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 ν>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

    (r_i/r_i+1)^d = ν  .





 

Such fixed ν strategies are often used for volume computation <cit.>. 
Given (<ref>), the radii are

r_i = r_0 ν^-i/d  .

Since we have already chosen the smallest radius, this gives a relationship between ν
and k as

        ν = (r_0/r_k)^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 ν.
If ν is too large then the ratios R_i (see (<ref>))
will not be estimated accurately – indeed they will be large so p_i in (<ref>)  will be small, increasing the relative standard deviation of the estimator. 
If ν 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 <ref>. 








§ INTEGRATION: EXAMPLES
 

This section presents some experiments with the integration algorithm of Section <ref>.
There are three sets of computations.
We first compute the surface area of the torus 𝕋^2 (section <ref>.)
This computation is cheap enough that we can repeat it many times and verify that the 
error bar formula (<ref>) is reasonably accurate.
We then calculate the volume of SO(n) (section <ref>), 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 <ref>), and make new predictions that could be tested experimentally.  



 §.§ Surface of a torus

 Consider the torus of Subsection <ref> with R=1
and r = .5.
The exact area is Z = 4 π^2 rR ≈ 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 𝕋^2⊂𝕋^2∩ 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∩ B_k →𝕋^2∩ 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, …, 20.
We repeated each calculation 50 times independently.
This allowed us to get a direct estimate of the variance of the volume estimator Z.
Figure <ref> shows that the standard deviation estimated from a single run using (<ref>)
is in reasonable agreement with the direct estimate from 50 independent runs.[We 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.





 §.§ The special orthogonal group
 
The volume of SO(n) as described in Subsection <ref> is (see, e.g., <cit.>)

    (SO(n)) = 2^n(n-1)/4∏_i=1^n-1(S^i),
        where  (S^i) = 2π^i+1/2/Γ(i+1/2)
 
is the volume of the ith dimensional sphere,. 
This formula may be understood as follows.
An element of SO(n) is a positively oriented orthonormal basis of ℝ^n.
The first basis element may be any element of the unit sphere in ℝ^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^n(n-1)/4 arises from the embedding of SO(n) into 
ℝ^n× n.
For example, SO(2) is a circle that would seem to have length 2π, but 
in ℝ^4 it is the set (cos(θ), sin(θ), -sin(θ), cos(θ)),
0 ≤θ < 2π.
This curve has length 2 √(2)π.
The difference is the ratio 2^2× 1/4 = √(2).

We applied the integration algorithm to compute the volume of SO(n), for 2≤ n ≤ 7. 
For each n, we ran the algorithm N = 50 independent times and computed the mean 
of the  estimates 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 s_v and the relative error ϵ_r via

    s_V = s_V/√(N) V_t,     ϵ_r = |V-V_t/V_t|.

We give the results in Table <ref>, 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.  





 §.§ 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 <cit.> and we briefly summarize the relevant ideas. 

Consider a system of N unit spheres in ℝ^3, represented as a vector 
x = (x_1, …, x_N)^t ∈^3N, where x_i ∈^3 is the center of the i-th sphere. 
We suppose there are m pairs of spheres in contact E = {(i_1,j_1), …, (i_m,j_m)}. 
For each pair in contact there is a constraint

    |x_i-x_j|^2-1 = 0,      (i,j)∈ E,
 
and for each pair not in contact there is an inequality which says the spheres cannot overlap:

    |x_i - x_j|> 1,      (i,j)∉ E.

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 

    ∑_i=1^N x_i = 0   .


We define M to be the set of points x satisfying the constraints (<ref>), 
(<ref>), and the inequalities (<ref>). 
It is possible that M has singular points where the equality constraints (<ref>) 
are degenerate (have gradients that are not linearly independent), but we ignore this possibility (note that the cone in section <ref> 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σ. 
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 <cit.>

    Z = κ^m z,     where   z = ∫_M f(x)  dσ(x).

Here κ 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




    f(x) = ∏_i=1^3N-m-3λ_i(x)^-1/2,



which is the product of the non-zero eigenvalues λ_i of R^tR, where R is half the Jacobian of the set of contact constraints (<ref>). In the mathematical theory of rigidity, the matrix R is called the rigidity matrix <cit.>. 

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 (<ref>) using n_t = 10^8 points and k = 4 steps. For each estimate we compute its error bar using formula (<ref>).   We compare these values with the values z_t calculated in <cit.>, which were obtained by  parameterizing and triangulating each manifold
and calculating the integrals using finite elements.[
In <cit.>, 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 π^2 z, where z is the value reported in <cit.>. The constant 8π^2 is the volume of the rotational space SO(3).
] 
Table <ref> 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 <ref>.  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 ∈{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 h = z/V, to that estimated by our sampling algorithm, h̃. The values of h̃, h̅ agree within statistical error as shown in Table <ref>. 

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 
 
    P(chain)/P(loop) = κ^-1n_C z_C/n_L z_L ,

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 calculate[
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. <cit.>.) 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_Cz_C)/(n_Lz_L), shown in Table <ref>, 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 (<ref>) still holds without calculating these additional entropies because the normalization factors cancel. 

 If we don't know the sticky parameter κ for the system, then (<ref>) can actually be used to measure it for a given system. Indeed, an estimate for κ using a fixed value of N would be 
 
    κ ≈ κ = # of loops of length N/# of chains of length N·n_C z_C/n_L z_L .

 
Such a measurement could be useful, for example, in systems of emulsion droplets whose surfaces have been coated by strands of sticky DNA <cit.>. 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, (<ref>) 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_Lz_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_Lz_L) is exactly the value of stickiness κ 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. 
 
 
 
 
 
 
 
 
 
 
 


§ MINIMIZING THE ERROR BARS: SOME TOY MODELS TO UNDERSTAND PARAMETER DEPENDENCE


Equation (<ref>) estimates the relative standard deviation of the estimator 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 <cit.>, 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 ∑_i (1-p_i)τ_i / (n_ip_i) at the right of (<ref>) which estimates the relative variance of a product of ratio estimates. This is because the estimator Z_k for ρ_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 ν_i = ν at each stage. The ratio ν and number of stages k are related by 

    ν^k = (r_0/r_k)^d = C,

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

We now make the approximation that R_i ∝ν. This would be exactly true if M were a hyperplane of dimension d. 
The relative variance we wish to minimize is

    g(ν) = C ν-1/log(ν)∑_i=0^k-1τ_i   .

 Recall that τ_i is the correlation time for the indicator function  1_B_i+1(X_j), assuming points are generated in M∩ B_i (see (<ref>), (<ref>).) 
 
 We now consider several models for the correlation times τ_i that will allow us to investigate the behavior of g(ν). Ultimately, we would like to find the value ν^* which minimizes this function. 




 §.§ Constant τ.


Suppose the τ_i= τ are equal and independent of k.
This would happen, for example, if we were able to make independent samples and achieve τ_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

    g(ν) = C ν-1/log(ν)^2  .

Optimizing this over ν gives ν^* ≈ 4.9, independent of dimension.
This result is interesting in that it suggests a more aggressive strategy than
the previously-suggested choice ν = 2 <cit.>, i.e. one should shrink the balls by a larger amount and hence have a smaller number of stages.




 §.§ Diffusive scaling


As a second scenario, we model τ_i as 

    τ_i ∝ r_i^2.

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 (<ref>) in section <ref> below.)
From (<ref>), we have τ_i ∝ν^-2i/d. 
Hence (<ref>) becomes, after some algebra and using (<ref>):

    g(ν) ∝ g_d(ν) ≡ν-1/log(ν)(1-ν^-2/d) .


Plots of g_d(ν) are given in Figure <ref>.
Examples of minimizers ν_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 ν^*_d →ν^* as d →∞, where ν^* is the optimal
ratio assuming constant τ above.
The same analysis shows that g_d(ν^*_d)≈d(ν^*-1)/2log^2(ν^*) 
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(ν) increases more slowly for ν > ν^* than for ν < ν^*.
This also suggests taking larger values of ν 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' balls. For the first model, this marginal cost is constant no matter what k' is, and equal to the cost of just the first ball. For the second model, the marginal cost can decrease as k' increases, since the contribution to the cost from the jth 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',Cr_k^d) increases less, when k'→ k'+1, than the area under a constant function connecting (0,C') to (k',C') (here C, C' are constants.) 




 §.§ A Brownian motion model in balls.


As a third model, we approximate the MCMC 
sequence X_0,X_1,… 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 τ_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 (<ref>)  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 |x|=R_0.
We consider the stochastic process 

    F(X_t) =  1_B_1(X_t),

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 (<ref>)) <cit.>

    τ = 2 ∫_0^+∞C(t)/C(0)  dt   .
 
This is expressed in terms of the 
stationary covariance function of F(X_t): 

    C(t) = 𝔼[ ( 1_B_1(X_0) - F)
                                  ( 1_B_1(X_t) - F)],     F = 𝔼[  1_B_1(X)] .

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 

    C(0) = (F(X_t)) = F(1-F)

 (as for a Bernoulli random variable), so we focus on computing the integral in (<ref>).  
 Notice that F = (R_1/R_0)^d. 

We start by writing (<ref>) 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 = (B_0).  
Using this expression to write
(<ref>) as an integral over dx,dy, integrating over time, and rearranging the integrals gives 

    ∫_0^∞   C(t) dt = Z_0^-1∫_B_0∫_B_0( 1_B_1(x) - F)
                                      ( 1_B_1(y) - F)∫_0^∞ p(t,x,y)  dtdydx   .

To calculate the time integral we must evaluate
p(x,y) = ∫_0^+∞ p(t,x,y)  dt .
Notice that p(t,x,y) satisfies the forward Kolmogorov  
equation with boundary and initial conditions (see, e.g., <cit.>):

    p_t = 1/2Δ_x p     (x ∈ B_0) ,   ∇ p · = 0     (x ∈∂ B_0) ,   
     p(0,x,y) = δ(x-y) ,
  
where  is the vector normal to the boundary. 
Therefore p(x,y) satisfies the equation

    Δ_x p(x,y) = ∫_0^+∞Δ_x p(t,x,y)  dt  
                                   =  2 ∫_0^+∞ p_t(t,x,y) dt 
                                   
                                   = -2 δ(x-y)  .

Therefore 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

    u(x)  = ∫_B_0 f(y)p(x,y) dy ,

where we have defined
f(x) =  1_B_1(x) - F to be the centered function whose covariance we wish to evaluate. 
One can check that u(x) solves the Poisson equation 

    Δ_x u = -2 f   for  x ∈ B_0,     ∇ u · = 0   for x ∈∂ B_0 .
 
We now solve this equation analytically for u(x). 
For d =1, (<ref>) reduces to an ODE, whose solution is 

    u(x) = {( R_1/R_0-1 ) x^2       |x| < R_1 
    R_1/R_0 x^2-2R_1|x|+R_1^2       R_1 < |x| < R_0
    .


For d ≥ 2, we solve (<ref>) by looking at radial solutions u = u(|x|) = u(r). 
In polar coordinates (<ref>) is

    u_rr + d-1/ru_r  = -2 f   for  r ∈ (0,R_0),      
     u_r(R_0) = 0.
 
We solve this equation using the elementary radial solutions  Φ(r) of the Laplace equation (Δ u = 0) in the disk,[These are given by 

    Φ(r) = {    C  log(r) + K     d = 2 
       C/2-dr^2-d+ K     d ≥ 3
    .
 
where C and K are constants <cit.>.
]
and requiring the solution to satisfy lim_r → 0 u(r) = 0 and to be continuous at r = R_1. 
For d = 2 we have

    u(r) = {   R_1^2-R_0^2/R_0^2 r^2       0≤ r≤ R_1 
        R_1^2(logR_1-logr)+R_1^2/2(r^2/R_0^2-1)     R_1 < r ≤ R_0
    .

and for d ≥ 3 we have

    u(r) = {   R_1^d-R_0^d/dR_0^d r^2     0 ≤ r ≤ R_1 
        -2R_1^d/d(2-d) r^2-d + R_1^d/dR_0^dr^2 + R_1^2/2-d    R_1 ≤ r ≤ R_0
    .


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

    ∫_0^+∞ C(t)  dt =Z_0^-1∫_B_0 f(x) u(x)   dx = (S_d-1)Z_0^-1∫_0^R_0 f(r)u(r)r^d-1dr


where in the last step we have written the integral using polar coordinates. 
Here (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 (<ref>) in (<ref>), we obtain

    τ = R_0^2 h_d(ν)  .

Here ν = (R_o/R_1)^d is the volume ratio as in (<ref>), and the functions h_d(ν) are given by

    h_d(ν) = {   4/3ν-1/ν^2    d = 1 
       ν^-1-1+log(ν)/ν-1    d = 2 
       4/d^2-4(d-2)ν^-1 -dν^2/d-1+2/ν^2/d(1-ν^-1)    d ≥ 3
    .



The functions h_d(ν) are proportional to the correlation time and are plotted in Figure <ref>  for d = 1,2,3,4,5. For fixed ν and d≥ 3 the correlation time decreases with the dimension, and for all d we have 

    lim_ν→ 1 h_d(ν) = lim_ν→ +∞ h_d(ν) = 0.

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



We use (<ref>) to approximate the correlation time τ_i in (<ref>) as 

    τ_i = r_i^2 h_d(ν).
 
 
After some algebra, we obtain

    g(ν) ∝ l_d(ν) ≡(ν-1)h_d(ν)/log(ν)(1-ν^-2/d)


One can verify from (<ref>), (<ref>) that l_d(ν) = g_d(ν)h_d(ν), so this model's only difference with that of section <ref> is that it additionally models how τ_i depends on ν. 

Plots of the function l_d(ν) for dimension up to d = 5 are given in Figure <ref>. The function l_1(ν) is strictly decreasing, while l_d(ν) is strictly increasing for d≥2. It follows that the ratio ν^* that minimizes (<ref>) is ν^* = 1 for d ≥ 2. 





 §.§ Comparing the toy models with data





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

Figure <ref> plots the relative standard deviation of Z obtained empirically from averaging over several simulations, as a function of log(ν), for several different manifolds of varying dimensions and co-dimensions. 
For this figure we calculated the volume using k = 1, …, 40 and computed ν using (<ref>), 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 ν. 
This correlation time was small for large ν, 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 ν; small enough to make the overall relative variance decrease monotonically as ν decreases. 
The fact that the correlation time approaches 0 as ν→ 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 ϵ≪ 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∩ M, B_0∩ 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 <cit.>, 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. 








§ 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 <cit.>. 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 <cit.>; one could calculate the entire free energy landscape of small clusters (building on results in <cit.>); or one could more accurately calculate the phase diagrams of colloidal crystals <cit.> 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 <cit.>. 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 <cit.>.   
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 <cit.> or origami-based models <cit.>, 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 <cit.> and <cit.>.
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 (<ref>) 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∩ 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 <cit.>. 
We could use additional ratios such as N_i,i+2 (see (<ref>)) 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 <ref>, 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 <cit.> and thermodynamic integration <cit.>, 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. 



unsrt








