Efficient Generation of Large Random Polygons

Clayton Shonkwiler

Colorado State University

shonkwiler.org

/crest25

this talk!

CREST Mini-Workshop

March 6, 2025

Collaborators

Funding

National Science Foundation (DMS–2107700), Simons Foundation (#709150)

Japan Science and Technology Agency CREST: Grant Number JPMJCR19T4

Henrik Schumacher

University of Georgia

Kandin Theis

Colorado State University

Jason Cantarella

University of Georgia

Key Point

Symplectic geometry, clever programming, and some combinatorics give surprisingly fast algorithms for sampling random polygons.

Random Knot Model: Equilateral Polygons

\(\operatorname{Pol}(n)\) is the space of equilateral \(n\)-gons in \(\mathbb{R}^3\); it consists of (congruence classes of) piecewise-linear maps \(S^1 \to \mathbb{R}^3\) with \(n\) unit-length segments.

Equilateral Polygons

\(\operatorname{Pol}(n)\) can be constructed as a symplectic reduction (see Kapovich–Millson and Hausmann–Knutson):

\operatorname{Pol}(n)=(S^2)^n/\!/_{\vec{0}}SO(3)

Continuous symmetry \(\Rightarrow\) conserved quantity

\(n-3\) commuting symmetries

Rotations around \(n-3\) chords \(d_i\) by \(n-3\) angles \(\theta_i\) commute.

More precisely, \(\operatorname{Pol}(n)\) is (almost) toric, and the \(d_i\) and \(\theta_i\) are action-angle coordinates.

Chord distributions

Theorem [with Cantarella]

The joint distribution of \(d_1,\ldots , d_{n-3}\) and \(\theta_1, \ldots , \theta_{n-3}\) are all uniform on their domains.

Therefore, sampling \(\operatorname{Pol}(n)\) is equivalent to sampling random points in the convex polytope of \(d_i\)’s and random angles \(\theta_i\).

A polytope

The \((n-3)\)-dimensional moment polytope \(\mathcal{P}_n \subset \mathbb{R}^{n-3}\) is defined by the triangle inequalities

0 \leq d_i \leq 2

1 \leq d_i + d_{i-1}

|d_i - d_{i-1}| \leq 1

0 \leq d_{n-3} \leq 2

Reconstruction

A sampling algorithm

Theorem [with Cantarella, Duplantier, Uehara, 2016]

A direct sampling algorithm (the Action-Angle Method) for equilateral \(n\)-gons with expected performance \(O(n^{5/2})\) per sample.

If we let \(s_i = d_i - d_{i-1}\) for \(i=1,\ldots , n-2\) and \(s_i \in [-1,1]\), then we satisfy the inequalities \(|d_i - d_{i-1}|\leq 1\).

Proposition [with Cantarella, Duplantier, Uehara]

If we build \(d_i\) from \(s_i\) sampled uniformly from the hypercube \([-1,1]^{n-3}\), the \(d_i\) obey the triangle inequalities with probability asymptotic to

6 \sqrt{\frac{6}{\pi}} n^{-3/2}.

Speeding This Up

Idea

Don’t generate all the \(s_i\) and then check inequalities at the end: progressively check inequalities as each \(s_i\) is generated.

Proposition [with Cantarella and Schumacher, 2024]

The expected number of \(s_i\) you generate before either failing an inequality (and restarting) or succeeding is \(\sim \sqrt{\frac{24}{\pi}n}\).

Proof sketch.

The probability \(p(k)\) that \(d_0 + d_1 \geq 1, \dots , d_{k-1} + d_k \geq 1\) is

\(p(k) = \frac{2}{\pi} \int_0^\infty \operatorname{sinc}^{k+1} t \, dt.\)

Expected number of inequalities checked is \(p(0) + \dots + p(n-3)\).

Proving this formula for \(p(k)\) uses surprisingly fancy math.

Theorem [with Cantarella and Schumacher, 2024]

A direct sampling algorithm (the Progressive Action-Angle Method) for equilateral \(n\)-gons with expected performance \(O(n^2)\) per sample.

But is it Fast in Practice?

Time to generate samples on an Apple M1 Max laptop (8 parallel cores)

Do It Yourself

git clone --recurse-submodules git@github.com:HenrikSchumacher/CoBarSLink.git

In terminal:

Then in Mathematica:

Get[FileNameJoin[{$HomeDirectory,"CobarSLink","CobarSLink.m"}]];

ActionAngleSample[17, 1]["VertexPositions"]

Probability of Knot Type \(K\)

Xiong, Taylor, Dennis, and Whittington (2021) propose the model

\(P_K(n) \simeq C_K n^{\nu_0 + n_p(K)} \exp(-n/n_0) [1+\beta_K n^{-1/2} + \gamma_K n^{-1}]\),

where \(n_p(K)\) is the number of prime summands in \(K\),

\(\nu_0 =-0.19 \pm 0.001\) and \(n_0 = 259.3 \pm 0.2\) are characteristic of the random polygon model, and \(\beta_K, \gamma_K\) are finite-size corrections.

Deguchi and Uehara (2020) propose

\(P_K(n) = \hat{C}_K(n/\hat{n}_K)^{\hat{\nu}_0 + n_p(K)}\exp(-n/\hat{n}_K)\exp(\hat{\beta}_K/n)\)

Unknot Probabilities

We generated 600 unknotted \(n\)-gons for \(n\) up to 3043, which took 58 hours on an Apple M1 Ultra desktop (16 parallel cores)

We got an excellent fit with \(C = 3.63 \pm 0.09\), \(\beta = 3.8 \pm 0.4\), \(\gamma = 7.2 \pm 1.7\).

Xiong et al. measured \(C \simeq 3.67\), \(\beta \simeq 3.8\), \(\gamma \simeq 8.3\).

Competing Models

Xiong et al. claim to rule out the “pure exponential” model

\(P_{0_1}(n) \simeq \exp(-n/n_0)[1+\beta n^{-1/2} + \gamma n^{-1}]\)

A Xiong et al 2021 J. Phys. A: Math. Theor. 54 405001

We Can’t Rule Out the Pure Exponential

We found that, with \(n_0 = 251.3 \pm 0.8\), \(\beta = -0.51 \pm 0.35\), \(\gamma = -1.2 \pm 1.9\), the pure exponential fits about as well as the power law!

Confined Polygons

Gelbart & Knobler Science 323 (2009) 1682–1683

If we want to sample polygons in rooted, spherical confinement of radius \(R\), then we simply add the constraints \(d_i \leq R\) for all \(i\).

Sampling

Algorithm [with Cantarella, 2016]

A uniformly ergodic Markov chain for sampling random equilateral \(n\)-gons in rooted spherical confinement (implemented in plCurve).

Algorithm [Diao–Ernst–Montemayor–Ziegler, 2012]

An algorithm for directly sampling successive marginals of Lebesgue measure on \(\mathcal{P}_n(R)\).

“The data shows that the runtime grows slower than an exponential function, but faster than a power function.”

Tightest Confinement

When \(R=1\), the inequalities reduce to

\(0\leq d_i \leq 1 \qquad 1 \leq d_i + d_{i+1}\)

Define \(\phi: (d_1, d_2, \dots ) \mapsto (d_1, 1-d_2, d_3, 1-d_4, \dots )\), which sends \(\mathcal{P}_n(1)\) to

which define a polytope \(\mathcal{P}_n(1) \subset [0,1]^{n-3}\).

\(\mathcal{O_{n-3}} := \{(s_1, \dots , s_{n-3}) \in [0,1]^{n-3} : s_1 \geq s_2 \leq s_3 \geq s_4 \leq \dots\}\),

The order polytope of the zig-zag (fence) poset.

Each possible total ordering \(s_{\sigma(i_1)} \leq \dots \leq s_{\sigma(i_{n-3})}\) corresponds to a linear extension of the zig-zag poset or, equivalently, an alternating permutation.

Each alternating permutation \(\sigma\) determines an orthoscheme \(\Delta_\sigma\) in a unimodular triangulation of \(\mathcal{O}_{n-3}\).