--- title: "Chapter C03: Total variation distances between multivariate normal densities" author: "Kjell Nygren" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true number_sections: true bibliography: REFERENCES.bib reference-section-title: References vignette: > %\VignetteIndexEntry{Chapter C03: Total variation distances between multivariate normal densities} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r setup} library(lmebayesCore) ``` This chapter reproduces the theory of @Nygren2020, *On the total variation distance between multivariate normal densities with applications to two-block Gibbs samplers*, which is the reference the package's convergence machinery is built on. Definitions, remarks, lemmas, claims, theorems and the corollary are stated in the manuscript's own numbering so that code comments and help pages citing "Remark 8" or "Theorem 3" can be read against this chapter directly. The chapter is deliberately model-free. Everything here concerns two multivariate normal densities and a two-block Gibbs sampler in the abstract. The application to the mixed model of `Chapter-C01` — the precision blocks, the eigenvalue computation, the sweep budget — is in `Chapter-C04`. The paper's argument runs in two stages. Section 1 bounds the total variation distance between two multivariate normal densities that can be ordered by their precision matrices, by splitting the problem with the triangle inequality into a same-precision comparison and a same-mean comparison, each expressible through multivariate generalizations of the error function. Section 2 applies those bounds to the \(l\)-step kernel of a two-block Gibbs sampler, whose mean and variance after \(l\) sweeps are available in closed form, and concludes geometric ergodicity. # Total variation distances between multivariate normal densities ## Definitions The following definition of a generalized \(n\)-dimensional error function is due to @Brown1963. > **Definition 1.** A generalized \(n\)-dimensional error function is a > function \(\mathrm{erf}_n(\cdot) : \mathbb{R}^{+} \to \mathbb{R}\) > defined by > \[ > \mathrm{erf}_n(x) \;=\; > \frac{\int_0^{x} e^{-u^2} u^{\,n-1}\,du} > {\int_0^{\infty} e^{-u^2} u^{\,n-1}\,du}. > \] Brown also provided the explicit expansions. > **Remark 1.** The error function \(\mathrm{erf}_n(x)\) can be shown to be > explicitly expressed as > \[ > \begin{aligned} > \mathrm{erf}_{2m}(x) &= 1 - e^{-x^2}\left[1 + \frac{x^2}{1!} + \frac{x^4}{2!} + \dots + \frac{x^{2(m-1)}}{(m-1)!}\right] \\[4pt] > \mathrm{erf}_{2m+1}(x) &= \mathrm{erf}_1(x) - \frac{e^{-x^2}}{\sqrt{\pi}}\left[\frac{(2x)\,0!}{1!} + \frac{(2x)^3\,1!}{3!} + \dots + \frac{(2x)^{2m-1}(m-1)!}{(2m-1)!}\right] > \end{aligned} > \] > where the first expression applies if \(n\) is even and the second if > \(n\) is odd. > **Remark 2.** The expression \(e^{-u^2}u^{\,n-1}\) is maximized when > \(u^{*} = \sqrt{(n-1)/2}\), which in turn implies that for \(x \ge 0\) > \[ > \mathrm{erf}_n(x) \;\le\; > x\,\frac{e^{-(n-1)/2}\left[\sqrt{(n-1)/2}\right]^{n-1}} > {\int_0^{\infty} e^{-u^2}u^{\,n-1}du}, > \] > and \(x_1 \ge x_2 \ge 0\) implies that > \[ > \mathrm{erf}_n(x_1) - \mathrm{erf}_n(x_2) \;\le\; > \left[x_1 - x_2\right] > \frac{e^{-(n-1)/2}\left[\sqrt{(n-1)/2}\right]^{n-1}} > {\int_0^{\infty} e^{-u^2}u^{\,n-1}du}. > \] Remark 2 is the source of every geometric envelope later in the paper: it converts a difference of error functions into a difference of their arguments times an explicit constant. > **Remark 3.** Let \(f(\cdot \mid \mu, \Sigma)\) be a multivariate normal > density of dimension \(n\). It can then be verified that > \[ > \mathrm{Prob}\!\left(\sqrt{(x-\mu)^{T}\Sigma^{-1}(x-\mu)} \le \alpha\right) > \;=\; \mathrm{erf}_n\!\left(\frac{\alpha}{\sqrt{2}}\right). > \] Remark 3 is what ties the error function to a probability statement, and it is the identity the implementation evaluates: the quantity on the left is \(P(\chi^2_n \le \alpha^2)\). > **Definition 2.** The total variation distance between two probability > measures \(P\) and \(Q\) on a sigma-algebra \(\mathcal{F}\) of subsets of > the sample space \(\Omega\) is defined via > \[ > \|P(\cdot), Q(\cdot)\|_{TV} \;=\; \sup_{A \in \mathcal{F}} |P(A) - Q(A)|. > \] > **Definition 3.** [@KullbackLeibler1951; @Kullback1959; @Kullback1987] > If \(P\) and \(Q\) are probability measures over a set > \(\mathcal{X}\), and \(P\) is absolutely continuous with respect to > \(Q\), then the Kullback–Leibler divergence from \(Q\) to \(P\) is > defined as > \[ > D_{KL}(P \| Q) \;=\; \int_{\mathcal{X}} \log\!\left(\frac{dP}{dQ}\right) dP, > \] > where \(dP/dQ\) is the Radon–Nikodym derivative of \(P\) with respect to > \(Q\), and provided the expression on the right-hand side exists. The following theorem, known as Pinsker's inequality, relates the total variation distance to the Kullback–Leibler divergence. > **Theorem 1.** If \(P\) and \(Q\) are two probability distributions on > the measurable space \((X, \Sigma)\), then > \[ > \|P(\cdot), Q(\cdot)\|_{TV} \;\le\; \sqrt{\tfrac{1}{2} D_{KL}(P\|Q)}. > \] > **Claim 1.** Suppose we have two multivariate normal distributions, with > means \(\mu_1, \mu_2\) and with (non-singular) covariance matrices > \(\Sigma_1, \Sigma_2\). If the two distributions have the same > dimensions, \(m\), then > \[ > D_{KL}(N_1 \| N_2) = \tfrac{1}{2}\left( > \mathrm{tr}(\Sigma_2^{-1}\Sigma_1) > + (\mu_2-\mu_1)^{T}\Sigma_2^{-1}(\mu_2-\mu_1) > - m + \ln\!\left(\frac{|\Sigma_2|}{|\Sigma_1|}\right)\right). > \] Theorem 1 and Claim 1 together give one route to a bound. The rest of Section 1 develops a sharper one, which does not pass through Pinsker's inequality at all. ## Distance between multivariate normal densities with the same variance-covariance matrix The following result provides an exact formula for the distance between multivariate normal densities sharing the same variance-covariance matrix. The result relies on just the univariate error function. > **Lemma 1.** Let \(N(\cdot \mid \mu_1, \Sigma)\) and > \(N(\cdot \mid \mu_2, \Sigma)\) be two multivariate normal densities of > dimension \(n\). Then > \[ > \|N(\cdot\mid\mu_1,\Sigma), N(\cdot\mid\mu_2,\Sigma)\|_{TV} > = \mathrm{erf}_1\!\left(\frac{0.5\sqrt{(\mu_2-\mu_1)^{T}\Sigma^{-1}(\mu_2-\mu_1)}}{\sqrt{2}}\right). > \] This is an equality, not a bound: the distance depends on the two means only through the Mahalanobis length of their difference. > **Remark 4.** The expression for the error function is bounded above as > follows > \[ > \mathrm{erf}_1\!\left(\frac{0.5\sqrt{(\mu_2-\mu_1)^{T}\Sigma^{-1}(\mu_2-\mu_1)}}{\sqrt{2}}\right) > \le > \left[\frac{0.5\sqrt{(\mu_2-\mu_1)^{T}\Sigma^{-1}(\mu_2-\mu_1)}}{\sqrt{2}}\right] > \frac{e^{-(n-1)/2}\left[\sqrt{(n-1)/2}\right]^{n-1}}{\int_0^{\infty}e^{-u^2}u^{\,n-1}du}. > \] ## Distance between multivariate normal densities with the same mean vectors The following result provides bounds for the distance between multivariate normal densities that can be ordered based on their inverse variance (precision) matrices. It makes explicit use of eigenvalues derived from the two variance-covariance matrices. Essentially, the result involves bounding the distance between the two multivariate densities using repeated applications of the triangle inequality. > **Lemma 2.** Let \(N(\cdot\mid\mu,\Sigma_1)\) and > \(N(\cdot\mid\mu,\Sigma_2)\) be two multivariate normal densities of > dimension \(n\) satisfying the property that > \(\Sigma_1^{-1} - \Sigma_2^{-1}\) is positive semi-definite and that the > matrix \(\Sigma_2^{1/2}\Sigma_1^{-1}\Sigma_2^{1/2}\) is a symmetric > positive definite matrix with eigenvalues greater than or equal to 1. > Denote by \(1 \le k_1 \le k_2 \le \dots \le k_n\) the eigenvalues of the > matrix \(\Sigma_2^{1/2}\Sigma_1^{-1}\Sigma_2^{1/2}\) and let \(k_0 = 1\). > Define ratios \(r_1, r_2, \dots, r_n\) by \(r_i = k_i / k_{i-1}\) and > \[ > d_i = > \begin{cases} > \mathrm{erf}_{n+1-i}\!\left(\sqrt{\dfrac{r_i}{r_i-1}\cdot\dfrac{(n+1-i)\ln[r_i]}{2}}\right) > - \mathrm{erf}_{n+1-i}\!\left(\sqrt{\dfrac{1}{r_i-1}\cdot\dfrac{(n+1-i)\ln[r_i]}{2}}\right) > & r_i > 1 \\[10pt] > 0 & \text{otherwise.} > \end{cases} > \] > Then > \[ > \|N(\cdot\mid\mu,\Sigma_1), N(\cdot\mid\mu,\Sigma_2)\|_{TV} \;\le\; \sum_{i=1}^{n} d_i. > \] The ordering requirement matters. The eigenvalues must be taken in **ascending** order, which is what makes the interpolating sequence used in the proof monotone. > **Remark 5.** Leveraging the known bounds on \(\ln(r_i)\) given by > \((1 - 1/r_i) \le \ln(r_i) \le (r_i - 1)\) and the definition of \(r_i\), > it follows that the expressions for the error functions in the above are > bounded as follows: > \[ > \begin{aligned} > &\mathrm{erf}_{n+1-i}\!\left(\sqrt{\tfrac{r_i}{r_i-1}\tfrac{(n+1-i)\ln[r_i]}{2}}\right) > - \mathrm{erf}_{n+1-i}\!\left(\sqrt{\tfrac{1}{r_i-1}\tfrac{(n+1-i)\ln[r_i]}{2}}\right) \\[4pt] > &\qquad \le\; > \mathrm{erf}_{n+1-i}\!\left(\sqrt{\tfrac{r_i}{1}\tfrac{(n+1-i)}{2}}\right) > - \mathrm{erf}_{n+1-i}\!\left(\sqrt{\tfrac{1}{r_i}\tfrac{(n+1-i)}{2}}\right) \\[4pt] > &\qquad \le\; > \left[\sqrt{r_i} - \sqrt{\tfrac{1}{r_i}}\right] > \sqrt{\frac{n+1-i}{2}} > \left[\frac{e^{-(n-i)/2}\left[\sqrt{(n-i)/2}\right]^{n-i}}{\int_0^{\infty}e^{-u^2}u^{\,n-i}du}\right]. > \end{aligned} > \] ## Distances between multivariate normal densities A simple application of the two previous results now implies the following bound for the distance between two ordered multivariate normal densities. > **Theorem 2.** Let \(N(\cdot\mid\mu_1,\Sigma_1)\) and > \(N(\cdot\mid\mu_2,\Sigma_2)\) be two multivariate normal densities of > dimension \(n\) satisfying the property that > \(\Sigma_1^{-1}-\Sigma_2^{-1}\) is positive semi-definite and that the > matrix \(\Sigma_2^{1/2}\Sigma_1^{-1}\Sigma_2^{1/2}\) is a symmetric > positive definite matrix with eigenvalues greater than or equal to 1. > Denote by \(1 \le k_1 \le \dots \le k_n\) the eigenvalues of > \(\Sigma_2^{1/2}\Sigma_1^{-1}\Sigma_2^{1/2}\), let \(k_0 = 1\), define > \(r_i = k_i/k_{i-1}\) and \(d_i\) as in Lemma 2. Then > \[ > \|N(\cdot\mid\mu_1,\Sigma_1) - N(\cdot\mid\mu_2,\Sigma_2)\|_{TV} > \;\le\; \left[\sum_{i=1}^{n} d_i\right] > + \mathrm{erf}_1\!\left(\frac{0.5\sqrt{(\mu_2-\mu_1)^{T}\Sigma_2^{-1}(\mu_2-\mu_1)}}{\sqrt{2}}\right). > \] *Proof.* By the triangle inequality and our two earlier claims, we have \[ \begin{aligned} \|N(\cdot\mid\mu_1,\Sigma_1), N(\cdot\mid\mu_2,\Sigma_2)\|_{TV} &\le \|N(\cdot\mid\mu_1,\Sigma_1), N(\cdot\mid\mu_1,\Sigma_2)\|_{TV} \\ &\quad + \|N(\cdot\mid\mu_1,\Sigma_2), N(\cdot\mid\mu_2,\Sigma_2)\|_{TV} \\ &\le \left[\sum_{i=1}^{n} d_i\right] + \mathrm{erf}_1\!\left(\frac{0.5\sqrt{(\mu_2-\mu_1)^{T}\Sigma_2^{-1}(\mu_2-\mu_1)}}{\sqrt{2}}\right). \end{aligned} \] \(\square\) The structure of the bound is now visible: **a variance term and a mean term**, obtained by passing through the intermediate density \(N(\cdot\mid\mu_1,\Sigma_2)\), which shares its mean with the first and its covariance with the second. > **Remark 6.** Combining the results from the earlier sections, we have > the following relationship > \[ > \begin{aligned} > &\mathrm{erf}_1\!\left(\frac{0.5\sqrt{(\mu_2-\mu_1)^{T}\Sigma_2^{-1}(\mu_2-\mu_1)}}{\sqrt{2}}\right) > + \left[\sum_{i=1}^{n} d_i\right] \\[4pt] > &\quad\le > \left[\frac{0.5\sqrt{(\mu_2-\mu_1)^{T}\Sigma^{-1}(\mu_2-\mu_1)}}{\sqrt{2}}\right] > \left[\frac{e^{-(n-i)/2}\left[\sqrt{(n-i)/2}\right]^{n-i}}{\int_0^{\infty}e^{-u^2}u^{\,n-i}du}\right] \\[4pt] > &\qquad + \sum_{i=1}^{n}\left[\left[\sqrt{r_i}-\sqrt{\tfrac{1}{r_i}}\right] > \sqrt{\frac{n+1-i}{2}} > \left[\frac{e^{-(n-i)/2}\left[\sqrt{(n-i)/2}\right]^{n-i}}{\int_0^{\infty}e^{-u^2}u^{\,n-i}du}\right]\right]. > \end{aligned} > \] # Convergence rates for two-block Gibbs samplers for multivariate normal densities ## Definitions > **Definition 4.** A two-block Gibbs sampler for a multivariate normal > density is a sampler that generates discrete time, time homogeneous > Markov chains \(X_1 = \{x_1^{(l)}, l = 0, 1, \dots\}\) and > \(X_2 = \{x_2^{(l)}, l = 0, 1, \dots\}\) for which the joint transition > kernel \(Q(\cdot,\cdot)\) is given by > \[ > \begin{aligned} > Q\!\left(\left(x_1^{(l-1)}, x_2^{(l-1)}\right), \left(x_1^{(l)}, x_2^{(l)}\right)\right) > =\;& N_1\!\left(x_1^{(l)} \,\middle|\, \mu_1 + \Sigma_{12}\Sigma_{22}^{-1}\left(x_2^{(l)}-\mu_2\right),\; \Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\right) \\ > &\times N_2\!\left(x_2^{(l)} \,\middle|\, \mu_2 + \Sigma_{12}\Sigma_{22}^{-1}\left(x_1^{(l-1)}-\mu_1\right),\; \Sigma_{22}-\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}\right). > \end{aligned} > \] > As the densities on the right above do not depend on \(x_2^{(l-1)}\), we > will simply write \(Q\!\left(x_1^{(l-1)}, \left(x_1^{(l)},x_2^{(l)}\right)\right)\) > instead of \(Q\!\left(\left(x_1^{(l-1)},x_2^{(l-1)}\right),\left(x_1^{(l)},x_2^{(l)}\right)\right)\). Note the order: **the second block is updated first, and the first block is updated second**. All the results below are statements about \(X_1\), the block updated second. > **Definition 5.** We will denote by \(Q_{X_1}(\cdot,\cdot)\) the marginal > transition kernel defined by > \(Q_{X_1}\!\left(x_1^{(l-1)}, x_1^{(l)}\right) = \int Q\!\left(x_1^{(l-1)},\left(x_1^{(l)},x_2^{(l)}\right)\right)dx_2^{(l)}\) > and we will use the notation \(Q\!\left(x_1^{(l-1)}, \cdot \mid x_1^{(l)}\right)\) > to denote the conditional transition kernel for \(X_2\) defined by > \[ > Q\!\left(x_1^{(l-1)}, x_2^{(l)} \,\middle|\, x_1^{(l)}\right) = > \frac{Q\!\left(x_1^{(l-1)},\left(x_1^{(l)},x_2^{(l)}\right)\right)} > {\int Q\!\left(x_1^{(l-1)},\left(x_1^{(l)},x_2^{(l)}\right)\right)dx_2^{(l)}}. > \] > **Definition 6.** For \(l = 1, 2, 3, \dots\), we will let > \(Q^{(l)}_{X_1}(\cdot,\cdot)\) denote the \(l\)-step transition kernel > for \(X_1\) defined by > \[ > Q^{(l)}_{X_1}\!\left(x_1^{(0)}, x_1^{(l)}\right) = \int > Q_{X_1}\!\left(x_1^{(0)},x_1^{(1)}\right) > Q_{X_1}\!\left(x_1^{(1)},x_1^{(2)}\right)\cdots > Q_{X_1}\!\left(x_1^{(l-1)},x_1^{(l)}\right) > dx_1^{(1)}dx_1^{(2)}\cdots dx_1^{(l-1)}. > \] > **Remark 7.** In the case of a two-block Gibbs sampler the \(l\)-step > transition kernel is itself multivariate normal, so the distance between > the \(l\)-step transition kernel and its target density satisfies the > property that > \[ > \left\|Q^{(l)}_{X_1}\!\left(x_1^{(0)},\cdot\right), N_1(\cdot\mid\mu_1,\Sigma_{11})\right\|_{TV} > = \left\|N_1\!\left(\cdot \,\middle|\, E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right], \Sigma^{(l)}_{11}\right), N_1(\cdot\mid\mu_1,\Sigma_{11})\right\|_{TV}. > \] Remark 7 is the hinge of the whole paper. It says the \(l\)-step kernel is *itself* a multivariate normal density, so the general bound of Theorem 2 applies to it directly — provided we can compute its mean and its covariance. Those are the subjects of the next two subsections. > **Definition 7.** The two-block Gibbs sampler above will be said to be > geometrically ergodic for the Markov chain \(X_1\) if there exists a > function \(g(\cdot)\) and a constant \(0 < t < 1\) such that for any > \(x_1^{(0)}\), > \[ > \left\|Q^{(l)}_{X_1}\!\left(x_1^{(0)},\cdot\right), N_1(\cdot\mid\mu_1,\Sigma_{11})\right\|_{TV} > \;\le\; g\!\left(x_1^{(0)}\right) t^{\,k} > \] > for \(l = 1, 2, 3, \dots\). ## Convergence of mean vectors > **Claim 2.** Let \(N(\cdot\mid\mu,\Sigma)\) be a multivariate normal > density. Consider a two-block Gibbs sampler applied to this density in > which the second block is updated first with variance \(P_{22}^{-1}\) and > the first block is updated second with variance \(P_{11}^{-1}\). The > expectation on the \(l\)-th observation is related to the previous > iteration through the iterative formula > \[ > E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right] - \mu_1 > = \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\Sigma_{11}^{-1}\left(E\!\left[x_1^{(l-1)}\middle|x_1^{(0)}\right]-\mu_1\right) > = P_{11}^{-1}P_{12}P_{22}^{-1}P_{21}\left(E\!\left[x_1^{(l-1)}\middle|x_1^{(0)}\right]-\mu_1\right) > \] > and hence > \[ > \begin{aligned} > &\left[E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right]-\mu_1\right]^{T}\Sigma_{11}^{-1}\left[E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right]-\mu_1\right] \\ > &\quad = \left[E\!\left[x_1^{(l-1)}\middle|x_1^{(0)}\right]-\mu_1\right]^{T} > \Sigma_{11}^{-1/2} > \left[\Sigma_{11}^{-1/2}\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\Sigma_{11}^{-1/2}\right] > \left[\Sigma_{11}^{-1/2}\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\Sigma_{11}^{-1/2}\right] > \Sigma_{11}^{-1/2}\left[\cdot\right] \\ > &\quad \le (\lambda^{*})^{2}\left[E\!\left[x_1^{(l-1)}\middle|x_1^{(0)}\right]-\mu_1\right]^{T}\Sigma_{11}^{-1}\left[E\!\left[x_1^{(l-1)}\middle|x_1^{(0)}\right]-\mu_1\right] \\ > &\quad \le (\lambda^{*})^{2l}\left[x_1^{(0)}-\mu_1\right]^{T}\Sigma_{11}^{-1}\left[x_1^{(0)}-\mu_1\right] > \end{aligned} > \] > where \(\lambda^{*}\) is the maximal eigenvalue of the matrix > \(\Sigma_{11}^{-1/2}\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\Sigma_{11}^{-1/2}\). > **Remark 8.** The eigenvalues for the matrix > \(\Sigma_{11}^{-1/2}\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\Sigma_{11}^{-1/2}\) > coincide with the eigenvalues for the matrix > \(P_{11}^{-1/2}P_{12}P_{22}^{-1}P_{21}P_{11}^{-1/2}\). Remark 8 is what makes the rate computable. The covariance form requires \(\Sigma\), which is the inverse of the joint precision; the precision form requires only the blocks of \(P\), which are read directly off the quadratic form of the log posterior. This is the identity `two_block_rate()` relies on. > **Remark 9.** The above also implies the following: > \[ > \begin{aligned} > E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right] - \mu_1 &= T^{\,l}\left[\left(x_1^{(0)}-\mu_1\right)\right] \\ > \mu_1 &= \left[I - T^{\,l}\right]^{-1}\left[E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right] - T^{\,l}x_1^{(0)}\right] > \end{aligned} > \] > and > \[ > \begin{aligned} > E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right] - \mu_1 > &= \left[I-T^{\,l}\right]^{-1}\left[\left[I-T^{\,l}\right]E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right] - \left[E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right]-T^{\,l}x_1^{(0)}\right]\right] \\ > &= \left[I-T^{\,l}\right]^{-1}T^{\,l}\left[x_1^{(0)} - E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right]\right]. > \end{aligned} > \] > **Remark 10.** The above result now implies the following relationship > (since \(\Sigma_{11}^{-1} = P_{11}-P_{12}P_{22}^{-1}P_{21}\)) > \[ > \begin{aligned} > &\mathrm{erf}_1\!\left(\frac{0.5\sqrt{\left(E[x_1^{(l)}|x_1^{(0)}]-\mu_1\right)^{T}\Sigma_{11}^{-1}\left(E[x_1^{(l)}|x_1^{(0)}]-\mu_1\right)}}{\sqrt{2}}\right) \\ > &\quad = \mathrm{erf}_1\!\left(\frac{0.5\sqrt{\left[x_1^{(0)}-E[x_1^{(l)}|x_1^{(0)}]\right]^{T}T^{\,l}\left[I-T^{\,l}\right]^{-1}\Sigma_{11}^{-1}\left[I-T^{\,l}\right]^{-1}T^{\,l}\left[x_1^{(0)}-E[x_1^{(l)}|x_1^{(0)}]\right]}}{\sqrt{2}}\right) \\ > &\quad \le \mathrm{erf}_1\!\left(\frac{0.5\sqrt{\left[x_1^{(0)}-E[x_1^{(l)}|x_1^{(0)}]\right]^{T}T^{\,l}\left[I-T^{\,l}\right]^{-1}P_{11}\left[I-T^{\,l}\right]^{-1}T^{\,l}\left[x_1^{(0)}-E[x_1^{(l)}|x_1^{(0)}]\right]}}{\sqrt{2}}\right). > \end{aligned} > \] > It is worth noting that any upper bound for \(P_{22}(\beta)\) can be used > in applications below. ## Convergence of variance-covariance matrices > **Claim 3.** Let \(N(\cdot\mid\mu,\Sigma)\) be a multivariate normal > density. Consider a two-block Gibbs sampler applied to this density in > which the second block is updated first with variance \(P_{22}^{-1}\) and > the first block is updated second with variance \(P_{11}^{-1}\). Let > \(A = P_{11}^{-1/2}P_{12}P_{22}^{-1}P_{21}P_{11}^{-1/2}\). After \(l\) > iterations, the variance for the first block is given by > \[ > \Sigma^{(l)}_{11} = P_{11}^{-1/2}\left[\sum_{i=1}^{2l} A^{\,i-1}\right]P_{11}^{-1/2} > \] > and > \[ > \Sigma_{11}^{-1/2}\Sigma^{(l)}_{11}\Sigma_{11}^{-1/2} = \left[I - A^{(2l)}\right]. > \] The exponent is \(2l\), not \(l\): **the variance converges at twice the rate of the mean**. This is an equality, which is what makes it usable as a diagnostic rather than merely a bound. > **Remark 11.** Denote by \(a_i\) the \(i\)th eigenvalue for the matrix > \(A\) and let \(a_0 = 0\). Then after \(l\) iterations, the constants > \(k_i^{(l)}, r_i^{(l)}\) are given by > \[ > k_i^{(l)} = \frac{1}{1-a_i^{2l}} > \qquad\text{and}\qquad > r_i^{(l)} = \frac{1-a_{i-1}^{2l}}{1-a_i^{2l}} > \] > respectively. Remark 11 is the bridge between the two halves of the paper. The \(k_i\) of Lemma 2 — eigenvalues of a ratio of covariance matrices — become explicit functions of the eigenvalues of \(A\) and the sweep count \(l\). ## Convergence bounds and geometric ergodicity for the two-block Gibbs sampler > **Theorem 3.** Let \(N(\cdot\mid\mu,\Sigma)\) be a multivariate normal > density. Consider a two-block Gibbs sampler applied to this density in > which the second block is updated first with variance \(P_{22}^{-1}\) and > the first block is updated second with variance \(P_{11}^{-1}\). Then > \[ > \begin{aligned} > &\left\|N\!\left(\cdot\middle| E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right],\Sigma^{(l)}_{11}\right), N(\cdot\mid\mu_1,\Sigma_{11})\right\|_{TV} \\ > &\quad\le \left\|N\!\left(\cdot\middle|E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right],\Sigma^{(l)}_{11}\right), N\!\left(\cdot\middle|E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right],\Sigma_{11}\right)\right\|_{TV} \\ > &\qquad + \left\|N\!\left(\cdot\middle|E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right],\Sigma_{11}\right), N(\cdot\mid\mu_1,\Sigma_{11})\right\|_{TV} \\[4pt] > &\quad\le \left[\sum_{i=1}^{n} d_i^{(l)}\right] > + \mathrm{erf}_1\!\left(\frac{0.5\sqrt{\left(E[x_1^{(l)}|x_1^{(0)}]-\mu_1\right)^{T}\Sigma_{11}^{-1}\left(E[x_1^{(l)}|x_1^{(0)}]-\mu_1\right)}}{\sqrt{2}}\right) \\[4pt] > &\quad = \left[\sum_{i=1}^{n} d_i^{(l)}\right] > + \mathrm{erf}_1\!\left(\frac{0.5\sqrt{\left(x_1^{(0)}-\mu_1\right)^{T}\left(T^{(l)}\right)^{T}\Sigma_{11}^{-1}\left(T^{(l)}\right)\left(x_1^{(0)}-\mu_1\right)}}{\sqrt{2}}\right) \\[4pt] > &\quad\le \left[\sum_{i=1}^{n} d_i^{(l)}\right] > + \mathrm{erf}_1\!\left(\frac{0.5(\lambda^{*})^{l}\sqrt{\left(x_1^{(0)}-\mu_1\right)^{T}\Sigma_{11}^{-1}\left(x_1^{(0)}-\mu_1\right)}}{\sqrt{2}}\right) > \end{aligned} > \] > where the constant \(r_i^{(l)}\) in each \(d_i^{(l)}\) term is given by > \(r_i^{(l)} = \dfrac{1-a_{i-1}^{2l}}{1-a_i^{2l}}\). The first inequality is Theorem 2's triangle-inequality split, applied to the \(l\)-step kernel via Remark 7. The variance term is Lemma 2 evaluated at the Remark 11 ratios; the mean term is Lemma 1 evaluated at the Claim 2 contraction. > **Corollary 1.** Let \(N(\cdot\mid\mu^{*},\Sigma)\) be a multivariate > normal density. Consider a two-block Gibbs sampler applied to this > density in which the second block is updated first with variance > \(P_{22}^{-1}\) and the first block is updated second with variance > \(P_{11}^{-1}\). Then > \[ > \begin{aligned} > &\left\|N\!\left(\cdot\middle|E\!\left[x_1^{(l)}\middle|x_1^{(0)}\right],\Sigma^{(l)}_{11}\right), N(\cdot\mid\mu_1,\Sigma_{11})\right\|_{TV} \\ > &\quad\le \frac{0.5(\lambda^{*})^{l}\sqrt{\left(x_1^{(0)}-\mu_1\right)^{T}\Sigma_{11}^{-1}\left(x_1^{(0)}-\mu_1\right)}}{\sqrt{2}} > \left[\frac{e^{-(n-i)/2}\left[\sqrt{(n-i)/2}\right]^{n-i}}{\int_0^{\infty}e^{-u^2}u^{\,n-i}du}\right] \\ > &\qquad + \sum_{i=1}^{n}\left[\frac{a_i^{2l}}{\sqrt{1-a_i^{2}}\sqrt{1-a_{i-1}^{2}}}\right] > \left[\frac{e^{-(n-i)/2}\left[\sqrt{(n-i)/2}\right]^{n-i}}{\int_0^{\infty}e^{-u^2}u^{\,n-i}du}\right] > \end{aligned} > \] > where the constants \(\lambda^{*}\) and \(a_i\) in each term are defined > as above. Hence the two-block Gibbs sampler is geometrically ergodic. Every term decays like \(a_i^{2l}\) (or \((\lambda^{*})^{l}\) for the mean term), with explicit constants, which is Definition 7's requirement with \(t = \lambda^{*}\). # Proofs ## Proof for densities with shared mean vector *Proof of Lemma 2.* Let \(C\Lambda C^{-1} = C\Lambda C^{T}\) be an eigenvalue decomposition for \(\Sigma_2^{1/2}\Sigma_1^{-1}\Sigma_2^{1/2}\) and assume without loss of generality that the eigenvalues are non-decreasing. Let \(A_0 = I_n\) and for \(0 < i \le n\) let \(A_i\) be the diagonal matrix with diagonal elements given by \[ a_{i,jj} = \begin{cases} k_j & j \le i \\ k_i & j \ge i. \end{cases} \] From the triangle inequality, it follows that \[ \begin{aligned} \|f(\cdot\mid\mu,\Sigma_1), f(\cdot\mid\mu,\Sigma_2)\|_{TV} &= \left\|f\!\left(\cdot\middle|\Sigma_2^{-1/2}\mu, \left[\Sigma_2^{1/2}\Sigma_1^{-1}\Sigma_2^{1/2}\right]^{-1}\right), f\!\left(\cdot\middle|\Sigma_2^{-1/2}\mu, I_n\right)\right\|_{TV} \\ &= \left\|f\!\left(\cdot\middle|\Sigma_2^{-1/2}\mu, \left[C\Lambda C^{T}\right]^{-1}\right), f\!\left(\cdot\middle|\Sigma_2^{-1/2}\mu, \left[CI_nC^{T}\right]^{-1}\right)\right\|_{TV} \\ &= \left\|f\!\left(\cdot\middle|C^{-1}\Sigma_2^{-1/2}\mu, \Lambda^{-1}\right), f\!\left(\cdot\middle|C^{-1}\Sigma_2^{-1/2}\mu, I_n\right)\right\|_{TV} \\ &= \left\|f\!\left(\cdot\middle|\tilde\mu, A_n^{-1}\right), f\!\left(\cdot\middle|\tilde\mu, A_0^{-1}\right)\right\|_{TV} \\ &\le \sum_{i=1}^{n}\left\|f\!\left(\cdot\middle|\tilde\mu, A_i^{-1}\right), f\!\left(\cdot\middle|\tilde\mu, A_{i-1}^{-1}\right)\right\|_{TV} \end{aligned} \] where \(\tilde\mu = C^{-1}\Sigma_2^{-1/2}\mu\). If \(r_i = 1\) then \(A_i = A_{i-1}\) and the corresponding distance is 0. If \(r_i > 1\), it is straightforward to verify that \[ \begin{aligned} f\!\left(x\middle|\tilde\mu, A_i^{-1}\right) &\ge f\!\left(x\middle|\tilde\mu, A_{i-1}^{-1}\right) \\ \iff\quad (x-\tilde\mu)^{T}\left[A_i - A_{i-1}\right](x-\tilde\mu) &\le \log\left[|A_i|\right]-\log\left[|A_{i-1}|\right] \\ \iff\quad \sum_{j=i}^{n}(x_j-\tilde\mu_j)^{T}\left[k_i-k_{i-1}\right](x_j-\tilde\mu_j) &\le \sum_{j=i}^{n}\left[\log[k_i]-\log[k_{i-1}]\right] \\ \iff\quad \sum_{j=i}^{n}(x_j-\tilde\mu_j)^{T}k_{i-1}(x_j-\tilde\mu_j) &\le \frac{1}{r_i-1}(n+1-i)\log[r_i] \\ \iff\quad \sum_{j=i}^{n}(x_j-\tilde\mu_j)^{T}k_i(x_j-\tilde\mu_j) &\le \frac{r_i}{r_i-1}(n+1-i)\log[r_i]. \end{aligned} \] Let \(d_i = \left\|f(\cdot\mid\tilde\mu, A_i^{-1}), f(\cdot\mid\tilde\mu, A_{i-1}^{-1})\right\|_{TV}\). From the last two inequalities above and Remark 3, it now follows that \[ d_i = \begin{cases} \mathrm{erf}_{n+1-i}\!\left(\sqrt{\tfrac{r_i}{r_i-1}\tfrac{(n+1-i)\ln[r_i]}{2}}\right) - \mathrm{erf}_{n+1-i}\!\left(\sqrt{\tfrac{1}{r_i-1}\tfrac{(n+1-i)\ln[r_i]}{2}}\right) & r_i > 1 \\ 0 & \text{otherwise.} \end{cases} \] No further standardization is required: in the first case the vector defined by \(\sqrt{k_{i-1}}(x_j-\tilde\mu_j)\) for \(j = i,\dots,n\) is a multivariate standard normal, and in the second case the vector defined by \(\sqrt{k_i}(x_j-\tilde\mu_j)\) for \(j = i,\dots,n\) is also multivariate standard normal. \(\square\) ## Proof for densities with shared variance-covariance matrix *Proof of Lemma 1.* Let \(P = \Sigma^{-1}\). We then have \[ \begin{aligned} f(x\mid\mu_1,\Sigma) &\ge f(x\mid\mu_2,\Sigma) \\ \iff\quad 0.5(x-\mu_1)^{T}P(x-\mu_1) &\le 0.5(x-\mu_2)^{T}P(x-\mu_2) \\ \iff\quad (x-\mu_1)^{T}P(\mu_2-\mu_1) &\le 0.5(\mu_2-\mu_1)^{T}P(\mu_2-\mu_1) \\ \iff\quad \frac{(x-\mu_1)^{T}P(\mu_2-\mu_1)}{\sqrt{(\mu_2-\mu_1)^{T}P(\mu_2-\mu_1)}} &\le 0.5\sqrt{(\mu_2-\mu_1)^{T}P(\mu_2-\mu_1)}. \end{aligned} \] It is straightforward to verify that \[ E\!\left[\frac{(x-\mu_1)^{T}P(\mu_2-\mu_1)}{\sqrt{(\mu_2-\mu_1)^{T}P(\mu_2-\mu_1)}}\,\middle|\,\mu_1,\Sigma\right] = 0 \] and that \[ E\!\left[\frac{(x-\mu_1)^{T}P(\mu_2-\mu_1)}{\sqrt{(\mu_2-\mu_1)^{T}P(\mu_2-\mu_1)}}\,\middle|\,\mu_2,\Sigma\right] = \sqrt{(\mu_2-\mu_1)^{T}P(\mu_2-\mu_1)}. \] From this it follows that the event above has probability \(\Phi\!\left(0.5\sqrt{(\mu_2-\mu_1)^{T}P(\mu_2-\mu_1)}\right)\) under \(\mu_1\) and \(\Phi\!\left(-0.5\sqrt{(\mu_2-\mu_1)^{T}P(\mu_2-\mu_1)}\right)\) under \(\mu_2\), which in turn implies that \[ \begin{aligned} \|f(\cdot\mid\mu_1,\Sigma), f(\cdot\mid\mu_2,\Sigma)\|_{TV} &= \Phi\!\left(0.5\sqrt{(\mu_2-\mu_1)^{T}P(\mu_2-\mu_1)}\right) - \Phi\!\left(-0.5\sqrt{(\mu_2-\mu_1)^{T}P(\mu_2-\mu_1)}\right) \\ &= \mathrm{erf}_1\!\left(\frac{0.5\sqrt{(\mu_2-\mu_1)^{T}P(\mu_2-\mu_1)}}{\sqrt{2}}\right). \qquad\square \end{aligned} \] ## Proof for convergence of variance-covariance matrices > **Remark 12.** Let \(f(\cdot\mid\mu^{*},P^{-1})\) be a multivariate > normal density. Then the block-diagonal of the variance-covariance matrix > satisfies: > \[ > \begin{aligned} > \Sigma_{11}^{-1} &= \left[P_{11}-P_{12}P_{22}^{-1}P_{21}\right] & > \Sigma_{22}^{-1} &= \left[P_{22}-P_{21}P_{11}^{-1}P_{12}\right] \\ > \Sigma_{12} &= -P_{11}^{-1}P_{12}\left(P_{22}-P_{21}P_{11}^{-1}P_{12}\right)^{-1} & > \Sigma_{12}\Sigma_{22}^{-1} &= -P_{11}^{-1}P_{12} \\ > \Sigma_{21}\Sigma_{11}^{-1} &= -P_{22}^{-1}P_{21}. & & > \end{aligned} > \] > **Remark 13.** Let \(f(\cdot\mid\mu^{*},\Sigma)\) be a multivariate > normal density. Consider a two-block Gibbs sampler applied to this > density in which the second block is updated first with variance > \(P_{22}^{-1}\) and the first block is updated second with variance > \(P_{11}^{-1}\). Then the variance of the first block after \(l\) > iterations is related to the variance after \(l-1\) iterations through > the iterative formula > \[ > \begin{aligned} > \Sigma^{(l)}_{11} > &= P_{11}^{-1} + \Sigma_{12}\Sigma_{22}^{-1}P_{22}^{-1}\Sigma_{22}^{-1}\Sigma_{21} > + \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\Sigma_{11}^{-1}\Sigma^{(l-1)}_{11}\Sigma_{11}^{-1}\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21} \\ > &= P_{11}^{-1} + P_{11}^{-1}P_{12}P_{22}^{-1}P_{21}P_{11}^{-1} > + P_{11}^{-1}P_{12}P_{22}^{-1}P_{21}\Sigma^{(l-1)}_{11}P_{12}P_{22}^{-1}P_{21}P_{11}^{-1}. > \end{aligned} > \] *Proof of Claim 3.* For \(l = 1\), we have \[ \begin{aligned} \Sigma^{(1)}_{11} &= P_{11}^{-1} + P_{11}^{-1}P_{12}P_{22}^{-1}P_{21}P_{11}^{-1} = P_{11}^{-1/2}\left[I + P_{11}^{-1/2}P_{12}P_{22}^{-1}P_{21}P_{11}^{-1/2}\right]P_{11}^{-1/2} \\ &= P_{11}^{-1/2}\left[I+A\right]P_{11}^{-1/2} = P_{11}^{-1/2}\left[\sum_{i=1}^{2l}A^{\,i-1}\right]P_{11}^{-1/2} \end{aligned} \] so the result holds. In order to demonstrate the result for \(l > 1\) assume that the result holds for \(l-1\). By induction, we have \[ \begin{aligned} \Sigma^{(l)}_{11} &= P_{11}^{-1}+P_{11}^{-1}P_{12}P_{22}^{-1}P_{21}P_{11}^{-1} + P_{11}^{-1}P_{12}P_{22}^{-1}P_{21}\Sigma^{(l-1)}_{11}P_{12}P_{22}^{-1}P_{21}P_{11}^{-1} \\ &= P_{11}^{-1/2}\left[I+A\right]P_{11}^{-1/2} + P_{11}^{-1/2}\left[A\right]\left[\sum_{i=1}^{2(l-1)}A^{\,i-1}\right]\left[A\right]P_{11}^{-1/2} \\ &= P_{11}^{-1/2}\left[I+A\right]P_{11}^{-1/2} + P_{11}^{-1/2}\left[\sum_{i=3}^{2l}A^{\,i-1}\right]P_{11}^{-1/2} = P_{11}^{-1/2}\left[\sum_{i=1}^{2l}A^{\,i-1}\right]P_{11}^{-1/2}. \end{aligned} \] To see the second part of the result, we note that \[ P_{11}^{-1/2}\sum_{i=1}^{2l}A^{\,i-1}P_{11}^{-1/2} = P_{11}^{-1/2}\left[I-A\right]^{-1/2}\left[I-A^{2l}\right]\left[I-A\right]^{-1/2}P_{11}^{-1/2}. \] Since \(A\) is a convergent matrix, it follows that \[ \Sigma_{11} = P_{11}^{-1/2}\left[I-A\right]^{-1/2}\left[I-A\right]^{-1/2}P_{11}^{-1/2}, \qquad \Sigma_{11}^{1/2} = P_{11}^{-1/2}\left[I-A\right]^{-1/2}. \] It follows that \(\Sigma^{(l)}_{11} = \Sigma_{11}^{1/2}\left[I-A^{2l}\right]\Sigma_{11}^{1/2}\) and hence \(\Sigma_{11}^{-1/2}\Sigma^{(l)}_{11}\Sigma_{11}^{-1/2} = \left[I-A^{2l}\right]\). \(\square\) ## Proof for convergence of mean vectors > **Remark 14.** Let \(f(\cdot\mid\mu^{*},P^{-1})\) be a multivariate > normal density. Then the conditional means are given by > \[ > E[\mu_1\mid\mu_2] = \mu_1^{*}+\Sigma_{12}\Sigma_{22}^{-1}(\mu_2-\mu_2^{*}), > \qquad > E[\mu_2\mid\mu_1] = \mu_2^{*}+\Sigma_{21}\Sigma_{11}^{-1}(\mu_1-\mu_1^{*}). > \] > **Claim 4.** Let \(f(\cdot\mid\mu^{*},\Sigma)\) be a multivariate normal > density. Consider a two-block Gibbs sampler applied to this density in > which the second block is updated first with variance \(P_{22}^{-1}\) and > the first block is updated second with variance \(P_{11}^{-1}\). The > expectation on the \(l\)-th observation is related to the previous > iteration through the iterative formula > \[ > E\!\left[\mu_1^{(l)}\middle|\mu_1^{(0)}\right]-\mu_1^{*} > = \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\Sigma_{11}^{-1}\left(E\!\left[\mu_1^{(l-1)}\middle|\mu_1^{(0)}\right]-\mu_1^{*}\right). > \] > **Remark 15.** In the multivariate normal case above, it follows that > \[ > \mu_1^{*} = \left[I-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\Sigma_{11}^{-1}\right]^{-1} > \left[E\!\left[\mu_1^{(l)}\middle|\mu_1^{(0)}\right]-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\Sigma_{11}^{-1}E\!\left[\mu_1^{(l-1)}\middle|\mu_1^{(0)}\right]\right]. > \] ## Proof of geometric ergodicity > **Remark 16.** Combining the results from the earlier sections, we have > the following relationship > \[ > \begin{aligned} > &\mathrm{erf}_1\!\left(\frac{0.5\sqrt{(\mu_2-\mu_1)^{T}\Sigma_2^{-1}(\mu_2-\mu_1)}}{\sqrt{2}}\right)+\left[\sum_{i=1}^{n}d_i\right] \\ > &\quad\le \left[\frac{0.5\sqrt{(\mu_2-\mu_1)^{T}\Sigma^{-1}(\mu_2-\mu_1)}}{\sqrt{2}}\right]\left[\frac{1}{\int_0^{\infty}e^{-u^2}u^{\,n-1}du}\right] \\ > &\qquad + \sum_{i=1}^{n}\left[\sqrt{r_i}-\sqrt{\tfrac{1}{r_i}}\right]\left[\frac{n+1-i}{2}\right]\left[\frac{e^{-(n-i)/2}\left[\sqrt{(n-i)/2}\right]^{n-i}}{\int_0^{\infty}e^{-u^2}u^{\,n-i}du}\right]. > \end{aligned} > \] > **Remark 17.** Denote by \(a_i\) the \(i\)th eigenvalue for the matrix > \(A\) and let \(a_0 = 0\). Then after \(l\) iterations, the constants > \(k_i^{(l)}, r_i^{(l)}\), and the expression > \(\sqrt{k_i^{(l)}}-\sqrt{1/k_i^{(l)}}\) are given by > \[ > k_i^{(l)} = \frac{1}{1-a_i^{2l}}, > \qquad > r_i^{(l)} = \frac{1-a_{i-1}^{2l}}{1-a_i^{2l}} > \] > and > \[ > \sqrt{r_i^{(l)}}-\sqrt{\frac{1}{r_i^{(l)}}} > = \frac{a_i^{2l}-a_{i-1}^{2l}}{\sqrt{1-a_i^{2l}}\sqrt{1-a_{i-1}^{2l}}} > \;\le\; \frac{a_i^{2l}}{\sqrt{1-a_i^{2}}\sqrt{1-a_{i-1}^{2}}}, > \qquad > \sqrt{k_i^{(l)}}-\sqrt{\frac{1}{k_i^{(l)}}} = \frac{a_i^{2l}}{\sqrt{1-a_i^{2l}}} \le \frac{a_i^{2l}}{\sqrt{1-a_i^{2}}} > \] > respectively. Remark 17 is what turns the \(d_i^{(l)}\) of Theorem 3 into the geometric terms of Corollary 1: the ratio-based quantity is bounded by \(a_i^{2l}\) times a constant that no longer depends on \(l\). > **Remark 18.** The above result now implies the following relationships > \[ > \begin{aligned} > \mathrm{erf}_1\!\left(\frac{0.5\sqrt{\left(E[\mu_1^{(l)}|\mu_1^{(0)}]-\mu_1^{*}\right)^{T}\Sigma_{11}^{-1}\left(E[\mu_1^{(l)}|\mu_1^{(0)}]-\mu_1^{*}\right)}}{\sqrt{2}}\right) > &\le \mathrm{erf}_1\!\left(\frac{0.5(\lambda^{*})^{l}\sqrt{\left(\mu_1^{(0)}-\mu_1^{*}\right)^{T}\Sigma_{11}^{-1}\left(\mu_1^{(0)}-\mu_1^{*}\right)}}{\sqrt{2}}\right) \\ > &\le \frac{0.5(\lambda^{*})^{l}\sqrt{\left(\mu_1^{(0)}-\mu_1^{*}\right)^{T}\Sigma_{11}^{-1}\left(\mu_1^{(0)}-\mu_1^{*}\right)}}{\sqrt{2}}\left[\frac{1}{\int_0^{\infty}e^{-u^2}u^{\,n-1}du}\right]. > \end{aligned} > \] ## Special case derivations In the multivariate normal case, \[ \left\|\pi(\cdot\mid\nu,\Sigma_{11}) - Q^{(l)}\!\left(\cdot\middle|\nu,\Sigma^{(l)}_{11}\right)\right\|_{TV} = \left\|\pi(\cdot\mid\nu,I) - Q^{(l)}\!\left(\cdot\middle|\nu,\Sigma_{11}^{-1/2}\Sigma^{(l)}_{11}\Sigma_{11}^{-1/2}\right)\right\|_{TV} = \left\|\pi(\cdot\mid\nu,I) - Q^{(l)}\!\left(\cdot\middle|\nu,\left[I-A^{(2l)}\right]\right)\right\|_{TV}. \] # How the package evaluates these results The bounds above are implemented in `two_block_ergodicity.R`. Three translation steps turn the manuscript's expressions into elementary computations. **The error function is a chi-squared probability.** By Remark 3, \(\mathrm{erf}_n(x) = P\!\left(\chi^2_n \le 2x^2\right)\), so the internal helper is one line: ```{r erfn} lmebayesCore:::.two_block_erfn ``` This agrees with the closed forms of Remark 1. Taking \(n = 2\) (\(m = 1\), so \(\mathrm{erf}_2(x) = 1 - e^{-x^2}\)) and \(n = 1\): ```{r erfn-check} x <- c(0.25, 0.5, 1, 2) rbind( erf2_remark1 = 1 - exp(-x^2), erf2_impl = lmebayesCore:::.two_block_erfn(x, 2L), erf1_classic = 2 * pnorm(x * sqrt(2)) - 1, erf1_impl = lmebayesCore:::.two_block_erfn(x, 1L) ) ``` **Each \(d_i^{(l)}\) is a difference of two chi-squared probabilities.** Writing \(m_i = n+1-i\) and \(r = r_i^{(l)}\), squaring and doubling the arguments of Lemma 2 gives \[ d_i^{(l)} = P\!\left(\chi^2_{m_i} \le \frac{m_i\, r \log r}{r-1}\right) - P\!\left(\chi^2_{m_i} \le \frac{m_i \log r}{r-1}\right), \] which is what `two_block_tv_bound(method = "theorem3")` sums, with the eigenvalues sorted ascending as Lemma 2 requires and \(a_0 = 0\) as in Remark 11. **The Corollary 1 envelope is assembled from Remarks 5, 17 and 18.** Writing \(c_m = 2e^{-m/2}(m/2)^{m/2}/\Gamma\!\left(\frac{m+1}{2}\right)\) for the bracketed constant — using \(\int_0^{\infty}e^{-u^2}u^{\,m}du = \Gamma\!\left(\frac{m+1}{2}\right)/2\) — the implemented variance term is \[ d_i^{(l)} \;\le\; \frac{a_i^{2l}}{\sqrt{1-a_i^{2}}\sqrt{1-a_{i-1}^{2}}} \sqrt{\frac{n+1-i}{2}}\; c_{\,n-i}, \] and the mean term is \((\lambda^{*})^{l}\sqrt{D_0/(2\pi)}\), which is Remark 18's final line with \(n = 1\). Two points where the implementation departs from the manuscript's display are worth recording. First, the Corollary 1 display omits the factor \(\sqrt{(n+1-i)/2}\) that appears in Remark 5; it is required for the chain of inequalities from Remark 5 through Remark 17 to close, and the implementation includes it. Second, the Corollary 1 display carries a free index \(i\) in the mean term's constant, which does not depend on \(i\); the implementation uses Remark 18's form, in which that constant is evaluated at \(n = 1\) and reduces to \(1/\sqrt{2\pi}\) after combining with the leading \(0.5/\sqrt{2}\). Finally, \(D_0 = \left(x_1^{(0)}-\mu_1\right)^{T}\Sigma_{11}^{-1}\left(x_1^{(0)}-\mu_1\right)\) is the only place the starting point enters either bound. Starting a chain at the posterior mean sets \(D_0 = 0\) and removes the mean term identically, leaving only the variance sum — which is why the package computes the exact posterior mean before sampling and starts every replicate chain there. # See also - `vignette("Chapter-C04")` — these results applied to the mixed model of `Chapter-C01`: the precision blocks, the eigenvalue spectrum, the sweep budget, and the Claim 3 diagnostic measured empirically. - `vignette("Chapter-B02")` — the same results stated in engine terms. - `vignette("Chapter-B03")` — turning \(\lambda^{*}\) into `m_convergence`. - `?two_block_rate`, `?plot_var_convergence` - `vignette("Chapter-C05")` — a different route past the exactly-Gaussian restriction: Doeblin minorization on a certified set, which yields a geometric total variation bound for the two-block chain restricted to that set on hierarchical GLMM targets, where Theorem 2's eigenvalue machinery does not apply. - `inst/SAFE_UNSAFE_TV_DECOMPOSITION.md` in the package sources — the safe/unsafe region decomposition that generalizes Theorem 2's comparison-density argument to precision structures beyond the exactly-Gaussian case here, including the heuristic that feeds Corollary 1's bound into it. - `inst/BLOCK_GIBBS_ERGODICITY.md`, `inst/BLOCK_GIBBS_ERGODICITY_ING.md` and `inst/ELLIPSOID_TV_BOUND.md` in the package sources, which extend this theory to rank conditions, sampled variance components, and non-Gaussian likelihoods respectively.