Rejection and Importance Sampling based Perfect Simulation for Gibbs hard-sphere models

TL;DR

This paper introduces an importance sampling rejection method for perfect simulation of Gibbs hard-sphere models, providing theoretical and numerical comparisons with existing algorithms.

Contribution

It proposes a novel importance sampling based rejection approach for perfect sampling of Gibbs hard-sphere models, analyzing its efficiency and asymptotic complexity.

Findings

The new method outperforms naive rejection algorithms in efficiency.

Asymptotic analysis shows favorable complexity as sphere radius decreases.

Numerical results confirm the theoretical advantages of the proposed approach.

Abstract

Coupling from the past (CFTP) methods have been used to generate perfect samples from finite Gibbs hard-sphere models, an important class of spatial point processes, which is a set of spheres with the centers on a bounded region that are distributed as a homogeneous Poisson point process (PPP) conditioned that spheres do not overlap with each other. We propose an alternative importance sampling based rejection methodology for the perfect sampling of these models. We analyze the asymptotic expected running time complexity of the proposed method when the intensity of the reference PPP increases to infinity while the (expected) sphere radius decreases to zero at varying rates. We further compare the performance of the proposed method analytically and numerically with a naive rejection algorithm and popular dominated CFTP algorithms. Our analysis relies upon identifying large deviations decay rates of the non-overlapping probability of spheres whose centers are distributed as a homogeneous PPP.

Figures

Equations

$$\mathscr{G}_{n}=\big{\{}\mathbf{x}=\{(y_{1},r_{1}),\dots,(y_{n},r_{n})\}:(y_{i},r_{i})\in[0,1]^{d}\times[0,\infty),\,\text{ for }\,i=1,\dots,n\big{\}}.$$

$$\frac{d\mu }{d{\mu }^0}(\mathbf{x})=\frac{I\left(\mathbf{x}\in \mathcal{A}\right)}{\mathcal{P}(\lambda )},\mathbf{x}\in \mathcal{G},$$

$$\mathcal{P}(\lambda )={\mathbb{P}}_{{\mu }^0}\left(\mathbf{X}\in \mathcal{A}\right).$$

$$S(x,r)=\left\{(y_1\mathsf{mod}1,\dots ,y_d\mathsf{mod}1):y=(y_1,\dots ,y_d)\in {\mathbb{R}}^d,\Vert x-y\Vert <r\right\},$$

$${\gamma }^{\prime }=\{\begin{array}{ll}\gamma ,& \text{ if [0,1]d is treated as the torus},\\ \gamma /2^d,& \text{otherwise}.\end{array}$$

$$\underset{\lambda \to \infty }{\lim }\mathcal{P}(\lambda )$$

$$\underset{\lambda \to \infty }{\lim }\left[\frac{1}{{\lambda }^{2-\eta d}}\log \mathcal{P}(\lambda )\right]$$

$$\text{ and }\underset{\lambda \to \infty }{\lim }\left[\frac{1}{\lambda }\log \mathcal{P}(\lambda )\right]$$

$$\underset{\lambda \to \infty }{\lim }\left[\mathcal{P}(\lambda )\exp \left(\frac{\gamma m_1}{2}{\lambda }^{2-\eta d}\right)\right]$$

$$\underset{\lambda \nearrow \infty }{\lim }\frac{\mathsf{VF}(\lambda )}{\gamma {\overline{r}}^d{\lambda }^{1-\eta d}}=1,\text{ if }\eta d>1,$$

$$\underset{\lambda \nearrow \infty }{\lim }\mathsf{VF}(\lambda )={\rho }^{\max }\gamma ,\text{ if }\eta d<1,$$

$$\underset{\lambda \nearrow \infty }{\lim }\mathsf{VF}(\lambda )<{\rho }^{\max }\gamma ,\text{ if }\eta d=1.$$

$${\mathcal{T}}_{\mathsf{NAR}}=\frac{C_{\mathsf{itr}}(\lambda )}{\mathcal{P}(\lambda )}.$$

$$C_{\mathsf{itr}}(\lambda )=\Theta \left({\lambda }^{\min \{\eta d,2\}}\right).$$

$$\mathcal{T}_{\mathsf{NAR}}=\begin{cases}\Theta\left(\lambda^{2}\right),&\quad\text{ if }\,\eta d\geq 2,\\ \Theta\left(\lambda^{\eta d}\,\exp\Big{(}\left(\gamma m_{1}/2+o(1)\right)\lambda^{2-\eta d}\Big{)}\right),&\quad\text{ if }\,1<\eta d<2,\\ \Theta\Big{(}\lambda^{\eta d}\,\exp\left(\delta\lambda\right)\Big{)},\text{ for some }\,0<\delta\leq 1,&\quad\text{ if }\,\eta d=1,\\ \Theta\left(\lambda^{\eta d}\,\exp\Big{(}(1+o(1))\lambda\Big{)}\right),&\quad\text{ if }\,0<\eta d<1.\end{cases}$$

$$L_{n,k}({\mathbf{x}}_n)\le {\sigma }_{n,k}\le 1,\text{ if }{\mathbf{x}}_n\in D_{n,k}\cap \mathcal{A},$$

$$\mu (\mathcal{B})$$

$$\displaystyle=\sum_{n\in\mathbb{N}_{0}}e^{-\lambda}\frac{\lambda^{n}}{n!}\Bigg{(}\sum_{k=1}^{K}\,\mathbb{E}_{\mu_{n,k}}\big{[}I\left(\mathbf{X}_{n}\in D_{n,k}\cap\mathscr{B}\cap\mathscr{A}\right)L_{n,k}(\mathbf{X}_{n})\big{]}\Bigg{)}$$

$$\displaystyle=\sum_{n\in\mathbb{N}_{0}}e^{-\lambda}\frac{\lambda^{n}\widetilde{\sigma}(n)}{n!}\Bigg{(}\sum_{k=1}^{K}\frac{\sigma_{n,k}}{\widetilde{\sigma}(n)}\,\mathbb{E}_{\mu_{n,k}}\Bigg{[}\frac{I(\mathbf{X}_{n}\in D_{n,k}\cap\mathscr{B}\cap\mathscr{A})L_{n,k}(\mathbf{X}_{n})}{\sigma_{n,k}}\Bigg{]}\Bigg{)}$$

$$\displaystyle=\sum_{n\in\mathbb{N}_{0}}\frac{\lambda^{n}\widetilde{\sigma}(n)}{n!}\Bigg{(}\sum_{k=1}^{K}\frac{\sigma_{n,k}}{\widetilde{\sigma}(n)}\,\mathbb{P}_{\mu_{n,k}}\big{(}J=1,\mathbf{X}_{n}\in D_{n,k}\cap\mathscr{B}\cap\mathscr{A}\big{)}\Bigg{)},$$

$$\mathbb{P}\left(M=m\right)=\frac{1}{C_{\lambda }}\frac{{\lambda }^m\tilde{\sigma }(m)}{m!},m\in {\mathbb{N}}_0,$$

$$P_{\mathsf{acc}}(\lambda )=\frac{\mathcal{P}(\lambda )}{\mathbb{E}[\tilde{\sigma }(N)]}.$$

$$\displaystyle P_{\mathsf{acc}}(\lambda)=\frac{1}{C_{\lambda}}\sum_{n\in\mathbb{N}_{0}}\frac{\lambda^{n}\widetilde{\sigma}(n)}{n!}\Bigg{(}\sum_{k=1}^{K}\frac{\sigma_{n,k}}{\widetilde{\sigma}(n)}\,\mathbb{E}_{\mu_{n,k}}\Bigg{[}\frac{L_{n,k}(\mathbf{X}_{n})}{\sigma_{n,k}};\mathbf{X}_{n}\in D_{n,k}\cap\mathscr{A}\Bigg{]}\Bigg{)}.$$

$${\mathcal{T}}_{\mathsf{ISAR}}=\frac{{\tilde{C}}_{\mathsf{itr}}(\lambda )\mathbb{E}[\tilde{\sigma }(N)]}{\mathcal{P}(\lambda )},$$

$$J_2\sim \mathsf{Bern}\left(\frac{L_{M,J_1}(\mathbf{X})\exp \left(-\beta V(\mathbf{X})\right)I\left(\mathbf{X}\in D_{M,J_1}\right)}{{\sigma }_{M,J_1}}\right).$$

$$\displaystyle\widetilde{L}_{n}(\mathbf{x}_{n})=\frac{d\mu^{0}}{d\widetilde{\mu}_{n}}(\mathbf{x}_{n})=\prod_{i=1}^{n}\Big{(}1-B_{i}\Big{)},$$

$$B_i\ge \frac{{\gamma }^{\prime }}{{\lambda }^{\eta d}}\sum _{j=1}^{i-1}{R}_j^d,$$

$$I({\mathbf{x}}_n\in \mathcal{A}){\tilde{L}}_n({\mathbf{x}}_n)\le \prod _{i=1}^n{\left(1-(i-1)\frac{{\gamma }^{\prime }}{{\lambda }^{\eta d}}{\overline{r}}^d\right)}^{+}=:{\delta }_n,$$

$$P_{\mathsf{acc}}(\lambda )=\frac{\mathcal{P}(\lambda )}{\mathbb{E}[{\tilde{\sigma }}_N]}=\frac{\mathcal{P}(\lambda )}{\mathbb{E}[{\delta }_N]}.$$

$${\hat{L}}_n({\mathbf{x}}_n):=\frac{d{\mu }^0}{d{\hat{\mu }}_n}({\mathbf{x}}_n)=\prod _{i=1}^n\left(1-{\hat{B}}_i\right),{\mathbf{x}}_n\in {\mathcal{G}}_n\cap \mathcal{A},$$

$$I({\mathbf{x}}_n\in \mathcal{A}){\hat{L}}_n({\mathbf{x}}_n)\le {\delta }_n=\prod _{i=1}^n{\left(1-(i-1)\frac{{\gamma }^{\prime }{\overline{r}}^d}{{\lambda }^{\eta d}}\right)}^{+},{\mathbf{x}}_n\in {\mathcal{G}}_n.$$

$${\mathcal{T}}_{\mathsf{ISAR}}\le c\mathbb{E}\left[{\delta }_N\right]\frac{{\lambda }^{\min \{\eta d,1\}}}{\mathcal{P}(\lambda )},$$

$$\underset{\lambda \nearrow \infty }{\mathrm{limsup}}\left[\frac{1}{{\lambda }^{2-\eta d}}\log \mathbb{E}\left[{\delta }_N\right]\right]$$

$$\underset{\lambda \nearrow \infty }{\mathrm{limsup}}\left[\frac{1}{\lambda }\log \mathbb{E}\left[{\delta }_N\right]\right]$$

Full text

Keywords: Exact Simulation, Dominated Coupling From The Past, Large Deviations, Non-overlapping Probability.

1 Introduction

Perfect sampling, that is, generating unbiased samples from a target distribution (also referred to as perfect simulation or exact sampling), is an important and exciting area of research in stochastic simulation. In this paper, we introduce and investigate a novel methodology for generating perfect samples of finite Gibbs hard-sphere models, which are an important family of Gibbs point processes. Roughly, a Gibbs hard-sphere model can be described as a set of spheres such that their centers constitute a Poisson point process on a bounded Euclidean space conditioned that no two spheres overlap with each other. The proposed methodology combines importance sampling (IS) and acceptance-rejection (AR) techniques to achieve substantial performance improvement in certain important regimes of interest. In statistical physics, there is a large body of work related to the Gibbs hard-sphere models; see, e.g., [35, 30, 1, 2, 28, 37, 24, 6]. These models are important also in modelling adsorption of latexes or proteins on solid surfaces [40, 38, and references therein]. For the analysis of wireless communication networks, it is common to use the Gibbs hard-sphere models to model base-stations in a cellular network because no two base-stations are to be normally placed closer than a certain distance from each other [39, 18]. Our results can be used to assess the stationary behaviour of Code Division Multiple Access (CDMA) wireless networks.

Literature Review: The existing literature offers several perfect sampling methods for Gibbs hard-sphere models. Among these, the dominated coupling from the past (dominated CFTP) methods are most prominent and they are based on the seminal paper by Propp & Wilson [36]; see [27, 25, 22, 26]. Another well-known perfect sampling method for the Gibbs hard-sphere models is called the backward-forward algorithm (BFA) by Ferrari et al. [12]; also see [23, 16]. To see some of the applications of perfect sampling for these models, refer to [8, 7, 33]. For other related literature on perfect sampling for spatial point processes, refer to [34, 19]. As mentioned in [16], all the existing methods are, in some sense, complementary to each other. They take advantage of an important property that the distribution of a Gibbs hard-sphere model can be realized as an invariant measure of a spatial birth-and-death process, call it the interaction process. For example, the main ingredient of the dominated CFTP method is to construct a birth-and-death process backward in time starting from its steady-state at time zero such that it dominates the interaction process, and then use thinning on the dominating process to construct coupled upper and lower bound processes forward in time such that the coalescence of these two bounding processes assures a perfect sample from the target measure, which is the invariant measure of the interaction process. The BFA is based on the construction of the clan of ancestors that uses thinning of a dominating process and extends the applicability to infinite-volume measures. A crucial drawback of the naive AR and the dominated CFTP methods is that they are guaranteed to be efficient only if the intensity of the Gibbs hard-sphere model is close to the intensity of the reference Poisson point process; see [23] for details. In addition, most of the dominated CFTP methods suffer from the so-called impatient-user bias (a bias that is induced when a user aborts long runs of the algorithm); see [13], [14] and [41].

Our Contributions: Acceptance-rejection methods are free of the impatient-user bias and involve neither thinning nor coupling (which are crucial for the other methods). Despite being an obvious alternative to the existing methods, to the best of our knowledge, in the context of Gibbs point processes, the use of AR methods is still largely unexplored (except brief discussions, e.g., in [15] and [23]). AR methods for Gibbs hard-sphere models are amenable to further algorithmic enhancements that may substantially decrease the expected running time of the algorithm. The proposed methodology provides one such enhancement. To highlight the significance of the proposed methodology, we compare its running time complexity with that of both the naive AR and the dominated CFTP methods. This effectiveness analysis is based on our large deviations analysis of the non-overlapping probability. A brief summary of our results is as follows.

• •

Our first key contribution is that we conduct a large deviations analysis of the probability of spheres not overlapping with each other when their centers constitute a homogeneous Poisson point process (PPP). More specifically, we consider a homogeneous marked PPP on $[0,1]^{d}$ with intensity $\lambda$ where the points are the center of spheres with independently and identically distributed (iid) radii as marks which are independent of the centers and identical in distribution to $R/\lambda^{\eta}$ for a positive bounded random variable $R$ and a constant ${\eta>0}$. We establish large deviations of the probability of spheres do not overlap with each other, as ${\lambda\nearrow\infty}$. This analysis is useful in the study of the asymptotic behavior of the expected running time complexities of the proposed and the existing perfect sampling methods for the Gibbs hard-sphere models. This analysis may also be of independent interest.

• •

Our second key contribution is that we propose a novel IS based AR algorithm for generating perfect samples of the Gibbs hard-sphere model obtained by considering the homogeneous marked PPP conditioned on no overlap of the spheres. This is achieved by partitioning the underlying configuration space and arriving at an appropriate change of measure on each partition. Applicability of the proposed algorithm is illustrated in two scenarios. In the first scenario, all the spheres are assumed to be of a fixed size (i.e., $R$ is a fixed positive constant). We develop a grid based IS technique under which spheres are generated sequentially such that the chance of spheres overlapping is small and the corresponding likelihood ratio has a better deterministic upper bound that improves the acceptance probability in each iteration of the algorithm. In the second scenario, we consider the general case where spheres have iid radii. In this scenario, we divide the underlying configuration space into two sets. On one set, the sum of the volumes of spheres is bounded from below and on the other set, the volume sum takes small values so that the set consists of rare configurations only. For the first set, we develop a grid based IS method that is similar to the one stated above, and for the second set, we use an exponential twisting on the sphere volume distribution. In both the scenarios, the new method provably substantially improves the performance of the algorithm compared to the naive AR method.

• •

We analytically and numerically compare the performance of the proposed IS based AR method with that of some of the dominated CFTP methods. The numerical results support our analytical conclusions that the proposed method is substantially efficient compared to the existing methods over the high density regime where $\eta d\leq 1$ and $\lambda$ is large.

Organization: Section 2 provides a definition of the hard-sphere model. The large deviations of the non-overlapping probability is presented in Section 3. In Section 4, we first review a naive AR method and analyze its expected running time complexity, and we then propose and analyze the IS based AR method. In Section 5, a review of the well-known dominated CFTP methods for the hard-sphere models is given. Section 6 illustrates the efficiency of the proposed methodology using numerical experiments. Section 7 is a brief conclusion of the paper. All proofs are presented in Appendix A.

2 Preliminaries

First we introduce some notation. $X\sim F$ denotes that the distribution of a random object $X$ is $F$. $\mathsf{Poi}(\lambda)$ and $\mathsf{Bern}(p)$ denote, respectively, Poisson distribution with mean $\lambda>0$ and Bernoulli distribution with success probability $p$. The uniform distribution on $[0,1]$ is denoted by $\mathsf{Unif}(0,1)$. For an event $A$, the indicator function $I(A)$ is equal to $1$ if $A$ occurs, otherwise it is equal to [math]. A measure $\mu_{1}$ is absolutely continuous with respect to a measure $\mu_{2}$ on a measurable set $A$ if $\mu_{1}(B\cap A)=0$ for any measurable $B$ such that $\mu_{2}(B\cap A)=0$. For any probability measure $\mu$, $\mathbb{P}_{\mu}(A)$ denotes the probability of an event $A$ under $\mu$, and $\mathbb{E}_{\mu}[\cdot]$ denotes the associated expectation. We drop the subscript $\mu$ when it is not relevant. For any non-negative real valued functions $f$ and $g$, write $f(x)=O(g(x))$ if $\limsup_{x\rightarrow\infty}f(x)/g(x)~{}\leq~{}c$ for some constant $c>0$, write $f(x)=\varOmega(g(x))$ if $g(x)=O(f(x))$, and write $f(x)=o(g(x))$ if $\limsup_{x\rightarrow\infty}f(x)/g(x)=0$. Write $f(x)=\Theta(g(x))$ if both $f(x)=O(g(x))$ and $f(x)=\varOmega(g(x))$ are true. For any real value $x$, the largest integer $n$ such that $n\leq x$ is denoted by $\lfloor x\rfloor$ and the smallest integer $n$ such that $n\geq x$ is denoted by $\lceil x\rceil$. The set of all the non-negative integers is denoted by $\mathbb{N}_{0}$.

A random finite subset $\mathbf{X}=\{X_{1},\dots,X_{N}\}$ of an observation window $W\subset\mathbb{R}^{d}$ is called a Poisson point process (PPP) with a finite intensity measure $\nu$ on $W$ if $N\sim\mathsf{Poi}(\nu(W))$ and for every $n\in\mathbb{N}_{0}$, conditioned on $N=n$, the points $X_{1},\dots,X_{n}$ are iid with distribution $\nu(\mathrm{d}x)/\nu(W)$. A PPP on $[0,1]^{d}$ is called $\lambda$-homogeneous PPP with intensity $\lambda>0$ if the intensity measure $\nu(\mathrm{d}x)=\lambda\,\mathrm{d}x$, where $\mathrm{d}x$ is Lebesgue measure on $W$. To each point $X_{i}$ of the $\lambda$-homogeneous PPP on $[0,1]^{d}$, we associate a mark which is a non-negative number interpreted as the radius of a sphere centered at $X_{i}$. In particular, a $\lambda$-homogeneous marked PPP on $[0,1]^{d}$ is a PPP on $W=[0,1]^{d}\times[0,\infty)$ with the intensity measure ${\nu(\mathrm{d}x\times\mathrm{d}r)=\lambda\mathrm{d}x\times F(\mathrm{d}r)}$ where $F$ is the distribution of each radius. That is, the centers constitute a $\lambda$-homogeneous PPP on $[0,1]^{d}$ which is independent of the radii, and the radii are iid with distribution $F$. A realization of the marked PPP with $n$ points is denoted by $\mathbf{x}=\{(y_{1},r_{1}),\dots,(y_{n},r_{n})\}$, where $r_{i}\geq 0$ is the radius of the sphere centered at $y_{i}\in[0,1]^{d}$. Define $\mathscr{G}=\cup_{n\in\mathbb{N}_{0}}\mathscr{G}_{n}$ where

[TABLE]

Now we define a Gibbs hard-sphere model. Suppose that $\mu^{0}$ is the distribution of a $\lambda$-homogeneous marked PPP as defined above with $F$ being the distribution of $R/\lambda^{\eta}$ for a constant $\eta>0$ and a non-negative random variable $R$. Let $\mathscr{A}\subset\mathscr{G}$ be the set of all configurations with no two spheres overlapping with each other. Then the distribution $\mu$ of the Gibbs hard-sphere model is absolutely continuous with respect to $\mu^{0}$ with the Radon-Nikodym derivative given by

[TABLE]

where the normalizing constant $\mathcal{P}(\lambda)$ is the non-overlapping probability given by

[TABLE]

We refer to the Gibbs hard-sphere model as a torus-hard-sphere model if the boundary of the underlying space $[0,1]^{d}$ is periodic, that is, a sphere $S(x,a)$ centered at $x\in[0,1]^{d}$ with radius $r$ is defined by

[TABLE]

where $\|\cdot\|$ is the $d$-dimensional Euclidean norm and ’$\mathsf{mod}$’ denotes the modulo operation [9]. If the boundary is not periodic, we refer to the model as a Euclidean-hard-sphere model.

From now onwards, the phrase ‘hard-sphere model’ refers to either of these two models and we assume that $R$ is bounded from above by a constant $\overline{r}>0$. In particular, if $R$ is a constant, we take $\overline{r}=R$. Furthermore, we assume that $2\overline{r}/\lambda^{\eta}<1$ to avoid certain trivial difficulties such as the possibility of a sphere on the torus overlapping with itself.

3 Large Deviations Results

In this section, we obtain large deviations results for the non-overlapping probability $\mathcal{P}(\lambda)$. We use these results for analyzing the running time complexity of both the naive and importance sampling based acceptance-rejection methods. Hereafter, $\gamma={\pi^{d/2}}/{\Gamma(d/2+1)}$, where $\Gamma(\cdot)$ is the gamma function. Note that the volume of a sphere with radius $r$ is given by $\gamma r^{d}$. Define $m_{1}:={\mathbb{E}}[(R+\widehat{R})^{d}]$, where $\widehat{R}$ is independent and identical in distribution to $R$, and let

[TABLE]

Theorem 1**.**

The non-overlapping probability $\mathcal{P}(\lambda)$ satisfies

[TABLE]

When $\eta d=1$, the limit $\delta:=\lim_{\lambda\rightarrow\infty}\left[\frac{1}{\lambda}\log\mathcal{P}(\lambda)\right]$ exists and $-1\leq\delta<0$. Furthermore, $\delta\nearrow 0$ if $\gamma m_{1}\searrow 0$, and $\delta\leq-\frac{1}{2}\left(1-\frac{1}{\gamma^{\prime}\overline{r}^{d}}\right)^{2}$ if $R\equiv\overline{r}$ and $\gamma^{\prime}\overline{r}^{d}>1$. In addition, for the torus-hard-sphere model,

[TABLE]

An important and fundamental characteristic of a Gibbs point process is its intensity; see, for example, [29, and references therein] and [5]. Roughly speaking, the intensity of a Gibbs point process is the expected number of points of the process per unit volume. There is an interesting connection between the regimes considered in Theorem 1 and the asymptotic intensity of the torus-hard-sphere model. To see this, assume that each sphere has a fixed radius $\overline{r}/\lambda^{\eta}$. Since the underlying space is $[0,1]^{d}$, the intensity $\rho(\lambda)$ of the model is exactly equal to the expected total number of points in a realization of the model. Equivalently, we may consider the fraction of the volume $\mathsf{VF}(\lambda)$ occupied by the spheres, given by $\mathsf{VF}(\lambda)=\rho(\lambda)\gamma\overline{r}^{d}\lambda^{-\eta d}$. For the torus-hard-sphere model, the volume fraction $\mathsf{VF}(\lambda)$ is bounded from above by $\rho^{\max}\gamma$, where $\rho^{\max}$ is the closest packing density defined by $\rho^{\max}=\lim_{n\rightarrow\infty}N_{n}/(n+1)^{d},$ with $N_{n}$ being the maximal number of mutually disjoint unit radius spheres which are included in the hypercube $[-(n+1/2),(n+1/2)]^{d}$; see [29]. Proposition 1 describes asymptotic behavior of $\mathsf{VF}(\lambda)$ as $\lambda\to\infty$ for different values of $\eta d$. In particular, the regime with $\eta d>1$ is a low density regime while the regime with $\eta d<1$ is a high density regime. In the high density regime, the intensity of the hard-sphere model is much smaller than the intensity $\lambda$ of the reference PPP.

Proposition 1**.**

For the torus-hard-sphere model with a fixed radius $R=\overline{r}$,

[TABLE]

4 Acceptance-Rejection Based Algorithms

In Section 4.1, we present a naive acceptance-rejection (AR) algorithm for generating perfect samples of the hard-sphere model and analyze its expected running time complexity. We then proceed to present and analyze our importance sampling (IS) based AR algorithm where the key idea is to partition the configuration space $\mathscr{G}$ so that a well chosen IS technique can be implemented on each partition. One such IS for the hard-sphere model is the reference IS presented in Section 4.3 where spheres are generated sequentially such that, whenever possible, the center of each sphere is selected uniformly over the region on $[0,1]^{d}$ that guarantees no overlap with the existing spheres. However, generating samples from this IS measure can be computationally challenging when $d\geq 2$. The grid based IS introduced in Sections 4.4 and 4.5 overcomes this difficulty by imitating the reference IS, and interestingly, it is more efficient than the reference IS.

In every algorithm presented in this paper, the running time complexity is calculated under the assumption that checking overlap of a newly generated sphere with an existing sphere is done in a sequential manner. That is, if there are $n$ existing spheres, the expected running time complexity of the overlap check is proportional to $n$. However, if enough computing resources are available, the overlap check can be done in parallel so that its running time complexity is a constant. We omit the discussion of this parallel overlap check because it is easy to modify the results to accommodate the parallel case, and also the key conclusions of the paper do not change.

4.1 Naive AR Algorithm

Algorithm 1 is a naive AR algorithm for generating perfect samples of the Gibbs hard-sphere model. The basic idea of the algorithm is standard [11], and its correctness is straightforward and hence omitted.

Let ${\mathcal{T}}_{\mathsf{NAR}}$ be the expected running time complexity of Algorithm 1, where the running time complexity denotes the number of elementary operations performed by the algorithm; every elementary operation takes at most a fixed amount of time. Note that the acceptance probability of each iteration is $\mathcal{P}(\lambda)$. Then the expected total number of iterations of the algorithm is $1/\mathcal{P}(\lambda)$. Suppose $C_{\mathsf{itr}}(\lambda)$ is the expected running time complexity of an iteration. Then,

[TABLE]

We now establish bounds on ${\mathcal{T}}_{\mathsf{NAR}}$, and then establish its asymptotic behavior as $\lambda\nearrow\infty$ using Theorem 1. In each iteration of Algorithm 1, spheres are generated in a sequential order until we see an overlap or a configuration with $N$ non-overlapping spheres. The key to prove Proposition 2 is to establish that the expected number of spheres generated per iteration is $\Theta\left(\lambda^{\min\{\eta d,2\}}\right)$.

Proposition 2**.**

The expected running time complexity $C_{\mathsf{itr}}(\lambda)$ of an iteration of the naive AR algorithm, Algorithm 1, satisfies

[TABLE]

Furthermore, the expected total running time $\mathcal{T}_{\mathsf{NAR}}$ satisfies:

[TABLE]

Remark 1**.**

From (5) and Theorem 1, we see that for large values of $\lambda$ and for $\eta d<2$, $\mathcal{T}_{\mathsf{NAR}}$ is mainly governed by $\mathcal{P}(\lambda)$, which can be very small for large $\lambda$. This suggests that any rejection based perfect sampling algorithm with a significant improvement in the acceptance probability will have a significantly improved running time complexity.

4.2 Importance Sampling Based Acceptance-Rejection Algorithm

A sequence of tuples ${\left\{\left(D_{n,k},\mu_{n,k},\sigma_{n,k}\right)_{k=1}^{K}\right\}_{n\in\mathbb{N}_{0}}}$ with some $K\leq\infty$ is called stable IS sequence if for each ${n\in\mathbb{N}_{0}}$, ${\left(D_{n,k}\right)_{k=1}^{K}}$ is a partition of $\mathscr{G}_{n}$, and $(\mu_{n,k})_{k=1}^{K}$ a sequence of probability measures such that $\mu^{0}$ is absolutely continuous with respect to $\mu_{n,k}$ on $D_{n,k}\cap\mathscr{A}$ and the corresponding likelihood ratio ${L_{n,k}(\mathbf{x}_{n}):=\frac{d\mu^{0}}{d\mu_{n,k}}(\mathbf{x}_{n})}$ satisfies

[TABLE]

for $k=1,\dots,K$. Under the stability condition, for every measurable subset ${\mathscr{B}\subseteq\mathscr{G}}$,

[TABLE]

where ${\widetilde{\sigma}(n):=\sum_{k=1}^{K}\sigma_{n,k}}$, ${U\sim\mathsf{Unif}(0,1)}$ and $J\sim\mathsf{Bern}\left(\frac{L_{n,k}(\mathbf{X}_{n})}{\sigma_{n,k}}\right)$. Let $M$ be a non-negative integer valued random variable with the pmf defined by,

[TABLE]

where $C_{\lambda}:=\sum_{n=0}^{\infty}\frac{\lambda^{n}\widetilde{\sigma}(n)}{n!}$. The pmf (7) is well defined because $\mathbb{E}\left[\widetilde{\sigma}(N)\right]$ is finite under the stability condition. Now consider Algorithm 2.

Proposition 3**.**

Algorithm 2 generates a perfect sample of the Gibbs hard-sphere model. Furthermore, let $N\sim\mathsf{Poi}(\lambda)$. Then the probability of accepting the configuration generated in an iteration of Algorithm 2 is given by

[TABLE]

We omit the proof of Proposition 3 because the correctness easily follows from (6), and (8) holds from the observation that

[TABLE]

Note that the expected number of iterations of Algorithm 2 is $1/P_{\mathsf{acc}}(\lambda)$. Corollary 1 is an important and trivial consequence of Proposition 3.

Corollary 1**.**

For all stable IS sequences $\left\{\left(D_{n,k},\mu_{n,k},\sigma_{n,k}\right)_{k=1}^{K}\right\}_{n\in\mathbb{N}_{0}}$ with the same $\mathbb{E}[\widetilde{\sigma}(N)]=\sum_{k=1}^{K}\mathbb{E}\left[\sigma_{N,k}\right]$, the expected number of iterations of Algorithm 2 is the same.

Suppose that $\widetilde{C}_{\mathsf{itr}}(\lambda)$ is the expected running time complexity of an iteration of Algorithm 2. Then the expected total running time of the algorithm is given by

[TABLE]

where ${N\sim\mathsf{Poi}(\lambda)}$. Recall that the acceptance probability of the naive AR method is $\mathcal{P}(\lambda)$. It is reasonable to seek a valid stable IS sequence ${\left\{\left(D_{n,k},\mu_{n,k},\sigma_{n,k}\right)_{k=1}^{K}\right\}_{n\in\mathbb{N}_{0}}}$ so that ${\widetilde{C}_{\mathsf{itr}}(\lambda)\mathbb{E}[\widetilde{\sigma}(N)]}$ is smaller than ${C_{\mathsf{itr}}(\lambda)}$. In Subsections 4.4 and 4.5, we present applications of Algorithm 2 where $\mathcal{T}_{\mathsf{ISAR}}$ is indeed much smaller than $\mathcal{T}_{\mathsf{NAR}}$.

Remark 2** (Extension of IS Based AR to General Gibbs Point Processes).**

Suppose that $\mu$ is the distribution of a Gibbs point process that is absolutely continuous with respect to $\mu^{0}$ with the corresponding Radon-Nikodym derivative given by ${\frac{d\mu}{d\mu^{0}}\left(\mathbf{x}\right)=\frac{\exp\left(-\beta\,V(\mathbf{x})\right)}{Z},\,\mathbf{x}\in\mathscr{G}},$ where the constant $\beta\in\mathbb{R}$ is known as inverse temperature, $V$ is called non-negative potential function, and the normalizing constant ${Z=\mathbb{E}_{\mu^{0}}\left[\exp\left(-\beta\,V(\mathbf{X})\right)\right]}$. If the stability condition holds true when $I(\mathbf{x}_{n}\in\mathscr{A})$ is replaced by $\exp\left(-\beta\,V(\mathbf{x}_{n})\right)$, then Algorithm 2 can generate perfect samples from $\mu$ if in line 5 of the algorithm,

[TABLE]

To see that the hard-sphere model is a special case of such a Gibbs point process, take $\beta>0$ and assume that ${V(\mathbf{x})=0}$ if $\mathbf{x}$ is a non-overlapping configuration of spheres, otherwise, ${V(\mathbf{x})=\infty}$.

4.3 Reference Importance Sampling

We now introduce an IS measure, called reference IS and denoted by $\widetilde{\mu}_{n}$ for each $n$, so that $\left\{\left(\mathscr{G}_{n},\widetilde{\mu}_{n},\sigma_{n}\right)\right\}_{n\in\mathbb{N}_{0}}$ is a stable IS sequence (with $K=1$) that can be used in Algorithm 2 for generating perfect samples of the hard-sphere model for an appropriate choice of the sequence $\left\{\sigma_{n}:n\in\mathbb{N}_{0}\right\}$. Under $\widetilde{\mu}_{n}$, first generate iid sequence $R_{1},\dots,R_{n}$ identical in distribution to $R$, and then $n$ spheres are generated sequentially as follows. Generate the center of the first sphere uniformly distributed on $[0,1]^{d}$. Suppose that $i-1$ spheres are already generated. For the $i^{th}$ sphere generation, a subset $\mathcal{B}_{i}\subseteq[0,1]^{d}$ is called blocking region if $\mathcal{B}_{i}$ is the largest set such that the center $Y_{i}$ of the $i^{th}$ sphere falling in this region (that is, $Y_{i}\in\mathcal{B}_{i}$) would result in an overlap of the $i^{th}$ sphere with one of the existing $i-1$ spheres. The center of the $i^{th}$ sphere is generated with uniform distribution over the non-blocking region $[0,1]^{d}\setminus\mathcal{B}_{i}$. If for any sphere $i\leq n$, the entire space is blocked (that is, $\mathcal{B}_{i}=[0,1]^{d}$), we select the centers of spheres $i,\dots,n$ arbitrarily. Figure 1 illustrates this for $d=2$ and $n=1,2$. In conclusion, $\widetilde{\mu}_{n}$ is the distribution of an output of Algorithm 3.

Observe that $\mu^{0}$ is absolutely continuous with respect to $\widetilde{\mu}_{n}$ on $\mathscr{G}_{n}\cap\mathscr{A}$, and the associated likelihood ratio satisfies

[TABLE]

for all $\mathbf{x}_{n}\in\mathscr{G}_{n}\cap\mathscr{A}$ and $n\in\mathbb{N}_{0}$, where $B_{i}$ is the volume of $\mathcal{B}_{i}$ and $\widetilde{L}_{0}=1$. Note that $\widetilde{L}_{n}(\mathbf{x}_{n})=0$ if and only if $\mathbf{x}_{n}\notin\mathscr{A}$ because for any $\mathbf{x}_{n}\notin\mathscr{A}$, there exists $i\leq n$ such that $B_{i}=1$.

Observe that the blocking volume added by the $i^{th}$ sphere is at least $\gamma^{\prime}\left(R_{i}/\lambda^{\eta}\right)^{d}$ when it does not overlap with any of the existing spheres. This is because, for the torus-hard-sphere model, the entire volume within an accepted sphere is added to blocking volume, and for the Euclidean-hard-sphere model, at least $1/2^{d}$ fraction of an accepted sphere is added to the blocking volume. Thus,

[TABLE]

for every configuration $\mathbf{x}_{i-1}\in\mathscr{G}_{i-1}\cap\mathscr{A}$. In particular, if all the spheres are of the same size with a fixed radius $\overline{r}$,

[TABLE]

for all $n\in\mathbb{N}_{0}$ and $\mathbf{x}_{n}\in\mathscr{G}_{n}$, where $x^{+}=\max(0,x)$ and $\delta_{0}=1$. Then the stability condition is satisfied with $K=1$, $D_{n,1}=\mathscr{G}_{n}$, $\mu_{n,1}=\widetilde{\mu}_{n}$ and $\sigma_{n,1}=\delta_{n}$ for $n\in\mathbb{N}_{0}$. Thus, Algorithm 2 generates perfect samples of the fixed radius hard-sphere model, and from Proposition 3, the corresponding acceptance probability

[TABLE]

Remark 3**.**

When the dimension $d=1$, spheres become line segments and thus it is easy to generate samples from the IS measure $\widetilde{\mu}_{n}$. However, for $d\geq 2$, generating samples under the reference IS is difficult because every time a new sphere is generated, we need to know the volume of the blocking region created by the existing spheres and then we need to generate a point uniformly on this non-blocking region; see line 11 in Algorithm 3. One possible way to implement the reference IS is by combining a well-known method called power tessellation and a simple rejection method in two steps: i) Using the power tessellation, we can compute the blocking volumes exactly; see, e.g, [4] and [32]. ii) Then, use a simple acceptance-rejection method where repeatedly a point is generated independently and uniformly on $[0,1]^{d}$ until it falls within the non-blocking region. Unfortunately, implementing the power tessellation method is computationally prohibitive. Besides, even if we have an efficient implementation of the power tessellation method, the above simple rejection step can be expensive when the non-blocking region is small. Fortunately, we can overcome both these difficulties by using a simple grid on $[0,1]^{d}$. From (9), it is evident that if there are two IS methods with the same $\mathbb{E}[\widetilde{\sigma}_{N}]$, it is computationally preferable to use the method that has smaller per iteration expected running time, $\widetilde{C}_{\mathsf{itr}}(\lambda)$. In Subsection 4.4, we introduce a hyper-cubic grid based IS method that continues to generate perfect samples while the blocking regions are closely approximated by grid cells. With a careful choice of the cell-edge length, we make sure that the inequality (12) holds for the grid IS as well (and thus, $\mathbb{E}[\widetilde{\sigma}_{N}]$ is same as that of the reference IS). As a consequence of Corollary 1, the expected iterations of Algorithm 2 is the same as that of the reference IS method. However, the grid method is easy to implement and has a much smaller expected iteration cost $\widetilde{C}_{\mathsf{itr}}(\lambda)$ compared to that of the reference IS. The choice of the hyper-cubic grid is just an option as it simplifies the implementation. However, the method can be implemented using other kinds of grids. In two dimensional case, for example, it is possible to use a hexagonal grid for implementing the IS method.

4.4 Grid Based Importance Sampling for Fixed Radius Case

Consider the hard-sphere model with a fixed radius $\overline{r}/\lambda^{\eta}$. Generation of $n$ spheres under the following grid based IS measure $\widehat{\mu}_{n}$ starts by partitioning the underlying space $[0,1]^{d}$ into a hyper-cubic grid with a cell-edge length $\varepsilon>0$ such that $1/\varepsilon$ is an integer. The centers of the spheres are generated in a sequential order: Suppose that $i-1$ spheres with centers $Y_{1},\dots,Y_{i-1}$ are already generated. At the time of $i^{th}$ sphere generation, a cell $C$ in the grid is labeled as fully-blocked if the cell is completely inside a sphere with radius $2\overline{r}/\lambda^{\eta}$ centered at an existing point, that is, $C\subseteq S(Y_{j},2\overline{r}/\lambda^{\eta})$ for some $j\leq i-1$; otherwise, the cell is labeled as non-fully-blocked. A non-fully-blocked cell $C$ is called partially-blocked if $C\cap S(Y_{j},2\overline{r}/\lambda^{\eta})\neq\varnothing$ for some $j\leq i-1$; otherwise, it is called non-blocked. The center $Y_{i}$ of the $i^{th}$ sphere is selected uniformly over the non-fully-blocked cells, because selecting $Y_{i}$ over a fully-blocked cell will certainly result in the $i^{th}$ sphere overlapping with an existing sphere. We then check for overlap only if $Y_{i}$ is generated over a partially-blocked cell, because the overlap is not possible if $Y_{i}$ is generated over a non-blocked cell. If either there is an overlap or all the cells are fully-blocked by the existing spheres, the centers $Y_{i},\dots,Y_{n}$ of the remaining spheres are selected arbitrarily (such a selection results in an overlapping configuration). Otherwise, for the next sphere $i+1$ generation, we repeat the same procedure by relabeling the non-fully-blocked cells by considering spheres $1,\dots,i$ as the existing spheres. Note that at the beginning of each iteration all the cells are labeled as non-blocked. Also note that since all the spheres have the same radius, for relabeling of the cells, we only need to focus on the cells that might interact with the last sphere generated. See Figure 2 for an illustration of this sequential procedure.

Suppose that $\widehat{\mu}_{n}$ is the probability measure under which $n$ spheres are generated by the above procedure. Then $\mu^{0}$ is absolutely continuous with respect to $\widehat{\mu}_{n}$ on ${\mathscr{G}}_{n}\cap\mathscr{A}$ and the corresponding likelihood ratio

[TABLE]

where $\widehat{B}_{i}$ is the volume of fully-blocking cells for the $i^{th}$ sphere generation, that is, $\widehat{B}_{i}$ equal to the product of the number of fully-blocked cells and $\varepsilon^{d}$. To apply Algorithm 2 for the fixed radius hard-sphere model, take $K=1$ and for each ${n\in\mathbb{N}_{0}}$, take ${D_{n,1}={\mathscr{G}}_{n}}$, ${\mu_{n,1}=\widehat{\mu}_{n}}$ and ${\sigma_{n,1}=\delta_{n}}$. Thus, ${\widetilde{\sigma}(n)=\delta_{n}}$ and ${L_{n,1}(\mathbf{x}_{n})=\widehat{L}_{n}(\mathbf{x}_{n})}$ for all ${\mathbf{x}_{n}\in{\mathscr{G}}_{n}\cap\mathscr{A}}$.

Selection of the cell-edge length $\boldsymbol{\varepsilon}$: Observe that the longest diagonal length of a cell is $\sqrt{d}\,\varepsilon$. Since we focus only on the non-overlapping configurations, in the implementation, we generate a sphere only if all the existing spheres are non-overlapping. Suppose that the cell-edge length $\varepsilon$ is selected so that $\displaystyle\sqrt{d}\,\varepsilon\leq\overline{r}/\lambda^{\eta}.$ Then for the $i^{th}$ sphere generation, every cell that has non-empty intersection with $S(Y_{j},\overline{r}/\lambda^{\eta})$, for any ${j=1,\dots,i-1}$, has to be fully-blocked, because such a cell is a subset of $S\left(Y_{j},2\overline{r}/\lambda^{\eta}\right)$. Thus, the non-overlapping condition of the existing spheres imply that $\cup_{j=1}^{i-1}S(Y_{j},\overline{r}/\lambda^{\eta})$ is a subset of the union of the fully-blocked cells, and hence ${\widehat{B}_{i}\geq\frac{\gamma^{\prime}(i-1)\overline{r}^{d}}{\lambda^{\eta d}}}.$ Thus, for $n\geq 1$,

[TABLE]

This upper bound is same as that we obtained in the case of the reference IS; see the inequality (12). Since the acceptance probability ${P_{\mathsf{acc}}(\lambda)=\mathcal{P}(\lambda)/\mathbb{E}[\delta_{N}]}$ is the same for both the grid IS and the reference IS methods, we need to choose the cell-edge length $\varepsilon\leq\overline{r}/\lambda^{\eta}$ so that the expected per iteration running time $\widetilde{C}_{\mathsf{itr}}(\lambda)$ is minimum. It is easy to see that the higher the value of $\varepsilon$, the smaller $\widetilde{C}_{\mathsf{itr}}(\lambda)$ due to the following reasons:

Labelling of the cells is faster if they are bigger in size; 2. 2.

Increase in the cell size increases the chances of overlap of the new sphere with the existing spheres, and hence on average each iteration generates fewer spheres;

In conclusion, we choose $\varepsilon=1/\lfloor\lambda^{\eta}/\overline{r}\rfloor$ for the implementation of the grid IS method.

To reduce the per iteration complexity of the algorithm, we make some changes to the steps 4 and 5 in Algorithm 2. Observe that a realization $\mathbf{X}_{n}$ generated under $\widehat{\mu}_{n}$ is accepted only if ${\mathbf{X}_{n}\in\mathscr{A}}$ and ${J=1}$, where ${J\sim\mathsf{Bern}\left(\widehat{L}_{n}(\mathbf{X}_{n})/\delta_{n}\right)}$. In the implementation, we generate an iid sequence $U_{1},\dots,U_{n}\sim\mathsf{Unif}(0,1)$ independent of everything else so far generated, and take $J_{i}=I\left(U_{i}\leq\frac{1-\widehat{B}_{i}}{(1-(i-1)\gamma^{\prime}\overline{r}^{d}\lambda^{-\eta d})}\right)$ for $i\leq n$. Since $J$ and the product $\prod_{i=1}^{n}J_{i}$ are Bernoulli random variables with the same success probability $L(\mathbf{X}_{n})/\delta_{n}$, to reduce the per iteration cost, we generate the $i^{th}$ sphere only if $J_{i}=1$ and the existing spheres do not overlap with each other.

Algorithm 4 implements the grid based IS for a given $n$ with the above mentioned enhancements. Algorithm 2 is restated as Algorithm 5.

Remark 4** (The pmf of $M$).**

Note that, for the current setup, the pmf of $M$, given by (7), becomes $\mathbb{P}\left(M=m\right)=\frac{1}{C_{\lambda}}\frac{\lambda^{m}\delta_{m}}{m!}$, $m\in\mathbb{N}_{0}$, where the normalizing constant $C_{\lambda}=\sum_{n\in\mathbb{N}_{0}}\frac{\lambda^{n}\delta_{n}}{n!}$. The support of the pmf is finite because $\delta_{m}=0$ for all $m\geq\lambda^{\eta d}/(\gamma^{\prime}\overline{r}^{d})+1$. To increase the performance of the algorithm, we can further truncate the support of the pmf. Using the maximum packing density, we can obtain an integer $m_{\max}$ such that $\mathbf{X}\notin\mathscr{A}$ for all $m\geq m_{\max}$ and configurations $\mathbf{X}$ with $|\mathbf{X}|=m$. In that case, we can take $\mathbb{P}\left(M=m\right)=\frac{1}{C_{\lambda}}\frac{\lambda^{m}\delta_{m}}{m!}$, $0\leq m\leq m_{\max}$, with $C_{\lambda}=\sum_{n=0}^{m_{\max}}\frac{\lambda^{n}\delta_{n}}{n!}$. For example, refer to [29] for finding maximum packing densities for $d=2$ and $d=3$.

We now focus on the expected running time analysis of Algorithm 5. By Proposition 3, the acceptance probability $P_{\mathsf{acc}}(\lambda)$ of Algorithm 5 is $\mathcal{P}(\lambda)/\mathbb{E}[\widetilde{\sigma}(N)]=\mathcal{P}(\lambda)/\mathbb{E}[\delta_{N}].$ A proof of Proposition 4 is given in Section A.4.

Proposition 4**.**

For the fixed radius hard-sphere model, there exists a constant $c>0$ such that

[TABLE]

where $N\sim\mathsf{Poi}(\lambda)$. Furthermore,

[TABLE]

The following result is a trivial consequence of Propositions 2 and 4.

Corollary 2**.**

For the fixed radius hard-sphere model, if $\eta d\geq 2$, both $\mathcal{T}_{\mathsf{ISAR}}$ and $\mathcal{T}_{\mathsf{NAR}}$ are of the same order, and if $0<\eta d<2$, there exists a constant $c>0$ such that $\displaystyle\mathcal{T}_{\mathsf{ISAR}}\leq c\,\mathbb{E}\left[\delta_{N}\right]\mathcal{T}_{\mathsf{NAR}}.$

Remark 5** (Better choice of $\delta_{n}$ for the Euclidean-hard-sphere model).**

If the spheres are Euclidean, further improvements in the choice of $\delta_{n}$ can be obtained by accounting for boundary effects. For instance, for $d=2$, the four corners of $[0,1]^{2}$ are covered by at most $4$ circles, each of which contributing a blocking area of at least $\gamma^{\prime}\overline{r}^{2}/\lambda^{2\eta}=\pi\overline{r}^{2}/4\lambda^{2\eta}$, while each of the remaining circles contributing a blocking area of at least $2\gamma^{\prime}\overline{r}^{2}/\lambda^{2\eta}=\pi\overline{r}^{2}/2\lambda^{2\eta}$. Let $b_{0}=0$, $b_{i}=(i-1)\frac{\pi\overline{r}^{2}}{4\lambda^{2\eta}}$ for $1\leq i\leq 5,$ and $b_{i}=\frac{\pi\overline{r}^{2}}{\lambda^{2\eta}}+(i-4)\frac{\pi\overline{r}^{2}}{2\lambda^{2\eta}}$ for $i\geq 6$. Then, for this particular scenario, a better choice of $\delta_{n}$ in (12) (as well as in (13)) is $\delta_{n}=\prod_{i=1}^{n}\left(1-b_{n}\right)^{+},\,n\in\mathbb{N}_{0}$.

4.5 Random Radii Case

We now consider another application of Algorithm 2 for the hard-sphere model when under the marked PPP the radii of the spheres are iid. For the fixed radius case presented in Section 4.4, the proposed IS method ensured a uniform bound $\delta_{n}$ on the likelihood ratio over ${\mathscr{G}}_{n}$ for every $n\in\mathbb{N}_{0}$, as shown in (13). Such upper bounds are possible for a random radii hard-sphere model if the radii are bounded below by a positive constant. Furthermore, a similar analysis can be established when the spheres are replaced with iid convex shapes such that each shape occupies a minimum positive volume. However, when the radii are not bounded from below almost surely, the associated blocking volumes can be arbitrarily small. We address this issue by partitioning ${\mathscr{G}}_{n}$ into two sets $D_{n,1}$ and $D_{n,2}$ for each $n$ so that the IS on $D_{n,1}$ is a grid based IS method that is similar to Algorithm 4 and the IS on $D_{n,2}$ is obtained by exponentially twisting the distribution of $R^{d}$ to put high probability mass on configurations with lower volume spheres.

We first introduce the exponential twisting of the distribution, say $G$, of $R^{d}$. Recall that $R$ is assumed to be a bounded non-negative random variable. Without loss of generality further assume that ${\alpha:=\mathbb{E}[R^{d}]>0}$. Thus the logarithmic moment generating function of $R^{d}$ defined by ${\Lambda(\theta):=\log\left(\mathbb{E}\left[\exp({\theta R^{d}})\right]\right)}$ is finite for every $\theta\in\mathbb{R}$. Furthermore, the derivative

[TABLE]

is finite and positive for all $\theta\in\mathbb{R}$ and in particular, ${\Lambda^{\prime}(0)=\alpha}$. In fact, using the results in Chapter 2 of [10], it can be seen that $\Lambda(\theta)$ is strictly convex. As a consequence, $\Lambda^{\prime}(\theta)$ is strictly increasing and hence

[TABLE]

Let $\widehat{\theta}$ be such that ${\Lambda^{\prime}(\widehat{\theta})=\varrho}$ for some ${\varrho\in(\alpha_{\min},\alpha)}$. Therefore, ${\widehat{\theta}<0}$. Now consider the distribution $\widetilde{G}$ obtained by exponentially twisting $G$ by the amount $\widehat{\theta}$, that is, ${\mathrm{d}\widetilde{G}(t)}=\exp\left(\widehat{\theta}t-\Lambda(\widehat{\theta})\right)\mathrm{d}G(t).$ Fix a constant $a\in(0,1)$ and for each integer $n\geq 1$, define

[TABLE]

We later see that $a=1/2$ is a good choice for increasing performance of the algorithm. Let $\Lambda^{*}(\cdot)$ be the Legendre-Fenchel transform of $\Lambda$, that is, ${\Lambda^{*}(t)=\sup_{\theta\in\mathbb{R}}\{\theta t-\Lambda(\theta)\}}$. This corresponds to the large deviations rate function associated with the empirical average of iid samples from $G$. From the definition of $\widehat{\theta}$ and the fact that $\Lambda(\theta)$ is strictly convex, ${\Lambda^{*}(\varrho)=\widehat{\theta}\varrho-\Lambda(\widehat{\theta})>0}$. Since $\widehat{\theta}<0$, for all ${(t_{1},t_{2},\dots,t_{\lceil na\rceil})\in H_{n}}$,

[TABLE]

and thus,

[TABLE]

Recall the definition of the distribution $\mu$ of the hard-sphere model given by (1). To apply Algorithm 2, select $K=2$ and define

[TABLE]

and $D_{n,2}={\mathscr{G}}_{n}\setminus D_{n,1},$ for each $n$ where $H^{\mathsf{c}}_{n}$ is the complement of $H_{n}$ within $[0,\overline{r}]^{\lceil na\rceil}$. To apply Algorithm 5, we are now left with identifying the IS measures $\mu_{n,1}$ and $\mu_{n,2}$, and the corresponding bounds $\sigma_{n,1}$ and $\sigma_{n,2}$ for each ${n\in\mathbb{N}_{0}}$.

The measure $\mu_{n,1}$ on $D_{n,1}$ is again a grid based IS method similar to the grid method introduced for the fixed radius case in Subsection 4.4. First, iid copies $\{R_{1},\dots,R_{n}\}$ of $R$ are generated. Then, we construct a new grid and label each cell every time a new sphere is generated as follows. For the generation of the $i^{th}$ sphere with radius $R_{i}/\lambda^{\eta}$, we take the cell-edge length ${\varepsilon=1/\lceil\lambda^{\eta}/R_{i}\rceil}$. A cell $C$ in the grid is labeled as fully-blocked if ${C\subseteq S(X_{j},(R_{j}+R_{i})/\lambda^{\eta})}$ for an existing sphere ${j\leq i-1}$ with the center $X_{j}$ and the radius $R_{j}/\lambda^{\eta}$; otherwise, the cell is labeled as non-fully-blocked. A non-fully-blocked cell $C$ is called partially blocked if ${C\cap S(X_{j},(R_{j}+R_{i})/\lambda^{\eta})\neq\varnothing}$ for some ${j\leq i-1}$; otherwise, it is called non-blocked. Then the next center $X_{i}$ is generated uniformly over the non-fully-blocking cells. Just like in the case of fixed radius, $X_{1}$ is generated uniformly over $[0,1]^{d}$ and we check the possibility of the overlap of $i^{th}$ sphere with an existing sphere only if $X_{i}$ falls on a partially-blocked cell.

The measure $\mu^{0}$ is absolutely continuous with respect to $\mu_{n,1}$ on ${D_{n,1}\cap\mathscr{A}}$ and the associated likelihood ratio $L_{n,1}$ is given by

[TABLE]

where $\widehat{B}_{i}$ is the volume of all the fully-blocked cells for the $i^{th}$ sphere generation. By (11) and the fact that the cell-edge length is $1/\lceil\lambda^{\eta}/R_{i}\rceil$, we have ${B_{i}\geq\min\left(1,\frac{\gamma^{\prime}\lceil na\rceil}{\lambda^{\eta d}}\varrho\right)}$ on ${D_{n,1}\cap\mathscr{A}}$ for all ${i\geq\lceil na\rceil+1}$ because $\frac{1}{\lceil na\rceil}\sum_{j=1}^{\lceil na\rceil}R_{j}^{d}\geq\varrho$ over the set $H^{\mathsf{c}}_{n}$. Consequently,

[TABLE]

The measure $\mu_{n,2}$ is induced by the following procedure: Generate iid samples $R^{d}_{1},\dots,R^{d}_{\lceil na\rceil}$ from $\widetilde{G}$, and independently of this, generate iid samples $R^{d}_{\lceil na\rceil+1},\dots,R^{d}_{n}$ from $G$. For $i=1,\dots,n$, the radius of the $i^{th}$ sphere is $R_{i}/\lambda^{d}$ and the center generated uniformly distributed over the non-blocking region created by the existing $i-1$ spheres. Since $R^{d}_{1},\dots,R^{d}_{\lceil na\rceil}$ are sampled from $\widetilde{G}$, by (15),

[TABLE]

In summary, $\left\{\left(D_{n,k},\mu_{n,k},\sigma_{n,k}\right)_{k=1}^{2}\right\}_{n\in\mathbb{N}_{0}}$ is a stable IS sequence, and hence Algorithm 2 generates perfect samples from $\mu$. However, to reduce the per iteration complexity (as in the fixed radius case), we make some modification to the algorithm. Algorithm 6 is similar to Algorithm 4 and Algorithm 2 is restated as Algorithm 7.

We now focus on the running time complexity of Algorithm 7. Notice that $\widetilde{\sigma}(n)=\sigma_{n,1}+\sigma_{n,2}$ for each $n\in\mathbb{N}_{0}$. By Proposition 3, $P_{\mathsf{acc}}(\lambda)=\mathcal{P}(\lambda)/\mathbb{E}\left[\widetilde{\sigma}(N)\right]$ with $N\sim Poi(\lambda)$. Observe that $\sigma_{n,1}\leq\exp\left(-\frac{\gamma^{\prime}n^{2}\,a(1-a)}{\lambda^{\eta d}}\varrho\right)$. The proof of Proposition 4 can be extended to the current scenario to show that

[TABLE]

It is now clear that a good choice for $a$ is $1/2$ because it maximizes $a(1-a)$. Furthermore, using the moment generating function of Poisson random variables, we have

[TABLE]

Recall that ${\mathcal{T}_{\mathsf{ISAR}}\leq\mathbb{E}\left[\widetilde{\sigma}(N)\right]\widetilde{C}_{\mathsf{itr}}(\lambda)/\mathcal{P}(\lambda)}$. The per iteration complexity $\widetilde{C}_{\mathsf{itr}}(\lambda)$ mainly determined by relabelling of cells in the new grid for each sphere generation. The grid size for the $i^{th}$ sphere generation is an order of $\lambda^{\eta d}/R_{i}^{d}$ and the total number of spheres generated in each iteration is at most an order of $\lambda^{\min\{\eta d,1\}}$. Therefore, for $\eta d\leq 2$, we can show that $\widetilde{C}_{\mathsf{itr}}(\lambda)$ is of order $\lambda^{\min\{\eta d,1\}}\mathbb{E}[1/R^{d}]$.

Remark 6**.**

If $\varrho$ is selected to equal $\operatornamewithlimits{argmin}_{\varrho\in(\alpha_{\min},\alpha)}\,\,(\sigma_{n,1}+\sigma_{n,2})$ for each $n=1,2,\dots$, then $\mathbb{E}\left[\widetilde{\sigma}(N)\right]$ is minimum. Note that $\sigma_{n,1}$ decreases and $\sigma_{n,2}$ increases as functions of $\varrho$. The above decompositions were chosen to illustrate ideas simply. More complex decompositions are easily created for further performance improvement. For instance, we could have defined $H_{n}$ above as

[TABLE]

and then arrived at appropriate $\{\varrho_{m}\}_{m\leq n}$ and appropriate changes of measures for configurations in $H_{n}$ and $H_{n}^{\mathsf{c}}$. While this should lead to substantial performance improvement, it also significantly complicates the analysis.

5 Dominated CFTP Methods

In this section, we review some of the well-known dominated CFTP algorithms for the hard-sphere models. We refer to [27] for a general description of the dominated CFTP for Gibbs point processes (this method was first proposed for area-interaction processes by Kendall [25]).

Let $\mathbf{D}=\{\mathbf{D}(t):t\in\mathbb{R}\}$ be the so-called dominating birth-and-death process on $[0,1]^{d}$ with births arriving as a Poisson process with rate $\lambda$, where each birth is a uniformly and independently generated marked point on $[0,1]^{d}$ that denoting the center of a sphere with the mark being its radius. Each birth is alive for an independent random time exponentially distributed with mean one. It is well known that the steady-state distribution of $\mathbf{D}$ is $\mu^{0}$. Furthermore, it is easy to generate the dominating process $\mathbf{D}$ both forward and backward in time so that ${\mathbf{D}(t)\sim\mu^{0}}$ for all $t$. To see this, let ${\dots<t_{-2}<t_{-1}<0<t_{1}<t_{2}<\dots}$ be the event instants of the process $\mathbf{D}$, where an event can be either a birth or a death. Assume that with each birth there is an additional mean one exponentially distributed independent mark to determine its life time. Since the births are arriving as a Poisson process, the interarrival times are exponential with mean $1/\lambda$. Generate ${\mathbf{D}(0)\sim\mu^{0}}$, determine the next event instant $t_{1}$ and take ${\mathbf{D}(t)=\mathbf{D}(0)}$ for ${0\leq t<t_{1}}$. If the next event is a birth, generate a new independent (marked) point; otherwise, remove the existing point with the smallest lifetime. Continue the same procedure starting with $\mathbf{D}(t_{1})$ to generate the process over $[t_{1},t_{2})$, and so on. For generating the dominating process $\mathbf{D}$ backward in time, observe that $\mathbf{D}$ is time-reversible, and hence we can generate $\{\mathbf{D}(t):-T\leq t\leq 0\}$ for any finite ${T>0}$ just by generating an independent copy $\{\widetilde{\mathbf{D}}(t):\,0\leq t\leq T\}$ of the dominating process $\left\{\mathbf{D}(t):\,0\leq t\leq T\right\}$ and taking ${\mathbf{D}(-t)=\widetilde{\mathbf{D}}(t)}$ for ${0\leq t\leq T}$.

Since the distribution $\mu$ of the hard-sphere model is absolutely continuous with respect to $\mu^{0}$, using coupling, it is possible to construct a spatial birth-and-death process ${\mathbf{Z}=\{\mathbf{Z}(t):t\in\mathbb{R}\}}$, called the interaction process, such that ${\mathbf{Z}(t)\subseteq\mathbf{D}(t)}$ and ${\mathbf{Z}(t)\sim\mu}$ for all ${t\in\mathbb{R}}$; see [27]. Each iteration of any dominated CFTP method essentially involves the following two steps:

Fix $n>0$ and construct ${\{\mathbf{D}(t):t_{-n}\leq t\leq 0\}}$ backward in time starting at time zero with ${\mathbf{D}(0)\sim\mu^{0}}$ 2. 2.

Then, as detailed in Sections 5.1-5.3, use thinning on the dominating process $\{\mathbf{D}(t):t_{-n}\leq t\leq 0\}$ to obtain an upper bounding process ${\{\mathbf{U}_{n}(t):t\geq t_{-n}\}}$ with ${\mathbf{U}_{n}(t_{-n})=D(t_{-n})}$ and a lower bounding process ${\{\mathbf{L}_{n}(t_{-n}):t\geq t_{-n}\}}$ with ${\mathbf{L}_{n}(t_{-n})=\varnothing}$ forward in time such that for ${t\geq t_{-n}}$, ${\mathbf{L}_{n}(t)\subseteq\mathbf{Z}(t)\subseteq\mathbf{U}_{n}(t)\subseteq\mathbf{D}(t)}$ and ${\mathbf{L}_{m}(t)\subseteq\mathbf{L}_{n}(t)\subseteq\mathbf{U}_{n}(t)\subseteq\mathbf{U}_{m}(t)}$ for $m\leq n$.

If $\mathbf{U}_{n}$ and $\mathbf{L}_{n}$ coalescence at time [math], that is, ${\mathbf{U}_{n}(0)=\mathbf{L}_{n}(0)}$, then $\mathbf{U}_{n}(0)$ is a perfect sample from the target distribution $\mu$. If there is no coalescence, then repeat the steps by increasing $n$ and extending the dominating process ${\{\mathbf{D}(t):-t_{-n}\leq t\leq 0\}}$ further backward to time $t_{-n}$ and repeat the same procedure. It is well known that a good strategy for increasing $n$ is doubling it after every iteration. The criteria for thinning depends on the coupling used for constructing $\mathbf{Z}$. However, the dominating process $\mathbf{D}$ depends only on $\lambda$. In summary, a dominated CFTP algorithm is described by Algorithm 8.

Consider the backward coalescence time ${N^{*}=\min\left\{n\in\mathbb{N}_{0}:L_{n}(0)=U_{n}(0)\right\}}$. The average running time complexity of Algorithm 8 depends on the number of operations involved within $N^{*}$, which further depends on the construction of the interaction process and the bounding processes. At each iteration, the length of the dominating process $\mathbf{D}$ is doubled on average backwards in time. Hence, on average the running time complexity doubles at each iteration. From the definition of $N^{*}$, the length of the last iteration is $2^{\lceil\log_{2}N^{*}\rceil}\geq N^{*}$. Let ${N^{f}=\min\left\{n\in\mathbb{N}_{0}:\mathbf{L}_{0}(t_{n})=\mathbf{U}_{0}(t_{n})\right\}}$ be the forward coalescence time. Due to the reversibility of the dominating process, it can be shown that $N^{*}$ and $N^{f}$ are identical in distribution [3], and hence the expected computational effort for constructing the dominating, upper bounding and lower bounding processes up to the forward coalescence time $N^{f}$, starting from time [math], is a lower bound on the expected running time of the algorithm.

Below we consider three dominated CFTP methods applicable for the hard-sphere models.

5.1 Method 1

This method is based on [27]. Note that the Papangelou conditional intensity of the hard-sphere model is given by

[TABLE]

with the convention that $0/0=0$. The interaction process $\mathbf{Z}=\{\mathbf{Z}(t):t\in(-\infty,\infty)\}$ is constructed as follows: Suppose $x$ is a birth to $\mathbf{D}$ that sees $\mathbf{Z}$ in a state $\mathbf{x}\in{\mathscr{G}}$. Then $x$ is added to $\mathbf{Z}$ if and only if $\ell(\mathbf{x},x)=1$. Every death in $\mathbf{D}$ reflects in $\mathbf{Z}$, that is, if there is a death of a point $y$ in $\mathbf{D}$, then $y$ is removed from the process $\mathbf{Z}$ as well if it is present. It can be shown that $\mu$ is the unique invariant probability measure of $\mathbf{Z}$; see, e.g., [16] or [12].

For each $n\geq 1$, the bounding processes are constructed as follows: As mentioned earlier take ${\mathbf{L}_{n}(t_{-n})=\varnothing}$ and ${\mathbf{U}_{n}(t_{-n})=\mathbf{D}(t_{-n})}$. Suppose that ${\mathbf{x}^{l}=\mathbf{L}_{n}(t_{i})}$ and ${\mathbf{x}^{u}=\mathbf{U}_{n}(t_{i})}$ for ${-n\leq i<0},$ then assign $\mathbf{L}_{n}(t)=\mathbf{x}^{l}$ and $\mathbf{U}_{n}(t)=\mathbf{x}^{u}$ for $t_{i}<t<t_{i+1}$. In case it is a birth $x$ in the dominating process $\mathbf{D}$ at time $t_{i+1}$, set $\mathbf{L}_{n}(t_{i+1})=\mathbf{x}^{l}\cup\{x\}$ if $\ell(\mathbf{x}^{u},x)=1$; otherwise, it will remain unchanged, that is, $\mathbf{L}_{n}(t_{i+1})=\mathbf{x}^{l}$. Similarly, set $\mathbf{U}_{n}(t_{i+1})=\mathbf{x}^{u}\cup\{x\}$ if $\ell(\mathbf{x}^{l},x)=1$; otherwise, set $\mathbf{U}_{n}(t_{i+1})=\mathbf{x}^{u}$. Every death in the dominating process reflects in both lower and upper bound processes. Note that a birth is accepted by the lower bounding process if the resulting state of the upper bounding process is in $\mathscr{A}$. Similarly, a birth in $\mathbf{D}$ is accepted in the upper bounding process if the resulting state of the lower bounding process is in $\mathscr{A}$.

Theorem 2**.**

The expected running time complexity $\mathcal{T}_{\mathsf{DC1}}$ of the above dominated CFTP algorithm satisfies

[TABLE]

for some constant $c>0$. In particular,

[TABLE]

As highlighted by the numerical results in Section 6, the lower bound (17) is a loose bound, because the bound is established by considering the running time complexity only up to the time at which the lower bounding process receives its first arrival. This can be much smaller than the running time complexity until the coalescence of the upper and lower bound processes.

5.2 Method 2

This method is an improved version of Method 1, again based on [27]. Observe that at any given time $t\in\mathbb{R}$, the interaction process $\mathbf{Z}(t)$ can have only non-overlapping spheres. This suggests a better way of constructing the bounding processes that we describe now. For each $n\geq 1$, just like in Method 1, start with ${\mathbf{L}_{n}(t_{-n})=\varnothing}$ and ${\mathbf{U}_{n}(t_{-n})=\mathbf{D}(t_{-n})}$ to guarantee that ${\mathbf{L}_{n}(t_{-n})\subseteq\mathbf{Z}(t_{-n})\subseteq\mathbf{U}_{n}(t_{-n})}$. Suppose that the event at $t_{i+1}$ is an arrival of sphere $x$. Irrespective of whether $\mathbf{U}_{n}(t_{i})\in\mathscr{A}$ or not, if $x$ is not overlapping with any sphere in the upper bounding process $\mathbf{U}_{n}(t_{i})$, then it can not overlap with any sphere in $\mathbf{Z}(t_{i})$ and hence is accepted to $\mathbf{Z}$. Thus, we add $x$ to both the bounding processes. (Observe that in Method 1, such an $x$ is added to both the bounding processes only if $\mathbf{U}_{n}(t_{i})\in\mathscr{A}$ because of the Papangelou conditional intensity (16).) If $x$ overlaps with any sphere in the lower bounding process $\mathbf{L}_{n}(t_{i})$, then it must overlap with a sphere in $\mathbf{Z}(t_{i})$ as well, and hence it is not added to any of the bounding processes $\mathbf{L}$ and $\mathbf{U}$. If $x$ does not overlap with any sphere in $\mathbf{L}_{n}(t_{i})$, but does overlap with a sphere in $\mathbf{U}_{n}(t_{i})$, its presence in the process $\mathbf{Z}$ cannot be ruled out, and hence we keep it in the upper bounding process, but not in the lower bounding process. Finally, every death in $\mathbf{D}$ is reflected in both the bounding processes $\mathbf{L}_{n}$ and $\mathbf{U}_{n}$. Under this construction, the lower bounding process accepts births more often and hence the upper bounding process accepts births less often when compared with the construction in Section 5.1. As a result the running time of Method 2 is shorter than that of Method 1.

5.3 Method 3

A different approach for dominated CFTP for repulsive pairwise interaction processes has been proposed by Huber [22]. Here, we discuss main ingredients of the method for hard-sphere model; refer to [22] and [23] for more details. In this method, the interaction process $\mathbf{Z}$ is different from $\mathbf{Z}$ in Sections 5.1-5.2 and is known as spatial birth-death swap process whose invariant distribution is again the distribution $\mu$ of the hard-sphere model. In addition to births and deaths of spheres, this process also allows swap moves; here a swap move is an event where an existing sphere is replaced by an arrival if it is the only sphere that is overlapping with the arrival. The lower and upper bounding processes are constructed as follows: As usual let ${\mathbf{U}_{n}(t_{-n})=\mathbf{D}(t_{-n})}$ and ${\mathbf{L}_{n}(t_{-n})=\varnothing}$. For any ${0<k<n}$, if $t_{-k}$ is an instant of a death in the dominating process $\mathbf{D}(t)$ then the death is reflected in both the upper and lower bound processes. Now suppose that $x\in\mathbf{D}(t_{-k})$ is born at $t_{-k}$.

• Case 1: If no sphere in $\mathbf{U}_{n}(t_{-k})$ is overlapping with $x$, then the arrival sphere $x$ is added to both $\mathbf{U}_{n}(t_{-k})$ and $\mathbf{L}_{n}(t_{-k})$. If only one sphere $y$ in $\mathbf{U}_{n}(t_{-k})$ is overlapping with $x$, then $y$ is removed from $\mathbf{U}_{n}(t_{-k})$ (from $\mathbf{L}_{n}(t_{-k})$ if it is present) and $x$ is added to both $\mathbf{U}_{n}(t_{-k})$ and $\mathbf{L}_{n}(t_{-k})$.
• Case 2: There are at least two spheres in $\mathbf{L}_{n}(t_{-k})$ overlapping with $x$. Then $x$ is rejected by both $\mathbf{U}_{n}(t_{-k})$ and $\mathbf{L}_{n}(t_{-k})$.
• Case 3: There is at most one sphere in $\mathbf{L}_{n}(t_{-k})$ and at least two spheres in $\mathbf{U}_{n}(t_{-k})$ overlapping with $x$. Then $x$ is added to $\mathbf{U}_{n}(t_{-k})$ (but not to $\mathbf{L}_{n}(t_{-k})$). If $y\in\mathbf{L}_{n}(t_{-k})$ is the one that is overlapping with $x$, then remove $y$ from $\mathbf{L}_{n}(t_{-k})$.

6 Simulations

We compare the performance of all the methods discussed in this paper using numerical experiments, and illustrate the effectiveness of the proposed IS based AR method over certain regimes where the other methods fail to work. For this, we consider the torus-hard-sphere model with a fixed radius $\overline{r}/\lambda^{\eta}$ on $2$-dimensional square $[0,1]^{2}$. Thus, $\eta d=2\eta$. In the first two experiments, by fixing values of $\eta$ and $\overline{r}$, we estimate the complexities of the algorithms as a function of the intensity $\lambda$ of the reference PPP by computing a sample average of the number of spheres (or, circles in this case) generated per generation of a perfect sample of the hard-sphere model. Instead of estimating the expected running time complexities, we take this approach to keep the discussion independent of the underlying data structures and programming language used in the implementation of the algorithms. In addition, we estimate the non-overlapping probability $\mathcal{P}(\lambda)$ using the conditional Monte Carlo rare event estimation for Gilbert graphs proposed by [20]. The Gilbert graph under consideration is a random graph where the nodes constitute a $\lambda$-homogeneous PPP on $[0,1]^{2}$ and there is an edge between two points if they are within a distance of $2\overline{r}/\lambda^{\eta}$. Therefore, $\mathcal{P}(\lambda)$ is the probability that there are no edges in the graph. The codes for all the methods discussed in this paper are available at https://github.com/saratmoka/PerfectSampling_HardSpheres.

For the implementation of the proposed IS based AR method, the grid is constructed using the cell-edge length $\epsilon=1/\lfloor\lambda^{\eta}/\overline{r}\rfloor$; see Section 4.4 for more details on the cell-edge length selection. The complexities of the algorithms are estimated using $1000$ samples. In the simulation results presented below, $\widehat{S}_{\mathsf{NAR}}$ and $\widehat{S}_{\mathsf{ISAR}}$ denotes the sample means of complexities of the naive AR and the IS based AR algorithms. Likewise, $\widehat{S}_{\mathsf{DCM1}}$, $\widehat{S}_{\mathsf{DCM2}}$ and $\widehat{S}_{\mathsf{DCM3}}$ are the corresponding estimates for the three dominated CFTP methods $1$, $2$, and $3$ presented in Section 5, respectively.

A standard software used for generating perfect samples of the hard-sphere model using the dominated CFTP is rHardcore(), which is a part of R package Spatstat that is available at https://spatstat.org/. Experiment 3 provides a perspective on the performance of the proposed method with respect to rHardcore() by comparing their expected running times as a function of $\overline{r}$ for a fixed $\lambda$. We note that rHardcore() does not support the torus-hard-sphere model. However, when selected ”expand=TRUE”, it reduces the boundary effects by generating a perfect sample on a larger window, and then clipping the result to the original window $[0,1]^{2}$.

Experiment 1: In this experiment, we consider the high density regime. Figure 3 compares the performance of all the algorithms for $\eta=0.5$ (that is $2\eta=1$) and $\overline{r}=1$ (this is identical to the regime where the underlying space is $[0,\sqrt{\lambda}]^{2}$ and the radius of each sphere is $\overline{r}$). This experiment suggests that the proposed IS based AR method can perform significantly better than every other method. To comprehend the rarity of the samples of the hard-sphere configurations under $\mu^{0}$, we plot $\log\mathcal{P}(\lambda)$ in Figure 4 and the expected intensity of the hard-sphere model in Figure 5.

Significance of the proposed IS method in the high density regime is more evident when $\eta=0.25$ (that is, $2\eta=0.5$) and $\overline{r}=1$. In this case, for values of $\lambda$ greater than $50$, almost all the times all the dominated CFTP algorithms terminated without giving an output. In particular, the rHardcore function terminated by producing the error: memory exhausted (limit reached?). On the other hand, the time taken (in secs) for generating $1000$ samples using the proposed method are $0.13,0.21,68.94$ and $271.72$, when $\lambda$ values are $50,100,200,300$ and $400$, respectively.

Experiment 2: In this experiment, we consider the low density regime. Figure 6 compares the performances of all the methods for $2\eta=1.5$ and $\overline{r}=0.5$ to illustrate the case where $1<\eta d<2$. As we can see, for large values of $\lambda$, the dominated CFTP methods $2$ and $3$ perform better than the other methods, including the proposed method. Figure 7 is a plot of $\log\mathcal{P}(\lambda)$ against $\lambda$ while Figure 8 is a plot of the intensity of the hard-sphere model against $\lambda$.

Experiment 3 Figure 9 compares the running times of the proposed IS based AR method and rHardcore() for generating $1000$ samples. The same computer is used for running both the softwares. Here, we vary $\overline{r}$ while fixing $\lambda=50$ and $2\eta=1$. Observe that for large values of $\overline{r}$ the density is higher, and the proposed method performs far better than the dominated CFTP. As we expect for this regime, as the radius $\overline{r}$ increasing, the intensity of the hard-sphere model is decreasing (Figure 11) while the rarity of the hard-sphere configurations under $\mu^{0}$ is increasing (Figure 10).

7 Conclusion

In this paper we considered the problem of perfect sampling for Gibbs hard-sphere models on $[0,1]^{d}$. We discussed the performance of the naive acceptance-rejection method and introduced importance sampling based enhancements to it. We also compared these methods to some of the popular coupling from the past based techniques prevalent in the existing literature. For the performance analysis and comparison (of expected running time complexity), we developed an asymptotic regime where the intensity $\lambda$ of the reference Poisson point process increased to infinity, while the (random) volume of each sphere is an order of $\lambda^{-\eta d}$ decreased to zero, for different regimes of $\eta d>0$. One main conclusion is that while the dominated coupling from the past methods perform better for $1<\eta d<2$ for large $\lambda$, our importance sampling based methods provide a significant improvement for $\eta d\leq 1$. Enroute, we established large deviations results for the probability that spheres do not overlap with each other when their centers constitute a Poisson point process. We also conducted numerical experiments to validate our asymptotic results.

The proposed importance sampling based acceptance-rejection methods rely on clever partitioning of the underlying configuration space and arriving at an appropriate change of measure on each partition. While we showed how this may be effectively conducted for hard-sphere models, further research is needed to develop effective implementations for perfect sampling from a broad class of Gibbs point processes.

Acknowledgement

SM acknowledges support of the Australian Research Council Centre of Excellence for Mathematical and Statistical Frontiers (ACEMS), under grant number CE140100049. SJ acknowledges support of the Department of Atomic Energy, Government of India, under project no. RTI4001. MM’s research is partly funded by the NWO Gravitation project NETWORKS, grant number 024.002.003.

Appendix A Proofs

The following lemmas are useful for proving Theorem 1, Proposition 2 and Proposition 4. Lemma 1 is a standard Chernoff bound for Poisson random variables and Lemma 2 is Hoeffding’s inequality for $U$-statistics [21].

Lemma 1** (Chernoff bound for Poisson).**

Let $N\sim Poi(\lambda)$. Then, for any $0<\epsilon<1$,

[TABLE]

Lemma 2** (Hoeffding, 1963).**

Suppose that $\xi_{1},\xi_{2},\dots,\xi_{n}$ are iid random variables and $g:\mathbb{R}^{k}\rightarrow[0,1]$ is a measurable function. Set

[TABLE]

for a positive integer $k\leq n$ (this is known as a $U$-statistics of order $k$). Then, for any $\epsilon>0$,

[TABLE]

The same estimate holds for \mathbb{P}\left(Z_{n}\leq{n\choose k}\Big{(}\mathbb{E}[g(\xi_{1},\xi_{2},\dots,\xi_{k})]-\epsilon\right)\Big{)}.

A.1 Proof of Theorem 1

Recall that $\lambda>0$, $\eta>0$, and $\frac{R_{1}}{\lambda^{\eta}},\dots,\frac{R_{n}}{\lambda^{\eta}}$ are the radii of $n$ spheres whose respective centers $Y_{1},\dots,Y_{n}$ are independently and uniformly generated on the $d$-dimensional unit cube $[0,1]^{d}$, where $R_{1},\dots,R_{n}$ are iid positive random variables bounded from above by a constant $\overline{r}$ and are independent of $Y_{1},\dots,Y_{n}$. Define $m_{i}:=\mathbb{E}\left[(R_{1}+R_{2})^{id}\right]$, for all $i=1,2,\dots$. Let $\mathcal{P}_{n}(\lambda)$ be the probability that these $n$ spheres do not overlap with each other. Since the number of spheres in a $\lambda$-homogeneous marked PPP on $[0,1]^{d}$ is a Poisson random variable with mean $\lambda$, the non-overlapping probability

[TABLE]

where $N\sim\mathsf{Poi}(\lambda)$.

In proving the theorem, we use the following lemmas which exploits the reference IS measure $\widetilde{\mu}_{n}$ introduced in Section 4.2. By (10) and the definition of $\mathcal{P}_{n}(\lambda)$,

[TABLE]

The following bound holds trivially,

[TABLE]

where the sum is taken to be zero when $i=1$. Let $\theta_{n,\lambda}=\frac{\gamma(2\overline{r})^{d}\,n}{\lambda^{\eta d}}$. We have the following upper and lower bounds on $\mathcal{P}_{n}(\lambda)$.

Lemma 3**.**

Under the above set-up,

[TABLE]

Furthermore, for any $\epsilon>0$,

[TABLE]

for any $n$ and $\lambda$ such that $\theta_{n,\lambda}<1$, where $N_{n,\lambda}$ is a function of $n,\lambda$ and $\overline{r}$, defined by (25) in the proof below, such that

[TABLE]

In particular, for the torus-hard-sphere model,

[TABLE]

Proof.

Lower Bound: To prove (21) notice that, by (A.1),

[TABLE]

using the Taylor’s expansion $\log(1-x)=-\sum_{j=1}^{\infty}x^{j}/j$ for $0\leq x\leq 1$. By Jensen’s inequality and (20),

[TABLE]

Again by Jensen’s inequality, $\left(\frac{1}{i-1}\sum_{l=1}^{i-1}\left(R_{l}+R_{i}\right)^{d}\right)^{j}\leq\frac{1}{i-1}\left(\sum_{l=1}^{i-1}\left(R_{l}+R_{i}\right)^{jd}\right)$, and thus

[TABLE]

We establish (21) using $\sum_{i=1}^{n}(i-1)^{j}\leq\int_{x=0}^{n}x^{j}\,dx=\frac{n^{j+1}}{j+1}$.

Upper Bound: Let $R_{(1)},R_{(2)},\dots,R_{(n)}$ be the order statistics of $R_{1},R_{2},\dots,R_{n}$. Since the non-overlapping probability $\mathcal{P}_{n}(\lambda)$ is independent of the order in which the spheres are generated, without loss of generality assume that the $i^{th}$ sphere has radius $R_{(i)}$. Let, for each $1\leq j\leq i-1$, $\widetilde{B}_{i}(j)$ be the volume of the blocked region for the $i^{th}$ sphere generation when the $(j+1)^{th},(j+2)^{th},\dots,(i-1)^{th}$ spheres are ignored, where $\widetilde{B}_{i}(0)=0$. We can think of $\widetilde{B}_{i}(j)-\widetilde{B}_{i}({j-1})$ as the blocking volume contributed by the $j^{th}$ sphere for the $i^{th}$ sphere. Under the new measure $\widetilde{\mu}_{n}$, the blocking volume seen by the $i^{th}$ sphere, $\displaystyle B_{i}=\sum_{j=1}^{i-1}\left(\widetilde{B}_{i}(j)-\widetilde{B}_{i}({j-1})\right).$ Consider the sets

[TABLE]

for $i\leq n$ and take $\bar{\mathscr{N}}^{(i)}:=\{1,2,\dots,i-1\}\setminus{\mathscr{N}}^{(i)}$. Using the inequality $1-x\leq e^{-x}$ and (A.1),

[TABLE]

where $Z_{n}=\sum_{i=1}^{n}\sum_{j=1}^{i-1}\left(\frac{R_{(j)}+R_{(i)}}{2\overline{r}}\right)^{d},$ and the last inequality holds from the assumption that each $R_{i}\leq\overline{r}$. Since $R_{(i)}$ is non-decreasing with $i$, from the definition of ${\mathscr{N}}^{(i)}$, it is easy to see that $|\bar{\mathscr{N}}^{(i)}|$ is non-decreasing with $i$. Therefore,

[TABLE]

We now show that $|\bar{\mathscr{N}}^{(n)}|$ is stochastically bounded by a binomial random variable, and as a consequence, the conditional expectation \mathbb{E}_{\widetilde{\mu}_{n}}\left[\exp\left(\theta_{n,\lambda}|\bar{\mathscr{N}}^{(n)}|\right)\Big{|}R_{1},\dots,R_{n}\right] is uniformly bounded from above by a constant, which is a function of $n,\lambda$ and $\overline{r}$. Let

[TABLE]

Clearly, $q_{j}$ is increasing with $j$, and therefore $q_{j}\leq q_{n-1}$ for all $j\leq n-1$. This implies that $|\bar{\mathscr{N}}^{(n)}|$ is stochastically bounded by a binomial random variable with parameters $n$ and $q_{n-1}$, and thus

[TABLE]

Due to the boundary effect, $q_{n-1}$ is not the same for the torus model (where boundary spheres loop over to the opposite boundaries) and the Euclidean model. Observe that, for the Euclidean model, $\widetilde{B}_{n}(n-1)-\widetilde{B}_{n}(n-2)<\frac{\gamma(R_{(j)}+R_{(n)})^{d}}{\lambda^{\eta d}}$ if either

• (1)

the center of the $(n-1)^{th}$ sphere is within $(R_{(j)}+R_{(n-1)}+2R_{(n)})/\lambda^{\eta}$ distance from the center of $j^{th}$ sphere for some $j\leq n-2$, or

• (2)

the center of $(n-1)^{th}$ sphere is within $(R_{(n-1)}+R_{(n)})/\lambda^{\eta}$ distance form the boundary of the unit cube.

Note that the boundary event (2) is irrelevant for the torus-hard-sphere model. The probability of the event (1) is maximized by $\displaystyle\frac{\gamma}{\lambda^{\eta d}}\sum_{j=1}^{n-2}\frac{\left(R_{(j)}+R_{(n-1)}+2R_{(n)}\right)^{d}}{1-B_{n-1}}$, while that for the event (2) is maximized by $\displaystyle\frac{1-\left(1-2(R_{(n-1)}+R_{(n)})/\lambda^{\eta}\right)^{d}}{1-B_{n-1}}$. Since the $R_{i}^{\prime}$s are bounded from above by $\overline{r}$ and $B_{n-1}\leq\theta_{n,\lambda}$ (from (20)), we have

[TABLE]

for any $n$ and $\lambda$ such that $\theta_{n,\lambda}<1$, and for some constant $c$. Let

[TABLE]

then $\mathcal{P}_{n}(\lambda)\leq N_{n,\lambda}\,\mathbb{E}_{\widetilde{\mu}_{n}}\left[\exp\left(-\frac{\gamma\,(2\overline{r})^{d}}{\lambda^{\eta d}}Z_{n}\right)\right].$ Using the definition of $Z_{n}$, for any $\epsilon>0$,

[TABLE]

By Lemma 2 (with $k=2$), $\mathbb{P}_{\widetilde{\mu}_{n}}\left(Z_{n}<\frac{n(n-1)}{2}\left(\frac{m_{1}-\epsilon}{(2\overline{r})^{d}}\right)\right)\leq 2\exp\left(-\frac{(n-1)\epsilon^{2}}{(2\overline{r})^{2d}}\right)$, and thus (22) is established.

It remains to prove (23) and (24) under the assumption that $\eta d>1$. For this case, we have that $\lim_{\lambda\nearrow\infty}\theta_{\lambda,\lambda}=0$ and hence $\lim_{\lambda\nearrow\infty}\bar{q}_{\lambda,\lambda}=0$. Since

[TABLE]

using Taylor’s expansion of the exponential function,

[TABLE]

Thus (23) holds. In particular, for the torus model with $\eta d>3/2$,

[TABLE]

goes to [math] as $\lambda\nearrow\infty$, and hence (24) holds. ∎

Lemma 4**.**

Suppose that $1<\eta d\leq 2$. Then, for any $0<a<0.5$,

[TABLE]

Furthermore, let $\bar{\lambda}=\lceil\lambda(1-\frac{1}{\lambda^{a}})\rceil$ for some constant $a$ such that $0<a<\frac{\eta d-1}{2}$. Then, for any $\epsilon>0$,

[TABLE]

where $N_{\lambda,\lambda}$ satisfies (23) and (24). In particular, (26) holds with $\epsilon=0$ if $\eta d>5/3$ and $2-\eta d<a<\frac{\eta d-1}{2}$.

Proof.

Lower Bound: Fix $a$ such that $0<a<0.5$. Since $\mathcal{P}_{n}(\lambda)$ is a decreasing function of $n$ for any fixed $\lambda$, by Lemma 3, we can say that for all $n<\lambda\left(1+\frac{1}{\lambda^{a}}\right)$,

[TABLE]

and from (18) and the Chernoff bound for the Poisson variable $N$ (see Lemma 1),

[TABLE]

Now (26) easily because $\exp\left(-\frac{1}{3}\lambda^{1-2a}\right)=o(1)$ as a function of $\lambda$.

Upper Bound: From (18),

[TABLE]

where the last inequality holds due the fact that $\mathcal{P}_{n}(\lambda)$ is a decreasing function of $n$ for fixed $\lambda$. We now analyze $\mathcal{P}_{\bar{\lambda}}\left(\lambda\right)$ and $\mathbb{P}(N<\bar{\lambda})$ separately.

By Lemma 3, for any $\epsilon>0$,

[TABLE]

where we used the fact that $\bar{\lambda}\leq\lambda$. We rewrite the above expression as follows:

[TABLE]

Note that $\frac{\gamma\bar{\lambda}^{2}\,m_{1}}{2\lambda^{\eta d}}=O\left(\lambda^{2-\eta d}\right)$ and $\frac{\bar{\lambda}\,\epsilon^{2}}{(2\overline{r})^{2d}}=\Omega\left(\lambda\right)$. Since $\eta d>1$,

[TABLE]

and since $\lim_{\lambda\rightarrow\infty}\exp\left(\frac{\gamma\,m_{1}}{2\lambda^{\eta d-1}}\right)=1$, we can say that the first term $\mathcal{P}_{\bar{\lambda}}\left(\lambda\right)$ in (28) satisfies the following inequality,

[TABLE]

By Lemma 1,

[TABLE]

Since $2a<1$ (because $\eta d\leq 2$), we have $1-2a>2-\eta d,$ and hence using (30) and the fact that $N_{\lambda,\lambda}\geq 1$,

[TABLE]

and hence (27) follows from (28) and (31).

In particular if $\eta d>5/3$, we choose $a$ such that $2-\eta d<a<\frac{\eta d-1}{2}$. Let $\epsilon=1/\lambda^{a}$. Then (29) and (31) holds. We complete the proof using the fact that $\lim_{\lambda\rightarrow\infty}\exp\left(\frac{\gamma\bar{\lambda}^{2}\,\epsilon}{2\lambda^{\eta d}}\right)=1$. ∎

Proof of Theorem 1.

The following upper and lower bounds together complete the proof.

Lower Bounds: Consider the inequality (26).

Case: $\boldsymbol{\eta d>2}$. Since $R\leq\overline{r}$, we have $m_{j}\leq(2\overline{r})^{jd}$. Thus, for $\eta d>2$, all the terms in the exponent of the right-hand side of (26) go to zero asymptotically. In other words, for any $0<a<0.5$,

[TABLE]

That means, $\lim_{\lambda\rightarrow\infty}\mathcal{P}(\lambda)=1,\,\,\text{ for }\,\eta d>2.$

Case: $\boldsymbol{3/2<\eta d\leq 2}$. Using (26), $\mathcal{P}(\lambda)\exp\left(\frac{\gamma m_{1}}{2}\lambda^{2-\eta d}\right)$ is bounded from below by

[TABLE]

By fixing $a>2-\eta d$, we see that the right-hand side of the above expression goes to one as $\lambda\nearrow\infty$. Thus, $\liminf_{\lambda\rightarrow\infty}\left[\mathcal{P}(\lambda)\exp\left(\frac{\gamma m_{1}}{2}\lambda^{2-\eta d}\right)\right]\geq 1.$

Case: $\boldsymbol{1<\eta d\leq 3/2}$. By applying $\log$ on both the sides of (26), we have for any $0<a<0.5$ that

[TABLE]

and see that

[TABLE]

Thus $\liminf_{\lambda\rightarrow\infty}\frac{1}{\lambda^{2-\eta d}}\log\mathcal{P}(\lambda)\geq-\frac{\gamma m_{1}}{2}$ for $\eta d>1$.

Case: $\boldsymbol{0<\eta d\leq 1}$. Configurations with one sphere or no sphere is always accepted, that is, $\mathcal{P}_{1}(\lambda)=\mathcal{P}_{0}(\lambda)=1$. The probability of generating no sphere is $e^{-\lambda}$. Consequently, $\mathcal{P}(\lambda)>e^{-\lambda}$ and for any $\eta d>0$,

[TABLE]

In particular, assume that $\eta d=1$. For this case, first we show that the limit $\delta:=\lim_{\lambda\rightarrow\infty}\left[\frac{1}{\lambda}\log\mathcal{P}(\lambda)\right]$ exists. To prove this, partition the cube $[0,1]^{d}$ into a cubic grid of cell-edge length $x^{1/d}\in(0,1)$. Ignore the cells at the boundary that have the edge length strictly smaller than $x^{1/d}$. So, the total intensity of the underlying PPP over a cell is $\lambda x$.

When $\eta d=1$, the radius of each sphere is identical in distribution to $R/\lambda$. Observe that the non-overlapping probability of the spheres restricting to a cell is $\mathcal{P}(\lambda x)$ (see the definition of the non-overlapping probability). Since the total number of cells is at least $1/x$, the non-overlapping probability $\mathcal{P}(\lambda)\leq\left(\mathcal{P}(\lambda x)\right)^{\frac{1}{x}},$ and thus $\frac{1}{\lambda}\log\mathcal{P}(\lambda)\leq\frac{1}{\lambda x}\log\mathcal{P}(\lambda x).$ We can increase $\lambda$ and decrease the cell-edge length $x^{1/d}$ such that $y:=\lambda x$ is fixed. Then the following inequality holds

[TABLE]

Now the existence of the required limit is established by applying limit on $y$:

[TABLE]

To show that $\delta\nearrow 0$ as $\gamma m_{1}\searrow 0$, assume that $\gamma m_{1}<\epsilon$ for a constant $\epsilon\in(0,1)$. By (A.1) and (20),

[TABLE]

Consider the following partial order on $\mathbb{R}_{+}^{n}$: for any $y,\,y^{\prime}\in\mathbb{R}_{+}^{n}$, we say that $y\preceq y^{\prime}$ if $y_{i}\leq y_{i}^{\prime}$ for all $i=1,\dots,n$. A function $f:\mathbb{R}_{+}^{n}\rightarrow\mathbb{R}$ is called increasing (or decreasing) if $f(y)\leq f(y^{\prime})$ (or $f(y)\geq f(y^{\prime})$) for all $y,\,y^{\prime}\in\mathbb{R}_{+}^{n}$ such that $y\preceq y^{\prime}$. If $f$ and $g$ are either increasing or decreasing functions then Theorem 2.4 of [17] (FKG inequality) can be trivially extended to show that $\mathbb{E}[f(Y)g(Y)]\geq\mathbb{E}[f(Y)]\mathbb{E}[g(Y)]$. Clearly the following function $f_{i}$ is a decreasing function on $\mathbb{R}_{+}^{n}$:

[TABLE]

Therefore,

[TABLE]

Using the convexity of the function $x^{+}$ and Jensen’s inequality, for each $i$,

[TABLE]

and thus $\mathcal{P}_{n}(\lambda)\geq\prod_{i=1}^{n-1}\left(1-i\,\frac{\gamma m_{1}}{\lambda}\right)^{+}.$ With $\underline{\lambda}=\lfloor\lambda+\lambda^{0.75}\rfloor$ and $N\sim Poi(\lambda)$,

[TABLE]

By applying $\log$ on both the sides of the above inequality and scaling with $1/\lambda$,

[TABLE]

From the definition of $\underline{\lambda}$ and Lemma 1, the second term, $\frac{1}{\lambda}\log\mathbb{P}\left(N\leq\underline{\lambda}\right)$, goes to zero as $\lambda\nearrow\infty$. We now focus on the first term, $\frac{1}{\lambda}\log\mathcal{P}_{\underline{\lambda}}(\lambda)$. Since $\gamma m_{1}<\epsilon<1$, for all $i\leq\underline{\lambda}$,

[TABLE]

for large values of $\lambda$. Thus, we can write using Bernoulli’s inequality that

[TABLE]

for large values of $\lambda$. Therefore, by combining the trivial bound (32) and the above conclusions,

[TABLE]

Upper Bounds: We have a complete proof of the large deviation of $\mathcal{P}(\lambda)$ for the case $\eta d>2$. So, it remains to prove the theorem for $0<\eta d\leq 2$. We first prove the required upper bounds for the case $1<\eta d\leq 2$. If $0<a<0.5$ and $\bar{\lambda}=\lceil\lambda(1-\frac{1}{\lambda^{a}})\rceil$, then from Lemma 4, for any $\epsilon>0$,

[TABLE]

Case: $\boldsymbol{\eta d>1}$. By applying $\log$ on both the sides of (33) and then dividing by $\lambda^{2-\eta d}$, we see that

[TABLE]

As a consequence of Lemma 3, $\limsup_{\lambda\rightarrow\infty}\frac{1}{\lambda^{2-\eta d}}\log\mathcal{P}(\lambda)\leq-\frac{\gamma(m_{1}-\epsilon)}{2}.$ Now take $\epsilon\searrow 0$.

In particular, consider torus-hard-sphere model with $5/3<\eta d\leq 2$. We can fix $a$ such that $2-\eta d<a<\frac{\eta d-1}{2}$. From Lemma 4, (33) holds with $\epsilon=0$. Therefore,

[TABLE]

and hence $\limsup_{\lambda\rightarrow\infty}\left[\mathcal{P}(\lambda)\exp\left(\frac{\gamma m_{1}}{2}\lambda^{2-\eta d}\right)\right]\leq 1$ from Lemma 3.

Case: $\boldsymbol{0<\eta d<1}$. Let $\underline{\lambda}=\lfloor\lambda^{\frac{1+\eta d}{2}}\rfloor$ and $N\sim Poi(\lambda)$. From the definition

[TABLE]

For any $\epsilon>0$, let $H_{n}(\epsilon):=\left\{\frac{1}{n}\sum_{i=1}^{n}R_{i}^{d}>\epsilon\right\}$. From (11),

[TABLE]

where the second inequality holds because $\left(1-\frac{\gamma^{\prime}}{\lambda^{\eta d}}\sum_{j=1}^{n}R_{j}^{d}\right)^{+}=0$ on $H_{n}\left(\frac{\lambda^{\eta d}}{\gamma^{\prime}n}\right)$. Since $\eta d<(\eta d+1)/2<1$, see that $\frac{\lambda^{\eta d}}{\gamma^{\prime}\underline{\lambda}}\searrow 0$ as $\lambda\nearrow\infty$, and thus for every $\epsilon>0$ there exists $\lambda_{\epsilon}$ such that

[TABLE]

for all $\lambda>\lambda_{\epsilon}$ and $n>\underline{\lambda}$.

Suppose there is a constant $c>0$ such that $R\geq c$. Then for all sufficiently small values of $\epsilon$, $\mathbb{P}\left(H^{c}_{n}\left(\epsilon\right)\right)=0$ for all $n>\underline{\lambda}$. Thus for large values of $\lambda$, $\mathcal{P}(\lambda)\leq\mathbb{P}\left(N\leq\underline{\lambda}\right)\leq e^{-\lambda}\lambda^{\underline{\lambda}}$, and from the definition of $\underline{\lambda}$,

[TABLE]

So we can assume that $\mathbb{P}(R<\epsilon)>0$ for every $\epsilon>0$. Recall that $\mathbb{P}(R>0)=1$ and $\Lambda(\theta)$ is the logarithmic moment generating function of $R^{d}$. As a consequence of positivity of $R$, we see that $\Lambda(\theta)\searrow-\infty$ as $\theta\searrow-\infty$. Let $\Lambda^{*}(x)=\sup_{\theta\in\mathbb{R}}\left\{\theta x-\Lambda(\theta)\right\}$. As a consequence of the assumption that $\mathbb{P}(R<\epsilon)>0$ for every $\epsilon>0$, we can show $\Lambda^{*}(x)\nearrow\infty$ as $x\searrow 0$. From Theorem 2.2.3 of [10],

[TABLE]

for all $n>\underline{\lambda}$ and $\epsilon<\mathbb{E}[R_{1}^{d}]$, where the last inequality holds because $\Lambda^{*}(x)$ is non-decreasing over $0<x\leq\mathbb{E}[R_{1}^{d}]$. By (35),

[TABLE]

for all $\lambda\geq\lambda_{\epsilon}$. To conclude that $\limsup_{\lambda\rightarrow\infty}\frac{1}{\lambda}\log\mathcal{P}(\lambda)\leq-1$, see from the definition of Poisson distribution and $\underline{\lambda}$ that

[TABLE]

where we used the fact that $\lambda^{n-1}/(n-1)!<\lambda^{n}/n!$ for all $n<\lambda$. Hence,

[TABLE]

From the definition of $\underline{\lambda}$, see that $\frac{\underline{\lambda}}{\lambda}\log\lambda\searrow 0$ as $\lambda\nearrow\infty$. As a consequence, as $\lambda\nearrow\infty$,

[TABLE]

goes to zero. Therefore,

[TABLE]

We have the required result by taking $\epsilon\searrow 0$.

Case: $\boldsymbol{\eta d=1}$. It remains to show that $\delta\leq-\frac{1}{2}\left(1-\frac{1}{\gamma^{\prime}\overline{r}^{d}}\right)^{2}$ if $R\equiv\overline{r}$ and $\gamma^{\prime}\overline{r}^{d}>1$. Since, from (11), $\mathcal{P}_{n+1}(\lambda)\leq\prod_{i=1}^{n}\left(1-\frac{\gamma^{\prime}}{\lambda}i\overline{r}^{d}\right)^{+}=0,$ for all $n>\lambda\frac{\lambda}{\gamma^{\prime}\overline{r}^{d}},$ we have $\mathcal{P}(\lambda)\leq\mathbb{P}\left(N\leq\lambda\frac{\lambda}{\gamma^{\prime}\overline{r}^{d}}\right).$ Now the proof is complete by Lemma 1. ∎

A.2 Proof of Proposition 1

From [29], the intensity of the torus-hard-sphere model is given by

[TABLE]

where $\mathcal{P}_{n}(\lambda)$ is the non-overlapping probability of $n$ uniformly and independently generated spheres with radius $\overline{r}/\lambda^{\eta}$.

Case $\boldsymbol{\eta d>1}$.** In this regime, we show that $\rho(\lambda)$ is of the order of $\gamma\overline{r}^{d}\lambda^{1-\eta d}$. Using inequalities (11) and (20),**

[TABLE]

Therefore,

[TABLE]

and

[TABLE]

Consequently,

[TABLE]

**and thus $\lim_{\lambda\nearrow\infty}\frac{\mathsf{VF}(\lambda)}{\gamma\overline{r}^{d}\lambda^{1-\eta d}}=1$.

**

Case $\boldsymbol{\eta d\leq 1}$.** We know show that $\lim_{\lambda\nearrow\infty}\mathsf{VF}(\lambda)\leq\rho^{\max}\gamma$ with equality if and only if $\eta d<1$. Towards this end, we first consider another torus-hard-sphere model on $[0,\lambda^{\eta}/\overline{r}]^{d}$ with unit radius spheres and absolutely continuous with respect to a $\kappa$-homogeneous Poisson point process for some intensity $\kappa>0$. Let $\rho(\kappa,\lambda)$ be the intensity of this new hard-sphere model. We can easily see that when $\kappa=\overline{r}^{d}\lambda^{1-\eta d}$, the fraction of the volume occupied by the spheres in the new hard-sphere model is also $\mathsf{VF}(\lambda)$.**

The proof of Proposition 1 and 2 of [29] can be easily modified to show that $\rho(\kappa,\lambda)$ is strictly increasing in $\kappa$ for any fixed $\lambda>0$, and

[TABLE]

where $\rho^{max}$ is the closest packing density. On the other hand, by fixing $\kappa$,we can further using [29] show that the limit $\lim_{\lambda\to\infty}\rho(\kappa,\lambda)$ exists and is equal to the intensity of the stationary hard-sphere model on $\mathbb{R}^{d}$ with unit radius spheres and the reference PPP being $\kappa$-homogeneous. (In fact, the difference between $\rho(\kappa,\lambda)$ and the limit $\lim_{\lambda\to\infty}\rho(\kappa,\lambda)$ is known to be insignificant for large values of $\lambda$; see, for example, [5].)

From the above discussion, when $\eta d<1$, for sufficiently small $\epsilon>0$, there exist constants $\kappa_{\epsilon}$ and $\lambda_{\epsilon}$ such that $\rho(\kappa,\lambda)>\rho^{max}-\epsilon,$ for all $\lambda>\lambda_{\epsilon}$ and $\kappa>\kappa_{\epsilon}$. If we take $\kappa=\overline{r}^{d}\,\lambda^{1-\eta d}$, since $\eta d<1$,

[TABLE]

which is the maximum packing intensity.

Finally, if $\eta d=1$ and $\kappa=\overline{r}^{d}$, the limit $\lim_{\lambda\to\infty}\rho(\overline{r}^{d},\lambda)$ is strictly less than $\rho^{max}$. Hence, we have $\lim_{\lambda\to\infty}\mathsf{VF}(\lambda)<\rho^{max}\gamma$.

A.3 Proof of Proposition 2

Let $N^{\prime}$ be the number of spheres generated sequentially, independently and identically before seeing an overlap. Let $N\sim\mathsf{Poi}(\lambda)$ independently of $N^{\prime}$. Then from the construction of Algorithm 1,

[TABLE]

for some positive constants $c$ and $c^{\prime}$. In the above expression, $\log(\lambda)$ appears because the cost to generate a sample of $N$ is at most an order of $\log(\lambda)$ (see, e.g., ****[11]****). Observe that

[TABLE]

where the last equality follows from the fact that $\mathbb{P}(N^{\prime}>n)=\mathcal{P}_{n}(\lambda)$.

Upper bound: For $\eta d\geq 2$, since $\mathcal{P}_{n}(\lambda)\leq 1$, we can upper bound (37) by a constant times $\mathbb{E}[N^{2}]$, which is further bounded from above by a constant times $\lambda^{2}$. Therefore we just need to consider the case $\eta d<2$. From (A.1), {\mathcal{P}_{n}(\lambda)={\mathbb{E}}_{\widetilde{\mu}_{n}}\left[\prod_{i=1}^{n}\Big{(}1-B_{i}\Big{)}\right]}. As a consequence of (11),

[TABLE]

**Let $\alpha=\mathbb{E}[R_{1}^{d}]$. Since $\overline{r}$ is an upper bound on the **$R_{i}^{\prime}$s, by Hoeffding’s inequality (Lemma 2) on the sequence $\{R^{d}_{1}/\overline{r}^{d},\dots,R^{d}_{n}/\overline{r}^{d}\}$ with $\epsilon=\alpha/2\overline{r}^{d}$, $k=2$ and $g(x,y)=x$,

[TABLE]

Let $a=\sqrt{\frac{2\lambda^{\eta d}}{\gamma^{\prime}\alpha}}$. Then from the above expression,

[TABLE]

Select $\lambda$ large enough so that $b>0$. Then with $p=1-\exp(-\alpha^{2}/(4\overline{r}^{2d}))$, the second term on the right side of (38) is $1/p$ times $\mathbb{E}\left[Z\right]$ for a geometric random variable $Z$ with success probability $p$ and support $\{1,2,3,\dots\}$. Since $\mathbb{E}\left[Z\right]=1/p$, the term $\sum_{n=1}^{\infty}n\exp\left(-n\alpha^{2}/(4\overline{r}^{2d})\right)$ bounded from above by a constant.

On the other hand, since $n\exp\left(-(n-1)^{2}/(2a^{2})\right)\leq\int_{n-2}^{n-1}(x+2)\exp\left(-\frac{x^{2}}{2a^{2}}\right)$ for any non-negative integer $n$, we can write that

[TABLE]

where $Z$ is a Gaussian random variable with mean [math] and variance $a^{2}$. Since $\mathbb{E}[|Z|]=a\sqrt{2/\pi}$, using the definition of $a$, the first term in (38) is bounded from above by a constant times $\lambda^{\eta d}$. Thus, the required upper bound on (36) established as a consequence of (37) and (38).

Lower bound: Let $\epsilon^{\prime}=\min(1,\eta d/2)$. Then from (37),

[TABLE]

for a constant $c>0$. From (21),

[TABLE]

where $m_{j}=\mathbb{E}\left[(R_{1}+R_{2})^{jd}\right]$. Note that $m_{j}\leq(2\overline{r})^{jd}$ since $R\leq\overline{r}$. Therefore, $\displaystyle\mathcal{P}_{\lfloor\lambda^{\epsilon^{\prime}}/2\rfloor}(\lambda)\geq\exp\left(-\frac{\lambda^{\epsilon^{\prime}}}{2}\sum_{j=1}^{\infty}\frac{1}{j}\left(\frac{\gamma(2\overline{r}^{d})}{2\lambda^{\eta d-\epsilon^{\prime}}}\right)^{j}\right)$. Using Taylor’s expansion of $\log(1-x)$ for $0<x<1$, and the fact that $\frac{\gamma(2\overline{r}^{d})}{2\lambda^{\eta d-\epsilon^{\prime}}}<1$ for sufficiently large values of $\lambda$,

[TABLE]

From the definition of $\epsilon^{\prime}$,

[TABLE]

In addition, from Lemma 1, ${\lim_{\lambda\rightarrow\infty}\mathbb{P}\left(N>\lambda^{\epsilon^{\prime}}/2\right)=1}$. Therefore, there exists a constant $c$ such that ${C_{\mathsf{itr}}(\lambda)\geq c\,\lambda^{2\epsilon^{\prime}}}$. The proof of Proposition 2 is complete using Theorem 1 and (4).

A.4 Proof of Proposition 4

First note that the sphere volume is at most a constant time the cell volume for all $\lambda$. Thus, after generating a sphere, the complexity of relabelling cells around the center of the new sphere plus the complexity of overlap check is a constant. For $\eta d>1$, the number of spheres generated in an iteration of Algorithm 5 is stochastically dominated by a Poisson random variable with rate $\lambda$. Therefore, there exists a constant $c$ such that $\widetilde{C}_{\mathsf{itr}}(\lambda)\leq c\,\lambda$. On the other hand, if $0<\eta d\leq 1$, the expected number of spheres generated per iteration is of order $\lambda^{\eta d}$ because the expected volume of each sphere is an order of $1/\lambda^{\eta d}$. It is clear that there exists a constant $c>0$ such that $\widetilde{C}_{\mathsf{itr}}(\lambda)\leq c\,\lambda^{\eta d}$. Thus, by (9) and Proposition 2,

[TABLE]

Thus, (14) holds, because $\widetilde{\sigma}(n)=\delta_{n}$ for each $n\in\mathbb{N}_{0}$. Furthermore, from the definition of $\widetilde{\sigma}(\cdot)$ and $N$,

[TABLE]

By the Chernoff bound (Lemma 1), for any $0<\epsilon<1$,

[TABLE]

If $\eta d>1$, then the second term on the right-hand side of the above expression decreases faster than the first term, and thus the claim holds true. For $\eta d=1$, take $\epsilon=1/2$, then we have the required result with $b=\min\left(1/8,\gamma\overline{r}^{d}/4\right)$. Furthermore, if $0<\eta d<1$ then the first term decreases faster than the second one, and hence the proof is completed by taking $b=1/2$.

A.5 Proof of Theorem 2

To derive the lower bound on $\mathcal{T}_{\mathsf{DC1}}$, we view the entire dominating process $\mathbf{D}$ as a Poisson Boolean model on a higher dimensional space and use an extension of the FKG inequality ****[31]**** (alternatively, see Theorem 2.2 in ****[31]****). Let $s_{0}=0$ and $s_{i}$ be the instant of the $i^{th}$ arrival in the dominating process after time zero. Let $C(\mathbf{x},\mathbf{x}^{u},\mathbf{x}^{l})$ be the running time complexity of updating the dominating, upper bound and lower bound processes at the instant of an arrival when their respective states are $\mathbf{x},\mathbf{x}^{u}$ and $\mathbf{x}^{l}$.

Since $U_{0}(0)=D(0)$ and $L_{0}(0)=\varnothing$, on $\bigcap_{j=0}^{i}\{D(s_{j})\notin\mathscr{A}\}$, for all $t\leq s_{i},$

[TABLE]

Thus, $L_{0}(t)\neq U_{0}(t)$ for all $t\leq s_{i}$ on $\bigcap_{j=0}^{i}\{D(s_{j})\notin\mathscr{A}\}$. Now take

[TABLE]

From the above conclusion, it is clear that $N^{f}\geq\tau$. Then,

[TABLE]

where $I\left(\cap_{j=0}^{-1}\left\{D(s_{j})\notin\mathscr{A}\right\}\right)=1$ and the last equality follows from (39).

Suppose that $\mathscr{D}$ is the state space of the entire process $\left\{D(t):t\in\mathbb{R}\right\}$. Then we can define a simple partial order on $\mathscr{D}$ as follows: For any $\omega,\omega^{\prime}\in\mathscr{D}$, we say $\omega\preceq\omega^{\prime}$ if and only if every sphere present in $\omega$ is also present in $\omega^{\prime}$, that is, either $\omega^{\prime}=\omega$ or $\omega^{\prime}$ is obtained by adding spheres to $\omega$. Define the following notion of increasing functions: A real valued function on $\mathscr{D}$ is increasing if $f(\omega)\leq f(\omega^{\prime})$ for all $\omega,\omega^{\prime}\in\mathscr{D}$ such that $\omega\preceq\omega^{\prime}$.

At each arrival, if there are $n$ points in the upper bounding process, the cost to decide whether to accept the new point is at least the the cost to check overlap condition in the upper bounding process and that cost is an order of $n$. Therefore, C\Big{(}D(s_{i}),U_{0}(s_{i}),\varnothing\Big{)}=c|U_{0}(s_{i})| for some constant. Since $|U_{0}(s_{i})|$ is a non-decreasing function on $\mathscr{D}$ as per the partial order stated above, by FKG inequality (Theorem 2.2 in ****[31]****),

[TABLE]

is bounded from below by

[TABLE]

Thus,

[TABLE]

for some constant $c>0$. Then (17) follows from (4) and (5). The proof is completed using Theorem 1.