Blocking Correlated Variables in Gibbs Sampling

Why Update Correlated Variables Together?

Gibbs sampling updates variables by drawing from their conditional distributions while holding the other coordinates fixed. This is convenient when those conditionals are easy to sample, but it can be inefficient when the target distribution has strong dependence. For a correlated Gaussian, each coordinate conditional is concentrated around a line determined by the other coordinate. Updating one variable at a time can therefore keep the chain close to its current location, even though the joint target extends along a much longer diagonal direction.

Blocking addresses this geometry by updating a group of variables jointly. We will use a bivariate Gaussian to derive the conditional distributions, prove that the updates preserve the target, and calculate the movement of a systematic single-component Gibbs sampler. These calculations describe dependence between successive states; they do not by themselves establish a general convergence rate.

The Bivariate Gaussian and Its Conditionals

Work on a probability space $(\Omega,\mathcal{F},\mathbb{P})$. The state space is $S=\mathbb{R}^2$, and a discrete-time Gibbs chain is written $Z_n=(X_n,Y_n)$, $n\geq 0$, with filtration $\mathcal{F}_n=\sigma(Z_0,\ldots,Z_n)$. Consider a target distribution for $(X,Y)$ that is bivariate normal, with means $\mu_x,\mu_y$, positive standard deviations $\sigma_x,\sigma_y$, and correlation $\rho$, where $|\rho|<1$. Its covariance matrix is positive definite under this condition.

Definition: A single-component Gibbs update replaces one coordinate by a draw from its full conditional distribution given the current values of the other coordinates. A block Gibbs update replaces a group of coordinates by a joint draw from their conditional distribution given the coordinates outside the block. Every draw described here is assumed to be generated from the specified conditional kernel.

For the bivariate Gaussian, standardize by writing $u=(x-\mu_x)/\sigma_x$ and $v=(y-\mu_y)/\sigma_y$. The joint density has exponent proportional to $-[u^2-2\rho uv+v^2]/[2(1-\rho^2)]$. Completing the square gives $u^2-2\rho uv+v^2=(u-\rho v)^2+(1-\rho^2)v^2$. Thus, conditional on $Y=y$, the part of the density depending on $x$ is a normal density with mean $\mu_x+\rho(\sigma_x/\sigma_y)(y-\mu_y)$ and variance $\sigma_x^2(1-\rho^2)$. Completing the square in $v$ gives the other conditional.

Conditional distributions: For every real $y$ and $x$, the Gaussian density gives positive marginal density at $y$ and $x$, so these conditional distributions are defined: $$ X\mid Y=y\sim N\!\left(\mu_x+\rho\frac{\sigma_x}{\sigma_y}(y-\mu_y),\ \sigma_x^2(1-\rho^2)\right), \qquad Y\mid X=x\sim N\!\left(\mu_y+\rho\frac{\sigma_y}{\sigma_x}(x-\mu_x),\ \sigma_y^2(1-\rho^2)\right). $$

The condition $|\rho|<1$ matters: it makes both conditional variances positive and the joint distribution have a density on $\mathbb{R}^2$. At $\rho=1$ or $\rho=-1$, the Gaussian is degenerate, so these density calculations do not apply.

Why Exact Gibbs Updates Preserve the Target

Theorem: If $(X,Y)$ has the bivariate Gaussian target above, replacing $X$ by a draw from the conditional distribution of $X$ given the current $Y$ leaves the joint target distribution unchanged. Replacing $Y$ by a draw from its conditional distribution given the current $X$ also leaves it unchanged. Consequently, a sequence of these exact updates, in either fixed order, preserves the target when the chain starts with that target distribution.

Proof. Let $f(x,y)$ be the joint density, $f_Y(y)$ the marginal density of $Y$, and $f_{X\mid Y}(x\mid y)=f(x,y)/f_Y(y)$ the conditional density. In an update of $X$, the old $X$ is discarded and a new value $x'$ is drawn using $f_{X\mid Y}(\cdot\mid y)$, while $Y=y$ is retained. The joint density of the output $(x',y)$ is therefore $f_Y(y)f_{X\mid Y}(x'\mid y)=f(x',y)$, exactly the original joint density. For an update of $Y$, the output density is $f_X(x)f_{Y\mid X}(y'\mid x)=f(x,y')$, again the target density. If a first update preserves the target and a second update preserves it, their composition does too. Repeating this argument proves the claim for any finite sequence of the stated updates. $\square$

This result concerns invariance, not speed: it says the target remains a distribution of the chain, not how quickly the chain explores it from other initial distributions. The same density argument applies to a block: resampling a group from its joint conditional distribution preserves the target, because the marginal density of the variables outside the block multiplied by the block conditional density equals the joint density.

Movement Under a Systematic Single-Component Sweep

To see the effect of correlation directly, take the standardized target with means zero, variances one, and correlation $\rho$. In each sweep, first update $X$ given the current $Y$, then update $Y$ given the newly drawn $X$. Let $Y_n$ denote the value after sweep $n$. Using independent standard normal random variables $\varepsilon_{n+1}$ and $\eta_{n+1}$ for the two conditional draws gives $X_{n+1}=\rho Y_n+\sqrt{1-\rho^2}\,\varepsilon_{n+1}$ and $Y_{n+1}=\rho X_{n+1}+\sqrt{1-\rho^2}\,\eta_{n+1}$. Substitution yields the one-coordinate recurrence below.

Theorem: Under the systematic update just specified, when the chain starts in the standardized target distribution, $\operatorname{Var}(Y_n)=1$ and $\operatorname{Cov}(Y_n,Y_{n+1})=\rho^2$. Also, $\mathbb{E}[(Y_{n+1}-Y_n)^2]=2(1-\rho^2)$. A joint block draw of $(X,Y)$ from the entire bivariate target, in contrast, gives a new $Y$ independent of the previous $Y$, with $\mathbb{E}[(Y_{n+1}-Y_n)^2]=2$ when the chain starts in the target.

Proof. Substituting the first update into the second gives $$ Y_{n+1}=\rho^2Y_n+\rho\sqrt{1-\rho^2}\,\varepsilon_{n+1} +\sqrt{1-\rho^2}\,\eta_{n+1}. $$ The two new normal variables are independent of the current state, have mean zero and variance one, and are independent of each other. The added noise consequently has variance $\rho^2(1-\rho^2)+(1-\rho^2)=(1-\rho^2)(1+\rho^2)=1-\rho^4$. When $Y_n$ has variance one and mean zero, the recurrence gives $\operatorname{Var}(Y_{n+1})=\rho^4+(1-\rho^4)=1$, so the marginal variance is preserved. Its covariance with the previous value is $\operatorname{Cov}(Y_n,Y_{n+1})=\rho^2\operatorname{Var}(Y_n)=\rho^2$, because the added noise is independent of $Y_n$. Therefore $\mathbb{E}[(Y_{n+1}-Y_n)^2]=\operatorname{Var}(Y_{n+1})+\operatorname{Var}(Y_n)-2\operatorname{Cov}(Y_n,Y_{n+1})=2-2\rho^2$. For a joint block draw, $Y_{n+1}$ and $Y_n$ are independent standard normal variables when the chain starts in the target. Their difference has mean zero and variance $1+1=2$, proving the final claim. $\square$

The coefficient $\rho^2$ describes the one-sweep linear dependence of $Y_{n+1}$ on $Y_n$ for this systematic order. When $|\rho|$ is close to one, $\rho^2$ is close to one and the single-component sweep changes $Y$ by little on average. A block update of both variables avoids that particular coordinate-by-coordinate persistence. The comparison is specific to this Gaussian target and update scheme; it is not a claim that every possible block is computationally cheaper or mixes faster.

Worked Example: Deriving the Conditional Means and Variances

Let the target have $\mu_x=2$, $\mu_y=-1$, $\sigma_x=2$, $\sigma_y=\tfrac12$, and $\rho=\tfrac35$. Since $|\rho|=\tfrac35<1$, the conditional distributions above apply. First, $1-\rho^2=1-\tfrac9{25}=\tfrac{16}{25}$. For $X\mid Y=y$, the conditional mean is $2+\tfrac35(2/(\tfrac12))(y+1)=2+\tfrac{12}{5}(y+1)$, and the conditional variance is $2^2(\tfrac{16}{25})=\tfrac{64}{25}$. For $Y\mid X=x$, the conditional mean is $-1+\tfrac35((\tfrac12)/2)(x-2)=-1+\tfrac3{20}(x-2)$, and the conditional variance is $(\tfrac12)^2(\tfrac{16}{25})=\tfrac4{25}$.

A joint block draw can also be constructed explicitly. Let $Z_1,Z_2$ be independent standard normal variables and set $X=2+2Z_1$, $Y=-1+\tfrac3{10}Z_1+\tfrac25Z_2$. Then $\mathbb{E}[X]=2$ and $\mathbb{E}[Y]=-1$. The variance of $X$ is $2^2=4$; the variance of $Y$ is $(\tfrac3{10})^2+(\tfrac25)^2=\tfrac9{100}+\tfrac{16}{100}=\tfrac14$; and the covariance is $2(\tfrac3{10})=\tfrac35$. Thus the standard deviations are $2$ and $\tfrac12$, and the correlation is $(\tfrac35)/(2\cdot\tfrac12)=\tfrac35$, as required.

Check: The calculated conditional variances are positive, and the constructed joint draw has covariance $\tfrac35$. Both verify the target parameters used in the conditional and block calculations.

Worked Example: A Moderate-Correlation Comparison

For the standardized target with $\rho=\tfrac35$, a systematic sweep has conditional updates $X_{n+1}=\tfrac35Y_n+\tfrac45\varepsilon_{n+1}$ and $Y_{n+1}=\tfrac35X_{n+1}+\tfrac45\eta_{n+1}$, since $\sqrt{1-\rho^2}=\sqrt{\tfrac{16}{25}}=\tfrac45$. Substituting the first equation into the second gives $Y_{n+1}=\tfrac9{25}Y_n+\tfrac{12}{25}\varepsilon_{n+1}+\tfrac45\eta_{n+1}$. The innovation variance is $(\tfrac{12}{25})^2+(\tfrac45)^2=\tfrac{144}{625}+\tfrac{400}{625}=\tfrac{544}{625}=1-\tfrac{81}{625}=1-\rho^4$. Thus a unit-variance $Y_n$ produces a unit-variance $Y_{n+1}$, and the lag-one covariance is $\tfrac9{25}$.

The expected squared change in $Y$ for this single-component sweep is $2(1-\tfrac9{25})=\tfrac{32}{25}$. For a joint block draw, the expected squared change is $2$. Their ratio, block movement divided by single-component movement, is $2/(\tfrac{32}{25})=\tfrac{25}{16}$. These are second-moment comparisons under the target distribution, not guarantees about every realized step.

Check: The recurrence coefficient $\tfrac9{25}$ and innovation variance $\tfrac{544}{625}$ give stationary variance $(\tfrac{544}{625})/(1-(\tfrac9{25})^2)=(\tfrac{544}{625})/(\tfrac{544}{625})=1$, agreeing with the target.

Worked Example: Strong Correlation Makes Single Updates Small

Now take $\rho=\tfrac{99}{100}$, with both marginal variances equal to one. The single-component recurrence coefficient is $\rho^2=(\tfrac{99}{100})^2=\tfrac{9801}{10000}$, while the conditional variance for each coordinate update is $1-\rho^2=1-\tfrac{9801}{10000}=\tfrac{199}{10000}$. The innovation variance in the sweep recurrence is $1-\rho^4=1-(\tfrac{9801}{10000})^2 =1-\tfrac{96059601}{100000000}=\tfrac{3940399}{100000000}$. The expected squared change in $Y$ per sweep is $2(1-\rho^2)=\tfrac{398}{10000}=\tfrac{199}{5000}$.

A joint block draw of both coordinates has expected squared change in $Y$ equal to $2$, since its new $Y$ is an independent standard normal draw when the current state has the target distribution. The ratio of these expected squared changes is $2/(\tfrac{199}{5000})=\tfrac{10000}{199}$. This large ratio illustrates the small one-coordinate movement in this example; it should not be interpreted as a general computational speed-up, because a joint conditional draw may be harder to implement in a larger model.

Check: The sweep innovation variance equals $1-(\tfrac{9801}{10000})^2=\tfrac{3940399}{100000000}$, and $1-(\tfrac{9801}{10000})^2$ is also the denominator in the stationary-variance calculation. Their ratio is $1$, as required for the standardized target variance.

Choosing a Block and Avoiding a Common Error

In a model with more than two variables, a block need not include the entire state. One can jointly update a strongly dependent pair while conditioning on the remaining coordinates, then update other blocks or coordinates in turn. The relevant conditional is the joint distribution of the selected block given everything held fixed. Its usefulness depends on whether that conditional can be sampled accurately and at an acceptable computational cost. Blocking is most directly motivated when the variables move together under the target and separate conditional updates tend to undo one another.

Common mistake: using the one-coordinate correlation as the full-sweep coefficient. With $\rho=\tfrac35$, it is wrong to say that a full $X$-then-$Y$ sweep gives coefficient $\tfrac35$ on the old $Y$. The first update gives $X_{n+1}=\tfrac35Y_n+\tfrac45\varepsilon_{n+1}$; the second multiplies this conditional mean by another $\tfrac35$. Hence the coefficient on $Y_n$ is $\tfrac35\cdot\tfrac35=\tfrac9{25}$, not $\tfrac35$. A related error is to use conditional variance $1-\rho$; the correct standardized conditional variance is $1-\rho^2=\tfrac{16}{25}$, not $1-\tfrac35=\tfrac25$.

The calculations also distinguish two claims that should not be conflated. Exact Gibbs updates preserve the target by conditional resampling. The movement calculation explains why, for this Gaussian and this systematic order, high absolute correlation makes the successive $Y$-values strongly dependent. Invariance alone does not imply rapid exploration, and the block update is advantageous only when the joint conditional is available and computationally manageable.

Check Your Understanding

Use the conditional formulas and recurrence derived above. State any assumptions needed for your answers.

For standardized coordinates with correlation $\rho=-\tfrac45$, what are the two conditional variances, and what is the coefficient of $Y_n$ in the full $X$-then-$Y$ sweep recurrence?
For $\rho=\tfrac12$, calculate the expected squared change in $Y$ under a single-component sweep and under a joint block draw, assuming the current state has the target distribution.
Explain why resampling $X$ from its conditional distribution given $Y$ preserves the joint density.
What fails in the density-based conditional formulas when $|\rho|=1$?
In a model with additional coordinates, what distribution must be sampled to update a selected pair as a block?