This is the chapter the rest of the engine rests on. A two-block Gibbs sampler is only useful here if we can say, for a specific model and dataset, how many sweeps separate a stored draw from the target — not asymptotically, but with a number. For the Gaussian model with fixed variance components that number is computable, and this chapter derives it.
The results are stated for the exactly-Gaussian case. Section “Beyond the Gaussian case” explains precisely what survives when the likelihood is Poisson, binomial, or Gamma, and what does not.
Start with one scalar hyper-mean \(\mu\) and one group coefficient \(\beta_b\), variance components fixed:
| Symbol | Role |
|---|---|
| \(\mu\) | population mean (Block 2) |
| \(\beta_b\) | group coefficient (Block 1) |
| \(P_{11}\) | prior precision on \(\mu\) |
| \(A = 1/\tau^2\) | coupling precision, \(\beta_b \mid \mu\) |
| \(C_b = x_b^2/\sigma^2\) | likelihood precision from group \(b\) |
Expanding the joint log-density as a quadratic form in \((\mu, \beta)\) gives the joint precision
\[ P = \begin{pmatrix} P_{11} + A & -A \\ -A & A + C_b \end{pmatrix}. \]
The off-diagonal is \(-A\): the two blocks are coupled only through the hierarchical prior, never through the data. For a two-dimensional Gaussian Gibbs sampler the contraction factor — the spectral radius of the Markov operator on the centered subspace — is the squared canonical correlation,
\[ r \;=\; \frac{P_{12}^2}{P_{11}^{\text{tot}} P_{22}} \;=\; \frac{A^2}{(P_{11}+A)(A+C_b)} . \]
Reading off the regimes:
| Regime | \(r\) | Behaviour |
|---|---|---|
| \(P_{11} \to \infty\) (tight prior on \(\mu\)) | \(\to 0\) | instant mixing |
| \(P_{11} > 0\), \(C_b = 0\) | \(0 < r < 1\) | geometrically ergodic |
| \(P_{11} \to 0\), \(C_b = 0\) | \(\to 1\) | mixing time \(\to \infty\) |
| \(P_{11} = 0\), \(C_b = 0\) | \(\det P = 0\) | improper joint, null recurrent |
The last row is the failure mode. When a group contributes no likelihood precision in some direction — a rank-deficient design — and the prior on \(\mu\) is flat in that direction, the joint posterior is improper and the chain does not converge at all.
Data rescues the chain. With \(C_b > 0\), even a single observation,
\[ r = \frac{A^2}{(P_{11}+A)(A+C_b)} \;<\; \frac{A}{P_{11}+A} \]
for any \(P_{11} \ge 0\). Under a flat prior on \(\mu\) this is still \(r \approx A/(A + C_b) < 1\): the data pins down \(\beta\), breaking the deterministic dependence on \(\mu\).
Notice what the formula says about design, not just propriety. \(r\) is close to 1 whenever \(C_b \ll A\), i.e. when each group carries little information relative to the shrinkage prior — small groups, or a small \(\tau^2\). Large, well-separated groups mix fast. This intuition is exactly what the multivariate rate below generalizes.
In the general model of Chapter-B00, order the
parameters so that the paper’s \(x_1\)
is the block updated second in each sweep — the Block 2
hyper vector \(\gamma\), of dimension
\(q\) — and \(x_2\) is the Block 1 stack \(b\), of dimension \(J p_{re}\). The joint posterior precision
blocks are
\[ \begin{aligned} P_{11} &= \sum_j H_j' P_b H_j + \mathrm{blockdiag}_k(V_k^{-1}) && (q \times q)\\ P_{22} &= \mathrm{blockdiag}_j\big(Z_j' W_j Z_j + P_b\big) && (Jp_{re} \times Jp_{re})\\ P_{12} &= \begin{bmatrix} -H_1'P_b & \cdots & -H_J'P_b \end{bmatrix} && (q \times Jp_{re}) \end{aligned} \]
where \(P_b = \Psi^{-1}\) is the
Block 1 prior precision, \(V_k\) is the
Block 2 prior covariance for coefficient \(k\), \(W_j\) collects the per-observation weights
(\(1/\sigma^2\) in the Gaussian case),
and \(H_j\) (\(p_{re} \times q\)) maps \(\gamma\) to group \(j\)’s Block 1 prior mean — row \(k\) of \(H_j\) is X_hyper[[k]][j, ]
placed in \(\gamma_k\)’s columns.
The Remark 8 rate matrix is
\[ A \;=\; P_{11}^{-1/2}\, P_{12}\, P_{22}^{-1}\, P_{21}\, P_{11}^{-1/2}, \]
a \(q \times q\) symmetric matrix with eigenvalues \(a_1 \ge a_2 \ge \dots \ge a_q\) in \([0, 1)\). The largest, \[ \lambda^* = a_1, \] is the geometric contraction rate of the two-block Gibbs sampler (Nygren 2020): the total-variation distance between the \(l\)-step kernel and the target decays like \((\lambda^*)^l\).
\(P_{22}\) is \(Jp_{re} \times Jp_{re}\) and is never built. Because \(P_{22}\) is block-diagonal by group and each \(H_j\) has exactly one non-zero segment per coefficient,
\[ P_{12} P_{22}^{-1} P_{21} \;=\; \sum_{j=1}^{J} H_j'\Big[P_b B_j^{-1} P_b\Big] H_j, \qquad B_j = Z_j' W_j Z_j + P_b, \]
accumulates over groups using only \(p_{re}
\times p_{re}\) solves. One \(q \times
q\) symmetric eigendecomposition finishes the job. Nothing larger
than \(q \times q\) is ever decomposed,
so the cost is linear in \(J\) and the
diagnostic is cheap enough to run before sampling — which is what
lmerb() does on every fit.
two_block_rate() implements exactly this.
For a balanced random-intercept model — \(J\) groups, \(n_j\) observations each, \(p_{re} = q = 1\), and a flat Block 2 prior — the algebra collapses to \[ \lambda^* = \frac{1/\tau^2}{1/\tau^2 + n_j/\sigma^2} = \frac{\sigma^2}{\sigma^2 + n_j \tau^2}, \] which is the scalar formula \(A/(A + C_b)\) with \(C_b = n_j/\sigma^2\).
J <- 12L # groups
n_j <- 8L # observations per group
sigma2 <- 1.0 # observation variance
tau2 <- 0.5 # random-intercept variance
group <- factor(rep(seq_len(J), each = n_j))
Z <- matrix(1, nrow = J * n_j, ncol = 1L,
dimnames = list(NULL, "(Intercept)"))
x_hyper <- list("(Intercept)" = matrix(1, nrow = J, ncol = 1L,
dimnames = list(levels(group),
"(Intercept)")))
rate <- two_block_rate(
x = Z,
group = group,
x_hyper = x_hyper,
prior_list_block1 = list(Sigma = matrix(tau2), dispersion = sigma2),
prior_list_block2 = list("(Intercept)" = list(mu = 0, Sigma = matrix(100)))
)
c(computed = rate$lambda_star,
closed_form = sigma2 / (sigma2 + n_j * tau2))
#> computed closed_form
#> 0.1999167 0.2000000The small discrepancy is the Block 2 prior: the closed form above assumes \(V \to \infty\), while the call used \(V = 100\). Tighten the prior and \(\lambda^*\) falls, exactly as the scalar table predicted.
The rate depends on \(n_j\) but not on \(J\): adding groups adds information about \(\gamma\) and adds parameters to \(b\) in equal measure. Making groups larger, on the other hand, helps a great deal.
lam <- function(n_j) {
g <- factor(rep(seq_len(J), each = n_j))
two_block_rate(
x = matrix(1, J * n_j, 1L, dimnames = list(NULL, "(Intercept)")),
group = g,
x_hyper = list("(Intercept)" = matrix(1, J, 1L,
dimnames = list(levels(g), "(Intercept)"))),
prior_list_block1 = list(Sigma = matrix(tau2), dispersion = sigma2),
prior_list_block2 = list("(Intercept)" = list(mu = 0, Sigma = matrix(100)))
)$lambda_star
}
round(vapply(c(2L, 4L, 8L, 16L, 32L), lam, numeric(1)), 4)
#> [1] 0.4998 0.3332 0.1999 0.1111 0.0588The printed form of a rate object adds the implied sweep counts for
common tolerances, which is the subject of Chapter-B03:
rate
#> Two-block Gibbs convergence rate (Nygren 2020, Remark 8)
#> groups J = 12, p_re = 1, q = 1, n_obs = 96 [family: gaussian]
#>
#> Eigenvalues of A (q x q):
#> [1] 0.199917
#>
#> lambda* = 0.199917
#>
#> Implied sweeps per tolerance:
#> tol (lambda*)^m proxy Theorem 3 bound Corollary 1 bound
#> 0.01 3 1 2
#> 0.001 5 2 3
#> 1e-06 9 4 5
#> (TV bounds assume start at the posterior mean: D0 = 0)Knowing that TV decays geometrically is not the same as knowing the constant. The bound actually used is Theorem 3 of (Nygren 2020), which gives the distance as an explicit sum over eigendirections rather than a single geometric envelope.
Order the eigenvalues ascending, \(a_0 = 0 \le a_1 \le \dots \le a_q\) (the Lemma 2 ordering). After \(l\) sweeps the bound is
\[ \mathrm{TV}\big(K^l(x^{(0)}, \cdot),\, \pi\big) \;\le\; \underbrace{\sum_{i=1}^{q} d_i^{(l)}}_{\text{variance convergence}} \;+\; \underbrace{(\lambda^*)^{l}\sqrt{D_0/(2\pi)}}_{\text{mean convergence}}, \]
capped at 1, where each \(d_i^{(l)}\) is a closed-form expression in the ratio \[ r_i^{(l)} = \frac{1 - a_{i-1}^{2l}}{1 - a_i^{2l}}, \] evaluated through the generalized error function \[ \mathrm{erf}_n(x) = P\big(\chi^2_n \le 2x^2\big), \] which reduces to the classical error function at \(n = 1\) (Nygren 2020, Remark 3; Brown 1963). The looser Corollary 1 form replaces the exact terms with a geometric envelope decaying like \(a_i^{2l}\) with explicit constants.
\(D_0\) is the squared standardized distance of the starting point from the posterior mean, \[ D_0 = (x^{(0)} - \mu)^{\top}\,\Sigma_{11}^{-1}\,(x^{(0)} - \mu). \]
Starting at the posterior mean sets \(D_0 = 0\) and the mean term vanishes
identically. This is the single most consequential
implementation choice in the package: lmerb() computes the
exact joint posterior mean by iterated conditional modes before sampling
and starts every replicate chain there, so only the variance-convergence
sum has to be certified. Chapter B03 shows what this buys in sweeps.
Two internal helpers evaluate and invert the bound:
lmebayesCore:::two_block_tv_bound(rate, l = 1:6) # theorem 3, D0 = 0
#> [1] 9.868989e-03 3.868176e-04 1.544796e-05 6.173849e-07 2.467480e-08
#> [6] 9.861698e-10
lmebayesCore:::two_block_l_for_tv(rate, tol = 1e-3) # sweeps for TV <= 1e-3
#> [1] 2The exact Theorem 3 terms are considerably sharper than the naive \((\lambda^*)^m\) proxy — often by a factor of two in sweeps — because the proxy ignores that all but the leading eigendirection have already converged.
Both helpers are \keyword{internal} topics: reach them
with lmebayesCore:::, or read them off the
print() output of a two_block_rate object,
which is the supported route.
A final bookkeeping point. The bound applies to the block updated
second in each sweep, i.e. \(\gamma\). The stored Block 1 draw lags by a
half-step, so the bound at \(l-1\)
applies to \(b\). The sampler adds one
sweep to compensate; that is why a certified m_convergence
is always one more than two_block_l_for_tv() returns.
Geometric ergodicity is not automatic. It requires the joint posterior to be proper in every direction, and that imposes two rank conditions — one at each level of the hierarchy.
A rank-deficient \(D_j\) contributes zero likelihood precision in the deficient direction, \(C_b = 0\), which is the degenerate scalar case above. In the \(p\)-dimensional Gaussian model with fixed \(\Sigma\), full column rank of \(D_j\) means \(C_j = D_j'D_j/\sigma^2\) is positive definite, so every eigendirection of the data block contributes positively to the Schur complement. Rank deficiency puts at least one direction back into the null-recurrent regime.
Full column rank of every \(D_j\) is therefore necessary for the joint posterior to be proper in all directions, and necessary for geometric ergodicity when the prior on \(\gamma\) is weak or improper.
The second condition is subtler and is the caller’s responsibility. Take the hyper-regression \(\beta_b \mid \gamma \sim N(x_b^{\top}\gamma, \tau^2)\) with group-level covariates \(x_b\). Integrating out a group with no data (\(C_b = 0\)) removes its contribution exactly: the Schur-complement correction cancels the coupling term it introduced. A data-free group is informationally equivalent to a group that never existed, and its covariate vector \(x_b\) is irrelevant.
What is left, after integrating out every non-identified group, is the marginal precision \[ P_{\gamma|\mathrm{marg}} \;=\; P_{11} \;+\; \sum_{b:\, C_b > 0} x_b x_b^{\top} \cdot \frac{A C_b}{A + C_b}. \]
With \(P_{11} = 0\) this is positive definite iff the group-level covariate vectors of the data-bearing groups span \(\mathbb{R}^q\): \[ \mathrm{rank}\big(X_{\mathrm{hyper}}[\text{data-bearing rows}, ]\big) = q . \]
| Level | Condition | Ensures |
|---|---|---|
| Data (per group) | qr(D_j)$rank == p_re |
\(C_j \succ 0\): \(\beta_j\) is identified by its own data |
| Hyper (across groups) | qr(X_hyper[kept rows, ])$rank == q |
\(\gamma\) is identified through the data-bearing groups |
Both must hold. If either fails there is a direction in \((\beta, \gamma)\) space along which the chain is null recurrent — it will still produce numbers, and they will be meaningless.
The condition is not vacuous: adding a group-level covariate that is
constant across all retained groups introduces a degenerate direction in
\(\gamma\) in exactly the way rank
deficiency in \(D_j\) introduces one in
\(\beta_j\).
check_identifiability() runs both checks; the
intercept-only special case (\(x_b =
1\) for all \(b\)) reduces level
2 to “at least one full-rank data-bearing group exists”.
The rate function also guards the boundary directly: if the computed
\(\lambda^*\) is numerically
indistinguishable from 1, two_block_rate() stops and prints
the offending precision blocks rather than returning a rate that would
certify nothing.
For Poisson, binomial, and Gamma responses the joint posterior is not normal and no theorem above applies as stated. What carries over, and what does not:
Carries over — the structural argument. For log-concave GLM likelihoods, \(C_b\) is replaced by the expected Fisher information at the current draw. A rank-deficient \(D_j\) still gives zero information in some direction, still makes that direction’s conditional improper under a weak prior on \(\gamma\), and still breaks ergodicity. Full-rank designs remain necessary. Sufficiency for the coupled chain needs additional conditions and is not claimed.
Does not carry over — the numerical rate. \(\lambda^*\) is only defined for a Gaussian
target. What the package does instead is evaluate the Gaussian rate at
IRLS-style curvature weights taken at the posterior mode: \[
w_i = \frac{\mathrm{wt}_i\,\mu'(\eta_i)^2}{V(\mu_i)\,\phi},
\qquad \eta_i = \mathrm{offset}_i + z_i^{\top} b_{\mathrm{mode},\,j(i)}
.
\] For canonical links this is the exact observed Hessian;
otherwise it is the expected (Fisher) information. Supplying these as
weights to two_block_rate() yields a
local-Gaussian heuristic rate — a good practical guide,
not a certificate. The print method labels such objects accordingly, and
non-Gaussian fits run a pilot stage rather than trusting the number
outright.
The same caveat applies, for a different reason, when variance
components are sampled rather than fixed: the target is then not
Gaussian in \((\gamma, \tau^2)\)
jointly. The package’s response is to evaluate the rate at the
conservative end of the truncated \(\tau^2\) support — smaller \(\tau^2\) means tighter coupling and a
larger \(\lambda^*\), so a bound
computed at disp_lower holds throughout the support.
Chapter-B05 covers the truncation window, and
Chapter-B06 covers the extended local rate diagnostic that
tracks \(\tau^2\) explicitly.
m_convergence, and the pilot/main cost trade-off.two_block_rate_ing(),
the local rate when \(\tau^2\) is
sampled.check_identifiability()
and the rank conditions in practice.?two_block_rate — the function reference, which also
documents two_block_rate_from_pfamily_list() and
two_block_rate_ing().inst/BLOCK_GIBBS_ERGODICITY.md in the package sources,
the maintainer note this chapter is drawn from.