Conclique-based Gibbs Sampling

• Conclique-based Gibbs Sampling is a method that partitions MRF nodes into non-neighboring groups (concliques) to allow efficient, parallel updates.
• By using graph-coloring techniques, the approach reduces the computational cost per iteration, achieving speedups of up to two orders of magnitude over single-site Gibbs sampling.
• Its effectiveness is model-dependent, showing limited gains in complete graphs or highly dependent systems, making it ideal for large-scale spatial and network models with bounded local dependencies.

Conclique-based Gibbs sampling is a methodological innovation in Markov random field (MRF) simulation designed to address the computational bottlenecks of standard single-site Gibbs sampling. By exploiting the independence structure in MRFs, the conclique-based approach simultaneously updates disjoint sets of non-neighboring sites—termed “concliques”—resulting in considerable computational acceleration, particularly for large or high-dimensional spatial and network models with bounded local dependencies (Kaplan et al., 2018).

1. Definition and Identification of Concliques

Let $G = (V, E)$ denote a finite undirected graph, where $V = \{1, \dots, n\}$ indexes sites or locations and $\mathcal{N}(i) = \{j \neq i : (i, j) \in E\}$ denotes the neighborhood of site $i$. A conclique $C \subset V$ is defined as a subset for which $i, j \in C$, $i \neq j$ implies $j \notin \mathcal{N}(i)$. No two sites in a conclique are neighbors, i.e., $\mathcal{N}(i) \cap C = \varnothing$ for each $i \in C$. This property ensures conditional independence within the conclique given the configuration outside.

To partition $V$ into concliques, one seeks a conclique cover $V = \bigcup_{k=1}^Q C_k$ with disjoint concliques $C_1, \dots, C_Q$. In graph-theoretic terms, this is equivalent to a proper vertex-coloring: the minimal number $Q$ is the chromatic number $\chi(G)$. Conclique covers can be found by:

• Greedy (Welch–Powell) algorithm: orders vertices by degree and assigns the smallest available color.
• DSatur algorithm: Brélaz's saturation-degree heuristic, often finding $\chi(G)$ in practice.
• Manual construction: for regular lattices, concliques can be described explicitly (e.g., $Q=2$ for two- or four-nearest neighbor, $Q=4$ for eight-nearest).

A brief pseudocode for the greedy algorithm is:

Input: Graph G=(V,E) Order vertices i₁,…,iₙ by decreasing degree. Initialize Q=0, color[i]=0 for all i. for t=1…n: let i = iₜ. let S = {1,2,…,Q} \ {color[j]: (i,j)∈E}. if S nonempty: assign color[i] = min(S) else: Q = Q + 1; color[i] = Q Output: concliques Cₖ = {i: color[i]=k}, k=1…Q.

2. Markov Random Field Model Structure

The MRF is defined on random variables $$X = (X_1, \dots, X_n) \in \mathcal{X} \subset \mathbb{R}^n$$, where the joint is implicitly specified by local conditionals:

$$f_i(x_i \mid x_{-i}) = f_i(x_i \mid x_{\mathcal{N}(i)}), \quad i = 1, \dots, n.$$

Assuming compatibility, there exists a joint density

$$\pi(x) = \frac{1}{Z} \prod_{i=1}^n f_i(x_i \mid x_{\mathcal{N}(i)})$$

with normalizing constant $Z$. Many relevant models use a one-parameter exponential family:

$$f_i(x_i \mid x_{\mathcal{N}(i)}) = \exp\left(A_i(x_{\mathcal{N}(i)}) T(x_i) - B_i(x_{\mathcal{N}(i)}) + C(x_i)\right).$$

This MRF setup underlies both spatial and network models, with structure enforced by $G$.

3. Conclique-based Gibbs Sampler: Algorithmic Description

With a conclique partition $C_1, \dots, C_Q$, the conclique-based Gibbs sampler (CGS) executes sequential “block” updates over concliques; within each, all $x_i$, $i \in C_k$, are updated independently and in parallel given current values outside $C_k$.

Precisely, for $t$ denoting the iteration and $\ell=1,\dots,Q$ indexing concliques,

• For each $i \in C_\ell$,

$$X_i^{(t,\ell)} \sim f_i\bigl(\cdot \mid X_{\mathcal{N}(i)}^{(t,\ell-1)}\bigr),$$

where

$$X_{\mathcal{N}(i)}^{(t,\ell-1)} = \{X_j^{(t+1)} : j \in \bigcup_{k < \ell} C_k\} \cup \{X_j^{(t)} : j \in \bigcup_{k > \ell} C_k\}.$$

• For $j \notin C_\ell$, set $X_j^{(t,\ell)} = X_j^{(t,\ell-1)}$.

After $Q$ conclique updates, let $X^{(t+1)} = X^{(t,Q)}$.

This mechanism ensures that for each conclique $C_k$:

$$\{X_i^{(t+1)}: i \in C_k\} \overset{\textrm{ind}}{\sim} \prod_{i \in C_k} f_i(x_i \mid x_{\mathcal{N}(i)})$$

holding the remaining variables fixed at their current or updated values.

4. Invariance and Ergodicity Properties

Each conclique update is a valid Gibbs step on the subvector $X_{C_k}$ conditional on the complement, preserving the target joint $\pi(\cdot)$. The entire composite step—sequentially applying all $Q$ conclique updates—also leaves $\pi$ invariant.

Ergodicity holds under the usual positivity condition: if every conditional $f_i$ has full support, the CGS chain $\{X^{(t)}\}$ is Harris ergodic with unique invariant $\pi$. When the number of concliques is $Q=2$, the chain is a two-component Gibbs sampler, which can be shown to be geometrically ergodic via standard drift/minorization theory.

5. Computational Complexity and Comparative Performance

Let $n = |V|$ and $Q$ denote the number of concliques. Single-site Gibbs sampling performs $n$ univariate updates per iteration ($O(n)$ per iteration). CGS decomposes each iteration into $Q$ block-updates, each processing $|C_k| \approx n/Q$ sites; the total per-iteration work remains $O(n)$, but the constant is reduced by a factor of roughly $Q$.

For models with bounded local neighborhoods (i.e., maximum neighborhood size does not increase with $n$), a classical result gives $Q \leq 1 + \max_i |\mathcal{N}(i)|$. For spatial lattices with neighborhood radius $R$ in $\mathbb{Z}^2$, $Q = O(R^2) \ll n$, leading to a practical speedup of $n/Q$ per full sweep. Empirical studies demonstrate speed gains of one to two orders of magnitude or more, depending on the model structure and $Q$ (Kaplan et al., 2018).

Model	$Q$ (No. concliques)	Single-site Gibbs Runtime	CGS Runtime	Reported Speedup
Gaussian MRF ($m \times m$ grid, $m=75$)	2	$\approx$ 3 hr	$\approx$ 15 sec	$\sim 700\times$
Binary autologistic ($40\times40$ grid)	2	$0.023$ s	$\sim 3\!\times\!10^{-4}$ s	$\sim 75\times$
Exponential-graph triad ($V=100$)	$Q=2\lceil V/2\rceil-1$	$\sim$ 9.6 hr (1000 sweeps)	$\sim$ 13 min	$\sim 45\times$

Empirical mixing (e.g., inverse IACT) remains essentially identical for both CGS and single-site Gibbs, confirming that gains are in computational cost per sweep, not in mixing rate per sweep.

6. Situations with Limited Efficacy

The speedup advantage of CGS is model-dependent. In the worst case $Q = n$, as for complete or nearly-complete graphs, so each conclique is a singleton and CGS collapses to single-site Gibbs with no computational gain. In models exhibiting long-range or global dependence (where $\max_i |\mathcal{N}(i)|$ grows with $n$), $Q$ may also grow rapidly, reducing $n/Q$. Additionally, CGS does not address intrinsic slow mixing of Gibbs samplers in near-critical or highly dependent regimes; it only lowers per-iteration cost, not the number of sweeps to reach stationarity.

7. Practical Implications and Applications

Conclique-based Gibbs sampling is a "plug-in" replacement for single-site updating in any conditionally specified MRF for which a conclique cover is available. The only additional computational requirement is a one-time graph-coloring (conclique-finding) step. For most spatial and network MRFs with bounded or slowly growing neighborhood size, $Q$ is small relative to $n$, enabling large speedups—often of two orders of magnitude or more—for applications requiring repeated simulation (e.g., bootstrap, null-distribution estimation). Available implementations (such as in the R package conclique) make the method accessible for practical use in high-performance MRF simulation environments (Kaplan et al., 2018).