{'fraction_non_alphanumeric': 0.056818181818181816, 'fraction_numerical': 0.013356164383561644, 'mean_word_length': 3.867479921200182, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast and perfect sampling of subgraphs and polymer systems

We give an efficient perfect sampling algorithm for weighted, connected induced subgraphs (or graphlets) of rooted, bounded degree graphs. Our algorithm utilizes a vertex-percolation process with a carefully chosen rejection filter and works under a percolation subcriticality condition. We show that this condition is optimal in the sense that the task of (approximately) sampling weighted rooted graphlets becomes impossible in finite expected time for infinite graphs and intractable for finite graphs when the condition does not hold. We apply our sampling algorithm as a subroutine to give near linear-time perfect sampling algorithms for polymer models and weighted non-rooted graphlets in finite graphs, two widely studied yet very different problems. This new perfect sampling algorithm for polymer models gives improved sampling algorithms for spin systems at low temperatures on expander graphs and unbalanced bipartite graphs, among other applications.


Introduction
Sampling is a fundamental computational task: given a specification of a probability distribution on a (large) set of combinatorial objects, output a random object with the specified distribution or with a distribution close to the specified distribution. This task becomes challenging when the specification of the distribution is much more succinct than the set of objects, and one wants to sample using time and space commensurate with the specification. Fundamental examples include sampling from Markov random fields and probabilistic graphical models and sampling substructures of graphs. We will address both of these examples here and connect them in a new way.
We consider a natural sampling problem: given a bounded-degree graph G, sample a graphlet (a connected, vertex-induced subgraph) of G containing a fixed vertex r with probability proportional to an exponential in the size of the subgraph. That is, sample a graphlet S containing vertex r with probability proportional to λ |S| , where λ > 0 is a distribution parameter and |S| denotes the number of vertices in S. In the paper we are concerned with small values of λ, where the expected size of a sampled graphlet is much smaller than the size of the graph.
Sampling graphlets is an important task in data science, network analysis, bioinformatics, and sociology, as it allows us to gain information about massive graphs from small sections of it; see, e.g., [38,27,42,4]. A number of variants of the problem have consequently been studied, including sampling graphlets of a given size uniformly at random or sampling weighted graphlets of all sizes [47,5,36,14,9,10,44,1,43,46,11,8]. The variant we consider here, i.e., sample a graphlet S with probability proportional to λ |S| , arises as a key subroutine in recent sampling algorithms for spin systems (hard-core model, Ising model, Potts model, etc.) in the regime of strong interactions via polymer models described below in Section 1.2; see [31,12,41,24,6,35,34,13,16].
One major limitation of previous sampling algorithms for graphlets and polymer models (those in, e.g., [31,41,46,15,23], among others) is the use of exhaustive enumeration of graphlets of a given size; this requires restrictive parameter regimes or large polynomial running times, with the maximum degree ∆ of the graph appearing in the exponent of the polynomial. Here we design a fast perfect sampling algorithm for weighted graphlets based on a vertex percolation process combined with a rejection filter. This method bypasses the enumeration barrier and allows us to design perfect sampling algorithms for a number of applications, substantially improving upon existing algorithms in three ways: 1) our algorithms have considerably faster running times, with no dependence on ∆ in the exponent; 2) our algorithms return perfect, rather than approximate, samples from the desired distributions; and 3) our algorithms are conceptually simple and practical to implement.
Our algorithm proceeds as follows. First, run a vertex percolation process on the graph G beginning at vertex r in a breadth-first search manner, repeatedly adding each adjacent vertex to the graphlet with a carefully-chosen probability p. Once the percolation process terminates, the graphlet is accepted as the random sample with a certain probability that depends on the graphlet and rejected otherwise; if the graphlet is rejected, the algorithm restarts another percolation process from r. Because of the careful way we choose the percolation and rejection probabilities, we can prove the final accepted sample is drawn exactly from the desired distribution and the expected running time is bounded by a constant that depends only on λ and the maximum degree ∆.

Sampling rooted graphlets
Our key contribution is a new algorithm for perfectly sampling weighted graphlets containing a given vertex r. Formally, let G = (V, E) be a finite or infinite graph of maximum degree ∆. For r ∈ V , let S(G, r) be the set of all connected, vertex-induced subgraphs of G containing r. (The subgraph induced by U ⊆ V has vertex set U and includes all the edges of G with both endpoints in U .) We call r the root of G and the elements of S(G, r) graphlets rooted at r. For λ > 0 define the probability distribution ν G,r,λ on S(G, r) by ν G,r,λ (S) = λ |S| Z G,r,λ , where Z G,r,λ = Ŝ ∈S(G,r) The distribution is well defined when the normalizing constant Z G,r,λ , known as the partition function, is finite. This is the case for every graph of maximum degree ∆ and every r when λ is below the critical threshold: (see Lemma 1.3 below). We give an efficient perfect sampling algorithm for ν G,r,λ for λ < λ * (∆).
Theorem 1.1. Fix ∆ ≥ 3 and let λ < λ * (∆). There is a randomized algorithm that for any graph G = (V, E) of maximum degree ∆ and any r ∈ V outputs a graphlet distributed according to ν G,r,λ with expected running time bounded by a constant that depends only on ∆ and λ.
We assume a model of computation that allows for querying the adjacency list of a given vertex in a bounded degree graph in constant time, the standard model used in the study of sublinear algorithms [26]. We also assume access to a stream of perfectly random real numbers in [0, 1]. The model of computation is chosen for consistency; in particular, our methods extend to other models, only requiring to adjust the running time to account for any additional computational overhead. Previous algorithms to generate ε-approximate samples from ν G,r,λ (e.g., those in [31,46,15,23]) exhaustively enumerate all graphlets of size ≤ k, for some k that depends on the error parameter ε that describes how accurate the sample must be. This results in algorithms with (1/ε) O(log ∆) running times. Applications such as sampling from polymer models require multiple samples from ν G,r,λ and, consequently, have small error tolerance per sample; in particular, they require ε ≪ 1/n, which results in inefficient algorithms with overall running time n O(log ∆) . The algorithm in Theorem 1.1, on the other hand, is an exact sampler whose expected running time depends only on ∆ and λ and thus provides a significant advantage in applications as we detail below.
We also show that Theorem 1.1 is sharp in two ways. First, we establish that there is no polynomial-time approximate sampling algorithm for ν G,r,λ when λ > λ * (∆) for the class of graphs of maximum degree at most ∆ unless RP=NP; see Definition 1 for the precise definition of a polynomial-time approximate sampler. Second, in the infinite setting, the normalizing constant Z G,r,λ may diverge (and consequently the distribution ν G,r,λ is not be well-defined) when λ > λ * (∆); conversely, we prove that Z G,r,λ is finite on every graph of maximum degree ∆ when λ ≤ λ * (∆). Lemma 1.2. If for every finite graph G = (V, E) of maximum degree ∆ ≥ 3 and every r ∈ V , there is a polynomial-time approximate sampler for ν G,r,λ when λ > λ * (∆), then RP=NP. Lemma 1.3. The partition function Z G,r,λ is finite for every (possibly infinite) graph G = (V, E) of maximum degree ∆ and every r ∈ V if and only if λ ≤ λ * (∆).
Finally, we mention that the algorithmic result in Theorem 1.1 cannot be extended even to the case λ = λ * (∆): for the infinite ∆-regular tree, we can show that the expected size of a graphlet sampled from ν G,r,λ is infinite when λ = λ * (∆), and so it is impossible to have sampling algorithms with finite expected running time. In summary, the algorithm in Theorem 1.1 for λ < λ * (∆), combined with the hardness/impossibility results in Lemmas 1.2 and 1.3 for λ > λ * (∆), provide a resolution to the computational problem of sampling from ν G,r,λ on graphs of maximum degree at most ∆.
As mentioned, our sampling algorithm is based on exploring the connected component of r in a vertex-percolation process. We carefully choose a specific percolation parameter p ∈ (0, 1) as a function of λ and ∆ (see Lemma 2.2). We then perform breadth-first search (BFS) from r, labeling each new vertex encountered 'active' with probability p and 'inactive' with probability 1 − p independently over all vertices; we continue the BFS exploring only unexplored neighbors of active vertices. In this way we uncover the 'active' component of r, call it γ. We then accept γ with a given probability depending on λ, ∆, |γ| and |∂γ|, where ∂γ denotes the set of vertices outside of γ that are adjacent to γ. If we reject γ, we begin again with a new percolation process. We note that only when λ < λ * (∆) there exists a suitable percolation probability p that results in a subcritical percolation process, so that the size of the active component has finite expectation and exponential tails. Sampling a graphlet by exploring a random component and performing a rejection step has been used in the past (most notably in the recent work of Bressan [8] to sample uniformly random graphlets of size k; see also [2]). The weighted model we sample from is particularly well suited to this type of exploration algorithm because of the direct connection to a subcritical percolation process.
We prove a more general version of Theorem 1.1 in Section 2, allowing for vertex-labeled graphlets and modifications of the weights by multiplication by a non-negative function bounded by 1. These generalizations are needed for the application to polymer models in Section 1.2.

Sampling from polymer models
We use our algorithm for sampling weighted rooted graphlets to design fast and perfect samplers for polymer models. Polymer models are systems of interacting geometric objects representing defects from pure ground states (i.e., most likely configurations) in spin systems on graphs in classical statistical physics [28,39,21]. These geometric objects are most often represented by vertex-labeled graphlets from a given host graph. Recently, polymer models have found application as an algorithmic tool to sample from spin systems on various classes of graphs in strong interaction regimes; see, e.g., [31,12,41,15,30,23,24,6,35,34,13,19,16]. In these applications, the problem of sampling weighted vertex-labeled rooted graphlets emerged as a significant computational barrier.
We will work with subset polymer models in which all polymers are vertex-labeled graphlets from a host graph G = (V, E). These models were defined in [28] and generalized in [39]. Such a polymer model consists of: • A set C = C(G) of polymers, each of which is a graphlet in G with the vertices of the graphlet labeled with colors from a set Σ of size q.
• An incompatibility relation ≁ defined by connectivity. We say two polymers γ, γ ′ ∈ C are incompatible and write γ ≁ γ ′ if the union of their corresponding vertices induces a connected subgraph in G. Otherwise they are compatible and we write γ ∼ γ ′ .
Let Ω(C) denote the set of all sets of pairwise compatible polymers from C. The polymer model is the Gibbs distribution µ on Ω(C) defined by is the polymer model partition function. We say the weights of a polymer model are computationally feasible if w γ can be computed in time polynomial in |γ|. The size |γ| of a polymer γ is the number of vertices in the corresponding graphlet. We will assume without loss of generality that all vertex-labeled graphlets of G, including each individual vertex v ∈ V , are elements of C, by setting w γ = 0 when necessary. We let C v be all polymers containing vertex v.
Algorithms for sampling polymer models fall into two classes: those based on truncating the so-called cluster expansion of a polymer model to approximate a partition function and using selfreducibility to sample, and those based on Markov chains on the set of collections of compatible polymers. The cluster expansion approach, while giving polynomial-time algorithms, generally is relatively inefficient, with the degree of the polynomial bound on the running time growing with the degree of the underlying graph; e.g., running time n O(log ∆) in n-vertex graph of maximum degree ∆. The Markov chain approach in principle can be much faster (near linear time in the size of the graph) but runs into one hitch: a stricter condition on the parameters of the model is needed to perform one update step of the Markov chain (the "polymer sampling condition" in [15,23]). We solve this problem by adapting our rooted graphlet sampler to sample polymers models, leading to a near linear-time perfect sampling algorithm for polymer models under the least restrictive conditions known (see Table 1.2). Theorem 1.4. Consider a subset polymer model on a family of n-vertex graphs of maximum degree ∆ with computationally feasible weights satisfying: Then, there is a perfect sampling algorithm for µ with expected running time O(n log n).
The threshold defined in (3) is the generalization of the critical threshold for rooted graphlet sampling to the labeled case (taking q = 1 recovers the definition in (2)). Theorem 1.4 improves upon the known results for sampling from polymer models in two ways. For a very general class of polymer models, our algorithm simultaneously provides perfect sampling with near-linear running time under weak conditions on the polymer weights. We now review previous works to illustrate these improvements; see the accompanying Table 1

.2.
A number of conditions on polymer weights have been used to provide efficient sampling algorithms. The first papers in this direction (including [31,34,12]) used the Koteckỳ-Preiss condition for convergence of the cluster expansion of the polymer model partition function [39]: γ ′ ≁γ w γ ′ e |γ ′ | ≤ |γ| ∀γ ∈ C. This condition is typically verified by ensuring that: γ≁v w γ e |γ| ≤ 1 ∀v ∈ V .
Since the number of vertex-labeled rooted graphlets of size k in a maximum degree ∆ graph grows roughly like (eq∆) k−1 (see [7]), weights of polymers of size k must decay roughly like (e 2 q∆) −k for the polymer model to satisfy (5), with the extra factor of e coming from the exponential in the left hand side of the condition (5). The major downside to the algorithms based on the cluster expansion, i.e., those using (5) or the Koteckỳ-Preiss condition, is that the running times obtained are of the form n O(log ∆) . Subsequent works, namely [15,23], addressed this downside but at the cost of a significantly stricter condition on the polymer weights.
In [15], the authors devised a new Markov chain algorithm for sampling from polymer models. The condition on the polymer weights for rapid mixing of this chain is somewhat less restrictive than the Koteckỳ-Preiss condition; it is the Polymer Mixing condition: This requires weights of polymers of size k to decay like (eq∆) −k , a savings of a factor e in the base of the exponent over (5). However, to implement a single step of this Markov chain in constant expected time, a stronger condition (the Polymer Sampling condition) was required: This is a significant loss of a factor e 3 ∆ 2 q 2 in the base of the exponent compared to (5), but the resulting sampling algorithm does run in near linear time.
In [22], the authors use a different Markov chain condition, the Clique Dynamics condition, similar to (6), which requires weights of polymers of size k to decay like (eq∆) −k , saving the same factor e over (5). Their running times, though, are again of the form n O(log ∆) since implementing one step of their Markov chain involves enumerating rooted polymers of size O(log n).
Our results are a "best-of-both-worlds" for polymer sampling: under the conditions (3) and (4) that both require polymer weights to decay like (eq∆) −k -the precise conditions are similar to but slightly less restrictive than the polymer mixing condition (6)-we obtain a near linear time algorithm. Moreover, unlike any of the previous results, our algorithm is a perfect sampler.
To conclude this section, we comment briefly on the algorithm we design to sample from µ. Our starting point is the polymer dynamics Markov chain from [15]. We use it to implement a Coupling from the Past (CFTP) algorithm (see [45]). To do so efficiently (in terms of the number of steps of the Markov chain), we design a new "bounding Markov chain" for the polymer dynamics, a method pioneered in [32,29], and to implement each step of the Markov chain efficiently, we turn to our sampler for weighted rooted graphlets from Theorem 1.1.

Applications to spin systems
Our new algorithm for sampling subset polymer models can be used as a subroutine in essentially all previous applications of polymer models for spin system sampling at low temperatures, including those in [34,12,41,15,30,23,24,13,19,16]. This results in faster sampling algorithms under less restrictive conditions on model parameters in all those settings. As examples, we fleshed out here the details in two of these applications; more details are provided in Section 4.
Hard-core model on bipartite graphs. The hard-core model on a graph G is the probability distribution µ hc G on I(G), the set of all independent sets of G, with The complexity of approximate counting and sampling from µ hc G on bounded-degree graphs is well understood: there is a computational threshold at some λ c (∆), with efficient algorithms for λ < λ c (∆) [50,3,17,18] and hardness above the threshold (no polynomial-time algorithms unless NP=RP) [48,25,49]. However, on bipartite graphs, the complexity of these problems is unresolved and is captured by the class #BIS (approximately counting independent sets in bipartite graphs) defined by Dyer, Goldberg, Greenhill, and Jerrum [20]. Theorem 1.4 implies the existence of a fast perfect sampling algorithm for the hard-core model in a certain class of bipartite graphs called unbalanced bipartite graphs, considered in [12,22].
There is a perfect sampling algorithm for µ hc G running in expected time O(n log n) for n-vertex bipartite graphs G with bipartition (L, R), with maximum degree ∆ L in L, maximum degree ∆ R in R, and minimum degree δ R in R if Approximate sampling algorithms with large polynomial run times were previously given for this problem when 6λ∆ L ∆ R < (1 + λ) δ R /∆ L in [12] and when 3.3353λ∆ L ∆ R < (1 + λ) δ R /∆ L in [22]. Our result applies to a comparable parameter range: inequality (9) holds, for instance, when (1 + e)λ∆ L ∆ R < (1 + λ) δ R /∆ L , or when 3λ∆ L ∆ R < (1 + λ) δ R /∆ L and ∆ L < 6. More importantly, our algorithm is the first to achieve perfect sampling and near-linear running time.
Potts model on expander graphs. The Q-color ferromagnetic Potts model on a graph G = (V, E) is the probability distribution µ potts G on the set of Q-colorings of the vertices of G; i.e., {1, . . . , Q} V . Each Q-coloring σ is assigned probability µ potts G (σ) ∝ e βm(G,σ) , where m(G, σ) is the number of monochromatic edges of G under the coloring σ and β > 0 is a model parameter. When the parameter β is large, and G has some structure (e.g., G is an expander graph), typical configurations drawn from µ potts G are dominated by one of the Q colors; that is, there is phase coexistence in the model. This enables sampling using subset polymer models.
Recall that an n-vertex graph G = (V, E) is an α-expander if for all subsets S ⊆ V with |S| ≤ n/2, the number of edges in E with exactly one endpoint in S is at least α|S|. Corollary 1.6. Consider the Q-color ferromagnetic Potts model on an α-expander n-vertex graph of maximum degree ∆. Suppose Then there is a sampling algorithm with expected running time O(n log n) that outputs a sample σ with distributionμ so that μ − µ potts G tv ≤ e −Ω(n) . Previously, [15] provided a ε-approximate sample for µ P otts , so our algorithm applies to a larger range of parameters and removes the dependence on ε from the running time. We do not achieve perfect sampling in this application only because the subset polymer models used give approximations of µ potts G rather than describing µ potts G exactly.

Sampling rooted graphlets in finite graphs
As another application of our algorithm for sampling weighted rooted graphlets, we consider next the problem of sampling weighted unrooted graphlets in a finite graph. Given a finite graph G, let S(G) be the set of all graphlets of G. Define the distribution ν G,λ on S(G) by Read-McFarland andŠtefankovič [46] gave a polynomial-time approximate sampling algorithm for ν G,λ for the class of maximum-degree ∆ graphs when λ < λ * (∆) and prove that there is no such algorithm for λ ∈ (λ * (∆), 1) unless NP=RP 1 . We give a new algorithm for this problem, covering the entire λ < λ * (∆) regime, and improving on the result of [46] in two ways: (i) our running time is constant in expectation (with no dependence on n), while the running time of the ε-approximate sampler in [46] is n · (1/ε) O(log ∆) ; and (ii) our algorithm outputs a perfect sample instead of an approximate one (and thus the running time has no dependence on any approximation parameter).
Theorem 1.7. Fix ∆ ≥ 3 and let λ < λ * (∆). Then for the class of finite graphs of maximum degree ∆ there is a randomized algorithm running in constant expected time that outputs a perfect sample from ν G,λ . The expected running time is bounded as a function of ∆ and λ.
The algorithm we use for this theorem is a modification of the one for sampling rooted graphlets.
We pick a uniformly random v ∈ V , run the same BFS percolation exploration, and accept the connected component of v with an adjusted probability (to account for the fact that a graphlet can be generated from any of its vertices). The acceptance probability is bounded away from 0 and so the algorithm runs in constant expected time. As mentioned earlier, the ε-approximate sampling algorithm from [46] is based on the exhaustive enumeration of all subgraphs of size ≤ k, for some k that depends on ε. Our new algorithm entirely bypasses this enumeration barrier.

Graphlet sampling: algorithms
In this section we present our efficient perfect sampling algorithm for weighted, vertex-labeled graphlets containing a fixed vertex r from a maximum degree ∆ graph; in particular, in Section 2.1, we prove a generalized version of Theorem 1.1 from the introduction. We also provide in Section 2.2 our algorithm for sampling weighted graphlets (i.e., the unrooted, unlabeled case) and establish Theorem 1.7. Our hardness results, that is Lemmas 1.2 and 1.3, are proved later in Section 5.
Define the probability distribution ν G,r,λ on S(G, r, q) ∪ {∅} by setting where Z(G, r, λ) = γ ′ ∈S(G,r,q)∪{∅} w γ ′ . We assume that G, f , q and λ are such that Z(G, r, λ) is finite, so that this distribution is well defined. When q = 1 and f (γ) = ½(γ = ∅), ν G,r,λ corresponds exactly to the distribution defined in (1) over the unlabeled graphlets of G rooted at r. We consider the problem of sampling from ν G,r,λ ; this more general version of the sampling problem is later used as a subroutine for sampling polymer systems in Section 3. Let cf., (2). Our main algorithmic result for sampling colored rooted graphlets is the following.
Theorem 2.1. Suppose ∆ ≥ 3, λ > 0, and q ≥ 1 are such that λ < λ * (∆, q) and let a > 0 be a fixed constant. There is a randomized algorithm to exactly sample from ν G,r,λ for graphs G of maximum degree ∆ and functions f : ; this randomized algorithm has expected running time bounded by C · Z −1 G,r,λ , where C > 0 is a constant that depends only on q, λ, ∆ and a. Theorem 1.1 from the introduction corresponds to the special case when q = 1 and f (γ) = ½(γ = ∅) (in this case Z G,r,λ ≥ λ). Other mild assumptions on the function f , e.g., f (∅) = 1 or f (r) = 1, ensure that Z G,r,λ is bounded away from 0 and, consequently, that the sampling algorithm in the theorem has constant expected running time.
As a warm-up, let us consider first our algorithm for sampling labeled rooted graphlets on a finite graph G = (V, E) with f = 1, and purposely omit certain non-essential implementation details for clarity. First, we find p ∈ (0, 1) such that p q (1 − p) ∆−2 = λ; this choice of p will be justified in what follows. The algorithm then repeats the following process until a vertex-labeled graphlet is accepted: 1. Each vertex of the graph is independently assigned with probability p a uniform random color from {1, . . . , q}, or it is marked as "not colored" with the probability 1 − p.
2. Let γ be the vertex-labeled graphlet from S(G, r, q)∪{∅} corresponding the colored connected component of r; i.e., the set of vertices connected to r by at least one path of colored vertices.
3. Observe that the probability of obtaining γ is (p/q) |γ| (1 − p) |∂γ| , where ∂γ denotes to set of vertices in G that are not in γ but are adjacent to a vertex in γ (with a slight abuse of notation, we let |∂∅| = 1). Our aim is to output γ with probability ∝ λ |γ| which has no dependence on ∂γ. Therefore, we use a "rejection filter" and only accept γ with probability (1 − p) (∆−2)|γ|+2−|∂γ| , so that the probability that γ is the output becomes: From (12), the choice of p such that p q (1 − p) ∆−2 = λ is apparent. We will prove that only when In the actual implementation of the algorithm, it will in fact suffice to find an approximation for p.
We comment briefly on the intuition for the rejection filter. The acceptance probability must include a factor of (1 − p) −|∂γ| , so that the final acceptance probability depends on |γ| but not on |∂γ|. However, (1 − p) −|∂γ| > 1 is not a valid probability, so we use instead (1 − p) (∆−2)|γ|+2−|∂γ| , which is at most 1 since (∆ − 2)|γ| + 2 ≥ |∂γ|. This bound on |∂γ| is best possible since it is tight for the ∆-regular tree. We note that using looser bounds for |∂γ| affects the range of the parameter λ for which we can find p ∈ (0, 1) so that p q (1 − p) ∆−2 = λ. Finally, we mention that the algorithm as described requires Ω(|V |) time per iteration and cannot be extended to infinite graphs. This is easily corrected by assigning colors starting from r and revealing only the colored component of r in a breadth-first fashion. The threshold λ * (∆, q) is sharp in the sense that only when λ < λ * (∆, q) is the value of p such that the revealing process is a subcritical process that creates a small component with high probability. This ensures the algorithm can be implemented efficiently. In particular, we stress that our algorithm avoids exhaustively enumerating labeled graphlets, as done in previous methods [15].
Before giving the implementation details of our algorithm and proving Theorem 2.1, we consider the problem of finding p ∈ (0, 1) such that p It can be readily checked that the function g is continuous and differentiable in [0, 1], has a unique maximum at This implies that only when λ < λ * (∆, q), there exists a value of p ∈ [0, 1 ∆−1 ) such that g(p) = λ. In particular, when λ > λ * (∆, q), there is no value of p for which g(p) = λ and when λ = λ * (∆, q), the only possible value is p = 1 ∆−1 . The latter case would result in a critical percolation process, corresponding to the fact that the expected size of a graphlet from ν G,r,λ has no uniform upper bound in the class of graphs of maximum degree ∆; in fact, it is infinite on the ∆-regular tree. We can find a suitable approximation for p when λ < λ * (∆, q) via a simple (binary search) procedure.
The proof of this lemma appears after the proof of Theorem 2.1. We now prove Theorem 2.1, including giving a more detailed version of the algorithm outlined above that includes the previously omitted implementation details and allows for general functions f : Proof of Theorem 2.1. For ease of notation, let λ * = λ * (∆, q). Our algorithm to sample from ν G,r,λ when λ < λ * explores from r in a breadth-first manner and stops once it has revealed the colored connected component of r. It proceeds as follows: 2. Let Q be a queue. With probability 1 −p do not add r to Q; otherwise, assign r a color uniformly at random from {1, . . . , q} and add r to Q. Mark r as explored.
3. While Q = ∅, repeat the following: 2) For each unexplored neighbor w of v, with probability 1−p do not add w to Q; otherwise, assign w a color uniformly at random from {1, . . . , q} and add w to Q. Mark w as explored (regardless of whether it was added to Q or not).
4. Let γ be the vertex-labeled graphlet from S(G, r, q)∪{∅} corresponding the colored connected component of r. Accept γ with probability: 5. If γ is rejected, go to Step 2 and repeat.
The probability of obtaining γ ∈ S(G, r, q) ∪ {∅} in an iteration of the algorithm is: and thus the overall acceptance probability in an iteration is: Then, We next bound the expected running time of the algorithm. We claim first that expected running per iteration is at most a constant that depends only on a, ∆ and q. If γ is the configuration generated in an iteration, it is discovered in O(|γ|+|∂γ|) = O(|γ|) time and, by assumption, f (γ) can be computed in at most O(|γ| a ) time, for suitable a constant a > 0. Letμ the output distribution of Step 3 of the algorithm. Then, there exists a constant C = C(q, ∆) > 0 such that the expected running time of each iteration is at most: We show next that |γ| (underμ) is stochastically dominated by a random variable W = X + Y (i.e., |γ| ≺ W ), where X and Y are i.i.d. random variables corresponding to the cluster size of a homogeneous Galton-Watson tree with offspring distribution Bin (∆ − 1,p). To see this, first note that |γ| ≺ L, where L is the cluster size of a heterogeneous Galton-Watson tree, in which the root vertex has offspring distribution Bin(∆,p) and every other vertex has offspring distribution Bin(∆ − 1,p). This is because the branching process generating γ includes the root only with probabilityp (the root is always present in the Galton-Watson tree), and, in addition, it considers at most ∆ (from the root) or ∆ − 1 (from any other vertex) potential branches. In turn, we have that L ≺ X + Y , since we can couple the first ∆ − 1 branches of the root with X (starting at the root) and the remaining branch with Y (starting at the child of the root not coupled with X).
It is well-known that X and Y have finite moments when (∆ − 1)p < 1 (see, e.g., [33]). In particular, there exists a constant A = A(a, ∆,p) > 0 such that This together with (13) shows that the expected running time in each iteration of the algorithm is bounded by C · A. Now, let R be the number of times Steps 2-5 are repeated, let T be the overall the running time of the algorithm. Then: and the result follows.
We conclude this section with the proof of Lemma 2.2.

Sampling unrooted graphlets
We consider next the problem of sampling weighted graphlets from a finite graph G = (V, E) of maximum degree ∆; specifically, in this variant of the sampling problem we consider unrooted, unlabeled, weighted graphlets of G. Let S(G) be the set of all graphlets of G. We define the probability distribution ν G,λ on S(G) by setting where Z G,λ = S ′ ∈S(G)∪{∅} λ |S ′ | . The problem of (approximately) sampling from ν G,λ is quite natural. In [46], it was established that this problem is computationally hard when λ > λ * (∆) = (∆−2) ∆−2 (∆−1) ∆−1 ; an ε-approximate sampling algorithm was also given in [46] for the case when λ < λ * (∆) with running time n · (1/ε) O(log ∆) . We establish the following: Theorem 2.3. Suppose ∆ ≥ 3 and λ > 0 are such that λ < λ * (∆). There is a randomized algorithm to exactly sample from ν G,λ with O(1) expected running time for finite graphs G of maximum degree ∆.
Proof. For ease of notation, we set λ * = λ * (∆) throughout this proof. Our algorithm to sample from ν G,λ is based on the algorithm to sample from ν G,r,λ (the rooted, vertex-labeled, weighted case). The idea is to pick a root uniformly at random and run the algorithm for the rooted case from this random vertex with the rejection filter adjusted to account for the fact that a graphlet can be generated from any of its vertices. It proceeds as follows: 1. Findp ∈ [0, 1 ∆−1 ) andλ ∈ [λ, λ * ) such that g(p) =λ using the method from Lemma 2.2. 2. Pick a vertex r ∈ V uniformly at random.
3. Let Q be a queue. With probabilityp add r to Q and mark it as colored. Mark r as explored.
4. While Q = ∅, repeat the following: 2) For each unexplored neighbor w of v, with probabilityp add w to Q and mark w as colored. Mark w as explored.
5. Let S ∈ S(G) be the graphlet corresponding to the colored connected component of v. Accept S with probability: 6. If S is rejected, go back to Step 2 and repeat.
The analysis of this algorithm is similar to that in the proof of Theorem 2.1. Let n = |V |. The probability that the algorithm outputs S in an iteration is: Hence, conditioned on acceptance, the probability of obtaining S ∈ S(G) is thus ν G,λ (S), and so the output distribution of the algorithm is ν G,λ . For the running time of the algorithm, we note that Step 4 of the algorithm is analogous to Step 3 of the algorithm in the proof of Theorem 2.1, and so the expected running time of each round is also bounded by a constant C = C(∆,p) > 0. Let T be the overall the running time of the algorithm. From (16), we have that the overall acceptance probability in a round is ρ = (1−p) 2 Z(G,λ) n . Then, as in (15), we deduce that

Applications to Polymer Models
In this section, we show how to use our algorithm for sampling rooted vertex-labeled graphlets from Section 2 to sample from subset polymer models and prove Theorem 1.4.
Consider a subset polymer model on an n-vertex graph G = (V, E); see Section 1.2 for the definition. Recall that we use C v for the set of all polymers containing vertex v ∈ V , and let γ ∅ denote the empty polymer. Define the distribution ν v on C v ∪ {γ ∅ } by where we assign w γ ∅ = 1. The following Markov chain on Ω(C), introduced in [15], has stationary distribution µ and mixes rapidly in O(n log n) steps under the polymer mixing condition (6).
Polymer dynamics. Given a configuration X t ∈ Ω(C), form X t+1 as follows: 1. Pick v ∈ V uniformly at random and let S v = {γ ∈ X t : v ∈ γ} (note that S v is either empty or contains 1 polymer).
To implement a single update step, one must sample from ν v in Step 3. To do so efficiently, in [15] the much stricter polymer sampling condition (7) was required. Here we give a fast perfect sampler for ν v under a much weaker condition.
Theorem 3.1. Consider a subset polymer model on a family of n-vertex graphs of maximum degree ∆ ≥ 3 with computationally feasible weights that satisfy w γ ≤ λ |γ| for some λ < λ * (∆, q). There is a randomized algorithm to sample perfectly from ν v for any v ∈ V with expected running time bounded by a function of λ, ∆, and q.
Proof. This follows from Theorem 2.1.
Using this theorem and the fast mixing result for the polymer dynamics of [15], one can approximately sample from µ whenever both the polymer mixing condition (6) and the assumptions of Theorem 3.1 hold. We further improve this by giving a perfect sampling algorithm that works whenever a new condition (4) is satisfied (our algorithm also requires the assumptions in Theorem 3.1). In all known examples, condition (4) is more permissive than the polymer mixing condition (6) from [15].

Perfect Sampling for polymer systems: Proof of Theorem 1.4
As mentioned, the polymer dynamics from [15] is not a perfect sampler. We propose here a different algorithm to output a perfect sample from µ. Our algorithm uses the coupling from the past method [45] and the notion of bounding Markov chains [32,29] to efficiently implement it.
We proceed with the proof of Theorem 1.4. We start with the description of a grand coupling for the polymer dynamics, which is then used to implement a coupling from the past algorithm. For an n-vertex graph G = (V, E), let {X Γ t } denote an instance of the polymer dynamics started from the polymer configuration Γ ∈ Ω(C). For all Γ ∈ Ω(C), the chains {X Γ t } are coupled by choosing the same uniform random random vertex v ∈ V , the same polymer γ sampled from ν v , and the same uniform random number in [0, 1] to decide whether to remove S v (Step 2) or to add γ (Step 3). A coupling from the past algorithm will find a time −T such the grand coupling started from all possible states at time −T coalesces to a single state by time 0. This guarantees that the output of the algorithm, that is the state at time 0, has distribution µ (see Theorem 1 from [45]). Such a T can be found with a binary search procedure. Unfortunately, implementing the coupling from the past algorithm in this manner for the polymer dynamics Markov chain is infeasible in our setting, since it requires simulating an exponential number of copies of the polymer dynamics, one from each Γ ∈ Ω(C).
To work around this, we consider a bounding Markov chain for the polymer dynamics rather than the polymer dynamics chain itself. Bounding Markov chains were pioneered in [32,29] as a method for efficiently implementing coupling from the past. The bounding chain for the polymer dynamics has state space Ω(C) × 2 C and will be denoted by {B t , D t }, where B t ∈ Ω(C) and D t ⊆ C are sets of polymers. The chain will maintain throughout that all polymers in B t are compatible and that every polymer in B t is compatible with every polymer in D t . The polymers in D t do not need to be compatible with each other. A step of the bounding Markov chain is defined next.

With probability 1/2, remove all polymers containing v by setting
3. With the remaining probability 1/2, draw a sample γ according to ν v and: Observe that polymers are only added to B t if they are compatible with all other polymers in B t ; hence, if B 0 is a valid polymer configuration, so is B t for all t ≥ 0.
To implement a step of the polymer dynamics bounding chain it suffices to pick a vertex v ∈ V uniformly at random, a uniform random number in [0, 1], and a polymer γ from ν v , just like for the polymer dynamics. Hence, we can couple the evolution of {B t , D t } with the grand coupling of the polymer dynamics described earlier. If we set B 0 = ∅ and D 0 = C, it can be checked that for all Γ ∈ Ω(C) and all t ≥ 0 : Indeed, this holds initially for t = 0, and the grand coupling ensures that whenever a polymer is removed from X Γ t it is also removed from B t , and whenever a polymer is added to X Γ t it is also added to B t or D t . Consequently, {B t , D t } is a bounding chain for the polymer dynamics. In particular, the first time B t = B t ∪ D t , all instances X Γ t have necessarily coalesced to the same configuration. This bounding chain allows us to implement the coupling from the past algorithm efficiently, as follows.
Coupling from the Past. Set k = 1.
(A) For t = −2 k , −2 k + 1, . . . , −2 k−1 generate ρ t = (v t , γ t , r t ) by choosing v t ∈ V uniformly at random, r t ∈ [0, 1] uniformly at random, and by sampling γ t ∈ C vt from ν vt . This implementation of the coupling from the past algorithm provides a perfect sample from µ; see [45]. It remains for us to show that it can be efficiently implemented. For this, we show first that the expected number of steps of the polymer dynamics bounding chain throughout the execution of the algorithm is O(n log n). Afterwards, we will show how to implement steps so that they can be executed in amortized constant expected time.

Lemma 3.2. Suppose the subset polymer model satisfies condition (4). Then, the expected number of steps of the polymer dynamics bounding chain in the coupling from past algorithm is O(n log n).
Proof. Let us first consider an evolution of the bounding chain {B t , D t } from initial state {B 0 , D 0 } = {∅, C}. Let T be the first time B t = B t ∪ D t ; that is, the time until D t = ∅. We bound first the expected value of T . Let φ t = |(∪ γ∈Dt γ) ∩ V |, the number of vertices of G that are included in at least one polymer of D t . Note that φ t ≤ n. We analyze the expected value of φ t+1 − φ t .
With probability φ t /(2n), a vertex in a polymer of D t is picked and all polymers in D t containing that vertex are removed, in which case φ t+1 ≤ φ t − 1. (A polymer may also be removed from D t in Step 3(a), but this case is omitted from our calculations as it never increases φ t .) The potential φ t may increase in Step 3(b), when a polymer γ may be added to D t . Such a polymer is only added if it is incompatible with D t . The probability this occurs is at most: For each such polymer, the expected increase in the size of φ t is at most |γ|. We can calculate: By Condition (4), there exists a constant β < 1 such that γ ∼{v} |γ|ν v (γ) ≤ β. It follows that This shows that {φ t } is a stochastic process with variable multiplicative drift. Since T is also the first time φ t < 1, using standard hitting time estimates (see, e.g., Theorem 10 in [40]), we get: It remains to consider how to efficiently implement the steps of the polymer dynamics bounding chain. This is subtle because D t may initially contain an exponentially large number of polymers, and care is thus needed in how D t is represented and stored. We describe and analyze efficient data structures for B t and D t in the following lemma. Proof. Data Structure for B t : Note B t is a collection of compatible polymers. The operations to be performed on B t include addition of a polymer (Step 3a), deletion of any polymer containing a particular vertex (Step 2), and comparison of a polymer to B t to determine whether it is compatible with B t (Step 3). B t will be stored as a length n array, and we use B t to denote the state of the array at time t. Initially, B 0 has only 0 entries. Each entry may also have a pointer to a node of a doubly linked list. Each polymer in B t is stored as a doubly linked list. This requires O(n) space.
When γ is added to B t , for each v ∈ γ we set B t (v) = 1. We also create a doubly linked list L γ , which has a node for each v ∈ γ that notes the color v is assigned in γ and contains a forward pointer, a backward pointer, and a pointer back to B t (v). Each entry of B t for v ∈ γ also points to the corresponding node of L γ . Doing these updates to the data structure for B t takes time O(|γ|).
The polymer removal step only occurs in Step 2 when any polymers containing v are removed from B t . There is at most one such polymer γ in B t containing v. The other vertices in γ can quickly be found using the doubly linked list L γ and its pointers back to B t . We can set B t (w) = 0 for each w ∈ γ and remove L γ . This removal of γ from B t occurs in time O(|γ|).
To compare a polymer γ to B t to determine whether it is compatible with B t , one looks up every vertex w that is in or adjacent to γ in B t to see if it is part of a polymer of B t , that is, if B t (w) = 1. If at least one of these vertices w has B t (w) = 1, γ is incompatible with B t ; if all these w have B t (w) = 0, then γ is compatible with B t . This takes time O(∆|γ|), which because we assume ∆ is constant is O(|γ|).
Each of these operations can be performed in time O(|γ|), where γ is a polymer that was at some (possibly earlier) time step drawn from ν v for some v. By the equation (3), this means γ has constant expected size, so each operation takes constant time in expectation. This data structure requires O(n) space.
Data Structure for D t : Much more care is needed with how D t is stored and accessed. Initially D 0 = C, which contains an exponential number of polymers. The operations we need to perform on D t include deleting all polymers containing a single vertex v (Step 2), checking whether a polymer is compatible with D t (Step 3b), adding a polymer to D t (Step 3b), and checking whether a polymer γ is compatible with D t \ {γ} and if so deleting γ from D t (Step 3b).
The key observation is that once a particular vertex v of G is selected in Step1 and the coin flip is such that Step 2 is performed, no polymers containing v remain in C and the size of D t has been reduced dramatically. Because of this, it makes sense to keep track of which vertices in D t have been the subject of a deletion in Step 2 at least once. We let D * t be (the state at time t) of an array of length n with an entry corresponding to each vertex v of G, where D * t (v) = 0 if Step 2, deleting C v from D t , has been performed for v at least once, and D * t (v) = 1 if it has not. Initially, D * 0 (v) = 1 for all v. At any time step, we know any polymer γ where D * t (v) = 1 for all v ∈ γ is in D t . Because of this, initially D * 0 completely describes D 0 . Amidst these deletions, polymers are also added to D t , and a polymer may contain both vertices where D * t (v) = 0 and vertices where D * t (v) = 1. All polymers added to D t in Step 3b will be stored separately from D * t . The data structure used here will have the same idea as that for B t . There will be an array D t with an entry for each vertex where D t (v) gives the number of polymers containing v that have been added to D t . Adding a polymer γ to D t involves both incrementing the corresponding entries in D t and adding a doubly linked list L γ connecting the vertices of the polymer together. However, because polymers in D t need not be compatible, there may be more than one polymer in D t containing a given vertex v. Because of this, instead of each vertex being in at most one doubly linked list L γ , it may be in many such doubly linked lists. To maintain all such pointers for vertex v (two for each L γ where v ∈ γ), these pointers are themselves stored in another doubly linked list L v .
The necessary operations that must be performed on D t can be implemented as follows. Note that initially D * 0 (v) = 1 for all v, D 0 (v) = 0 for all v, there are no doubly linked lists L γ , and each doubly linked list L v is empty. These steps are presented in a logical order for explaining the data structures for D t , rather than in the order in which they occur in the description of the Polymer Dynamics Bounding Chain.
• Adding polymer γ to D t (Step 3b): For each v in γ, increment D t (v). Create a doubly linked list L γ , which has a node for each v ∈ γ that notes the color v is assigned in γ and containing a forward pointer and a backward pointer. For each v, this node containing the pointers of L γ is inserted at the front of the doubly linked list L v . This takes time O(|γ|).
• Deleting all polymers containing v from D t (Step 2): Set D * t (v) = 0 and D t (v) = 0. For each node in L v , we explore the corresponding linked list L γ using the pointers in the node. We delete every other node of L γ , possible in time O(|γ|) because for list L w also containing a node of L γ , the nodes before and after this node in L w can easily be connected to each other because L w is doubly linked. For each node in L w that is deleted, we also decrement D t (w). This takes time O( γ∋v |γ|).
• Checking if polymer γ is compatible with D t (Step 3a: For each w that is in or adjacent to γ, check whether D * t (w) = 0 and D t (w) = 0. If both are 0 for all such w, then γ is compatible with D t . If at least one is nonzero, then γ is not compatible with D t . This takes time O(∆|γ|).
• Checking whether a polymer γ is compatible with D t \ {γ} and if so, deleting γ from D t (Step 3a): We only need to check this case when γ is compatible with B t . We first check whether γ is compatible with D t . It it is, we add γ to B t as described above. If γ is compatible with D t then it cannot be in D t , so D t and D t \ {γ} are the same and no further steps are needed.
If γ is incompatible with D t , the next step is to check if it is also incompatible with D t \ {γ}.
We check if D * (v) = 0 for all v in or adjacent to γ, D t (v) = 0 for all v adjacent to γ, and D t (v) = 1 for all v ∈ γ; if at least one of these does not hold, γ cannot be compatible with D t \ {γ} and so we are not in Step 3a and we proceed to Step 3b. If all these hold, we check further to see whether all D t (v) = 1 for v ∈ γ because of the presence of the single polymer γ or due to the presence of other polymers, which we can detect by looking at the doubly linked list(s) connecting the vertices in γ, including the colors assigned to each vertex. If there are polymers other than γ here, either multiple smaller polymers or a polymer with the same vertices as γ but differently assigned colors, then γ is incompatible with D t \ {γ} and we are not in Step 3a. If instead we find exactly polymer γ, then γ is compatible with D t \ {γ} and we add γ to B t and delete it from D t as previously described. This takes time O(∆|γ|).
Because any polymer added to D t was drawn from the distribution ν v for some v, its expected size is constant. All of the implementations above, with the possible exception of deleting all polymers containing a single vertex, takes expected constant time. Setting D * t (v) = 0 takes constant time, but it may take longer to remove any polymers that were added as linked lists L γ . However, we can amortize the cost of deletion if we pay for the cost of deleting a polymer when we add it, as all doubly linked lists L γ must be added before they are deleted. Doing this amortization makes the cost of adding a polymer to D t be O(|γ|) + O(|γ|) = O(|γ|), while the amortized cost of deleting all polymers containing a vertex v is now O(1). Thus all necessary operations for D t can be performed in amortized expected constant time.
The total space used for these data structures after t steps of the algorithm is in expectation at most O(n + t): O(n) for B t , D * t and D t , and O(|γ|) for each polymer γ added to D t , which is constant in expectation. There are at most t polymers in D t , one added each step, so this is O(t) space at most in expectation.

Applications to low-temperature spin systems
We present in this section the details of the applications spin system sampling given in Section 1.3.

The hard-core model on unbalanced bipartite graphs
We start by proving Corollary 1.5 for which we use Theorem 1.4.
Proof of Corollary 1.5. Let G be an n-vertex bipartite graph with partite sets L and R. Suppose that the vertices in L have maximum degree ∆ L , and that the vertices in R have maximum degree ∆ R and minimum degree δ R .
We can sample from the hardcore model distribution µ hc G (see (8) for its the definition) by considering an auxiliary distribution µ, which is a distribution on subsets S ⊆ R. For a subset S ⊆ R with neighbors N (S) in L, this distribution is given by where Z = S⊆R (1 + λ) |L\N (S)| · λ |S| . One can sample from µ hc G by first sampling S ⊆ R according to µ, adding S to the independent set, and then for each vertex v ∈ L \ N (S) include v in the independent set with probability λ/(1 + λ) independently. This results in exactly the desired distribution µ hc G . To sample from µ we use a polymer model, which we define next. Let R 2 be the graph whose vertices are the vertices of R, where two vertices are adjacent if they are at distance two in G. Note the maximum degree in R 2 is ∆ R (∆ L − 1). We define a polymer γ to be a connected subset of R 2 , its neighborhood N (γ) to be all vertices in L adjacent to a vertex in γ, and its weight as Polymers will not be labelled; that is, q = 1. Compatible polymer configurations in R 2 correspond exactly to subsets S ⊆ R, and the weights w γ mean a compatible polymer configuration corresponding to subset S has probability exactly µ(S).
To sample from this subset polymer model in R 2 , we show that when condition (9) holds, the conditions (3) and (4) of Theorem 1.4 hold, implying the existence of a perfect sampling algorithm with expected running time O(n log n).
From (9) it also follows that .
Using Claim 5.2 (which precisely counts the number of rooted subtrees of size k of the ∆-regular tree), the number of polymers of size k that are incompatible with v is at most (∆ R (∆ L − 1) + 1)(e∆ R (∆ L − 1)) k−1 /k; recall that ∆ R (∆ L − 1) is the maximum degree of R 2 . It follows that Thus, condition (4) also holds, and Theorem 1.4 supplies the perfect sampling algorithm with the desired running time.

Potts Model on Expander Graphs
The Q-color Potts model on G at inverse temperature β is a distribution over all (not necessarily proper) Q-colorings of V . Let Ω G,Q be all colorings ω : V → [Q]. For a coloring ω ∈ Ω G,Q with m(G, ω) monochromatic edges, this distribution has where Z = Ω G,q e βm(G,ω) . For j ∈ [Q], let Ω G,Q,j be all colorings in Ω G,Q such that strictly more than n/2 vertices are assigned color j. Let Ω G,Q := j∈[Q] Ω G,Q,j and consider the distribution µ potts given by: It follows from [34] that, under suitable conditions, (as we detailed next) a perfect sample from µ potts is an e −n -approximate sample from µ potts .
Let µ potts j be the distribution µ potts conditioned on being close to the ground state that is entirely color j, that is, conditioned on being in Ω G,Q,j . Because all Q ground states are symmetric, one can sample from µ potts by first picking a uniformly random j ∈ [Q], and then sampling from µ potts j . Sampling from µ potts j can be done using the polymer model we define next.
For ω ∈ Ω G,Q,j , let Γ(ω) = {v ∈ V : ω(v) = j}. Consider a subset polymer model whose host graph is G where a polymer is a graphlet of G whose vertices have colors other than j (there are thus Q − 1 colors available to color the vertices of a polymer). There is a bijection between Potts configurations in Ω G,Q,j and compatible polymer configurations consisting of at most n/2 total vertices, where vertices assigned a color other than j are identified. A sampling algorithm for this polymer model can give a sampling algorithm for µ potts j , where if the polymer configuration produced has more than n/2 total vertices, that configuration is rejected and resampling occurs. As we will see, the probability of needing to resample is small. Using this polymer model representation, [34] gives an efficient ε-approximate sampling algorithm with for µ potts G , derived from an approximate sampling algorithm for µ potts j . This algorithm applies to all α-expanding graphs G with maximum degree ∆ whenever α > 0, ∆ ≥ 3, Q ≥ 2, and β > 4 log((Q − 1)∆)/α. However, it involves a polymer enumeration step, and so its runtime is (omitting the dependence on other parameters) of the form n O(log ∆) .
In [15], the authors take a Markov chain approach and give an ε-approximate sampling algorithm for µ potts , again via approximate sampling from µ potts j , on α-expanding graphs G with maximum degree ∆ whenever β ≥ 5+3 log((Q−1)∆) α . This has running time O(n log(n/ε) log(1/ε)). In both of these prior works, it must hold that ε ≥ Qe −n .
Our algorithm gives an even larger range of β in which O(n log n) sampling is possible, and removes the dependence on ε. It produces an exact sample for µ potts rather than an approximate sample, although this sample is still only an e −n -approximate sample for µ potts G .
Proof of Corollary 1.6. From the discussion above, it suffices to generate a perfect sample from µ potts j for any j ∈ [Q] using the subset polymer model described above. We will show that conditions (3) and (4)  Using Claim 5.2, which implies the number of polymers incompatible with v of size k is at most Thus, condition (4) holds, and Theorem 1.4 provides a perfect sampling algorithm with the desired running time.

Graphlet sampling: hardness
In this section we establish the sharpness of the threshold λ * (∆, q). In particular, we prove Lemmas 1.2 and 1.3 from the introduction. We first prove the following more general variant of Lemma 1.3, which corresponds to the q = 1 case.
Lemma 5.1. The partition function Z G,r,λ is finite for every (possibly infinite) graph G = (V, E) of maximum degree ∆ and every r ∈ V if and only if λ ≤ λ * (∆, q).
Proof. Let C v,k be the number of graphlets of G with k vertices that contain v. Then, since f ≤ 1 by assumption, we have Z G,r,λ = γ∈S(G,r,q)∪{∅} λ |γ| f (γ) ≤ k≥0 γ∈S(G,r,q)∪{∅}:|γ|=k Consider the infinite ∆-regular tree rooted at ρ, and let T k be the number of subtrees of size k rooted at ρ so that C v,k ≤ T k . We have the following closed form expression for T k . Therefore, We apply the ratio test for the summation on the right hand side and consider the limit: We have that Using the inequality h(n)e 1/(12n+1) ≤ n! ≤ h(n)e 1/12n where h(n) = √ 2πn n+1/2 e −n , and setting we get From a direct calculation, it can be checked that Hence, From this and (18) we deduce that lim k→∞ A k = (∆−1) ∆−1 (∆−2) ∆−2 , which implies that From the ratio test we can then conclude that the series converges when λ < λ * (∆, q) and diverges when λ > λ * (∆, q).
It remains for us to consider the λ = λ * (∆, q) case, where . We also have for k ≥ 2 that Therefore, . Hence, and so Z G,r,λ is finite when λ = λ * (∆, q).
We conclude the proof of Lemma 5.1 by providing a proof for Claim 5.2.

Proof of Claim 5.2. For positive integers a and b, let
The numbers A k (a, b) are a generalization of the Catalan numbers and satisfy the recurrence see [37]. It is known (see, e.g., Lemma 1 in [46]) that the number of subtrees of size k containing the root of the infinite (∆ − 1)-ary tree is which is also equal to A k (1, ∆ − 1). Hence, for k ≥ 2, we have Using (19), we see that and using (19) repeatedly we get for k ≥ 2 The claim also holds trivially for k = 1 and the result follows.
We proceed next with the proof of Lemma 1.2. Let us first formally define the notion of a polynomial-time approximate sampler. Definition 1. An algorithm A is a polynomial-time approximate sampler for a distribution µ, if for a given error parameter ε > 0, it returns a sample from a distribution µ A such that µ A − µ tv ≤ ε, and the running time of A is polynomial in the size of the input and 1/ε.
Proof. Let u 1 , u 2 , . . . , u n−1 be the vertices in V \ {v} in a fixed order we define later. Let G i be the graph that results from removing vertices u 1 , . . . , u i from G and let Z i = Z G i ,v,λ . We let Z 0 = Z G,v,λ and note that Z n−1 = 1 + λ. Then: For each p i , we will findp i such that with probability at least 1 − δ/n: This implies that and so lettingZ G,v,λ = (1 + λ) n−2 i=0 1 p i , we get that with probability at least 1 − δ by a union bound.
To obtain (20), we define the sequence of vertices u 1 , u 2 , . . . , u n−1 such that u i is a leaf in the breadth-first search tree of G i−1 rooted at v. This way, we can conveniently claim that 1 ≥ p i ≥ 1 2 . To see this, note that S ∈ S(G i+1 ) if and only if S ∈ S(G i ) : u i+1 ∈ S, so Now, if S ∈ S(G i ) and S contains u i+1 , then, since u i+1 is a leaf in the breadth-first search tree rooted at v, S \ {u i+1 } is a connected subgraph that contains v. That is, every S ∈ S(G i ) that contains u i+1 can be mapped to a unique subgraph in S(G i ) that does not contain u i+1 of size |S| − 1. Hence, since λ < 1, and so Combined with (21), this gives p i = Z i+1 Z i ≥ 1 2 . Therefore, to deduce (20), it suffices to findp i such thatp .
For this, we draw L samples from ν G i ,v,λ and let X j be the indicator random variable for the event that the j-th sample contains u i+1 . Lettingp i = 1 L L i=1 X j , we get from a Chernoff bound that and setting ρ = ε 16np i and L = 384 n 2 ε 2 log(2n/δ): ≤ 2e − Lρε 24n ≤ 2e − Lε 2 384n 2 = δ n .
In summary, we have provided algorithm that computesp i such that (20) holds with probability at least 1 − δ/n. The algorithm has sample complexity and running time poly(n, 1/ε, log(1/δ)).
We can now provide the proof of Lemma 1.2.
Proof of Lemma 1.2. It suffices to prove the result for f = 1. Let G = (V, E) be a finite graph of maximum degree ∆, and for ease of notation let λ * = λ * (∆, 1). We show that if there is a polynomial-time approximate sampler for ν G,v,λ for each v ∈ V (see Definition 1), then there is also a polynomial-time approximate sampler for for ν G,λ . In [46], it was shown that there is no polynomial-time approximate sampler for for ν G,λ when λ ∈ (λ * , 1) unless NP=RP, and thus our result follows 2 . First, given exact sampling algorithms for ν G,v,λ for each v ∈ V with poly(n) running times, and a desired accuracy parameterε ∈ (0, 1), we provide a sampling algorithm whose output distribution is withinε total variation distance from ν G,λ . The running time of the algorithm will be poly(n, 1/ε). We then modify the algorithm to use instead the polynomial-time approximate samplers for ν G,v,λ (instead of the exact samplers) and show that this has a negligible effect on its running time and output distribution.
The algorithm is the following: 1. Let ε = min{ 1 8 ,ε 6 − 1 4n }. For each v ∈ V , using the exact sampler for ν G,v,λ and the algorithm from Lemma 5.3, obtainZ G,v,λ such that, with probability at least 1 − 1 n 2 2 n , 2. Set t = 0 and choose a vertex v ∈ V with probabilityZ G,v,λ wZ G,w,λ ; 3. Using the exact sampler for ν G,v,λ , draw a sample S from ν G,v,λ ; 4. Accept and output S with probability 1/|S|; 5. If S is rejected, increase t, and if t ≤ 2n 2 go to Step 2 and repeat; otherwise output ∅.
First, we note that by Lemma 5.3, Step 1 can be implemented in poly(n, 1/ε). Steps 2 to 5 of the algorithm are repeated at most O(n 2 ) times, and each of those step takes poly(n) time. Therefore, the overall running time of the algorithm is poly(n, 1/ε) = poly(n, 1/ε). Moreover, (22) holds for every v ∈ V with probability at least 1 − 1 n2 n by a union bound. Let µ alg be the output distribution of the algorithm, and let µ * alg be the output distribution when (i) Step 1 succeeds in finding the required approximations and (ii) a subgraph is accepted in Step 4 before t > 2n 2 . We show next that µ * alg − ν G,λ tv ≤ε. LetZ = v∈VZ G,v,λ and let S = S(G) ∪ {∅}. The probability that the algorithm outputs S ∈ S in Step 4 in a give iteration is: Then, µ * alg (S) = φ(S) Ŝ ∈S φ(Ŝ) , and µ * alg − ν G,λ tv = Plugging this into (24) we deduce that: Let E 1 be the event that Step 1 succeeds in finding the approximations, let E 2 the event that the algorithm accepts and outputs a graphlet in Step 4, and let E be the event that both E 1 and E 2 occur. Then, the triangle inequality and (25) imply: µ alg − ν G,λ tv ≤ µ * alg − ν G,λ tv + µ * alg − µ alg tv ≤ 4ε + µ * alg − µ alg tv .
So it remains for us to bound µ alg (¬E 2 | E 1 ). For this, note that the probability that a subgraph is accepted in Step 4 in an iteration is: we get the acceptance probability in Step 4 in an iteration is at least provided E 1 occurs. Therefore, if X is geometric random variable with parameter 1/(2n), we have and plugging this bound into (27) and (26) we obtain: µ alg − ν G,λ tv ≤ 4ε + 3 2n .
We have established so far the hardness of exactly sampling from ν G,v,λ for λ ∈ (λ * , 1). The same reduction (i.e., algorithm) with minor adjustments works for approximate sampling. Suppose that in Step 3 of the algorithm, we instead generate T = 2n 2 samples from a distribution µ v such that µ v − ν G,v,λ tv ≤ ε n 2 . Let µ ⊗T v and ν ⊗T G,v,λ be the product distributions corresponding to T independent samples from µ v and ν G,v,λ respectively. We have Letμ alg be the output distribution of the algorithm when using samples from µ v in Step 3. Then: μ alg − ν G,λ tv ≤ μ alg − µ alg tv + µ alg − ν G,λ tv ≤ μ alg − µ alg tv + 4ε + 3 2n .
That is, we obtain a polynomial-time approximate sampler for ν G,λ which completes the reduction.