m_convergenceChapter-B02 produced a computable bound on the
total-variation distance between the \(l\)-step kernel and the target. This
chapter spends it. The question is operational: given a total-variation
tolerance tv_tol, how many inner Gibbs sweeps should
separate the chain start from a stored draw?
The answer is m_convergence, and the user never has to
choose it. They choose a tolerance; the engine derives the sweeps.
pfamily_list + design
|
v
two_block_rate_from_pfamily_list() --> lambda*, eigenvalue spectrum
|
v
two_block_l_for_tv(rate, tv_tol) --> smallest l with TV(l) <= tv_tol
|
v
m_convergence = l + 1 (+1: Block 1 lags a half-step)Three properties make this practical.
The bound is decreasing in \(l\). So inverting it is exact: a doubling search followed by bisection finds the smallest sufficient \(l\), with no tuning constants.
The cost is logarithmic in \(1/\texttt{tv\_tol}\). Tightening the tolerance by three orders of magnitude does not triple the work; it adds a handful of sweeps. This is what makes the “certify every draw” contract affordable.
The chain starts at the posterior mean. \(D_0 = 0\), the mean term of the bound vanishes, and only the variance-convergence sum remains.
J <- 12L
n_j <- 8L
sigma2 <- 1.0
tau2 <- 0.5
group <- factor(rep(seq_len(J), each = n_j))
Z <- matrix(1, J * n_j, 1L, dimnames = list(NULL, "(Intercept)"))
x_hyper <- list("(Intercept)" = matrix(1, J, 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)))
)
rate$lambda_star
#> [1] 0.1999167The print method lays out the whole calibration table at once — the naive geometric proxy alongside the two certified bounds:
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)Read the columns as three increasingly good answers to the same
question. The (lambda*)^m proxy asks only when the
slowest eigendirection has decayed to tol,
ignoring that the others converged sweeps ago. Corollary 1 keeps the
geometric shape but supplies the correct constants. Theorem 3 evaluates
the exact per-direction terms and is the sharpest of the three; it is
what the sampler uses.
tv_tol meanstv_tol is a per-draw guarantee: each
stored draw is within tv_tol of the exact joint posterior
in total variation. It is not a statement about the sample as a
whole.
The default is tv_tol = 0.01, the conventional threshold
of the honest burn-in literature (Jones and
Hobert 2001). To certify an entire n-draw sample at
level \(\alpha\), pass
tv_tol = alpha / n — a union bound over draws. Because the
sweep cost is logarithmic, this is cheap:
tols <- c(1e-2, 1e-3, 1e-4, 1e-6, 1e-8)
data.frame(
tv_tol = tols,
sweeps = vapply(tols,
function(t) lmebayesCore:::two_block_l_for_tv(rate, t),
integer(1))
)
#> tv_tol sweeps
#> 1 1e-02 1
#> 2 1e-03 2
#> 3 1e-04 3
#> 4 1e-06 4
#> 5 1e-08 6Going from a per-draw tolerance of \(10^{-2}\) to a sample-wide certificate at \(n = 1000\) (i.e. \(10^{-5}\)) costs a few sweeps, not a few thousand.
The bound itself, evaluated directly:
The Theorem 3 bound applies to the block updated second in each sweep, which is Block 2 (\(\gamma\)). The stored Block 1 draw (\(b\)) is produced before that Block 2 update, so it lags by a half-step: the bound at \(l-1\) is what applies to it.
The sampler therefore sets
m_convergence = two_block_l_for_tv(rate, tv_tol) + 1so that both blocks of every stored draw satisfy the
tolerance. If you compare a hand-computed
two_block_l_for_tv() against a fit’s reported
m_convergence and find an off-by-one, this is it.
Everything so far assumed \(D_0 =
0\). That is exact for the Gaussian fixed-variance-component
case, where lmerb_posterior_mean() computes the joint
posterior mean in closed form by iterated conditional modes.
Two situations break it:
glmerb_posterior_mode()), but the mode is not
the mean, so a chain started there is at some unknown positive distance
from the target mean.In both cases the engine runs a pilot stage: a batch
of n_pilot short chains whose cross-chain mean is used as
the start for the main chains. The point of the pilot is not to produce
draws — they are discarded — but to produce a starting point with a
certified distance from the truth.
If the pilot mean is an average over n_pilot independent
chains, its sampling distribution has covariance \(\Sigma_{11}/n_{\mathrm{pilot}}\), so the
squared standardized distance obeys \[
D_0 \;=\; \frac{\chi^2_{p,\,\alpha}}{n_{\mathrm{pilot}}}
\] at confidence level \(\alpha\) (default 0.95), where \(p\) is the dimension of the Block 2 hyper
vector. Larger pilots give smaller \(D_0\), which in turn lowers the main-stage
m_convergence — but the pilot itself costs sweeps.
That is a genuine optimization. Total inner-sweep cost is \[ C(n_{\mathrm{pilot}}) \;=\; n_{\mathrm{pilot}}\, m_p \;+\; n \cdot m_{\mathrm{conv}}(n_{\mathrm{pilot}}), \] where \(m_p\) is sweeps per pilot draw and \(m_{\mathrm{conv}}(\cdot)\) is the main-stage sweep count implied by the \(D_0\) above. The first term grows linearly in \(n_{\mathrm{pilot}}\); the second decreases, and flattens out once \(D_0\) is small enough that the mean term no longer binds. There is an interior minimum.
Four internal helpers form the stack:
| Function | Role |
|---|---|
two_block_d0_pilot_start() |
\(D_0 = \chi^2_{p,\alpha}/n_{\mathrm{pilot}}\) |
two_block_m_convergence_for_pilot_start() |
main-stage sweeps at that \(D_0\) |
two_block_pilot_sampling_cost() |
total cost \(C\) |
two_block_optimize_pilot_cost() |
minimize \(C\) over \(n_{\mathrm{pilot}}\) |
opt <- lmebayesCore:::two_block_optimize_pilot_cost(
n = 1000L,
rate = rate,
tv_tol = 0.01,
m_convergence_pilot = 5L,
p = 1L
)
opt[c("n_pilot_opt", "m_convergence_opt", "total_cost_opt")]
#> $n_pilot_opt
#> [1] 11
#>
#> $m_convergence_opt
#> [1] 3
#>
#> $total_cost_opt
#> [1] 3055head(opt$cost_curve, 10)
#> n_pilot total_cost
#> 1 6 4030
#> 2 7 4035
#> 3 8 4040
#> 4 9 4045
#> 5 10 4050
#> 6 11 3055
#> 7 12 3060
#> 8 13 3065
#> 9 14 3070
#> 10 15 3075The curve is flat near its minimum, which is good news: the optimum
does not have to be found precisely for the plan to be near-optimal.
This is the Gibbs analogue of the build-versus-sample cost optimization
the iid envelope samplers perform when deciding how finely to partition
the envelope (glmbayesCore::EnvelopeOpt; see
vignette("Chapter-A08", package = "glmbayes")).
Note the asymmetry with the mode-start minimum. As \(n_{\mathrm{pilot}} \to \infty\), \(D_0 \to 0\) and \(m_{\mathrm{conv}}\) approaches the \(D_0 = 0\) sweep count from
Chapter-B02 — the pilot can never do better than starting
at the exact mean, only approach it.
pfamily_listIn practice the Block 2 priors arrive as pfamily objects
from pfamily_list(), not as raw precision matrices.
two_block_rate_from_pfamily_list() is the adapter:
rate <- two_block_rate_from_pfamily_list(
x = design$Z,
group = design$group,
x_hyper = design$X_hyper,
prior_list_block1 = list(Sigma = Sigma_ranef, dispersion = sigma2),
pfamily_list = pfamily_list(ps)
)It extracts the dispersion each component implies:
dNormal — the fixed
dispersion, i.e. the known \(\tau^2_k\).dIndependent_Normal_Gamma — the
conservative disp_lower plug-in, the lower end of the
truncation window.The second choice is the one that makes the calibration honest for
ING models. Smaller \(\tau^2\) means
tighter cross-block coupling and a larger \(\lambda^*\), so a bound computed at
disp_lower is valid for every \(\tau^2\) in the truncated support. It is
conservative rather than correct-on-average, which is the right
direction for a certificate. That conservatism is also why ING routes
carry a heavier eigenvalue bound than their fixed-\(\tau^2\) counterparts, and why
disp_lower must be set with some care — see
Chapter-B05.
two_block_l_for_tv() searches up to l_max
sweeps. If the bound has still not reached the tolerance — which happens
when \(\lambda^*\) is very close to 1 —
it issues a warning and returns the uncertified geometric
proxy \(l\) from
rate$m_for_tol(tol) rather than l_max. The
result is a usable sweep count that is explicitly not a certificate.
Upstream of that, two_block_rate() refuses to return at
all if \(\lambda^*\) is numerically
indistinguishable from 1: it prints the \(P_{11}\), \(P_{12}\), \(P_{21}\), \(P_{22}\) blocks and stops. A \(\lambda^*\) of 1 means the joint posterior
is improper in some direction, so the honest response is a failed rank
condition (Chapter-B08), not a larger sweep count.
Between those, a warning fires whenever \(\lambda^*\) exceeds 0.95: convergence is still geometric, but the coupling is strong enough that the model is worth a second look. The usual causes are small groups (small \(C_j\) relative to \(1/\tau^2\)) or a variance component driven toward zero.
disp_lower.plot_mean_convergence() and
plot_var_convergence().?two_block_rate, ?two_block_tv_bound,
?two_block_optimize_pilot_cost.