A faster direct sampling algorithm for equilateral closed polygons and the probability of knotting

We present a faster direct sampling algorithm for random equilateral closed polygons in three-dimensional space. This method improves on the moment polytope sampling algorithm of Cantarella et al (2016 J. Phys. A: Math. Theor. 49 275202) and has (expected) time per sample quadratic in the number of edges in the polygon. We use our new sampling method and a new code for computing invariants based on the Alexander polynomial to investigate the probability of finding unknots among equilateral closed polygons.


I. INTRODUCTION
The tendency of long, flexible structures to become entangled is familiar to anyone who has managed an electrical extension cord, carried headphones in their pockets, or brushed a child's hair.This phenomenon has been studied for several decades in statistical physics (cf.[33]).A standard ensemble for these investigations is the space of equilateral random polygons.These are polygonal walks in 3-space forming closed loops and consisting of unit-length steps.The closure condition imposes subtle global correlations between edge directions, which means it is not obvious how to generate random equilateral polygons.Indeed, algorithms have been proposed for at least 4 decades [2, 9, 16-18, 22, 27, 29, 30, 36, 37], though most are numerically unstable or have not been proved to sample from the correct probability distribution.
In previous work [8], we introduced the action-angle method (AAM), which is a numerically stable and provably correct algorithm for generating random equilateral n-gons in R 3 based on rejection sampling the hypercube.The action-angle method is the fastest extant method: it produces samples in expected time Θ(n 5/2 ).Xiong et al. [38] used large-scale computing resources to sample 1.6 × 10 9 polygons using AAM, classifying their knot type with three invariants based on the Alexander polynomial.These invariants can be computed extremely quickly using sparse matrix methods (we estimate O(n 1.18 ), see Figure 3) and AAM is faster than invariant computation only for n < 280.In this paper, we give an improved algorithm, the progressive action angle method (PAAM), which we prove produces samples in Θ(n 2 ) time.Empirically, PAAM is faster than AAM for n > 20 and faster than invariant computation for n < 1850.
One of the most interesting findings in [38] concerns the probability P 0 1 (n) of finding an unknotted equilateral polygon among random polygons.Diao [14] proved that P 0 1 (n) is at most exponential in n, and it had been widely assumed that (up to subdominant terms) that P 0 1 (n) ≃ C 0 1 exp(−n/n 0 1 ) for some constants C 0 1 and n 0 1 [28].(To be precise, expressions of the form P 0 1 (n) ≃ C 0 1 n v 0 exp(−n/n K ) appeared early on in the literature (cf.[13, eq. (1.4)], [26, eq. ( 2)], [11, eq.(1.3)]) but the power v 0 was thought to be close to zero based on the data sets available at the time.)This made the unknot an anomaly, since every other knot type studied fit to the general form P K (n) ≃ C K n v K exp(−n/n K ).Xiong et al. [38] presents strong evidence that the probability of a given knot type is in the form where v 0 ≃ −0.19 ± 0.001 and n 0 ≃ 259.3 ± 0.2 are constants that do not depend on the knot type, and n p (K) is the number of prime components of the knot type.Deguchi and Uehara [12] proposed a very similar model (with a different finite-size correction) and, as Xiong et al. point out, a model of a similar form was proposed for self-avoiding polygons on lattices by Orlandini et al. [32].
In this model, the unknot is better thought of as "a composite knot with zero components" instead of as an anomalous knot type.If this is true, then with free parameters C, β, and γ instead of with free parameters N, β and γ.Xiong et al. find that (1) gives a good fit, but (2) does not.
Using PAAM, we were able to generate data for unknot probabilities for n up to 3043 (see Table I) using a desktop computer.We confirm Xiong et al.'s finding that (1) fits the data very well and we obtain similar estimates for the free parameters.However, unlike those authors, we find that (2) fits the data equally well, though with very different values for β and γ.We conclude that more will be required to distinguish between these models.

II. THE PROGRESSIVE ACTION-ANGLE METHOD
For an n-gon in R 3 , let v 1 , . . ., v n ∈ R 3 be the coordinates of its vertices, and let e 1 , . . ., e n be the edge vectors, meaning that e i = v i+1 − v i for i = 1, . . ., n − 1 and e n = v 1 − v n .We will assume throughout that our polygons are equilateral, so that |e i | = 1 for all i; equivalently, e 1 , . . ., e n ∈ S 2 , the unit sphere in R 3 .The space Pol(n) consists of sets of edge vectors in (S 2 ) n which obey the closure condition n i=1 e i = 0.One can show that the set Pol(n) × ⃗ e ∈ (S 2 ) n : n i=1 e i = 0 and for all i j: e i e j is a (2n−3)-dimensional submanifold of (S 2 ) n and that the (2n−3)-dimensional Hausdorff measure of Pol(n) \ Pol(n) × vanishes.In this sense Pol(n) is almost everywhere a submanifold of (S 2 ) n .We may give it the submanifold metric and corresponding volume; it is equivalent to taking the (2n − 3)-dimensional Hausdorff measure on Pol(n) with respect to the metric on (S 2 ) n .
We focus on the quotient space Pol(n) = Pol(n)/ SO(3).This space has a Riemannian metricdefined by the condition that the quotient map Pol(n) → Pol(n) is a Riemannian submersion-and hence a natural probability measure after normalizing the Riemannian volume form.Now we introduce some new coordinates on this space.Connecting the vertices v 3 , . . ., v n−1 to v 1 , as in Figure 1 (far left), produces a collection of n − 3 triangles.The shape of the triangulated surface determined by these triangles (and hence also its boundary, which is the n-gon) is completely determined by the lengths d i of the diagonals joining v 1 and v i+2 and the dihedral angles between triangles meeting at each diagonal.Hence, we can reconstruct the surface (and hence the polygon) up to orientation from the data d 1 , . . ., d n−3 , θ 1 , . . ., θ n−3 , and so these give a system of coordinates for Pol(n).
Indeed, as we have shown [9], these coordinates are natural from the symplectic geometry point of view: in that context, they are called action-angle coordinates.Note that, while the dihedral angles can be chosen completely independently, the diagonal lengths cannot: they must obey the system of triangle inequalities ) n−3 is the (n − 3)-dimensional torus realized as the product of unit circles, then the action-angle coordinates are defined on P n−3 × T n−3 , and we have previously shown that the standard probability measure on this space-that is, the one coming from the product of Lebesgue measure on P n and the standard product measure on T n−3 -is measure-theoretically equivalent to Pol(n): Theorem 1 (Cantarella-Shonkwiler [9]).The reconstruction map α : P n × T n−3 → Pol(n) defining action-angle coordinates (i.e., the procedure illustrated in Figure 1) is measure-preserving.
Therefore, to sample points in Pol(n) (that is, equilateral n-gons), it suffices to sample ⃗ d from Lebesgue measure on P n and ⃗ θ uniformly from T n−3 .Of course, the only challenge is to produce the sample ⃗ d ∈ P n .
In [8], we showed how to do this efficiently.The key observation is that the consecutive differences s i d i+1 − d i lie in the hypercube [−1, 1] n−3 .Therefore, we can generate points in P n by rejection sampling: generate proposed differences (s 0 , . . ., s n−4 ) uniformly from [−1, 1] n−3 , and simply check whether the proposed diagonal lengths (d 1 , . . ., d n−3 ) given by . This is surprisingly efficient: Theorem 2 (Cantarella-Duplantier-Shonkwiler-Uehara [8]).The probability that a random point (s 0 , . . ., In the above and throughought the rest of the paper, we say that g(n) and h(n Since the time it takes to generate points in [−1, 1] n−3 is linear in n, rejection sampling the hypercube yields a valid point in P n in expected time Θ(n 5/2 ).The steps of generating dihedral angles and assembling the n-gon from (d 1 , . . ., d n−3 ) and (θ 1 , . . ., θ n−3 ) do not affect the time bound since they are both linear in n.Therefore, this gives a numerically stable algorithm for generating random equilateral n-gons in expected time Θ(n 5/2 ).We called it the action-angle method (AAM).
We now give an improved version of this method.By fixing d 0 = 1 and generating proposed consecutive differences s i uniformly from [−1, 1] n−3 , we guarantee that the inequalities 0 ≤ d 1 ≤ 2 and −1 ≤ d i+1 − d i ≤ 1 for i = 1, . . ., n − 4 are automatically satisfied.The final inequality 0 ≤ d n−3 ≤ 2 can only be checked at the very end, but the inequalities 1 ≤ d i + d i+1 can be checked one at a time as each s i is generated, and we can abort and start over as soon as one of these inequalities fails.We will show below that the expected number of coordinates that get generated before failure is Θ( √ n), yielding an overall time bound of Θ(n 2 ).Algorithm 1 gives an implementation of our approach, which we call the progressive action-angle method (PAAM).A reference C implementation of this algorithm is included in plCurve [5], where it is now the default algorithm for producing random equilateral polygons in R 3 .

III. COMPLEXITY
Let I(n) be expected number of coordinates s i that are generated in the innermost loop before either failing one of the d i + d i+1 ≥ 1 inequalities (and resetting i to 0) or passing all of them (when i = n − 3).Since we know from Theorem 2 that the overall acceptance probability is ∼ 6 , the expected number of times we will have to re-start the innermost loop is Θ(n 3/2 ).Multiplying these quantities yields Θ(n 3/2 I(n)) for the expected time to produce a valid list of diagonals.The postprocessing steps of generating dihedral angles and assembling the polygon are both linear in n, so they do not affect the overall time bound.Proof.Let p(k) be the probability that the proposed diagonal lengths generated by a random point in ⃗ s ∈ [−1, 1] n−3 satisfy each of the first k inequalities d 0 +d 1 ≥ 1, . . ., d k−1 +d k ≥ 1.In other words, p(k) (thought of as a function of k) is the survival function for the distribution of the index of the first inequality which fails.For ease of notation, we also declare p(0) 1.By a standard integration by parts argument (see, for example, [20, Exercise 1.7.2]), the expected value of a nonnegative random variate is the integral of the survival function, so I(n) = n−3 k=0 p(k).We will now show that for any nonnegative integer k, where sinc(t) = sin(t) t for t 0 and sinc(0) = 1 is the sinc function.For each k, we denote by S k ⊂ [−1, 1] k the subset of points ⃗ s that satisfy the inequalities Here we see both torus actions on the space APol(2).On the left, we may rotate each of the three edges (independently) around the z-axis, sweeping out three cones.Then we identify configurations which are the same under the diagonal circle action rotating the entire configuration around the z-axis (indicated by the dark circle in the xy-plane).On the right, we may rotate the first two edges around the diagonal joining vertices v 1 and v 3 or rotate the entire polygon around the diagonal joining v 1 and v 4 .In each case, this is a Hamiltonian 2-torus action on APol (2).
to ⃗ d given by d i = 1 + i−1 j=0 s j is volume-preserving.Let Q k be the image of S k under this map; i.e., Q k is the polytope of (d 1 , . . ., d k ) satisfying the inequalities 2 k , to prove the above formula for p(k) it suffices to prove that Notice that the defining inequalities for Q k are precisely the inequalities (3) except the last inequality 0 ≤ d k ≤ 2. This suggests that these inequalities may just be the diagonal lengths of an equilateral polygonal path in R 3 which is not required to close up.

APol(k)
(e 1 , . . ., e k+1 ) ∈ (S 2 ) k+1 : where e i = (x i , y i , z i ), so that the defining condition says that the path starts and ends in the xy-plane, and SO(2) acts by simultaneously rotating all edges around the z-axis; this is the diagonal subgroup of the T k+1 = SO(2) k+1 action on (S 2 ) k+1 which rotates edges around the z-axis.APol(k) is the space of abelian polygons introduced by Hausmann and Knutson [21], and we will see that it admits two different effective, Hamiltonian T k actions, as shown in Figure 2.
The first is the residual T k+1 / SO(2) ≃ T k action above.The moment map for this action simply records the z-coordinates of the edges; since the defining equation z 1 + • • • + z k+1 = 0 implies that z k+1 is determined by the remaining z i 's, the last coordinate can be dropped and we can think of the moment map as recording the vector (z 1 , . . ., z k ).Let H k be the image of this map; that is, the moment polytope for this torus action.Of course, −1 ≤ z i ≤ 1 for all i, and, since −1 ≤ z k+1 ≤ 1, we see that the defining inequalities of H k are In other words, H k is the central slab of the hypercube [−1, 1] k of points whose coordinates sum to between −1 and 1.This is a well-studied polytope, and its volume has been known at least since Pólya [34] to be (see also Borwein, Borwein, and Mares [7] for generalizations of the above formula).
On the other hand, we get a T k action on APol(k) which is analogous to the bending flows on Pol(n) described above and in more detail in [9].Specifically, the ith SO(2) factor acts by rotating the first i + 1 edges of the polygonal arm around the ith diagonal, which is the axis through v 1 and v i+2 .Just as in the case of Pol(n), the moment map records the lengths d 1 , . . ., d k of the diagonals, and the image of the moment map is precisely Q k .
But now we've realized APol(k) as a toric symplectic manifold in 2 ways, with moment polytopes H k and Q k .Since the Duistermaat-Heckman theorem [19] implies that the volume of APol(k) must be the product of the volume (2π) k of the torus T k and the volume of either of its moment polytopes, it follows that as desired.We can now estimate the integral using classical methods [23, pp. 172-173] (see also [25]) to get p(k) = 6 πk + O 1 k 3/2 and the result follows.□ Therefore, Θ(n 3/2 I(n)) = Θ(n 2 ) and we have proved a sharp complexity bound on the progressive action-angle method: Theorem 4. The progressive action-angle method generates uniform random samples of closed, equilateral n-gons in expected time Θ(n 2 ).

IV. EXPERIMENTS
As in [38], we classified each random polygon P as knotted or unknotted based on the three invariants ∆ 2 (P), ∆ 3 (P), and ∆ 4 (P), which are absolute values of the Alexander polynomial at roots of unity.The ∆ i are determinants of principal minors of an m × m Alexander matrix A(P) [1], where m is the number of crossings in a projection of P to a plane.The matrix A(P) is usually large, since ) from n = 16 to 8192 using the action-angle method (AAM) and the progressive action-angle method (PAAM).The time needed to generate samples scales as predicted by Theorem 2 and Theorem 4. We compare these to the time required for the computation of all three invariants ∆ 2 , ∆ 3 and ∆ 4 .Invariant computation shows some small-n effects.The most important one is that it is faster to use a dense matrix code to compute the determinants for n < 128.Dense matrix determinants scale as O(n 3 ), so this portion of the graph is steeper.For n < 64 the determinant computation is so fast that (constant time) computational overhead is the dominant factor in the computation, so this portion of the timing curve is flat.AAM is faster than invariant computation for n < 280 while PAAM is faster than invariant computation for n < 1850.
the expected value of m is ∼ 3 16 n log n [15], but sparse, with at most 3 nonzero entries in each row.We computed det(A T A) by using a sparse Cholesky factorization: First we determined a fill-in reducing reordering of the matrix A T A with the Approximate Minimum Degree algorithm [3,4] (routine amd order from the SuiteSparse library).Then we factorized the reordered matrix and extracted the diagonal entries of the factor matrix with our own code [35].Finally we computed the determinant as the product of these diagonal entries.
Our first experiment compared the empirical performance of AAM, PAAM, and invariant computation on an Apple M1 Max laptop (8 cores in parallel).Fitting the data shown in Figure 3 confirmed the expected runtime of O(n 5/2 ) for AAM and O(n 2 ) for PAAM and revealed that invariant computation seems to scale as O(n 1.18 ).
Our second experiment recomputed unknot probabilities using our invariant code.Since unknotted n-edge polygons are very uncommon for large n, we used inverse sampling [24].We sampled polygons with n = round(2 k/7 ) edges for integer k ∈ [28,81], covering the range n = 16 to n = 3043.For each n, we sampled until observing 600 polygons classified as unknots, computing between 608 and ∼ 10 8 samples to do so.The computations took a total of 58 hours on an Apple M1 Ultra personal computer with 16 CPU cores.Using Bennett's approximation for the negative binomial distribution [6], we computed 95% confidence intervals for our unknot probabilities of relative size about 8% of the probability (see Table I).
We then fit the data to models of the form (1) and (2).As shown in Figure 4, model (1) provides an excellent (R 2 > 0.9999) fit to the data, with best-fit parameters C ≃ 3.63±0.09,β ≃ 3.8±0.4,and γ ≃ 7.2±1.7.Xiong et al. measured C ≃ 3.67, β ≃ 3.8 and γ ≃ 8.3 for this model [38], all of which are within our confidence intervals, so we replicate their finding that this model explains the data.(We note that the confidence intervals for the parameters derive from our sampling uncertainty.)On the other hand, the same authors report a poor fit to model (2).In contrast, we find that setting N ≃ 251.3 ± 0.8, β ≃ −0.51 ± 0.35, and γ ≃ −1.2 ± 1.9 in model (2) fits our data equally well.We conclude that more will be needed to distinguish between (1) and (2).At left, we show a log plot of the probability that a given n-edge equilateral polygon has ∆ 2 , ∆ 3 , and ∆ 4 all equal to 1 together with 95% confidence intervals.Following [38], we use this as a proxy for the unknot probability.We also show our best fits to (1) and (2).At right, we show the relative fit errors together with the measurement uncertainty (dotted lines).We see that both models fit the data extremely well.

2 FIG. 1 .
FIG. 1. Constructing an equilateral pentagon from diagonals and dihedrals.The far left shows the fan triangulation of an abstract pentagon.Given diagonal lengths d 1 and d 2 of the pentagon which obey the triangle inequalities, build the three triangles in the triangulation from their side lengths (middle left).Given dihedral angles θ 1 and θ 2 , embed these triangles as a piecewise-linear surface in space (middle right).The far right shows the final space polygon, which is the (solid) boundary of this triangulated surface.

5 FIG. 3 .
FIG.3.This plot compares time per million samples to n for random equilateral n-gons with n = round(2 k/7 ) from n = 16 to 8192 using the action-angle method (AAM) and the progressive action-angle method (PAAM).The time needed to generate samples scales as predicted by Theorem 2 and Theorem 4. We compare these to the time required for the computation of all three invariants ∆ 2 , ∆ 3 and ∆ 4 .Invariant computation shows some small-n effects.The most important one is that it is faster to use a dense matrix code to compute the determinants for n < 128.Dense matrix determinants scale as O(n 3 ), so this portion of the graph is steeper.For n < 64 the determinant computation is so fast that (constant time) computational overhead is the dominant factor in the computation, so this portion of the timing curve is flat.AAM is faster than invariant computation for n < 280 while PAAM is faster than invariant computation for n < 1850.
FIG.4.At left, we show a log plot of the probability that a given n-edge equilateral polygon has ∆ 2 , ∆ 3 , and ∆ 4 all equal to 1 together with 95% confidence intervals.Following[38], we use this as a proxy for the unknot probability.We also show our best fits to (1) and(2).At right, we show the relative fit errors together with the measurement uncertainty (dotted lines).We see that both models fit the data extremely well.