Partial Rejection Sampling

• Partial rejection sampling is an algorithmic technique that generates perfect samples from conditional product measures by locally resampling only the variables involved in bad events.
• It leverages principles from the Lovász Local Lemma and the Moser–Tardos paradigm to ensure correctness and efficiency even in high-dimensional settings.
• PRS adapts to both discrete and continuous models, achieving polylogarithmic expected iterations under bounded dependency conditions.

Partial rejection sampling (PRS) is a family of algorithmic techniques for generating perfect samples from target distributions defined as product measures conditioned on the non-occurrence of combinatorially or spatially structured “bad events.” PRS generalizes classical rejection sampling by performing local, event-driven resampling of only those variables implicated in witnessed constraint violations, thereby achieving substantially higher algorithmic efficiency, especially in regimes where naive acceptance probabilities are exponentially small. PRS methods are rooted in combinatorial probability, with foundational links to the algorithmic Lovász Local Lemma and the Moser–Tardos resampling paradigm, but introduce additional requirements ensuring samples are drawn from the correct conditional (rather than merely any feasible) law.

1. Formal Problem Statement and Algorithmic Structure

Let $X = (X_1, \dots, X_n)$ be random variables, each $X_i$ taking values in a finite or measurable set $D_i$, independent under product measure $D = D_1 \times \cdots \times D_n$. The target distribution is the product law $D$ conditioned on a conjunction of constraints (encoded as a Boolean formula $\Phi(X) = \bigwedge_{k=1}^m \phi_k(X_{(\phi_k)})$), where each clause $\phi_k$ (the “bad event”) depends only on a specified subset of variables (the scope). The goal is to efficiently generate $X \sim D$ conditioned on $\Phi(X) = \text{true}$, i.e., to produce a sample from $\Pr_D(\cdot | \Phi(X))$ (Jerrum, 2021).

The core PRS algorithm proceeds as follows for discrete variables:

1. Initialization: Independently sample $X_i$0 for $X_i$1.
2. Event detection: Identify the current set $X_i$2 of violated constraints.
3. Resampling: While $X_i$3, choose any violated clause $X_i$4 and resample all variables in its scope.
4. Termination: Stop when all constraints are satisfied; return $X_i$5.

For product measures on continuous spaces or configurations (e.g., spatial point processes), resampling regions are defined by the union of the “scopes” of the bad events present, effectively localizing corrective updates (Guo et al., 2018, Moka et al., 2019).

2. Theoretical Guarantees and Correctness

PRS is rigorously correct for extremal constraint systems—where any two clauses sharing a variable cannot be simultaneously violated. It generates samples exactly according to the desired conditional law, regardless of the order in which violated clauses are selected for resampling (Jerrum, 2021). The confluence property of such extremal systems ensures that the final distribution is invariant to arbitrary resampling order.

For more general “quasi-extremal” systems (where interactions between bad events are governed by a bounded degree dependency graph), the expected number of resampling steps can be bounded using witness-tree analysis analogous to Moser–Tardos arguments. If the relative likelihood of each bad event $X_i$6 satisfies the LLL-like inequality $X_i$7 for some $X_i$8, then the expected number of resamplings of each clause $X_i$9 is at most $D_i$0 (Jerrum, 2021).

3. Methodological Variants and Continuous Extensions

Partial rejection sampling has been extended to infinite-dimensional and continuous settings, including spatial point processes and high-dimensional geometric models. In such contexts, PRS typically operates on partitioned state spaces (e.g., spatial cells in Gibbs processes), and bad events correspond to localized violations such as pairwise overlaps or forbidden configurations (Guo et al., 2018, Moka et al., 2019).

A general structure emerges: sample an initial configuration from a base (e.g., Poisson) process, identify localized constraint violations, define a minimal resampling region (union of the neighborhoods/scopes of all violated constraints), and resample the variables/configuration within this region only. The process iterates until all bad events are eliminated, maintaining correctness through a loop-invariant: the current state is always distributed as the product measure conditioned on the fixed set of failed constraints (Guo et al., 2018).

Correctness in continuous settings relies on monotonicity properties of the target density (e.g., if $D_i$1 and $D_i$2, then $D_i$3), as well as a lattice condition for incremental update probabilities. These ensure the coupling and conditional law preservation essential for perfect sampling (Guo et al., 2018).

4. Computational Complexity and Empirical Performance

In favorable parameter regimes—specifically, when the marginal probability that any bad event occurs is small and the corresponding dependency graph has bounded degree—PRS exhibits polylogarithmic (often $D_i$4 or $D_i$5) expected number of iterations, compared to the exponential scaling of naive rejection algorithms (Jerrum, 2021, Guo et al., 2018, Wellens, 2018, Moka et al., 2019).

For hard-disk and hard-sphere models, this translates to perfect samples in $D_i$6 rounds and $D_i$7 total time, provided the packing intensity is below an explicit threshold (e.g., $D_i$8 for the 2D hard-disk model) (Wellens, 2018, Guo et al., 2018). In finite-range Gibbs point processes, the expected running time is $D_i$9 iterations, where $D = D_1 \times \cdots \times D_n$0 is the interaction range, and each iteration can be parallelized due to localized dependency structure (Moka et al., 2019).

Empirical simulation studies support these theoretical findings, highlighting the competitiveness of PRS in domains like uniform proper graph coloring, statistical physics models (hard disks/spheres), and random constraint satisfaction, especially when implemented with parallelism (Moka et al., 5 Apr 2026, Moka et al., 2019).

5. Connections, Extensions, and Algorithmic Comparisons

PRS generalizes classical rejection sampling (which requires resampling all variables on any violation) and improves on its efficiency whenever high-dimensional or rare-event conditioning causes acceptance probabilities to become negligible. The main innovation is strictly local resampling, controlled by the dependency graph of bad events.

Comparisons with the Moser–Tardos algorithm highlight that PRS preserves the correct output distribution (i.e., draws a perfect sample from the target conditioned law), whereas Moser–Tardos guarantees only feasibility (returns some satisfying assignment but not with the right probability law) (Jerrum, 2021).

Extensions to sequential Monte Carlo and approximate Bayesian computation take the form of partial rejection control (PRC), which is not a direct generalization of PRS but rather applies its acceptance/rejection ideas within SMC frameworks to achieve variance reduction and increased effective sample size. Here, rejections are controlled adaptively to maintain a minimal particle quality, with provable variance reduction and unbiasedness properties (Peters et al., 2008, Kudlicka et al., 2019).

In Markov chain Monte Carlo, the philosophy has also inspired designs such as rejection-free jump chains and partial neighbor search, which subsample local neighborhoods for enhanced parallel speedup while controlling bias and maintaining detailed balance (Chen et al., 2022).

6. Practical Implementation and Algorithmic Applications

The modular structure of PRS allows hybrid and parallelized variants. For example, in uniform sampling of proper graph colorings, a soft relaxation approach creates a hierarchy of constraints, each enforced using PRS locally within subcritical percolation clusters of $D = D_1 \times \cdots \times D_n$1 size, yielding near-linear expected run time with inherent parallelism (Moka et al., 5 Apr 2026).

In continuum settings, domain decomposition enables scalable, distributive implementation: each spatial cell or resampling region is solved independently, often via a smaller-scale perfect sampling method (e.g., local dominating CFTP). Localized dependency management ensures that each region can be sampled in parallel without global synchronization, leveraging the bounded dependency structure induced by short-range interactions in the target law (Moka et al., 2019, Guo et al., 2018).

Applications span hard constraint models in statistical mechanics (hard disks, hard spheres), Gibbs random fields with finite-range interactions, random graph models, and combinatorial configurations such as sink-free orientations or proper colorings (Jerrum, 2021, Moka et al., 5 Apr 2026, Guo et al., 2018).

7. Limitations and Open Directions

The main bottleneck in PRS-based algorithms is the intensity of dependencies and the probability of bad events. Above certain model- or dimension-specific thresholds (e.g., packing density for geometric models), the expected number of rounds can become polynomial or worse, as resampling regions may grow to include a non-negligible fraction of the system (Guo et al., 2018, Wellens, 2018, Moka et al., 2019).

Other limitations include:

• Extension to non-monotone constraint families and infinite-dimensional non-lattice models may require substantial additional theory.
• In high-density regimes, overlap of constraint scopes undermines the efficacy of local resampling.
• Generalizations to heterogeneous or continuous constraints in high-dimensional or non-Euclidean spaces remain technically challenging.

Despite these challenges, PRS represents a unifying and widely applicable paradigm for exact sampling in structured probabilistic models. Continued development in coupling theory, dependency management, and parallel algorithm design are expected to further expand its range of effective applications (Moka et al., 5 Apr 2026, Moka et al., 2019, Guo et al., 2018, Jerrum, 2021).