Quasipolynomial-time algorithms for Gibbs point processes

We demonstrate a quasipolynomial-time deterministic approximation algorithm for the partition function of a Gibbs point process interacting via a finite-range stable potential. This result holds for all activities $\lambda$ for which the partition function satisfies a zero-free assumption in a neighborhood of the interval $[0,\lambda]$. As a corollary, for all finite-range stable potentials we obtain a quasipolynomial-time determinsitic algorithm for all $\lambda</(e^{B + 1} \hat C_\phi)$ where $\hat C_\phi$ is a temperedness parameter and $B$ is the stability constant of $\phi$. In the special case of a repulsive potential such as the hard-sphere gas we improve the range of activity by a factor of at least $e^2$ and obtain a quasipolynomial-time deterministic approximation algorithm for all $\lambda<e/\Delta_\phi$, where $\Delta_\phi$ is the potential-weighted connective constant of the potential $\phi$. Our algorithm approximates coefficients of the cluster expansion of the partition function and uses the interpolation method of Barvinok to extend this approximation throughout the zero-free region.


Introduction
Gibbs point processes are a fundamental model of random spatial phenomena in the continuum.Most classically, such processes are used to model a gas under local interactions (see Ruelle's [Rue99] text).Beyond that, Gibbs point processes are used to model various physical phenomenon that often exhibit local repulsion, such as the locations of galaxies in the universe, the time and place of earthquakes, and the growth of trees in a forest; see [MW07,DVJ07] for these applications and more.A simple and well-studied example of a Gibbs point process is the hard-sphere model, where one samples a Poisson point process in a set of finite volume in R d and conditions on no two points having distance less than some parameter r > 0. In order to better understand these models one often wants to approximately compute the partition function of the model, which may be understood as a weighted count of allowable configurations of points.The partition function grows exponentially in the volume, making exact computation intractable even for basic examples.Additionally, the rate of exponential growth is equal to the infinite-volume pressure, a central quantity in statistical physics.
Approximating the partition function and sampling-either approximately or exactly-are the two main algorithmic problems associated to Gibbs point processes.Under very mild assumptions, polynomial-time approximate sampling of the point process can be used to provide a randomized approximation to the partition function.Many techniques have been applied to Gibbs point processes for approximate and exact sampling in certain regimes; in fact, the seminal Metropolis-Hastings algorithm was developed to sample from the hard-sphere model in dimension 2 [MRR + 53].
On the other hand, deterministic approximation algorithms for partition functions of Gibbs point processes are less well-understood.For Gibbs point processes, to our knowledge the only rigorous result giving a deterministic algorithm is that of Friedrich, Göbel, Katzmann, Krejca and Pappik [FGK + 22], which shows that for the special case of hard-spheres, one may approximate the partition function in quasipolynomial time for a certain range of parameters.The authors show that the hard-sphere model can be well approximated by its discrete analogue (the hard core model) allowing them to apply known algorithmic results from the discrete setting (see Section 1.3 for more detail).
In this paper, we provide a quasipolynomial time deterministic approximation algorithm for the partition function for a general class of stable Gibbs point processes.Our main result is stated only under the assumption of zero-freeness of the partition function; from there, we deduce two main corollaries using existing zero-freeness results from the literature, one which applies for all stable potentials and a stronger result that applies for the more restricted class of repulsive potentials.We defer formal statements to Section 1.2.
Our approach is via Barvinok's interpolation method [Bar16] combined with use of the cluster expansion for Gibbs point processes.This allows us to work with the Gibbs process directly, rather than a discrete approximation of the process.By combining our result with the zero-free region for stable potentials guaranteed by the cluster expansion, we obtain what appears to be the first algorithmic result for stable, non-repulsive potentials.In the special case of a repulsive potential, combining our result with previous work of the second author and Perkins on zero-freeness [MP21] yields the first quasipolynomial time deterministic approximation algorithm for a large class of repulsive potentials (which includes the hard-sphere model) and range of parameters.For the hardsphere potential φ, we note that [MP21,Lemma 12] gives an explicit bound of ∆ φ < C φ (1 − 8 −d−1 ) where ∆ φ , C φ denote the potential-weighted connective constant and temperedness constant of φ respectively (defined in the next sections).This demonstrates that our algorithm works for a wider range of parameters than the previous deterministic algorithms of [FGK + 22].Additionally, [FGKP21] argues that the connective constant of the discretization used in [FGKP21, FGK + 22] would not provide an improvement to their results; as such, working in the continuum and using the zero-freeness result of [MP21] gets around this issue.
1.1.Formal definition of the model.The point processes we consider are defined by three parameters: • a measurable set S ⊂ R d of finite volume, • a parameter λ 0 referred to as the activity or fugacity, The temperedness constant of a potential is defined as (1) The temperedness constant may be understood as a measure of the strength of the potential.We say that a potential φ is tempered if C φ < ∞, and always assume that φ is tempered.The energy of a configuration of points {x 1 , . . ., x N } ⊂ R d is defined by We will always assume that φ is stable, meaning that there is a constant B 0 so that for all N and x 1 , . . ., x N we have The infimum over such B is called the stability constant of φ.The assumption of stability is used to show that the partition function-and thus the Gibbs point process-itself is well-defined.Under certain conditions on φ, the assumption of stability is a necessary condition for the point process to be well-defined (see [Rue99, Section 3.2] for a detailed discussion and many examples).A potential 0 for all x.In particular, repulsive potentials are stable with stability constant B = 0.
The Gibbs point process in S with potential φ at activity λ is the point process in S whose density against the Poisson point process of activity λ is proportional to e −H(x 1 ,...,x N ) .The grand canonical partition function at activity λ is defined by Throughout, we work with the case of S = Λ n := [−n, n] d ⊂ R d , i.e. the axis parallel box of side-length 2n.One of the most studied examples of a Gibbs point process is the hard sphere model which is defined by the potential for fixed r > 0. The hard sphere model is supported on configurations {x 1 , . . ., x N } such that H(x 1 , . . ., x N ) = 0 i.e. configurations consisting of the centers of a packing of spheres of radius r/2.Another similar example is the Strauss potential, in which the +∞ in the definition of the hard sphere potential is replaced with a parameter a > 0.
Among the most common examples of a stable potential that is not repulsive is a Lennard-Jones potential (see [Rue99, Section 3.2.10]).While there are many examples of potentials that are called Lennard-Jones potentials, they are characterized by being strongly repulsive at short distances and weakly attractive at far distances.A large family of widely used potentials is of the form for parameters A, B > 0 and α > d, where we recall d is the underlying dimension.One may also truncate this potential to be of finite range by simply setting it equal to 0 if x 2 T for some parameter T , and this still yields a stable potential.
1.2.Statement of results.Our results will require only two additional assumptions on the potential: first a basic assumption on its form so that we may approximately compute certain volumes; and second a zero-freeness assumption.The zero-freeness assumption we use for stable (nonrepulsive) potentials is a classical result that follows from the cluster expansion [Rue99, PY17], while in the repulsive case we will use the work of the second author and Perkins [MP21].We begin with the assumption required for computational purposes.
Assumption 1.1.There are compact centrally symmetric convex sets We work in a real-valued model of computation, and assume unit cost for elementary operations; we also assume that evaluation of e −φ as well as checking if a point lies in a given K j have unit cost.Assumption 1.1 essentially has two elements: it assumes that the potential φ has a certain amount of (piecewise) regularity, and also assumes that φ is of finite range.We are not aware of any natural potential φ in the literature that fails to satisfy the piecewise regularity portion of the assumption.There are however many natural potentials that are not of finite range.
We remark that the hard sphere potential (5), the Strauss potential, and the truncated Lennard-Jones potential are easily seen to satisfy Assumption 1.1.
A note on asymptotic notation.Throughout the paper we think of the parameters d, ℓ, L and R from Assumption 1.1 as fixed and allow the implicit constants in our asymptotic notation to depend on these parameters.
Our main theorem is that we obtain deterministic approximation algorithms for Z S (λ) for potentials satisfying Assumption 1.1 and the following zero-freeness assumption on Z S (λ).
Assumption 1.2.We say a stable potential φ satisfies Assumption 1.2 at λ 0 > 0 if there exist constants δ, C > 0 so that the following holds.For all bounded, measurable S ⊂ R d we have where Assumption 1.2 may be understood as saying that the Gibbs point process exhibits no phase transition in the regime [0, λ 0 ] in the Lee-Yang sense (see Section 1.3).
Given a number Z, an ε-approximation to Z is a value Ẑ so that e −ε Ẑ Z e ε Ẑ.Our main theorem asserts a quasipolynomial-time approximation for the partition function under these two assumptions.Recall that Λ Theorem 1.3.Let φ be a stable pair potential that satisfies the regularity Assumption 1.1.Suppose φ satisfies the zero-freeness Assumption 1.2 at λ 0. Let ε ∈ (0, 1), n ∈ N. Then there is a deterministic algorithm to compute an ε-approximation to Z Λn (λ) with runtime at most exp(O(log We note that the implicit constant in the exponent depends on the potential φ as well as the activity λ. We will deduce two main corollaries from Theorem 1.3.The first of which deduces an algorithm for all stable potentials in the regime of cluster expansion convergence.The fact that for λ < (e 1+2B C φ ) −1 a stable potential φ satisfies Assumption 1.2 is a direct consequence of the cluster expansion (see, e.g., (11) or [Rue99, Thm 4.2.3]).
This classical result was recently improved by Procacci and Yuhjtman [PY17] who showed that any stable tempered potential φ satisfies Assumption 1.2 at all λ < (e 1+B Ĉφ ) −1 , where We note that Ĉφ C φ for all φ with equality if φ is repulsive.
Our second corollary applies for repulsive potentials-i.e. the case of B = 0-for a much wider range of parameters.To apply Theorem 1.3 in the repulsive case, we will use pre-existing work on zero-freeness of Gibbs point process partition functions.In [MP21], the potential-weighted connective constant ∆ φ was introduced, which captures a relationship between the strength of the potential and the geometry of the underlying space.
First, define V k via where we write dv = dv 1 dv 2 . . .dv k and interpret v 0 = 0, and in the case of j = 1 interpret the empty sum as equal to 0. Since the potential φ is repulsive, the sequence {V k } k 1 is submultiplicative and so we may define the potential-weighted connective constant ∆ φ via For any repulsive potential φ we have ∆ φ C φ where C φ is the temperedness constant defined at (1).So long as φ is non-trivial we in fact have ∆ φ < C φ .In [MP21] it was shown that Assumption 1.2 is satisfied for each tempered repulsive potential for all λ 0 < e/∆ φ .We therefore obtain the following immediate corollary.

Context and related work. Much of the classical work on
Gibbs point processes has consisted of showing the absence of a phase transition when λ is small.There are many definitions and notions of a phase transition.One of the most robust definitions is due to Lee and Yang [LY52,YL52], which states that a phase transition is a point at which the pressure fails to be analytic.Often, various other definitions of phase transitions, e.g. in terms of infinite-volume Gibbs measures, can be shown to coincide with the Lee-Yang definition (see, e.g., [DS85] for some rigorous equivalences in the discrete case).
It remains a major open problem to demonstrate-or rule out-a phase transition for even a single non-trivial pair potential φ 1 .Groenveld [Gro62] used the cluster expansion to prove that no phase transition occurs for λ < 1/(eC φ ) for repulsive potentials; this was extended to the broader class of stable potentials by Penrose [Pen63] and Ruelle [Rue63].For repulsive potentials, works of Meeron [Mee62,Mee70] extends this regime to 1/C φ .Taking inspiration from Weitz's groundbreaking work [Wei06] on the hard-core model-a discrete repulsive spin system-the work of the second author and Perkins [MP20] proved that there is no phase transition for λ < e/C φ .This was extended further in [MP21] which showed that C φ may be replaced with the potential-weighted connective constant ∆ φ .
In practice, various Markov chain Monte Carlo algorithms are used to sample from Gibbs point processes, both approximately and exactly.There is a large literature dedicated to this study, see for example [MRR + 53,AW57,BK11,PW96,HVLM99,Gar00,Møl01,Hub16,Pre75,Møl89,FGKP21, MP22].We refer to [MP22] and the references therein for more context and background on these 1 The first example of a continuous system for which a phase transition was proven was the Widom-Rowlinson model, a result due to Ruelle [Rue71].However this does not fit into the framework of indistinguishable particles that we consider in this paper.
results; we also note that [MP22] demonstrates a fast mixing approximate sampling algorithm and polynomial-time random approximation algorithm for a repulsive Gibbs point process in the regime λ < e/∆ φ , the same regime in which the results of this work hold.
When it comes to deterministic algorithms, there are several results for discrete systems.In this setting, three approaches for obtaining approximation algorithms have emerged in recent years.The first is due to Weitz [Wei06] who pioneered an approach based on a notion of correlation decay (strong spatial mixing) related to the absence of phase transitions.The second is the interpolation method introduced by Barvinok [Bar16]-a framework for showing that under a zero-freeness assumption, one can approximate the logarithm of a polynomial using a small number of Taylor coefficients (see also [PR17] for an important extension of this method).The third is a method based on the cluster expansion in statistical physics pioneered by Helmuth, Perkins and Regts [HPR20] that is closely related to Barvinok's method.
For Gibbs point processes, to our knowledge the only rigorous result giving a deterministic algorithm is that of Friedrich, Göbel, Katzmann, Krejca and Pappik [FGK + 22].They show that for the hard-spheres model, one may approximately compute the partition function in quasipolynomial time for λ < e/C φ .The approach of [FGK + 22] works by showing that one may approximate the partition function of the hard-sphere model with the partition function for the hard-core model on a graph given by discretizing Euclidean space with a small mesh.After this approximation is in place, an application of Weitz's method [Wei06] provides an algorithm.We note also that [FGK + 22] applies not only to the hard-sphere model, but also to a class of multi-type Gibbs point processes with hard constraints, an example of this more general class being the Widom-Rowlinson model.Additionally, while both the algorithm of [FGK + 22] and Theorem 1.3 have quasipolynomial runtime, the algorithm of [FGK + 22] runs in time exp(O(log 2 (|Λ n |/ε)) rather than our runtime of exp(O(log 3 (|Λ n |/ε)).
To our knowledge, Theorem 1.3 marks the first approximation algorithm of any kind for the partition function for stable, non-repulsive potentials.
1.4.Our approach.Our approach is via Barvinok's interpolation method combined with use of the cluster expansion.The cluster expansion-also called the Meyer series-is a combinatorial description of the Taylor coefficients of log Z S (λ).In particular for bounded, measurable S ⊂ R d and |λ| < (e 1+B Ĉφ ) −1 we have where G k is the set of connected labeled graphs with k vertices and E(G) is the set of edges in a graph G (see, e.g., [PY17] for this and similar expansions).The algorithmic significance of the cluster expansion stems from two observations.The first is that in order to compute 1 |S| log Z S (λ) up to an additive error of ε, it suffices to compute the first O(log(|S|/ε)) coefficients C k (S).This simple but crucial observation is at the heart of Barvinok's interpolation method and algorithmic applications of the cluster expansion.The second is that the description of the Taylor coefficients as a sum over connected graphs allows us to approximate them efficiently.
At first sight, it seems that the use of the cluster expansion restricts the range of λ for which we can obtain algorithms to the interval (0, (e 1+B Ĉφ ) −1 ) where the cluster expansion is known to converge.Moreover, the work of Groeneveld [Gro62] and Penrose [Pen63] implies that for a repulsive potential the radius of convergence of the cluster expansion is at most 1/C φ (see Remark 3.7 in [NF20]).In order to use the cluster expansion to produce an algorithm throughout the zero-free region (which for repulsive potentials is known to include values of λ outside of the radius of convergence), we use an idea of Barvinok and apply a well-chosen polynomial to map between the zero-free region and the disk.This is handled in Section 3.
The main technical contribution of this work is to approximate the coefficients of the cluster expansion (11) in quasipolynomial time.We note that the problem of approximating coefficients of the cluster expansion was pointed out in [FGK + 22] as a central obstacle to obtaining an efficient deterministic algorithm that works entirely in the continuum.To approximate the coefficients C k (S), the main challenge is that for each graph G ∈ G k , the integrand in (11) need not be wellbehaved; in particular, it need not be Lipschitz.To handle this, we break each such term into many subsequent terms, each of which will give us an integral of a function that is Lipschitz over its support.We then approximate each of these integrals by taking a sufficiently fine mesh and comparing the integral to a weighted sum over this mesh.Our assumption that the potential φ has finite range allows us to restrict our attention to approximating integrals over regions whose volume is independent of n.It would be interesting to extend the results of this paper to include stable potentials of infinite range.

Approximating coefficients in the cluster expansion
Throughout this section, we fix a stable potential φ with with stability constant B that satisfies Assumption 1.1.Recall that Λ n := [−n, n] d ⊂ R d .Our goal in this section is to provide an approximation algorithm for C k (Λ n ) (as defined at (11)).In the next section we show how to use these approximate coefficients to arrive at an approximation of the partition function Z Λn (λ).
Proposition 2.1.Let φ be a stable potential satisfying Assumption 1.1.There are constants C, c > 0 depending only on φ and the dimension d so that the following holds.For each ε ∈ (0, 1) and k, n ∈ N we may approximate C k (Λ n )/|Λ n | up to an additive error of ε in time at most C|Λ n |ε −dk e ck 3 .Note that Proposition 2.1 makes no use of a zero-freeness or correlation decay type assumption.In order to prove Proposition 2.1, we will approximate each summand in (11).Further, we will break up the term for each graph into the terms that are Lipschitz on their support.The main work of this section is Lemma 2.3 below, which provides an approximation for each of these Lipschitz terms.Before stating this precisely, we make a few definitions.Definition 2.2.Let ∆ 1 , . . ., ∆ ℓ be defined as in Assumption 1.1.
(2) For each labelling σ and edge {i, j} ∈ E(G), define the function and note that supp(f σ ) ⊆ A σ .
Recall that G k is the set of connected labeled graphs with k vertices.Given a graph G ∈ G k we have where the sum ranges over all edge-labellings σ of G.We note that there are ℓ |E(G)| ℓ k 2 /2 such labellings, and there are at most 2 k 2 /2 connected graphs on k vertices.With this setup in mind, in order to approximate C k (Λ n ), it is sufficient to approximate the integral of f σ for each labelling σ of each graph G ∈ G k .In the next subsection, we will demonstrate such an approximation: Further, the set S σ,δ may be computed in time Proposition 2.1 follows quickly from here.
Proof of Proposition 2.1.Note that by (12) Taking δ = ε exp(−Ck 2 ) in Lemma 2.3 with C sufficiently large as a function of R, d, ℓ and B we may approximate ) where c is a constant depending only on R, d and ℓ and B.
2.1.Approximating the integral over a labelling.In this section, we prove Lemma 2.3.Fix G ∈ G k and set S = Λ n .Recall from Assumption 1.1 that the sets ∆ j are defined via ∆ j = K j \ K j−1 .For a given γ 0, define ∆ (γ) Before continuing, let us motivate the definition of U γ .We will show that if x ∈ U γ , then f σ is Lipschitz on a small neighbourhood of x (see Lemmas 2.8 and 2.7 below).This will allow us to approximate the integral Uγ f σ (x) dx by a sum over a mesh.Indeed, we define S σ,δ := U γ ∩ ((δZ) d ) k , and prove the following lemma.
To complete the proof of Lemma 2.3 we show that Uγ f σ (x) dx closely approximates Λ k n f σ (x) dx and that the sum in Lemma 2.4 can be computed efficiently.We turn our attention first to the latter task.
Lemma 2.5.The set S σ,δ may be enumerated in time Proof.Fix a spanning tree T of G and let ) and that T may be enumerated in time O(|S|(2R) d(k−1) δ −dk ) as well.Further, for each point in T , we may check membership in S σ,δ in time O(k 2 ), completing the lemma.
To show that Uγ f σ (x) dx closely approximates S k f σ (x) dx it will be enough to show that the measure of To see (13), fix a spanning tree T of G so that {a, b} ∈ E(T ).We then have that This provides a bound of σ({a,b}) | .
Write j = σ({a, b}) for notational simplicity.By convexity of the sets K j and K j−1 we have (1 − γ)K j ⊂ K j and K j−1 ⊂ (1 + γ)K j−1 .We may therefore bound Combining the two previous displayed equations shows (13) and completes the lemma.Assuming Lemma 2.4, the proof of Lemma 2.3 now follows quickly.
Proof of Lemma 2.3.
Since |e −φ(x) − 1| e 2B , we have that |f σ | e Bk 2 .The first term on the RHS of (15) may thus be bounded by e Bk 2 |W γ | which we bound using Lemma 2.6, and the latter term may be bounded by Lemma 2.4.It remains to prove Lemma 2.4.As discussed, a key step will be to show that if x ∈ U γ , then f σ is Lipschitz on a small neighbourhood of x.This is carried out in the following two lemmas.
To bound (16) we use the inequality where we used that x → exp(−φ(x)) is L-Lipschitz on each ∆ j by Assumption 1.1.Set Since we are working in the context of Lemma 2.3, we may also assume that γ 1/d by taking δ small enough as a function of d and R.
We let D = {z ∈ C : |z| < 1} denote the open unit disk in C. Our first step is finding a polynomial to map D into the region N γ so that some point in D is mapped to 1.While there are many ways to find such a polynomial (for example, [Bar16, Lemma 2.2.3] gives an explicit construction) the form of the polynomial is not important for us here, and so we state the relevant properties: Lemma 3.2.For each γ > 0 we may find a polynomial Φ = Φ γ so that Φ(D) ⊂ N γ , Φ(0) = 0 and there is a point z 1 ∈ D so that Φ(z 1 ) = 1.
The point now will be to work with the function g : D → C defined by g = f • Φ.Under the assumptions of Theorem 3.1, such a function g is analytic in D and |g(z)| 1 for all z ∈ D. The Cauchy integral formula will imply that we can approximate g by its Taylor polynomial provided we are not near the boundary: Fact 3.3.Let g be analytic in D with |g(z)| 1 for all z ∈ D. Then for all z ∈ D and k ∈ N we have and note that by Cauchy's integral formula we have Summing over j k completes the proof.Let z 1 be as in Lemma 3.2.It follows that in order to approximate f (1) = g(z 1 ) up to an additive error of ε, it is sufficient to approximate the first k = C γ log(1/ε) terms g (j) (0)/j! up to an error of ε/(2k) each.To do so, we will relate the derivatives of g to those of f .Iterating the chain rule gives the classical Faà di Bruno formula: Fact 3.4 (Faà di Bruno's formula).Let Φ and f be analytic at 0 with Φ(0) = 0 and set g = f • Φ.Then for each n ∈ N we have where B n,k (x 1 , x 2 , . . ., x n−k+1 ) are the Bell polynomials given by B n,k (x 1 , x 2 , . . ., where the sum is over all sequences of non-negative integers j 1 , j 2 , . . ., j n−k+1 satisfying To make use of this fact, we will need to evaluate the Bell polynomials B n,k at derivatives of our polynomial Φ.A simple term-by-term bound will be good enough for our purposes.Lemma 3.5.For each polynomial Φ there is a constant C = C Φ > 0 so that the following holds.For all n, k ∈ N we may compute B n,k (Φ ′ (0), Φ ′′ (0), . . ., Φ (n−k+1) (0)) in time O(k C ) and we have the bound k! n! B n,k (Φ ′ (0), Φ ′′ (0), . . ., Φ (n−k+1) (0)) e Ck .
Proof.First compute all derivatives of Φ at 0, which takes O Φ (1) time since Φ is a polynomial.Set x j = Φ (j) (0) and note that x j = 0 for j > d where d is the degree of Φ.We may thus reduce the sum in (25) to those for which j i = 0 for all i > d as all other summands are 0. Since the values j 1 , . . ., j d are non-negative integers of size at most k, we have (k + 1) d choices for the sequence (j 1 , . . ., j d ).Moreover the constraints j i = k and ij i = n imply that n dk and so each summand in (25) can be computed in time O d (k C ′ ) for some C ′ = C ′ d > 0. We may therefore calculate the sum (25) in time O Φ (k C ) for some C = C d > 0.
To see the claimed bound, apply the triangle inequality to see Relaxing the sum to consist of sequences satisfying j i = k, the binomial theorem gives Since g(z 1 ) = f (1), it is thus sufficient to approximate g (j) (0)/j! for j = 0, 1, . . ., k − 1 up to an additive error of at most δ := ε/(2k).By Fact 3.4 and Lemma 3.5 we may expand where the coefficients B i,j satisfy |B i,j | e Ci e Ck = O(ε −C ) and may be computed in time O(k C ) = O(log C (1/ε)).Thus, if we take C ′ large enough (as a function of γ) and approximate f (i) (0)/i! up to additive error ε C ′ , we obtain an approximation to g(z 1 ) = f (1) up to additive error ε.