Uniform Sampling through the Lovász Local Lemma

The Lovász Local Lemma [9] is a classical gem in combinatorics that guarantees the existence of a perfect object that avoids all events deemed to be “bad”. The original proof is non-constructive but there has been great progress in the algorithmic aspects of the local lemma. After a long line of research [3, 2, 30, 8, 34, 37], the celebrated result by Moser and Tardos [31] gives efficient algorithms to find such a perfect object under conditions that match the Lovász Local Lemma in the so-called variable framework. However, it is natural to ask whether, under the same condition, we can also sample a perfect object uniformly at random instead of merely finding one.

Roughly speaking, the resampling algorithm by Moser and Tardos [31] works as follows. We initialize all variables randomly. If bad events occur, then we arbitrarily choose a bad event and resample all the involved variables. Unfortunately, it is not hard to see that this algorithm can produce biased samples. This seems inevitable. As Bezáková et al. showed [4], sampling can be NP-hard even under conditions that are stronger than those of the local lemma. On the one hand, the symmetric Lovász Local Lemma only requires $ep\Delta\leq 1$, where $p$ is the probability of bad events and $\Delta$ is the maximum degree of the dependency graph. On the other hand, translating the result of [4] to this setting, one sees that as soon as $p\Delta^{2}\geq C$ for some constant $C$, then even approximately sampling perfect objects in the variable framework becomes NP-hard.

The starting point of our work is a new condition (see Condition 5) under which we show that the output of the Moser-Tardos algorithm is in fact uniform (see Theorem 8). Intuitively, the condition requires any two dependent bad events to be disjoint. Indeed, instances satisfying this condition are called “extremal” in the study of Lovász Local Lemma. For these extremal instances, we can in fact resample in a parallel fashion, since the occurring bad events form an independent set in the dependency graph. We call this algorithm “partial rejection sampling”,¹¹1Despite the apparent similarity in names, our algorithm is different from “partial resampling” in [20, 21]. We resample all variables in certain sets of events whereas “partial resampling” only resamples a subset of variables from some bad event. in the sense that it is like rejection sampling, but only resamples an appropriate subset of variables.

Our result puts some classical sampling algorithms under a unified framework, including the “cycle-popping” algorithm by Wilson [39] for sampling rooted spanning trees, and the “sink-popping” algorithm by Cohn, Pemantle, and Propp [7] for sampling sink-free orientations of an undirected graph. Indeed, Cohn et al. [7] coined the term “partial rejection sampling” and asked for a general theory, and we believe that extremal instances under the variable framework is a satisfactory answer. With our techniques, we are able to give a new algorithm to sample solutions for a special class of $k$-CNF formulas, under conditions matching the Lovász Local Lemma (see Corollary 19), which is an NP-hard task for general $k$-CNF formulas. Furthermore, we provide explicit formulas for the expected running time of these algorithms (see Theorem 13), which matches the running time upper bound given by Kolipaka and Szegedy [26] under Shearer’s condition [35].

The next natural question is thus whether we can go beyond extremal instances. Indeed, our main technical contribution is a general uniform sampler (Algorithm 6) that applies to any problem under the variable framework. The main idea is that, instead of only resampling occurring bad events, we resample a larger set of events so that the choices made do not block any perfect assignments in the end, in order to make sure of uniformity in the final output.

As a simple example, we describe how our algorithm samples independent sets. The algorithm starts by choosing each vertex with probability $1/2$ independently. At each subsequent round, in the induced subgraph on the currently chosen vertices, the algorithm finds all the connected components of size $\geq 2$. Then it resamples all these vertices and their boundaries (which are unoccupied). And it repeats this process until there is no edge with both endpoints occupied. What seems surprising is that this simple process does yield a uniformly random independent set when it stops. Indeed, as we will show in Theorem 34, this simple process is an exact sampler for weighted independent sets (also known as the hard-core model in statistical physics). In addition, it runs in expected linear time under a condition that matches, up to a constant factor, the uniqueness threshold of the model (beyond which the problem of approximate sampling becomes NP-hard).

In the more general setting, we will choose the set of events to be resampled, denoted by ${\sf Res}$, iteratively. We start from the set of occurring bad events. Then we include all neighbouring events of the current set ${\sf Res}$, until there is no event $A$ on the boundary of ${\sf Res}$ such that the current assignment, projected on the common variables of $A$ and ${\sf Res}$, can be extended so that $A$ may happen. In the worst case, we will resample all events (there is no event in the boundary at all). In that scenario the algorithm is the same as a naive rejection sampling, but typically we resample fewer variables in every step. We show that this is a uniform sampler on assignments that avoid all bad events once it stops (see Theorem 24).

One interesting feature of our algorithm is that, unlike Markov chain based algorithms, ours does not require the solution space (or any augmented space) to be connected. Moreover, our sampler is exact; that is, when the algorithm halts, the final distribution is precisely the desired distribution. Prior to our work, most exact sampling algorithms were obtained by coupling from the past [32]. We also note that previous work on the Moser-Tardos output distribution, such as [19], is not strong enough to guarantee a uniform sample (or $\varepsilon$-close to uniform in terms of total variation distances).

We give sufficient conditions that guarantee a linear expected running time of our algorithm in the general setting (see Theorem 25). The first condition is that $p\Delta^{2}$ is bounded above by a constant. This is optimal up to constants in observance of the NP-hardness result in [4]. Unfortunately, the condition on $p\Delta^{2}$ alone does not make the algorithm efficient. In addition, we also need to bound the expansion from bad events to resampling events, which leads to an extra condition on intersections of bad events. Removing this extra condition seems to require substantial changes to our current algorithm.

To illustrate the result, we apply our algorithm to sample satisfying assignments of $k$-CNF formulas in which the degree of each variable (the number of clauses containing it) is at most $d$. We say that a $k$-CNF formula has intersection $s$ if any two dependent clauses share at least $s$ variables. The extra condition from our analysis naturally leads to a lower bound on $s$. Let $n$ be the number of variables. We (informally) summarize our results on $k$-CNF formulas as follows (see Corollary 30 and Theorem 32):

• •

  If $d\leq\frac{1}{6e}\cdot 2^{k/2}$, $dk\geq 2^{3e}$ and $s\geq\min\{\log_{2}dk,k/2\}$, then the general partial rejection resampling algorithm outputs a uniformly random solution to a $k$-CNF formula with degree $d$ and intersection $s$ in expected running time $O(n)$.

• •

  If $d\geq 4\cdot 2^{k/2}$ (for an even $k$), then even if $s=k/2$, it is NP-hard even to approximately sample a solution to a $k$-CNF formula with degree $d$ and intersection $s$.

As shown in the hardness result, the intersection bound does not render the problem trivial.

Previously, sampling/counting satisfying assignments of $k$-CNF formulas required the formula to be monotone and $d\leq k$ to be large enough [4] (see also [5, 28]). Although our result requires an additional lower bound on intersections, not only does it improve the dependency of $k$ and $d$ exponentially, but also achieves a matching constant $1/2$ in the exponent. Furthermore the samples produced are exactly uniform. Thus, if the extra condition on intersections can be removed, we will have a sharp phase transition at around $d=O(2^{k/2})$ in the computational complexity of sampling solutions to $k$-CNF formulas with bounded variable occurrences. A similar sharp transition has been recently established for, e.g., sampling configurations in the hard-core model [38, 36, 12].

Simultaneous to our work, Hermon, Sly, and Zhang [24] showed that Markov chains for monotone $k$-CNF formulas are rapidly mixing, if $d\leq c2^{k/2}$ for a constant $c$. In another parallel work, Moitra [29] gave a novel algorithm to sample solutions for general $k$-CNF when $d\lesssim 2^{k/60}$. We note that neither results are directly comparable to ours and the techniques are very different. Both of these two samplers are approximate while ours is exact. Moreover, ours does not require monotonicity (unlike [24]), and allows larger $d$ than [29] but at the cost of an extra intersection lower bound. Unfortunately, our algorithm can be exponentially slow when the intersection $s$ is not large enough. In sharp contrast, as shown by Hermon et al. [24], Markov chains mix rapidly for $d\leq c2^{k}/k^{2}$ when $s=1$.

While the study of algorithmic Lovász Local Lemma has progressed beyond the variable framework [22, 1, 23], we restrict our focus to the variable framework in this work. It is also an interesting future direction to investigate and extend our techniques of uniform sampling beyond the variable framework. For example, one may want to sample a permutation that avoids certain patterns. The classical sampling problem of perfect matchings in a bipartite graph can be formulated in this way.

Since the conference version of this paper appeared [18], a number of applications of the partial rejection sampling method have been found [17, 16, 15]. One highlight is the first fully polynomial-time randomised approximation scheme (FPRAS) for all-terminal network reliability [17]. For the extremal instances, tight running time bounds have also been obtained [14]. Moreover, partial rejection sampling is adapted to dynamic and distributed settings as well [10].

We first describe the “variable” framework. Let $\{X_{1},\dots,X_{n}\}$ be a set of random variables. Each $X_{i}$ can have its own distribution and range $D_{i}$. Let $\{A_{1},\dots,A_{m}\}$ be a set of “bad” events that depend on $X_{i}$’s. For example, for a constraint satisfaction problem (CSP) with variables $X_{i}$ ($1\leq i\leq n$) and constraints $C_{j}$ ($1\leq j\leq m$), each $A_{j}$ is the set of unsatisfying assignments of $C_{j}$ for $1\leq j\leq m$. Let ${\sf var}(A_{i})$ be the (index) set of variables that $A_{i}$ depends on.

The dependency graph $G=(V,E)$ has $m$ vertices, identified with the integers $\{1,2,\ldots,m\}$, corresponding to the events $A_{i}$, and $(i,j)$ is an edge if $A_{i}$ and $A_{j}$ depend on one or more common variables, and $i\neq j$. In other words, for any distinct $i,j$, $(i,j)\in E$ if ${\sf var}(A_{i})\cap{\sf var}(A_{j})\neq\emptyset$. We write $A_{i}\sim A_{j}$ if the vertices $i$ and $j$ are adjacent in $G$. The asymmetric Lovász Local Lemma [9] states the following.

Theorem 1 has a condition that is easy to verify, but not necessarily optimal. Shearer [35] gave the optimal condition for the local lemma to hold for a fixed dependency graph $G$. To state Shearer’s condition, we will need the following definitions. Let $p_{i}:=\mathop{\mathrm{Pr}}\nolimits(A_{i})$ for all $1\leq i\leq m$. Let $\mathcal{I}$ be the collection of independent sets of $G$. Define the following quantity:

where $\textbf{p}=(p_{1},\ldots,p_{m})$. When there is no confusion we also simply write $q_{I}$ instead of the more cumbersome $q_{I}(\textbf{p})$. Moreover, to simplify the notation, let $q_{i}:=q_{\{i\}}$ for $1\leq i\leq m$. Note that if $I\notin\mathcal{I}$, $q_{I}=0$.

In particular, if the condition holds with $q_{\emptyset}>0$, then $\mathop{\mathrm{Pr}}\nolimits\left(\bigwedge_{i=1}^{m}\overline{A_{i}}\,\right)>0$.

Neither Theorem 1 nor Theorem 2 yields an efficient algorithm to find the assignment avoiding all bad events, since they only guarantee an exponentially small probability. There has been a long line of research devoted to an algorithmic version of LLL, culminating in Moser and Tardos [31] with essentially the same condition as in Theorem 1. The Resample algorithm of Moser and Tardos is very simple, described in Algorithm 1.

In [31], Moser and Tardos showed that Algorithm 1 finds a good assignment very efficiently.

Unfortunately, the final output of Algorithm 1 is not distributed as we would like, namely as a product distribution conditioned on avoiding all bad events.

In addition, Kolipaka and Szegedy [26] showed that up to Shearer’s condition, Algorithm 1 is efficient. Recall that $q_{i}:=q_{\{i\}}$ for $1\leq i\leq m$.

On the other hand, Wilson’s cycle-popping algorithm [39] is very similar to the Resample algorithm but it outputs a uniformly random rooted spanning tree. Another similar algorithm is the sink-popping algorithm by Cohn, Pemantle, and Propp [7] to generate a sink-free orientation uniformly at random. Upon close examination of these two algorithms, we found a common feature of both problems.

In other words, any two events $A_{i}$ and $A_{j}$ are either independent or disjoint. These instances have been noticed in the study of Lovász Local Lemma. They are the ones that minimize $\mathop{\mathrm{Pr}}\nolimits\left(\bigwedge_{i=1}^{m}\overline{A_{i}}\,\right)$ given Shearer’s condition (namely $\mathop{\mathrm{Pr}}\nolimits\left(\bigwedge_{i=1}^{m}\overline{A_{i}}\,\right)=q_{\emptyset}$). Instances satisfying Condition 5 have been named extremal [26].

We will show that, given Condition 5, the final output of the Resample algorithm is a sample from a conditional product distribution (Theorem 8). Moreover, we will show that under Condition 5, the running time upper bound $\sum_{i=1}^{m}\frac{q_{i}}{q_{\emptyset}}$ given by Kolipaka and Szegedy (Theorem 4) is indeed the exact expected running time. See Theorem 13.

In fact, when Condition 5 holds, at each step of Algorithm 1, the occurring events form an independent set of the dependency graph $G$. Think of the execution of Algorithm 1 as going in rounds. In each round we find the set $I$ of bad events that occur. Due to Condition 5, ${\sf var}(A_{i})\cap{\sf var}(A_{j})=\emptyset$ for any $i,j\in I$, i.e., $I$ is an independent set in the dependency graph. Therefore, we can resample all variables involved in the occurring bad events without interfering with each other. This motivates Algorithm 2.

We call Algorithm 2 the Partial Rejection Sampling algorithm. This name was coined by Cohn, Pemantle, and Propp [7]. Indeed, they ask as an open problem how to generalize their sink-popping algorithm and Wilson’s cycle popping algorithm. We answer this question under the variable framework. Partial Rejection Sampling differs from the normal rejection sampling algorithm by only resampling “bad” events. Moreover, Algorithm 2 is uniform only on extremal instances, and is a special case of Algorithm 6 given in Section 5, which is a uniform sampler for all instances.

In fact, Algorithm 2 is the same as the parallel version of Algorithm 1 by Moser and Tardos [31] for extremal instances. Suppose each event is assigned to a processor, which determines whether the event holds by looking at the variables associated with the event. If the event holds then all associated variables are resampled. No conflict will be created due to Condition 5.

In the following analysis, we will use the resampling table idea, which has appeared in both the analysis of Moser and Tardos [31] and Wilson [39]. Note that we only use this idea to analyze the algorithm rather than to really create the table in the execution. Associate each variable $X_{i}$ with an infinite stack of random values $\{X_{i,1},X_{i,2},\dots\}$. This forms the resampling table where each row represents a variable and there are infinitely many columns, as shown in Table 1. In the execution of the algorithm, when a variable needs to be resampled, the algorithm draws the top value from the stack, or equivalently moves from the current entry in the resampling table to its right.

Let $t$ be a positive integer to denote the round of Algorithm 2. Let $j_{i,t}$ be the index of the variable $X_{i}$ in the resampling table at round $t$. In other words, at the $t$-th round, $X_{i}$ takes value $X_{i,j_{i,t}}$. Thus, the set $\sigma_{t}=\{X_{i,j_{i,t}}\mid 1\leq i\leq n\}$ is the current assignment of variables at round $t$. This $\sigma_{t}$ determines which events happen. Call the set of occurring events, viewed as a subset of the vertex set of the dependency graph, $I_{t}$. (For convenience, we shall sometimes identify the event $A_{i}$ with its index $i$; thus, we shall refer to the “events in $S$” rather than the “events indexed by $S$”.) As explained above, $I_{t}$ is an independent set of $G$ due to Condition 5. Then variables involved in any of the events in $I_{t}$ are resampled. In other words,

Therefore, any event that happens in round $t+1$, must share at least one variable with some event in $I_{t}$ (possibly itself). In other words, $I_{t+1}\subseteq\Gamma^{+}(I_{t})$ where $\Gamma^{+}(\cdot)$ denotes the set of all neighbours of $I$ unioned with $I$ itself. This inspires the notion of independent set sequences (first introduced in [26]).

We adopt the convention that the empty list is an independent set sequence with $\ell=0$. Note that we allow $S_{\ell}$ to be $\emptyset$.

Let $M$ be a resampling table. Suppose running Algorithm 2 on $M$ does not terminate up to some integer $\ell\geq 1$ rounds. Define the log of running Algorithm 2 on $M$ up to round $\ell$ as the sequence of independent sets $I_{1},I_{2},\dots,I_{\ell}$ created by this run. Thus, for any $M$ and $\ell\geq 1$, the log $I_{1},I_{2},\dots,I_{\ell}$ must be an independent set sequence. Moreover, if Algorithm 2 terminates at round $T$, let $\sigma_{t}:=\sigma_{T}$ if $t>T$. Denote by $\mu(\cdot)$ the product distribution of all random variables.

We remark that Lemma 7 is not true for non-extremal instances (that is, if Condition 5 fails). In particular, Lemma 7 says that given any log, every valid assignment is not only reachable, but also with the correct probability. This is no longer the case for non-extremal instances — some valid assignments from the desired conditional product distribution could be “blocked” under the log $\mathcal{S}$. In Section 5 we show how to instead achieve uniformity by resampling an “unblocking” set of bad events.

In other words, let $\Sigma$ be the set of assignments that avoid all bad events. Let $U$ be the output of Algorithm 2. In the case that all variables are sampled from the uniform distribution, we have $\mathop{\mathrm{Pr}}\nolimits(U=\sigma)=\frac{1}{|\Sigma|}$, for all $\sigma\in\Sigma$.

In this section we give an explicit formula for the running time of Algorithm 2. We assume that Condition 5 holds throughout the section.

We first give a combinatorial explanation of $q_{I}$ for any independent set $I$ of the dependency graph $G$. To simplify the notation, we denote the event $\bigwedge_{i\in S}A_{i}$, i.e., the conjunction of all events in $S$, by $A(S)$.

For any set $I$ in the dependency graph, we denote by $p_{I}$ the probability $\mathop{\mathrm{Pr}}\nolimits_{\mu}(A(I))$ that all events in $I$ happen (and possibly some other events too). If $I$ is an independent set in the dependency graph, any two events in $I$ are independent and

Moreover, for any set $J$ of events that is not an independent set, we have $p_{J}=0$ due to Condition 5.

On the other hand, the quantity $q_{I}$ is in fact the probability that exactly the events in $I$ happen and no others do. This can be verified using inclusion-exclusion, together with Condition 5:

where $\mathcal{I}$ denotes the collection of all independent sets of $G$. Since the events $(\bigwedge_{i\in I}A_{i}\wedge\bigwedge_{i\notin I}\overline{A_{i}})$ are mutually exclusive for different $I$’s,

Moreover, since the event $A(I)$ is the union over $J\supseteq I$ of the events $(\bigwedge_{i\in J}A_{i}\wedge\bigwedge_{i\notin J}\overline{A_{i}})$, we have

Essentially the proof above is a delayed revelation argument. At each round $1\leq t\leq\ell-1$, we only reveal variables that are involved in $S_{t}$. Thus, at round $\ell$, we have revealed all variables that are involved in $\mathcal{S}$. With respect to these variables, the sequence $\mathcal{S}$ happens with probability $p_{\mathcal{S}}$. Condition 5 guarantees that what we have revealed so far does not interfere with the final output (cf. Lemma 7). Hence the final state happens with probability $q_{S_{\ell}}$.

We write $p_{\mathcal{S}}=\prod_{i=1}^{\ell}p_{S_{i}}$ for an independent set sequence $\mathcal{S}$ of length $\ell\geq 0$. Note the convention that $p_{\mathcal{S}}=1$ if $\mathcal{S}$ is empty and $\ell=0$. Then, Lemma 9 implies the following equality, which is first shown by Kolipaka and Szegedy [26] in the more general (not necessarily extremal) setting of the local lemma.

As a sanity check, the probability of all possible logs should sum to $1$ when $q_{\emptyset}>0$ and the algorithm halts with probability $1$. Indeed, by Corollary 10,

where $\mathcal{S}$ is an independent set sequence. In other words,

where $\mathcal{S}$ is an independent set sequence. This fact is also observed by Knuth [25, Page 86, Theorem F] and Harvey and Vondrák [23, Corollary 5.28] in the more general settings. Our proof here gives a combinatorial explanation of this equality.

Equation (10) holds whenever $q_{\emptyset}>0$. Recall that $q_{\emptyset}$ is the shorthand of $q_{\emptyset}(\textbf{p})$, which is

where $\mathcal{I}$ is the collection of independent sets of the dependency graph $G$.

Let $T_{i}$ be the number of resamplings of event $A_{i}$ and $T$ be the total number of resampling events. Then $T=\sum_{i=1}^{m}T_{i}$.

The quantity $\sum_{i=1}^{m}\frac{q_{i}}{q_{\emptyset}}$ is not always easy to bound. Kolipaka and Szegedy [26] have shown that when the probability vector p satisfies Shearer’s condition with a constant “slack”, the running time is in fact linear in the number of events in the more general setting. We rewrite it in our notations.