FORS: First-Order Rejection Sampling

• FORS is a variant of classical rejection sampling that employs first-order information to adapt proposals and tighten empirical envelope bounds.
• It uses gradient-refined techniques by optimizing a smooth surrogate loss to directly minimize the envelope ratio, resulting in improved acceptance rates.
• FORS also incorporates structural conditional methods that reduce rejection costs and offer integration with advanced MCMC and divide-and-conquer approaches.

Searching arXiv for the cited papers and the phrase “First-Order Rejection Sampling” to ground the article in the relevant literature. First-Order Rejection Sampling (FORS) is not a standardized label in the cited papers. In the interpretation suggested by gradient-refined rejection sampling and by probabilistic divide-and-conquer, it denotes rejection-sampling procedures that preserve the classical accept/reject mechanism while using first-order information to improve proposal construction or to relocate rejection to a more favorable stage of sampling (Raff et al., 2023, Arratia et al., 2011). In one formulation, “first order” refers to gradients with respect to proposal parameters $\theta$, with the proposal learned by minimizing a smooth surrogate of the envelope ratio $f(x)/g(x;\theta)$. In another, more structural formulation, it refers to a first-stage conditional correction based on $p_a=\mathbb E[h(a,B)]$, after which the remaining coordinates are completed exactly. This suggests that FORS is best understood as an umbrella concept spanning both parameter-gradient and structural conditional refinements of classical rejection sampling.

1. Classical rejection sampling and the FORS viewpoint

Classical rejection sampling assumes a target density $f(x)$ on a domain $\mathcal X$, a proposal density $g(x;\theta)$ that is easy to sample from, and a finite constant $C$ such that

$$f(x)\le C\,g(x;\theta)\qquad \forall x\in\mathcal X.$$

The algorithm draws $x\sim g(\cdot;\theta)$ and $u\sim \mathcal U(0,1)$, then accepts $f(x)/g(x;\theta)$0 if

$f(x)/g(x;\theta)$1

Its acceptance rate is $f(x)/g(x;\theta)$2, and accepted samples are i.i.d. from $f(x)/g(x;\theta)$3.

The difficulty is not the accept/reject rule itself but the construction of a good $f(x)/g(x;\theta)$4 and a near-minimal valid $f(x)/g(x;\theta)$5. The gradient-refined approach of “An Easy Rejection Sampling Baseline via Gradient Refined Proposals” is explicitly motivated by the observation that rejection sampling is a common tool for low dimensional problems $f(x)/g(x;\theta)$6, yet is often non-trivial to apply because devising $f(x)/g(x;\theta)$7 and selecting $f(x)/g(x;\theta)$8 requires considerable mathematical effort; more advanced samplers may require additional mathematical derivations, limitations on $f(x)/g(x;\theta)$9, or even cross-validation (Raff et al., 2023).

Within a FORS interpretation, the essential point is that the classical rejection-sampling law is retained, but the envelope is no longer treated as a purely analytic object fixed in advance. Instead, either the proposal is adapted by first-order optimization, or the rejection event is reorganized through conditional factorization so that the effective acceptance probability becomes substantially larger.

2. Gradient-refined proposals as parameter-first-order rejection sampling

The most explicit FORS-like construction in the supplied literature is ERS, which keeps the rejection-sampling structure but makes both $p_a=\mathbb E[h(a,B)]$0 and $p_a=\mathbb E[h(a,B)]$1 learned objects. ERS chooses a parametric proposal family

$p_a=\mathbb E[h(a,B)]$2

with $p_a=\mathbb E[h(a,B)]$3, where $p_a=\mathbb E[h(a,B)]$4 denotes a Gaussian truncated to $p_a=\mathbb E[h(a,B)]$5. The key design choice is that the loss is derived directly from the acceptance threshold rather than from an $p_a=\mathbb E[h(a,B)]$6-, KL-, or MSE-style approximation criterion.

For a set of points $p_a=\mathbb E[h(a,B)]$7 from the current accept/reject history, ERS defines

$p_a=\mathbb E[h(a,B)]$8

and the sample-induced envelope constant

$p_a=\mathbb E[h(a,B)]$9

Direct minimization of $f(x)$0 is undesirable because the gradient passes only through the argmax. ERS therefore uses the smooth surrogate

$f(x)$1

which spreads gradient mass across multiple high-ratio regions. Its gradient takes the form

$f(x)$2

since $f(x)$3 is independent of $f(x)$4.

A common misconception is that “first order” here means the method explicitly uses $f(x)$5. It does not. ERS requires only the ability to evaluate a differentiable target density $f(x)$6 or $f(x)$7, but it does not explicitly use $f(x)$8. The first-order information is with respect to proposal parameters, obtained by automatic differentiation of $f(x)$9 and of the cached log-ratio $\mathcal X$0. The implementation described in the paper uses JAX and the first-order optimizer AdaBelief, and it does not differentiate through the stochastic accept/reject step itself; it optimizes the underlying envelope bound (Raff et al., 2023).

3. Empirical supremum, refinement loops, and asymptotic correctness

ERS alternates between proposal sampling, empirical envelope tracking, mixture re-fitting, and gradient refinement. The high-level loop samples from the current proposal $\mathcal X$1, updates an empirical bound $\mathcal X$2 on $\mathcal X$3, periodically fits a new GMM to accumulated data, and refines $\mathcal X$4 by gradient descent on the log-ratio loss. After training, rejection sampling is run with the final $\mathcal X$5 and $\mathcal X$6.

Its correctness mechanism is adapted from the empirical supremum method of Caffo et al. For a fixed proposal $\mathcal X$7, the estimator is updated by

$\mathcal X$8

and the acceptance rule becomes

$\mathcal X$9

ERS uses this idea in batches. For a batch $g(x;\theta)$0,

$g(x;\theta)$1

The concrete algorithm initializes accepted and rejected sets $g(x;\theta)$2 and $g(x;\theta)$3, sets $g(x;\theta)$4, maintains a lower estimate $g(x;\theta)$5, and draws candidate batches of size $g(x;\theta)$6 with $g(x;\theta)$7. If $g(x;\theta)$8 or $g(x;\theta)$9, the refinement flag is triggered. The algorithm also updates

$C$0

so that degradation in batch supremum can trigger later refinements. Accepted and rejected points are then reused in two ways: an occasional EM re-fit of the GMM, and the refinement step

$C$1

optimized for up to 800 steps. After refinement, the empirical envelope

$C$2

is recomputed; the new parameters are accepted only if $C$3, otherwise the method reverts to the previous $C$4.

The theoretical guarantee is asymptotic rather than exact. For fixed $C$5, the paper states that Algorithm 1 converges to the same or better solution, in terms of fewer false samples, as Caffo et al. (2002), and thus retains the $C$6 convergence rate of correctness. When $C$7 changes after refinement, the argument is conceptually restarted using a new initial $C$8 computed from existing samples under the new proposal. ERS is therefore “correct with high probability,” not strictly exact, because $C$9 is optimized on a finite sample set and $$f(x)\le C\,g(x;\theta)\qquad \forall x\in\mathcal X.$$0 is only empirically approximated (Raff et al., 2023).

4. Proposal geometry and relation to adjacent samplers

ERS uses a Gaussian mixture model with diagonal covariances, nonnegative mixture weights summing to one, and an increasing number of components $$f(x)\le C\,g(x;\theta)\qquad \forall x\in\mathcal X.$$1 as more data arrive. In the re-fit stage, the number of components is

$$f(x)\le C\,g(x;\theta)\qquad \forall x\in\mathcal X.$$2

the covariances are diagonal, and both accepted and rejected samples are used with weights proportional to $$f(x)\le C\,g(x;\theta)\qquad \forall x\in\mathcal X.$$3, with accepted samples inflated further by a constant factor. The proposal family is therefore adaptive, but still deliberately simple.

The diagonal covariance choice is computational rather than theoretical. It enables efficient truncated sampling $$f(x)\le C\,g(x;\theta)\qquad \forall x\in\mathcal X.$$4; full covariance would require heavier truncated-Gaussian machinery and was reported as empirically $$f(x)\le C\,g(x;\theta)\qquad \forall x\in\mathcal X.$$5 slower in edge cases. The method does not enforce a global analytic inequality $$f(x)\le C\,g(x;\theta)\qquad \forall x\in\mathcal X.$$6 for all $$f(x)\le C\,g(x;\theta)\qquad \forall x\in\mathcal X.$$7; instead it uses empirical maxima and minimizes the empirical maximum of the log-ratio. Violations outside the observed region remain possible.

Relative to adjacent rejection samplers, ERS changes what is optimized. OS* and A* require a user-supplied decomposition

$$f(x)\le C\,g(x;\theta)\qquad \forall x\in\mathcal X.$$8

together with problem-specific bounds and proposal design. PRS and NNARS use kernel density estimators, optimize an $$f(x)\le C\,g(x;\theta)\qquad \forall x\in\mathcal X.$$9-type discrepancy

$x\sim g(\cdot;\theta)$0

and require bounded-support assumptions, hyper-parameters such as bandwidths or Hölder constants, and sometimes cross-validation or transformations to bounded domains. ERS instead optimizes the acceptance ratio $x\sim g(\cdot;\theta)$1 directly, requires no log-concavity assumptions, no analytic envelope derivation, no special factorization of $x\sim g(\cdot;\theta)$2, and no cross-validation. Relative to HMC and related gradient-based MCMC, the distinction is equally sharp: HMC uses $x\sim g(\cdot;\theta)$3 to evolve a Markov chain, whereas ERS is a pseudo-i.i.d. rejection sampler using gradients with respect to proposal parameters rather than the state (Raff et al., 2023).

5. Empirical behavior in low-dimensional regimes

The empirical scope of ERS is the regime where rejection sampling is still viable. The paper emphasizes the classical low-dimensional setting $x\sim g(\cdot;\theta)$4, while also reporting some experiments up to $x\sim g(\cdot;\theta)$5 for a synthetic multi-modal benchmark. Three benchmark families structure the evaluation: a one-dimensional peakiness problem on $x\sim g(\cdot;\theta)$6, a bounded multi-modal scaling problem on $x\sim g(\cdot;\theta)$7, and Minka’s clutter problem in one and two dimensions (Raff et al., 2023).

In the peakiness problem,

$x\sim g(\cdot;\theta)$8

increasing $x\sim g(\cdot;\theta)$9 makes the target increasingly peaked. For $u\sim \mathcal U(0,1)$0, NNARS and PRS slightly beat ERS, but all methods have high acceptance. As $u\sim \mathcal U(0,1)$1 increases to $u\sim \mathcal U(0,1)$2, ERS’s acceptance rate drops only about 1 percentage point per step, whereas NNARS and PRS deteriorate dramatically, down to approximately $u\sim \mathcal U(0,1)$3 and $u\sim \mathcal U(0,1)$4. In the multi-modal scaling problem,

$u\sim \mathcal U(0,1)$5

ERS has significantly higher acceptance rates than NNARS, PRS, and other baselines for $u\sim \mathcal U(0,1)$6. As $u\sim \mathcal U(0,1)$7 grows to $u\sim \mathcal U(0,1)$8, however, it suffers from the curse of dimensionality and becomes statistically indistinguishable from NNARS and PRS in acceptance rate.

The strongest numerical gains are reported on the clutter problem. In one dimension, ERS achieves approximately $u\sim \mathcal U(0,1)$9 acceptance with $f(x)/g(x;\theta)$00, compared with $f(x)/g(x;\theta)$01 for PRS, $f(x)/g(x;\theta)$02 for A*, and $f(x)/g(x;\theta)$03 for a simple rejection-sampling baseline. In two dimensions, ERS reaches approximately $f(x)/g(x;\theta)$04 with $f(x)/g(x;\theta)$05, while PRS achieves $f(x)/g(x;\theta)$06, A* $f(x)/g(x;\theta)$07, and simple rejection sampling is essentially zero. The reported gain is up to approximately $f(x)/g(x;\theta)$08 in acceptance rate relative to some baselines.

Distributional checks are used to support the approximate-correctness claim. For $f(x)/g(x;\theta)$09, the paper reports two-sample Kolmogorov–Smirnov tests between ERS samples and either A* samples or direct draws from the known ground-truth density, with no significant difference detected. For $f(x)/g(x;\theta)$10, a two-sample Cramér-type test is used, again with no significant discrepancies. Runtime is also favorable on the clutter problem with 100k accepted samples: ERS on CPU requires approximately $f(x)/g(x;\theta)$11 s, ERS on GPU approximately $f(x)/g(x;\theta)$12 s, A* approximately $f(x)/g(x;\theta)$13 s, and OS* approximately $f(x)/g(x;\theta)$14 s. Taking acceptance-rate differences into account, ERS is reported as roughly $f(x)/g(x;\theta)$15 faster than A* on CPU and approximately $f(x)/g(x;\theta)$16 faster with GPU acceleration.

6. Structural FORS through probabilistic divide-and-conquer

A different route to FORS-like behavior appears in “Probabilistic divide-and-conquer: a new exact simulation method, with integer partitions as an example.” The setup begins with two measurable spaces $f(x)/g(x;\theta)$17 and $f(x)/g(x;\theta)$18, independent random variables

$f(x)/g(x;\theta)$19

and an indicator $f(x)/g(x;\theta)$20 with

$f(x)/g(x;\theta)$21

The target object is

$f(x)/g(x;\theta)$22

Classical rejection sampling is recovered by drawing $f(x)/g(x;\theta)$23 from the product law and accepting if $f(x)/g(x;\theta)$24, at expected cost $f(x)/g(x;\theta)$25.

PDC changes where rejection occurs. If $f(x)/g(x;\theta)$26 is sampled from $f(x)/g(x;\theta)$27 and, conditional on $f(x)/g(x;\theta)$28, $f(x)/g(x;\theta)$29 is sampled from $f(x)/g(x;\theta)$30, then $f(x)/g(x;\theta)$31. In practice this is implemented by rejection on $f(x)/g(x;\theta)$32 alone. Defining

$f(x)/g(x;\theta)$33

and $f(x)/g(x;\theta)$34, the optimal acceptance function is

$f(x)/g(x;\theta)$35

The corresponding acceptance cost is

$f(x)/g(x;\theta)$36

and the speedup over naïve rejection is

$f(x)/g(x;\theta)$37

This is exact simulation, not approximate correction.

The integer-partition example is the paper’s central case study. Under Fristedt’s representation, with independent geometric coordinates $f(x)/g(x;\theta)$38 and total weight $f(x)/g(x;\theta)$39, basic rejection sampling for a uniform random partition of $f(x)/g(x;\theta)$40 has success probability

$f(x)/g(x;\theta)$41

so the expected number of proposals is asymptotically $f(x)/g(x;\theta)$42. A deterministic second-half decomposition takes $f(x)/g(x;\theta)$43, $f(x)/g(x;\theta)$44, and completes $f(x)/g(x;\theta)$45 deterministically after acceptance. The resulting acceptance cost drops from order $f(x)/g(x;\theta)$46 to order $f(x)/g(x;\theta)$47. A more elaborate self-similar recursive factorization based on

$f(x)/g(x;\theta)$48

yields constant asymptotic rejection cost per recursion layer: first $f(x)/g(x;\theta)$49, then $f(x)/g(x;\theta)$50 after a parity refinement. Theorem 3.7 further states that if the proposal-generation cost is $f(x)/g(x;\theta)$51 with $f(x)/g(x;\theta)$52, then

$f(x)/g(x;\theta)$53

The same deterministic-second-half pattern appears for partitions with bounded largest part, random set partitions, and plane partitions. For random set partitions, the proposal count is reduced from $f(x)/g(x;\theta)$54 to $f(x)/g(x;\theta)$55. For plane partitions, PDC yields a speedup of order $f(x)/g(x;\theta)$56 in the rejection stage, reducing that stage to $f(x)/g(x;\theta)$57 and leaving the overall complexity $f(x)/g(x;\theta)$58, dominated by Pak’s bijection. The paper also introduces a mix-and-match variant for generating many i.i.d. outputs under a “simple matching” condition

$f(x)/g(x;\theta)$59

where color maps $f(x)/g(x;\theta)$60 and $f(x)/g(x;\theta)$61 enable a coupon-collector effect in the $f(x)/g(x;\theta)$62-phase and can yield sublinear dependence on the number of desired samples (Arratia et al., 2011).

A plausible implication is that PDC represents a structural form of FORS. It does not use derivatives of a density, but it performs a first-stage correction using the local conditional success probability $f(x)/g(x;\theta)$63, thereby reweighting only part of the proposal and postponing exact completion to a conditional second stage.

7. Limits, exactness, and prospective extensions

The two FORS interpretations differ most sharply in their correctness status. ERS is approximate: it does not guarantee a global analytic envelope $f(x)/g(x;\theta)$64 outside the observed sample set, and rare false accepts are theoretically possible. The paper attributes its practical reliability to empirical-supremum convergence and reports no evidence of false acceptances in the final empirical checks. PDC, by contrast, is exact whenever the conditional laws are sampled correctly. A common misconception is therefore to treat all rejection-sampling refinements as equally exact; the supplied literature does not support that conclusion (Raff et al., 2023, Arratia et al., 2011).

ERS also inherits the usual dimensionality limitations of rejection sampling. Its efficiency collapses as dimension increases; empirically, beyond approximately $f(x)/g(x;\theta)$65–$f(x)/g(x;\theta)$66, acceptance rates become similar to KDE-based methods and relatively low. The refinement step is non-convex, so some updates worsen $f(x)/g(x;\theta)$67 and must be discarded. The GMM with diagonal covariance is flexible for $f(x)/g(x;\theta)$68 and some moderate-dimensional problems, but may be inadequate for more complex high-dimensional targets. PDC has different limitations: it requires an explicit decomposition $f(x)/g(x;\theta)$69, tractable access to $f(x)/g(x;\theta)$70, and, for the strongest recursive results, a self-similar factorization such as the partition identity $f(x)/g(x;\theta)$71.

The practical guidance in the ERS paper is correspondingly narrow and concrete: use ERS as a default baseline when a differentiable $f(x)/g(x;\theta)$72 is available and $f(x)/g(x;\theta)$73 is small; there is no need to specify $f(x)/g(x;\theta)$74, $f(x)/g(x;\theta)$75, or domain transformations, though $f(x)/g(x;\theta)$76 may optionally be supplied when bounded. The paper also suggests using ERS as a front-end to more advanced MCMC, for example by generating a modest number of high-quality i.i.d. samples to initialize chains. For a broader FORS program, the same source proposes three extensions: generalizing the proposal family beyond GMMs, for example to normalizing flows; incorporating explicit first-order information in $f(x)/g(x;\theta)$77, such as $f(x)/g(x;\theta)$78, to shape the proposal more aggressively; and adding regularization to balance envelope tightness against proposal variance. These directions remain programmatic rather than settled.