MacEachern-Berliner (1994) Subsampling the Gibbs Sampler

Summary

MacEachern and Berliner (1994) prove in three pages that systematically subsampling (thinning) the output of a stationary Markov chain never improves posterior mean estimates and generally degrades them. The argument requires only stationarity — not reversibility — and therefore covers the standard fixed-scan Gibbs sampler. The paper is the canonical reference for the practitioner rule "use the full chain, not a thinned subset."

Key Claims

• Thinning always loses. For any function with finite variance, the full-chain mean has weakly lower variance than any systematic subsample of the same total run length. Equality holds only for a degenerate (constant) distribution.
• Extension beyond Geyer (1992). Geyer's stronger result (asymptotic ratio var($\bar{x}$)/var($\bar{x}_1$) bounded below 1) required chain reversibility. MacEachern-Berliner need only stationarity, covering the typical non-reversible fixed-scan Gibbs sampler.
• Stationary time series trap. The spectral variance formula var$(\bar{x}) = (\sigma^2/nk)\{1 + 2\sum(1-i/nk)\rho_i\}$ can superficially suggest subsampling wins when $\rho_1 > 0.5$, but stationarity constraints on the autocorrelation sequence close off this possibility.
• Only computational rationale. The sole legitimate reason to subsample is computational: limited storage capacity before processing, or when processing cost per draw exceeds generation cost. In the common case where processing is cheap, subsampling wastes information.
• Preferred alternatives. Batching (block-mean variance estimator) and time-series spectral methods use the full chain and are uniformly preferred for variance estimation.

How It Works

The Inefficiency Proof (Section 2)

Let $x_1, \ldots, x_{nk}$ be a length-$nk$ run of a Gibbs sampler in its stationary distribution. Define $k$ systematic subsample means:

$$\bar{x}_i = \frac{1}{n}\sum_{j=0}^{n-1} x_{i+jk}, \quad i = 1, \ldots, k$$

By stationarity: $\text{var}(\bar{x}_i) = \text{var}(\bar{x}_1)$ and $|\text{cov}(\bar{x}_i, \bar{x}_j)| \leq \text{var}(\bar{x}_1)$ for all $i,j$. The overall mean $\bar{x} = \frac{1}{k}\sum_{i=1}^k \bar{x}_i$ satisfies:

$$\text{var}(\bar{x}) = \frac{\text{var}(\bar{x}_1)}{k} + \frac{1}{k^2}\sum_{i \neq j}\text{cov}(\bar{x}_i, \bar{x}_j) \leq \text{var}(\bar{x}_1)$$

Equality holds iff $\text{cov}(\bar{x}_i, \bar{x}_j) = \text{var}(\bar{x}_1)$ for all $i \neq j$, which forces $\Pr(\bar{x}_i = \bar{x}_1) = 1$ — a degenerate distribution. Hence the full chain always yields weakly lower variance.

The result extends to any function $f(x)$ estimated by $\sum_{j=0}^{n-1} f(x_{jk+i})/n$, provided $\text{var}(f(x)) < \infty$.

The Stationary Time Series Warning (Section 3)

The variance of the full-chain mean is:

$$\text{var}(\bar{x}) = \frac{\sigma^2}{nk}\left\{1 + 2\sum_{i=1}^{nk-1}\left(1-\frac{i}{nk}\right)\rho_i\right\} \tag{1}$$

The variance of the $k$-thinned mean is:

$$\text{var}(\bar{x}_1) = \frac{\sigma^2}{n}\left\{1 + 2\sum_{i=1}^{n-1}\left(1-\frac{i}{n}\right)\rho_{ik}\right\}$$

Example. Suppose $\rho_1 = \rho$ and $\rho_i = 0$ for $i > 1$, with $k = 2$. Then:

$$\text{var}(\bar{x}) = \frac{\sigma^2}{2n}\{1 + 2(1-1/2n)\rho\}, \qquad \text{var}(\bar{x}_1) = \frac{\sigma^2}{n}$$

For $\rho > 0.5$ and large $n$: $\text{var}(\bar{x}_1) < \text{var}(\bar{x})$ — subsampling appears to win. But stationarity with this autocorrelation structure requires $|\rho| \leq 0.5$, so the scenario never arises. The full sample is always better.

The key point: "The autocorrelations of a stationary process must satisfy a very large number of conditions" (Box-Jenkins 1976, p. 29). These constraints prevent the parameter combinations that would make thinning beneficial.

Concepts Introduced or Extended

• MCMC Convergence Diagnostics
• Gibbs Sampler

Entities Mentioned

• Steven N. MacEachern
• L. Mark Berliner

Quotes

"Hence the subsampled chain produces a poorer estimator than does the full chain."

"When, then, is subsampling worthwhile? Only when it is done for computational reasons: When storage space is at a premium and the output must be stored before being processed, or when the cost of processing the output outweighs the benefits of reduced variance of the estimator."

My Take

One of the most useful three-page papers in Bayesian computation. The proof is elementary — nothing beyond the variance of a sum — yet the implication is strong: the "thin to reduce autocorrelation" intuition is simply wrong from a variance standpoint. The stationarity-constraint argument in Section 3 is particularly instructive because it shows exactly why the naive time-series reasoning fails.

The paper's limitation is that it addresses only the mean estimator under a fixed total run budget of $nk$ draws. In practice, the thinning question usually arises under a fixed computation budget: the choice is between running $nk$ steps and storing $nk$ draws (expensive) vs. running $nk$ steps and storing $n$ draws. MacEachern-Berliner cover this case too (Section 4) and still recommend against thinning. The only genuine exception they acknowledge is when storage cost is the binding constraint before processing — a scenario that was more relevant in 1994 than today.