Tight Bounds for Sampling $q$-Colorings via Coupling from the Past

The coupling from the past (CFTP) technique, introduced by Propp and Wilson [43], provides a general framework for the perfect simulation of the stationary distribution of a Markov chain. The key idea is to consider an idealized Markov chain that has been running from time $-\mathrm{\infty }$ up to time $0$, so that its final state $X_0\sim \pi$ is an exact (perfect) sample from the stationary distribution $\pi$. To generate such a sample, the CFTP algorithm searches for a (possibly random) time $-T$ in the past such that, when the chain is simulated forward from every possible initial state at time $-T$ using the same random bits, all trajectories

coalesce into the same state at time $0$. Once such a coalescence time is identified, the final state $X_0\sim \pi$ is completely determined by the random bits used in the simulation, independent of the initial state at time $-T$. Consequently, $X_0$ is a perfect sample from the stationary distribution $\pi$.

At the heart of CFTP lies the notion of a grand coupling, which provides a unified probabilistic construction that simultaneously drives the evolution of the chain from all possible initial states. Formally, for a Markov chain with transition kernel $P(x,y)=\mathrm{Pr}{X}_{t+1}=y\mid X_t=x$ on state space $\mathrm{\Omega }$, a grand coupling is a mapping $g:\mathrm{\Omega }\times [0,1]\to \mathrm{\Omega }$, such that for $U\sim Uniform[0,1]$ and every $x,y\in \mathrm{\Omega }$,

In other words, $g(\cdot ,U)$ simulates a single Markov chain transition under shared randomness $U$, thereby coupling all trajectories in a common probability space. This notion extends naturally to the $T$-step transitions, represented by a grand coupling $F:\mathrm{\Omega }\times [0,1]\to \mathrm{\Omega }$ for the evolution $X_{-T}\to X_0$. The CFTP algorithm identifies a time $-T$ in the past and terminates once $F(\cdot ,U)$ becomes constant, that is, when all trajectories coalesce into a single state at time $0$.

Designing grand couplings that achieve fast coalescence is therefore the central problem in CFTP. Beyond CFTP, grand couplings have found applications in a variety of contexts, including derandomization of MCMC [15], analytic stability (absence of complex zeros) [41], and local algorithms for sampling [39]. Despite being a natural and fundamental concept, their systematic design and analysis remain poorly understood. While fast-coalescing grand couplings immediately imply rapid mixing of the underlying Markov chain, essentially nothing is known about the converse in general. This raises a fundamental question:

Does fast coalescence of grand couplings require stronger conditions than rapid mixing?

A key algorithmic challenge in implementing CFTP lies in detecting coalescence, that is, determining whether the induced random mapping on the state space $\mathrm{\Omega }$ has become a constant function. Formally, let ${\mathrm{\Omega }}_t\subseteq \mathrm{\Omega }$ denote the set of all possible states of $X_t$ that can be reached from different initial configurations $X_{-T}\in \mathrm{\Omega }$ under the same randomness. Specifically,

where each $g_t=g_t(\cdot ,U)$ denotes the grand coupling for the transition $X_t\to X_{t+1}$ with $U\sim Uniform[0,1]$. Initially, ${\mathrm{\Omega }}_{-T}=\mathrm{\Omega }$, and clearly $\left|{\mathrm{\Omega }}_{t+1}\right|\le \left|{\mathrm{\Omega }}_t\right|$ since $g_t$ is a function. Coalescence occurs once $\left|{\mathrm{\Omega }}_t\right|=1$ for some $t\le 0$. In certain special cases, most notably in monotone systems such as the ferromagnetic Ising model, where a natural partial order on $\mathrm{\Omega }$ is preserved by the grand coupling, checking coalescence is straightforward. In contrast, for general systems, verifying coalescence is substantially more intricate and can be computationally intractable.

The bounding chains method, introduced by Huber [30, 31] and independently by Häggström and Nelander [24], provides an efficient way to detect coalescence by tracking a “bounding box” of the subset ${\mathrm{\Omega }}_t$ in the space $\mathrm{\Omega }\subseteq {[q]}^V$. The bounding chain maintains, for each vertex $v\in V$, a bounding list ${ℒ}_t(v)\subseteq [q]$ of colors, which contains all colors that vertex $v$ may take under any configuration in ${\mathrm{\Omega }}_t$. The corresponding bounding box is given by

Initially, we set ${ℒ}_{-T}(v)=[q]$ for every vertex $v$. The bounding list configuration ${ℒ}_t={({ℒ}_t(v))}_{v\in V}$ evolves as a Markov chain: at each time step $t\ge -T$, the next bounding list ${ℒ}_{t+1}$ is constructed from the current list ${ℒ}_t$, under the conservative assumption that every configuration $\omega \in {ℬ}_t$ may appear at time $t$. By induction, we have ${\mathrm{\Omega }}_t\subseteq {ℬ}_t$ for any $t\ge -T$. Thus, coalescence at time $t$ occurs if the bounding box has collapsed to the size $\left|{ℬ}_t\right|=1$, or equivalently, when $\left|{ℒ}_t(v)\right|=1$ for all $v\in V$.

The fast coalescence of bounding chains relies on the following standard contraction property.

Intuitively, this condition requires the total size of the bounding lists to be monotonically non-increasing in expectation at each step, even in the worst case. All existing CFTP algorithms based on bounding chains satisfy this condition, including [30, 4, 35]. The seminal work of Huber [30] proposed a bounding chain algorithm for uniform sampling of $q$-coloring that is fast coalescing under the condition $q\ge (1+o(1)){\mathrm{\Delta }}^2$. For many years, this remained the best achievable bound for any CFTP algorithm for $q$-coloring, until recent breakthroughs [4, 35] improved the condition for fast-coalescing bounding chains first down to $q>3\mathrm{\Delta }$ and later to $q\ge (\frac{8}{3}+o(1))\mathrm{\Delta }$.

We establish an asymptotically optimal threshold of $q=2.5\mathrm{\Delta }$ for CFTP algorithms based on bounding chains for proper $q$-colorings. In particular, we design an efficient CFTP algorithm that samples uniform proper $q$-colorings under the improved condition $q\ge (2.5+o(1))\mathrm{\Delta }$.

The algorithm in Theorem 1.2 is based on a new construction of bounding chains for the Glauber dynamics on uniform proper $q$-colorings, which coalesce rapidly using substantially fewer colors.

This result advances the state of the art for CFTP algorithms and, more broadly, for all perfect sampling algorithms for $q$-colorings that achieve exponential convergence.

We show that the upper bound in Theorem 1.2 is asymptotically tight by establishing a matching lower bound for bounding-chain-based grand couplings satisfying the contraction property in ˜1.1.

Together, our results establish a tight bound for CFTP for $q$-colorings based on bounding chains. They show that contractive bounding chains require a significantly stronger condition on the number of colors $q$ than is needed for rapid mixing of the underlying Markov chain. In particular, the critical threshold $q>2.5\mathrm{\Delta }$ identified here is substantially stronger than the classical Dobrushin condition for $q$-colorings, which holds when $q>2\mathrm{\Delta }$ and corresponds to the worst-case contraction condition in path coupling [36, 45, 5], thereby guaranteeing rapid mixing.

In the CFTP framework, the random evolution $X_{-T}\to X_0$ of the Markov chain $(X_t)$ is captured by a random mapping $F:\mathrm{\Omega }\to \mathrm{\Omega }$, generated by a grand coupling $F=F(\cdot ,U)$ using a random seed $U\in [0,1]$. Equivalently, $F$ can be expressed as a composition of single-step couplings, $F=g_{-1}\circ g_{-2}\circ \mathrm{\cdots }\circ g_{-T}$, where each $g_t=g_t(\cdot ,U_t)$ corresponds to the grand coupling of the one-step transition at time $t$, determined by an independent random seed $U_t\sim Uniform[0,1]$. Conceptually, the global seed $U$ can thus be viewed as an infinite sequence of independent random variables $(U_{-1},U_{-2},\mathrm{\dots })$, each determining one step of the coupled evolution.

When the Markov chain $(X_t)$ is a single-site dynamics, each transition $g_t$ is specified by a pair $(v_t,f_{v_t})$, where $v_t\in V$ denotes the vertex selected for update at time $t$, and $f_{v_t}$ is a local random update function that recolors $v_t$ based on the current coloring of its neighbors. Once the sequence of update vertices $v_{-1},v_{-2},\mathrm{\dots }$ is fixed, the CFTP process is fully determined by these local updates.

Conceptually, the bounding-chain construction can then be viewed as a localization of this global coupling: rather than tracking the evolution of the entire configuration space $\mathrm{\Omega }\subseteq {[q]}^V$, it maintains for each vertex $v$ a local bounding list $ℒ(v)\subseteq [q]$ of possible colors, which compactly encodes all configurations consistent with the current stage of the coupled evolution.

The critical condition for rapid mixing of local Markov chains on proper $q$-colorings, such as Glauber dynamics, has been a central question in the study of sampling algorithms. A classic result of Jerrum [36] (and independently, Salas and Sokal [45]) established that when $q>2\mathrm{\Delta }$, the Glauber dynamics mixes in $O(n\log n)$ time. This foundational work laid the groundwork for coupling-based analyses, notably inspiring the path coupling technique of Bubley and Dyer [5] and the coupling criterion of Hayes [25]. Vigoda [46] introduced a coupling for the flip dynamics and proved an $O(n^2)$ mixing-time bound for the Glauber dynamics when $q>\frac{11}{6}\mathrm{\Delta }$. This threshold remained the state of the art for nearly two decades until Chen et al. [7] further refined it to $q>\left(\frac{11}{6}-{\epsilon }_0\right)\mathrm{\Delta }$ for an explicit ${\epsilon }_0=1/84000$. Most recently, Carlson and Vigoda [6] improved the bound to $q>1.809\mathrm{\Delta }$, which remains the best-known rapid-mixing threshold for general graphs.

Under additional structural assumptions, such as a lower bound on the girth, sharper thresholds are known. For instance, on triangle-free graphs, the Glauber dynamics mixes optimally whenever $q>{\alpha }^{\ast }\mathrm{\Delta }$ for some constant ${\alpha }^{\ast }\approx 1.763$; see [27, 9, 16, 11, 34]. Even stronger results hold for graphs of larger girth [26, 14, 10]. In particular, the recent breakthrough of Chen et al. [10] establishes an almost tight condition $q\ge \mathrm{\Delta }+3$ for an $O(n\log n)$ mixing time of the Glauber dynamics on graphs with constant maximum degree $\mathrm{\Delta }=O(1)$ and sufficiently large girth.

The conjectured optimal threshold for rapid mixing is $q>\mathrm{\Delta }$, corresponding to the uniqueness threshold of the Gibbs measure on the infinite $\mathrm{\Delta }$-regular tree [38]. In contrast, approximate sampling from $q$-colorings becomes computationally intractable when , unless $\mathrm{NP}=\mathrm{RP}$ [20].

In a separate line of work focused on deterministic approximate counting, Gamarnik and Katz [21] gave a polynomial-time approximate counting algorithm for proper $q$-colorings when $q>2.84\mathrm{\Delta }$ for bounded-degree graphs. Subsequent improvements were obtained by Lu and Yin [42] ($q>2.58\mathrm{\Delta }$) and by Liu, Sinclair, and Srivastava [40] ($q>2\mathrm{\Delta }$). This $2\mathrm{\Delta }$ threshold has recently been improved by Bencs, Berrekkal, and Regts [3], who obtained a deterministic approximation algorithm for $q\ge (2-\eta )\mathrm{\Delta }$ for some $\eta \ge 0.002$ using zero-freeness and Barvinok interpolation. A different recent approach of Chen, Feng, Guo, Zhang, and Zou [8] derives deterministic counting algorithms from coupling independence, giving FPTASes for colorings in regimes matching the best known randomized approximate-counting bounds.

Let $G=(V,E)$ be an undirected graph and $[q]=\{1,2,\mathrm{\dots },q\}$ a set of $q\ge 2$ colors. The configuration space of proper $q$-colorings of $G$ is

For each vertex $v$, let $\mathrm{\Gamma }(v):=\{u\in V\mid \{u,v\}\in E\}$ denote $v$’s neighborhood.

We consider the (heat-bath) Glauber dynamics for proper $q$-colorings. Fix a predetermined, possibly random, update schedule $(v_t)$, chosen independently of the current coloring. At time $t$, the dynamics updates the color of $v_t$ and leaves all other vertices unchanged. For a vertex $v\in V$ and a neighborhood coloring $\omega \in {[q]}^{\mathrm{\Gamma }(v)}$, the local update distribution is

where $Uniform(S)$ denotes the uniform distribution over a set $S$. Here, by slight abuse of notation, we treat the tuple $\omega$ equivalently as a set, so that $[q]\setminus \omega =[q]\setminus \{\omega (u):u\in \mathrm{\Gamma }(v)\}$.

We define bounding chains for the above Glauber dynamics. Let $g_t:\mathrm{\Omega }\times [0,1]\to \mathrm{\Omega }$ denote a grand coupling for the Glauber update at time $t$, as defined in (1). The randomness in the update schedule is identically coupled across all trajectories: at each time $t$, every coupled trajectory updates the same vertex $v_t$ using the same random seed $U_t$. Since a Glauber update at $v_t$ depends only on the colors in $\mathrm{\Gamma }(v_t)$, the outcome at $v_t$ is determined by the neighborhood configuration and the shared seed.

A bounding chain $({ℒ}_t)$ maintains, for every vertex $v\in V$, a bounding list ${ℒ}_t(v)\subseteq [q]$ containing all colors that $v$ may take at time $t$ under the coupled evolution. Given the shared random seed $U_t$ and the update vertex $v_t$, the bounding chain evolves as

In words, when vertex $v_t$ is updated, the new bounding list ${ℒ}_{t+1}(v_t)$ consists of all possible colors that $v_t$ can take, given the current bounding lists of its neighbors and the shared randomness $U_t$. For all other vertices, the bounding lists remain unchanged. In practice, one may instead maintain a superset of the set defined in (4); this preserves correctness, though it may affect the efficiency of coalescence.

In coupling from the past (CFTP), the bounding chain is initialized with full uncertainty: ${ℒ}_{-T}(v)=[q]$ for all $v\in V$. As the chain $({ℒ}_t)$ evolves forward from time $-T$ to time $0$, the bounding lists ${ℒ}_t{(v)}_{v\in V}$ gradually shrink as uncertainty about each vertex’s color decreases. When all bounding lists collapse to singletons, i.e., $|{ℒ}_t(v)|=1$ for every $v\in V$, the coupled Glauber dynamics have coalesced to a single configuration. At this point, the bounding-chain-based CFTP procedure outputs a perfect sample from the stationary distribution of the Glauber dynamics.

We introduce a general formulation for grand couplings constructed using bounding chains, which applies to single-site dynamics. Let $G=(V,E)$ be a graph and $[q]=\{1,2,\mathrm{\dots },q\}$ a set of $q\ge 2$ colors. Recall that a single-site dynamics on a state space $\mathrm{\Omega }\subseteq {[q]}^V$ is specified via a collection of local update distributions: $\{p_v^{\omega }:v\in V,\omega \in {[q]}^{\mathrm{\Gamma }(v)}\}$.

In the bounding-chain framework, we maintain for each vertex $v\in V$ a bounding list $ℒ(v)\subseteq [q]$, representing the set of possible colors that $v$ may take across coupled trajectories. The grand coupling is then constructed based on these bounding lists.

Previous bounding-chain based algorithms (e.g., [4, 35]) implicitly employ two stages. These stages can be unified in our framework for grand couplings via bounding chains.

• $\mathrm{■}$

  Generating stage: Once a random seed $U\sim Uniform[0,1]$ is generated, the function $f_v$ determines, for every $p\in {𝒫}_v$, a possible update outcome of $v$. The collection of all outcomes $\{f_v(p,U):p\in {𝒫}_v\}$ forms the new bounding list ${ℒ}^{\prime }(v)$ of vertex $v$, representing all colors that could be assigned under the current randomness.

• $\mathrm{■}$

  Decoding stage: Once the actual neighbor configuration $\omega \in {[q]}^{\mathrm{\Gamma }(v)}$ (and thus the corresponding update distribution $p_v^{\omega }\in {𝒫}_v$) is revealed, the updated color of $v$ is determined by $f_v(p_v^{\omega },U)$.

This formulation unifies the generating–decoding paradigm into a single mathematical framework, enabling an analysis of the optimality of grand coupling designs.

In this subsection we prove Theorem 1.4 by establishing a fundamental lower bound for Problem 3.3. Our analysis reveals an intrinsic obstacle that limits the effectiveness of all known bounding-chain constructions [30, 4, 35]:

The difficulty of reducing every bounding list from size $2$ to $1$.

This step represents the critical bottleneck in coalescence of bounding chains. To understand this limitation, we focus on the configuration where each neighbor $u\in \mathrm{\Gamma }(v)$ has a bounding list $ℒ(u)$ of size $2$. This setting captures the worst-case local uncertainty at vertex $v$, and hence determines the strongest contraction that any grand coupling can possibly attain. Our lower-bound argument is therefore developed under this extremal configuration.

Theorem 3.4 is established by explicitly constructing a bounding-list configuration exhibiting a structural bottleneck, formalized as follows.

We now present a global version of the above lower bound, which extends the local obstruction in Theorem 3.4 to the entire graph, thereby proving Theorem 1.4.

The proofs of Theorems 3.4 and 3.6 are deferred to the full version.

As noted in Section 2.3, any bounding-chain–based grand coupling selects a vertex $v$ at random in each step and then applies a local grand coupling at $v$. Thus Theorem 3.6 is in fact a formal restatement of Theorem 1.6, and the latter can be derived as a straightforward corollary of the former.

Based on the bounding-chain framework, we formulated in ˜3.3 an optimization problem for designing grand couplings.

To obtain a more tractable and intuitive formulation, we approximate ${𝒫}_v$ by considering only the set of colors $𝒮={\bigcup }_{u\in \mathrm{\Gamma }(v)}ℒ(u)$ appearing in the neighborhood bounding list of $v$, which coarsely tracks the local uncertainty around $v$. This motivates the following simplified grand coupling design problem, which serves as an approximation of ˜3.3.

It is straightforward to verify that any grand coupling $f:𝒫\times [0,1]\to [q]$ that is a feasible solution to ˜4.1 on an instance $(𝒮,\mathrm{\Delta })$ is also a feasible solution to ˜3.3 for any vertex $v\in V$ with neighborhood bounding list configuration ${(ℒ(u))}_{u\in \mathrm{\Gamma }(v)}$ satisfy $|\mathrm{\Gamma }(v)|=\mathrm{\Delta }$ and $𝒮={\bigcup }_{u\in \mathrm{\Gamma }(v)}ℒ(u)$. Indeed, for any $\omega \in {⨂}_{u\in \mathrm{\Gamma }(v)}ℒ(u)$, we have

Hence ${𝒫}_v\subseteq 𝒫({𝒮}_{\mathrm{\Delta }})$, where ${𝒫}_v$ denotes the family of local update distributions at $v$ defined in (6). Consequently, any grand coupling $f$ constructed on $𝒫({𝒮}_{\mathrm{\Delta }})$ in ˜4.1 remains a valid grand coupling on ${𝒫}_v$ in ˜3.3.

The simplified formulation in ˜4.1 is not only conceptually cleaner, but also captures the key structure underlying previous designs. In particular, the core grand couplings employed in the algorithms of Bhandari et al. [4] (Algorithm 3, Contract) and Jain et al. [35] (Algorithm 3, Seeding) can both be viewed within the framework of ˜4.1. Building on this connection, we present a new construction that optimally solves ˜4.1.

In this subsection, we present an optimal construction of the grand coupling for ˜4.1, focusing on the nontrivial case where .

For a grand coupling $f:𝒫\times [0,1]\to [q]$ that is feasible for ˜4.1, define, for each $U\in [0,1]$, the bounding set of possible outcomes

The grand coupling $f$ naturally induces a probability distribution ${\{r_k\}}_{k\ge 1}$ over the sizes of ${ℒ}_U$, where

We next present a system of linear constraints on the probabilities $\{r_k\}$ which correspond to the grand couplings that are feasible for ˜4.1.

Since the objective of ˜4.1 is to minimize $𝔼[|{ℒ}_U|]={\sum }_{k}k\cdot r_k$ where $r_k=\mathrm{Pr}|{ℒ}_U|=k$, we can equivalently reformulate ˜4.1 as a linear program over $\{r_k\}$ subject to (7).

Theorem 4.4 follows directly from Lemma 4.2, since Lemma 4.2 provides both necessary and sufficient conditions for any sequence ${\{r_k\}}_{k=1}^{\mathrm{\Delta }}$ that corresponds to a feasible solution to ˜4.1.

The proof of Lemma 4.2 and the discussion on solving the linear program in (8) are deferred to the full version.

In this subsection, we present several concrete constructions of grand couplings that serve as building blocks for the CFTP algorithm stated in Theorem 1.2.