Technical Summary: A Coupling-Based Approach to $f$-Divergences Diagnostics for Markov Chain Monte Carlo

Problem Statement

Markov Chain Monte Carlo (MCMC) methods are standard for approximating expectations under a target distribution $\pi$. While theoretical convergence is often analyzed using statistical divergences (e.g., Total Variation, Wasserstein, $f$-divergences), practical diagnostics have historically relied on metrics like the Gelman–Rubin statistic, which assess agreement on specific functionals rather than the fundamental convergence of the chain to stationarity.

Existing coupling-based diagnostics, such as those by Biswas et al. (2019), provide computable upper bounds for Total Variation (TV) and Wasserstein distances. However, these methods have limitations:

1. Scope: They are generally restricted to TV and Wasserstein metrics, lacking a unified framework for general $f$-divergences (e.g., Kullback–Leibler, $\chi^2$).
2. Utility: They serve primarily as alert tools (warning of non-convergence) but do not provide a mechanism to correct the discrepancy between the current chain distribution $\mu_t$ and the target $\pi$. Specifically, knowing $\|\mu_t - \pi\|_{TV} < \epsilon$ does not allow one to translate knowledge of $\mu_t$ into knowledge of $\pi$.
3. Warm-up Requirements: Some coupling-based methods require a pre-determined lag or warm-up period to ensure the chains have mixed sufficiently before diagnostics become valid.

The authors identify a gap: there is no general, online, computable diagnostic for arbitrary $f$-divergences that simultaneously provides consistent estimates of the Radon–Nikodym derivative $d\pi/d\mu_t$, enabling importance-weighted inference.

Methodology: Weight Harmonization via Couplings

The paper proposes a novel algorithm, Weight Harmonization, which runs parallel interacting Markov chains to produce consistent importance weights and convergence diagnostics.

Core Mechanism

The method utilizes $2N$ particles (chains) evolving under a $\pi$-invariant kernel $K$. Unlike independent parallel chains, these particles interact through couplings and weight averaging:

1. Coupling: At each step $t$, pairs of particles $(X^n_t, X^{m}_t)$ are evolved using a coupled kernel $\bar{K}$. This coupling is designed such that the marginal dynamics of each chain remain $K$, but the joint probability of the chains meeting ($X^n_{t+1} = X^m_{t+1}$) is positive.
2. Weight Averaging (Harmonization): When two coupled particles meet (i.e., their states coincide), their accumulated importance weights are averaged. If $X^n_{t+1} = X^{m}_{t+1}$, the new weights are set to $W^n_{t+1} = W^m_{t+1} = (W^n_t + W^m_t)/2$. If they do not meet, weights remain unchanged.
3. Reshuffling (Exchange): To prevent the system from stagnating in isolated pairs, the algorithm periodically reshuffles the pairing of particles. Pairs that have already coupled are randomized to exchange information with other chains, ensuring global mixing of weight information.

Theoretical Foundation

The method relies on the concept that the weighted empirical measure $\sum W^m_t \delta_{X^m_t}$ converges to the target $\pi$.

• Consistency: Theorem 2 establishes that for any bounded function $\varphi$, the weighted estimator $\sum W^m_t \varphi(X^m_t)$ converges in probability to $\int \varphi d\pi$ as the number of particles $N \to \infty$.
• Divergence Bounds: By applying Lemma 1 to the discrete weighted measure, the authors derive that the quantity $\frac{1}{2N} \sum f(2N W^m_t)$ serves as a probabilistic upper bound for the $f$-divergence $D_f(\pi \| \mu_t)$.
• Convergence: Theorem 3 proves that under ergodicity and coupling assumptions, the weights converge to uniformity ($1/2N$), causing the divergence bounds to converge to zero almost surely. Theorem 4 further shows that under a uniform coupling assumption, this convergence occurs exponentially fast.

Key Contributions

1. General $f$-Divergence Diagnostics: The paper introduces the first general convergence diagnostic for MCMC based on any $f$-divergence (including KL, $\chi^2$, Hellinger, and TV), rather than being limited to TV or Wasserstein distances.
2. Online Validity: Unlike lag-based methods (e.g., Biswas et al., 2019), the proposed method is valid from the first iteration ($t=0$) and requires no pre-defined warm-up period or lag parameter.
3. Radon–Nikodym Estimation: Beyond diagnostics, the method produces consistent estimates of the Radon–Nikodym derivative $d\pi/d\mu_t$. This allows the weights to be used for importance-weighted inference, correcting expectations under the current chain distribution $\mu_t$ to estimate expectations under $\pi$.
4. Theoretical Guarantees: The authors prove that the proposed divergence bounds are non-increasing over time and converge to zero as the chains approach stationarity.

Empirical Results

The paper evaluates the method on three benchmarks:

1. Tractable Gaussian System: A 100-dimensional Gaussian target where theoretical $f$-divergences are computable. The method consistently provides conservative upper bounds for the $\chi^2$-based effective sample size (ESS). The bounds are shown to be robust across different mixing rates but exhibit a "worst-case" degradation as the number of particles increases, attributed to the pairwise nature of the coupling.
2. Logistic Regression (Pólya–Gamma Gibbs): Applied to the German credit dataset. The method produces more conservative Total Variation bounds than the state-of-the-art Biswas et al. (2019) approach but does so without requiring the 350-iteration warm-up period necessary for the latter.
3. Stochastic Volatility (MALA): A high-dimensional ($D=2,500$) stochastic volatility model. The method again yields conservative TV bounds compared to Biswas et al. (2019), particularly during the "warm" phase where the lag-based method's bound has already begun to decrease. The authors note that the lag-based method required a 500-step warm-up (which was sufficient), whereas the harmonization method required no such tuning.

In all cases, the method is described as conservative: it tends to overestimate the divergence (or underestimate the ESS) compared to theoretical optima or other estimators. However, the authors argue this conservatism is a feature that ensures the bounds are reliable upper limits.

Significance and Claims

The paper claims to bridge the gap between theoretical convergence analysis and practical diagnostics by providing a unified, computable framework for $f$-divergences.

• Diagnostic Utility: It offers a practical tool to monitor convergence for a broader class of metrics than currently available.
• Inference Utility: The primary novelty is the dual use of the method: it is not just a diagnostic but a mechanism to generate consistent importance weights for inference, a capability lacking in previous coupling-based diagnostics.
• Limitations: The authors acknowledge that the method is more conservative than existing TV estimators and that the efficiency of the bounds degrades with the number of particles in the current pairwise implementation. They identify Rao–Blackwellization (integrating over all possible permutations rather than sampling one) and backward-correcting (offline smoothing) as potential avenues for future improvement to reduce this conservatism.

The work positions itself as a foundational step toward "importance-weighted MCMC," where the coupling mechanism is used not just to detect convergence but to actively correct the distribution of the chain samples.