Coupling-based convergence assessment of some Gibbs samplers for high-dimensional Bayesian regression with shrinkage priors

We consider the setting of high-dimensional regression where the number of observations $n$ is smaller than the number of covariates $p$ and the true signal is sparse. This problem formulation is ubiquitous in modern applications ranging from genomics to the social sciences. Optimization-based methods, such as the LASSO (Tibshirani, 1994) or Elastic Net (Zou and Hastie, 2005), allow sparse point estimates to be obtained even when the number of covariates is on the order of hundreds of thousands. More specifically, iterative optimization procedures to obtain these estimates are practical because the following conditions are met: 1) the cost per iteration scales favorably with the size of the input ($n$ and $p$), 2) the number of iterations to convergence also scales favorably, and 3) there are reliable stopping criteria to detect convergence.

In the Bayesian paradigm, Markov chain Monte Carlo (MCMC) methods are commonly used to sample from the posterior distribution. Default MCMC implementations, for example using general-purpose software such as Stan (Carpenter et al., 2017), can lead to high computational costs in high-dimensional settings, both per iteration and in terms of the number of iterations required to reach convergence. In the setting of high-dimensional regressions, tailored algorithms can provide substantial improvements on both fronts (e.g. Bhattacharya et al., 2016; Nishimura and Suchard, 2018; Johndrow et al., 2020; Bhattacharya and Johndrow, 2021). Comparatively less attention has been put on the design of reliable stopping criteria. Stopping criteria for MCMC, such as univariate effective sample size (ESS, e.g Vats et al., 2019), or the $\hat{R}$ convergence diagnostic (e.g. Vats and Knudson, 2020), rely on the asymptotic behavior of the chains as time goes to infinity, and thus effectively ignore the non-asymptotic “burn-in” bias, which vanishes as time progresses. This is acceptable in situations where we have solid a priori estimates of the burn-in period; otherwise the lack of non-asymptotic stopping criteria poses an important practical problem. With Bayesian analysis of high-dimensional data sets, MCMC algorithms can require seconds or minutes per iteration, and thus any non-asymptotic insight on the number of iterations to perform is highly valuable.

Diagnostics of convergence are particularly useful when considering the number of factors that affect the performance of MCMC algorithms for Bayesian regressions. Beyond the size of the data set, the specification of the prior and the (much less explicit) signal-to-noise ratio in the data all have an impact on the convergence of MCMC algorithms. Across applications, the number of iterations required for the chain to converge varies by multiple orders of magnitude. This could lead users to either waste computational resources by running overly long chains, or, worrisomely, to base their analysis on chains that have not converged. This manuscript proposes a concrete method to avoid these pitfalls in the case of a Gibbs sampler for Bayesian regression with heavy-tailed shrinkage priors. Specifically, we follow the approach of Glynn and Rhee (2014); Jacob et al. (2020); Biswas et al. (2019) and use coupled lagged chains to monitor convergence.

Consider Gaussian linear regression with $n$ observations and $p$ covariates, with a focus on the high-dimensional setting where $n\ll p$. The likelihood is given by

where $\|\cdot\|$ denotes the $L_{2}$ norm, $X\in\mathbb{R}^{n\times p}$ is the observed design matrix, $y\in\mathbb{R}^{n}$ is the observed response vector, $\sigma^{2}>0$ is the unknown Gaussian noise variance, and $\beta\in\mathbb{R}^{p}$ is the unknown signal vector that is assumed to be sparse. We study hierarchical Gaussian scale-mixture priors on $(\beta,\sigma^{2})$ given by

where $a_{0},b_{0}>0$, and $\pi_{\xi}(\cdot),\pi_{\eta}(\cdot)$ are continuous densities on $\mathbb{R}_{>0}$. Such global-local mixture priors induce approximate sparsity, where the components of $\beta$ can be arbitrarily close to zero, but never exactly zero. This is in contrast to point-mass mixture priors (e.g Johnson and Rossell, 2012), where some components of $\beta$ can be exactly zero a posteriori. The global precision parameter $\xi$ relates to the number of signals, and we use $\xi^{-1/2}\sim\text{Cauchy}_{+}(0,1)$ throughout. The local precision parameters $\eta_{j}$ determine which components of $\beta$ are null.

We focus on the Half-$t(\nu)$ prior family for the local scale parameter, $\eta_{j}^{-1/2}\sim t_{+}(\nu)$ for $\nu\geq 1$, where a $t_{+}(\nu)$ distribution has density proportional to $(1+x^{2}/\nu)^{-(\nu+1)/2}\,\mathbbm{1}_{(0,\infty)}(x)$. The induced prior on the local precision $\eta_{j}$ has density

The case $\nu=1$ corresponds to the Horseshoe prior (Carvalho et al., 2009, 2010) and has received overwhelming attention in the literature among this prior class; see Bhadra et al. (2019) for a recent overview. The use of Half-$t(\nu)$ priors on scale parameters in hierarchical Gaussian models was popularized by Gelman (2006). However, the subsequent literature has largely gravitated towards the default choice (Polson and Scott, 2012) of $\nu=1$, in part because a convincing argument for preferring a different value of $\nu$ has been absent to date. Here we give empirical evidence that Half-$t(\nu)$ priors yield statistical estimation properties similar to the Horseshoe prior, whilst leading to significant computational gains when $\nu>1$. We note that Ghosh and Chakrabarti (2017) established optimal posterior concentration in the Normal means problem for a broad class of priors on the local scale that includes the Half-$t(\nu)$ priors, extending earlier work by van der Pas et al. (2014) for the Horseshoe, providing theoretical support for the comparable performance we observe with $\nu=1$ and $\nu>1$ in simulations.

Our stopping criterion relies on couplings. Couplings have long been used to analyze the convergence of MCMC algorithms, while methodological implementations have been rare. The work of Johnson (1996) pioneered the use of couplings for practical diagnostics of convergence, but relied on the assumption of some similarity between the initial distribution and the target, which can be hard to check. Here we instead follow the approach of Glynn and Rhee (2014); Jacob et al. (2020) based on couplings of lagged chains, applied to the question of convergence in Biswas et al. (2019). That approach makes no assumptions on the closeness of the initial distribution to the target, and in our experiments we initialize all chains independently from the prior distribution.

The approach requires the ability to sample pairs of Markov chains such that 1) each chain marginally follows a prescribed MCMC algorithm and 2) the two chains “meet” after a random –but almost surely finite– number of iterations, called the meeting time. When the chains meet, i.e. when all of their components become identical, the user can stop the procedure. From the output the user can obtain unbiased estimates of expectations with respect to the target distribution, which can be generated in parallel and averaged, as well as estimates of the finite-time bias of the underlying MCMC algorithm. Crucially, this information is retrieved from the distribution of the meeting time, which is an integer-valued random variable irrespective of the dimension of the problem, and thus provides a convenient summary of the performance of the algorithm. A difficulty inherent to the approach is the design of an effective coupling strategy for the algorithm under consideration. High-dimensional parameter spaces add a substantial layer of complication as some of the simpler coupling strategies, that are solely based on maximal couplings, lead to prohibitively long meeting times.

In this paper we argue that couplings offer a practical way of assessing the number of iterations required by MCMC algorithms in high-dimensional regression settings. Our specific contributions are summarized below.

We introduce blocked Gibbs samplers (Section 2) for Bayesian linear regression focusing on Half-$t(\nu)$ local shrinkage priors, extending the algorithm of Johndrow et al. (2020) for the Horseshoe. Our algorithm has an overall computation cost of $\mathcal{O}(n^{2}p)$ per iteration, which is state-of-the-art for $n\ll p$. We then design a one-scale coupling strategy for these Gibbs samplers (Section 2.3) and show that it results in meeting times that have exponential tails, which in turn validates their use in convergence diagnostics (Biswas et al., 2019). Our proof of exponentially-tailed meeting times is based on identifying a geometric drift condition, as a corollary to which we also establish geometric ergodicity of the marginal Gibbs sampler. Despite a decade of widespread use, MCMC algorithms for the Horseshoe and Half-$t(\nu)$ priors generally have not previously been shown to be geometrically ergodic, in contrast to MCMC algorithms for global-local models with exponentially-tailed local scale priors (Khare and Hobert, 2013; Pal and Khare, 2014). Our proofs utilize a uniform bound on generalized ridge regression estimates, which enables us to avoid a technical assumption of Johndrow et al. (2020), and auxiliary results on moments of distributions related to the confluent hypergeometric function of the second kind. Shortly after an earlier version of this manuscript was posted, Bhattacharya et al. (2021) have also established geometric ergodicity of related Gibbs samplers for the Horseshoe prior.

The meeting times resulting from the one-scale coupling increase exponentially as dimension increases. To address that issue we design improved coupling strategies (Section 3). We develop a two-scale strategy, based on a carefully chosen metric, which informs the choice between synchronous and maximal couplings. In numerical experiments in high dimensions this strategy leads to orders of magnitude shorter meeting times compared to the naïve approach. We establish some preliminary results in a stylized setting, providing a partial understanding of why the two-scale strategy leads to the observed drastic reduction in meeting times. We also describe other couplings that do not require explicit calculation of a metric between the chains, have exponentially-tailed meeting times, and empirically result in similar meeting times as the two-scale coupling. For our two-scale coupling, we further note that priors with stronger shrinkage towards zero can give significantly shorter meeting times in high dimensions. This motivates usage of Half-$t(\nu)$ priors with greater degrees of freedom $\nu>1$, which can give similar statistical performance to the Horseshoe ($\nu=1$) whilst enjoying orders of magnitude computational improvements.

Finally, we demonstrate (Section 4) that the proposed coupled MCMC algorithm is applicable to big data settings. Figure 1 displays an application of our coupling strategy to monitor the convergence of the Gibbs sampler for Bayesian regression on a genome-wide association study (GWAS) dataset with $(n,p)\approx(2000,100000)$. Our method suggests that a burn-in of just $700$ iterations can suffice even for such a large problem, and confirms the applicability of coupled MCMC algorithms for practical high-dimensional inference.

Scripts in R (R Core Team, 2013) are available from https://github.com/niloyb/CoupledHalfT to reproduce the figures of the article.

Blocked Gibbs samplers are popularly used for Bayesian regression with global-local shrinkage priors, due to the analytical tractability of full conditional distributions (Park and Casella, 2008; Carvalho et al., 2010; Polson et al., 2014; Bhattacharya et al., 2015; Makalic and Schmidt, 2016). Several convenient blocking and marginalization strategies are possible, leading to conditionals that are easy to sample from. For the case of the Horseshoe prior ($\nu=1$), Johndrow et al. (2020) have developed exact and approximate MCMC algorithms for high-dimensional settings, building on the algorithms of Polson et al. (2014) and Bhattacharya et al. (2016). We extend that sampler to general Half-$t(\nu)$ priors, and summarize it in (LABEL:eq:blocked_gibbs_ergodicity_version). Given some initial state $(\beta_{0},\eta_{0},\sigma_{0}^{2},\xi_{0})\in\mathbb{R}^{p}\times\mathbb{R}^{p}_{>0}\times\mathbb{R}_{>0}\times\mathbb{R}_{>0}$, it generates a Markov chain $(\beta_{t},\eta_{t},\sigma^{2}_{t},\xi_{t})_{t\geq 0}$ which targets the posterior corresponding to Half-$t(\nu)$ priors.

Above $L(y|\xi_{t+1},\eta_{t+1})$ represents the marginal likelihood of $y$ given $\xi_{t+1}$ and $\eta_{t+1}$, and we use the notation $M=I_{n}+\xi_{t+1}^{-1}X\,\text{Diag}(\eta_{t+1}^{-1})\,X^{T}$, and $\Sigma=X^{T}X+\xi_{t+1}\text{Diag}(\eta_{t+1})$. We sample $\eta_{t+1}|\beta_{t},\sigma_{t}^{2},\xi_{t}$ component-wise independently using slice sampling and sample $\xi_{t+1}|\eta_{t+1}$ using Metropolis–Rosenbluth–Teller–Hastings with step-size $\sigma_{\text{MRTH}}$ (Hastings, 1970). Full details of our sampler are in Algorithms 3 and 4 of Appendix D, and derivations are in Appendix E.

Sampling $\eta_{t+1}|\beta_{t},\sigma_{t}^{2},\xi_{t}$ component-wise independently is an $\mathcal{O}(p)$ cost operation. Sampling $\xi_{t+1}|\eta_{t+1}$ requires evaluating $L(y|\xi_{t+1},\eta_{t+1})$, which involves an $\mathcal{O}(n^{2}p)$ cost operation (calculating the weighted matrix cross-product $X\text{Diag}(\eta_{t+1})^{-1}X^{T}$) and an $\mathcal{O}(n^{3})$ cost operation (calculating $y^{T}M^{-1}y$ and $|M|$ using Cholesky decomposition). Sampling $\sigma_{t+1}^{2}|\eta_{t+1},\xi_{t+1}$ is only an $\mathcal{O}(1)$ cost operation, as $y^{T}M^{-1}y$ is pre-calculated. Similarly, as $M^{-1}$ is pre-calculated, sampling $\beta_{t+1}|\eta_{t+1},\xi_{t+1},\sigma_{t+1}^{2}$ using the algorithm of Bhattacharya et al. (2016) for structured multivariate Gaussians only involves an $\mathcal{O}(np)$ and an $\mathcal{O}(n^{2})$ cost operation. In high-dimensional settings with $p\gg n$, which is our focus, the weighted matrix cross-product calculation is the most expensive and overall the sampler described in (LABEL:eq:blocked_gibbs_ergodicity_version) has $\mathcal{O}(n^{2}p)$ computation cost.

We develop couplings for the blocked Gibbs sampler presented in Section 2.1. We will use these couplings to generate coupled chains with a time lag $L\geq 1$, in order to implement the diagnostics of convergence proposed in Biswas et al. (2019). Consider an $L$-lag coupled chain $(C_{t},\tilde{C}_{t-L})_{t\geq L}$ with meeting time $\tau:=\inf\{t\geq L:C_{t}=\tilde{C}_{t-L}\}$. Here $C_{t}$ and $\tilde{C}_{t}$ denote the full states $(\beta_{t},\eta_{t},\sigma_{t}^{2},\xi_{t})$ and $(\tilde{\beta}_{t},\tilde{\eta}_{t},\tilde{\sigma}_{t}^{2},\tilde{\xi}_{t})$ respectively. For the coupling-based methods proposed in Biswas et al. (2019); Jacob et al. (2020) to be most practical, the meeting times should occur as early as possible, under the constraint that both marginal chains $(C_{t})_{t\geq 0}$ and $(\tilde{C}_{t})_{t\geq 0}$ start from the same initial distribution, and that they both marginally evolve according to the blocked Gibbs sampler in (LABEL:eq:blocked_gibbs_ergodicity_version).

Briefly, we recall how $L$-lag coupled chains enable the estimation of upper bounds on $\text{TV}(\pi_{t},\pi)$. Since $C_{t}=\tilde{C}_{t-L}$ for all $t\geq\tau$, then the indicator $\mathbbm{1}(C_{t}\neq\tilde{C}_{t-L})$ can be evaluated for all $t$ as soon as we observe $\tau$. Such indicators are estimators of $\mathbb{P}(C_{t}\neq\tilde{C}_{t-L})$, which themselves are upper bounds on $\text{TV}(\pi_{t},\pi_{t-L})$, since $C_{t}\sim\pi_{t}$ and $\tilde{C}_{t-L}\sim\pi_{t-L}$. We can then obtain upper bounds on $\text{TV}(\pi_{t},\pi)$ by the triangle inequality: $\text{TV}(\pi_{t},\pi)\leq\mathbb{E}[\sum_{j\geq 0}\mathbbm{1}(C_{t+(j+1)L}\neq\tilde{C}_{t+jL})]$. This can be justified more formally if the meeting times are exponentially-tailed (Biswas et al., 2019). Such tail condition also ensures that the expected computation time of unbiased estimators as in Jacob et al. (2020), generated independently in parallel, scales at a logarithmic rate in the number of processors. Thus we will be particularly interested in couplings that result in exponentially-tailed meeting times.

In Section 2.3, we consider an algorithm based on maximal couplings. We show that the ensuing meeting times have exponential tails, and discuss how our analysis directly implies that the marginal chain introduced in Section 2.1 is geometrically ergodic. We then illustrate difficulties encountered by such scheme in high dimensions, which will motivate alternate strategies described in Section 3. For simplicity, we mostly omit the lag $L$ from the notation.

For the blocked Gibbs sampler in Section 2.1, we first consider a coupled MCMC algorithm that attempts exact meetings at every step. We will apply a maximal coupling algorithm with independent residuals (see Thorisson (2000, Chapter 1 Section 4.4) and Johnson (1998)). It is included in Algorithm 5 of Appendix D, and has an expected computation cost of two units (Johnson, 1998).

Our initial coupled MCMC kernel is given in Algorithm 1, which we refer to as a one-scale coupling. That is, before the chains have met (when $C_{t}\neq\tilde{C}_{t}$), the coupled kernel on Steps $(1)$ and $(2)(a-c)$ does not explicitly depend on the distance between states $C_{t}$ and $\tilde{C}_{t}$. After meeting, the coupled chains remain together by construction, such that $C_{t}=\tilde{C}_{t}$ implies $C_{t+1}=\tilde{C}_{t+1}$. When $C_{t}\neq\tilde{C}_{t}$, Step $(1)$ uses the coupled slice sampler of Algorithm 2 component-wise for $(\eta_{t+1},\tilde{\eta}_{t+1})$, which allows each pair of components $(\eta_{t+1,j},\tilde{\eta}_{t+1,j})$ to meet exactly with positive probability. In Algorithm 2, we use a common random numbers or a “synchronous” coupling of the auxiliary random variables $(U_{j,*},\tilde{U}_{j,*})$, and alternative couplings could be considered. Steps $(2)(a-b)$ are maximal couplings of the conditional sampling steps for $(\xi,\tilde{\xi})$ and $(\sigma^{2},\tilde{\sigma}^{2})$, such that $\big{(}\sigma_{t+1},\xi_{t+1},\eta_{t+1}\big{)}=\big{(}\tilde{\sigma}_{t+1},\tilde{\xi}_{t+1},\tilde{\eta}_{t+1}\big{)}$ occurs with positive probability for all $t\geq 0$. Step $(2)(c)$ uses common random numbers, such that $\big{(}\sigma^{2}_{t+1},\xi_{t+1},\eta_{t+1}\big{)}=\big{(}\tilde{\sigma}^{2}_{t+1},\tilde{\xi}_{t+1},\tilde{\eta}_{t+1}\big{)}$ implies $\beta_{t+1}=\tilde{\beta}_{t+1}$. This allows the full chains to meet exactly with positive probability at every step.

Under the one-scale coupling, we prove that the meetings times have exponential tails and hence finite expectation. Our analysis is based on identifying a suitable drift function for a variant of the marginal chain in Section 2.1 and an application of Jacob et al. (2020, Proposition 4). We assume that the global shrinkage prior $\pi_{\xi}(\cdot)$ has a compact support on $(0,\infty)$. Such compactly supported priors for $\xi$ have been recommended by van der Pas et al. (2017) to achieve optimal posterior concentration for the Horseshoe in the Normal means model. For simplicity, we also assume that each component $\eta_{t+1,j}|\beta_{t,j},\xi_{t},\sigma_{t}^{2}$ and $\xi_{t+1}|\eta_{t+1}$ in (LABEL:eq:blocked_gibbs_ergodicity_version) are sampled perfectly, instead of slice sampling and Metropolis–Rosenbluth–Teller–Hastings steps respectively. Perfect sampling algorithms for each component $\eta_{t+1,j}|\beta_{t,j},\xi_{t},\sigma_{t}^{2}$ and $\xi_{t+1}|\eta_{t+1}$ are provided in Algorithms 10 and 11 of Appendix D for completeness; see also Bhattacharya and Johndrow (2021) for perfect sampling of the components of $\eta_{t+1,j}|\beta_{t,j},\xi_{t},\sigma_{t}^{2}$ for the Horseshoe. Perfectly sampling $\xi_{t+1}|\eta_{t+1}$ and $\eta_{t+1,j}|\beta_{t,j},\xi_{t},\sigma_{t}^{2}$ retains the $\mathcal{O}(n^{2}p)$ computational cost for the full blocked Gibbs sampler, though in practice this implementation is more expensive per iteration due to the computation of eigenvalue decompositions in lieu of Cholesky decompositions, and the computation of inverse of the confluent hypergeometric function of the second kind.

Proposition 2.1 gives a geometric drift condition for this variant of the blocked Gibbs sampler.

The drift (or “Lyapunov”) function in (5), approaches infinity when any $\beta_{j}$ approaches the origin or infinity. This ensures that the corresponding sub-level sets, defined by $S(R)=\{(\beta,\xi,\sigma^{2})\in\mathbb{R}^{p}\times\mathbb{R}_{>0}\times\mathbb{R}_{>0}:V(\beta,\xi,\sigma)<R\}$ for $R>0$, exclude the pole at the origin and are bounded. Such geometric drift condition, in conjunction with Jacob et al. (2020, Proposition 4), helps verify that the meeting times under the one-scale coupling have exponential tails and hence finite expectation.

Exponentially-tailed meeting times immediately imply that the blocked Gibbs sampler in Section 2.1 is geometrically ergodic. This follows from noting that Jacob et al. (2020, Proposition 4) is closely linked to a minorization condition (Rosenthal, 1995), and details are included in Proposition C.10 of Appendix C. For a class of Gibbs samplers targeting shrinkage priors including the Bayesian LASSO, the Normal-Gamma prior (Griffin and Brown, 2010), and Dirichlet-Laplace prior (Bhattacharya et al., 2015), Khare and Hobert (2013) and Pal and Khare (2014) have proven geometric ergodicity based on drift and minorization arguments (Meyn and Tweedie, 1993; Roberts and Rosenthal, 2004). For the Horseshoe prior, Johndrow et al. (2020) has recently established geometric ergodicity. In the high-dimensional setting with $p>n$ for the Horseshoe, the proof of Johndrow et al. (2020) required truncating the prior on each $\eta_{j}$ below by a small $\ell>0$ to guarantee the uniform bound $\eta_{j}\geq\ell$. Our proof of geometric ergodicity for Half-$t(\nu)$ priors including the Horseshoe ($\nu=1$) works in both low and high-dimensional settings, without any modification of the priors on the $\eta_{j}$. In parallel work, Bhattacharya et al. (2021) have also established geometric ergodicity for the Horseshoe prior without requiring such truncation.

We now consider the performance of Algorithm 1 on synthetic datasets. In particular, we empirically illustrate that the rate $\kappa_{0}$ in Proposition 2.2 can tend to $1$ exponentially fast as dimension $p$ increases. The design matrix $X\in\mathbb{R}^{n\times p}$ is generated with $[X]_{i,j}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,1)$, and a response vector as $y\sim\mathcal{N}(X\beta_{*},\sigma^{2}_{*}I_{n})$, where $\beta_{*}\in\mathbb{R}^{p}$ is the true signal and $\sigma_{*}\geq 0$ is the true error standard deviation. We choose $\beta_{*}$ to be sparse such that given sparsity parameter $s$, $\beta_{*,j}=2^{(9-j)/4}$ for $1\leq j\leq s$ and $\beta_{*,j}=0$ for all $s>j$. Figure 2 shows the meeting times $\tau$ of coupled Markov chains (with a lag $L=1$) targeting the Horseshoe posterior ($\nu=1$) with $a_{0}=b_{0}=1$ under Algorithm 1. We consider $n=100,s=10,\sigma_{*}=0.5$ and vary the dimension $p$. For each dimension $p$, we simulate $100$ meeting times under the one-scale coupling for independently generated synthetic datasets. We initialize chains independently from the prior distribution. For the $\xi$ updates, we use a proposal step-size of $\sigma_{\text{MRTH}}=0.8$ as in Johndrow et al. (2020). Figure 2, with a vertical axis on a logarithmic scale, shows that the meeting times increase exponentially and have heavier tails with increasing dimension.

We develop a two-scale coupling algorithm, which slightly differs from the terminology used in stochastic dynamics and in analysis of pre-conditioned HMC (Bou-Rabee and Eberle, 2020). It is also closely linked to the delayed one-shot coupling of Roberts and Rosenthal (2002), and refers to coupled kernels that attempt exact meetings when the chains are close, and aim for a contraction when the chains are far. The motivation for this construction is that in Algorithm 1, the components of $(\eta_{t},\tilde{\eta}_{t})$ that fail to meet instead evolve independently. As dimension grows, the number of components of $(\eta_{t},\tilde{\eta}_{t})$ that fail to meet also tends to grow, making it increasingly difficult for all components to exactly meet on subsequent iterations. In the two-scale coupling, we only attempt to obtain exact meetings when the associated probability is large enough. This is done by constructing a metric $d$ and a corresponding threshold parameter $d_{0}\geq 0$ such that when the current states are $d$-close with $d(C_{t},\tilde{C}_{t})\leq d_{0}$, we apply Algorithm 2 and try to obtain exact meetings. When $d(C_{t},\tilde{C}_{t})>d_{0}$, we instead employ common random numbers to sample $(\eta_{t},\tilde{\eta}_{t})$.

Our chosen metric $d$ on $\mathbb{R}^{p}\times\mathbb{R}^{p}_{>0}\times\mathbb{R}_{>0}\times\mathbb{R}_{>0}$ is

where recall that $C_{t}=(\beta_{t},\eta_{t},\sigma^{2}_{t},\xi_{t})$, $\tilde{C}_{t}=(\tilde{\beta}_{t},\tilde{\eta}_{t},\tilde{\sigma}^{2}_{t},\tilde{\xi}_{t})$, and $(\eta_{t+1},\tilde{\eta}_{t+1})|C_{t},\tilde{C}_{t}$ is sampled using the coupled slice sampler in Algorithm 2. As the coupling in Algorithm 2 is independent component-wise, we can obtain a simplified expression involving only univariate probabilities as

where $m_{t,j}={\xi_{t}\beta_{t,j}^{2}}/(2\sigma_{t}^{2})$ and $m_{t,j}={\tilde{\xi}_{t}\tilde{\beta}_{t,j}^{2}}/(2\tilde{\sigma}_{t}^{2})$. Under this metric $d$ and a threshold $d_{0}\in[0,1]$, an exact meeting is attempted when $\mathbb{P}(\eta_{t+1}=\tilde{\eta}_{t+1}|C_{t},\tilde{C}_{t})\geq 1-d_{0}$. Once $\eta_{t+1}$ and $\tilde{\eta}_{t+1}$ coincide, and if the scalars $\xi_{t+1}$ and $\tilde{\xi}_{t+1}$ also subsequently coincide, then the entire chains meet. When $d_{0}=1$, the two-scale coupling reduces to the one-scale coupling of Algorithm 1 since the maximal coupling slice sampler (Algorithm 2) is invoked at each step. On the other hand, when threshold $d_{0}=0$, we always apply the common random numbers. Then the chains $C_{t},\tilde{C}_{t}$ may come arbitrarily close but will never exactly meet. Empirically we find that different thresholds $d_{0}$ values away from $0$ and $1$ give similar meeting times. Simulations concerning these choices can be found in Appendix B.2.

We now discuss computation of the metric $d(C_{t},\tilde{C}_{t})$. Since the probability in (7) is unavailable in closed form, we can resort to various approximations. One option is to evaluate the one-dimensional integrals in (8) using numerical integration. Alternatively, one can form the Monte Carlo based estimate

where $(\eta^{(r)}_{t+1},\tilde{\eta}^{(r)}_{t+1})$ are sampled independently for $r=1,\ldots,R$ using Algorithm 2. We recommend a Rao–Blackwellized estimate by combining Monte Carlo with analytical calculations of integrals as

where $\big{(}U^{(r)}_{j},\tilde{U}^{(r)}_{j}\big{)}$ are sampled independently for $r=1,\ldots,R$ using Algorithm 2 and each $\mathbb{P}(\eta_{t+1,j}=\tilde{\eta}_{t+1,j}\Big{|}U^{(r)}_{j},\tilde{U}^{(r)}_{j},m_{t,j},\tilde{m}_{t,j})$ is calculated analytically based on meeting probabilities from Algorithm 2. Further metric calculation details are in Appendix B.3. For a number of samples $R$, estimates (9) and (10) both have computation cost of order $pR$. Compared to the estimate in (9), the estimate in (10) has lower variance and faster numerical run-times as it only involves sampling uniformly distributed random numbers. It suffices to choose a small number of samples $R$, and often we take $R=1$. Indeed it appears unnecessary to estimate $d(C_{t},\tilde{C}_{t})$ accurately, as we are only interested in comparing that distance to a fixed threshold $d_{0}\in[0,1]$. Often the trajectories $(d(C_{t},\tilde{C}_{t}))_{t\geq 0}$ initially take values close to $1$, and then sharply drop to values close to $0$. This leads to the estimates in (9) and (10) having low variance. Henceforth, we will use the estimate in (10) in our experiments with two-scale couplings, unless specified otherwise.

Overall the two-scale coupling kernel has twice the $\mathcal{O}(n^{2}p)$ cost of the single chain kernel in (LABEL:eq:blocked_gibbs_ergodicity_version). When $C_{t}\neq\tilde{C}_{t}$, we calculate a distance estimate and sample from a coupled slice sampling kernel. Calculating the distance estimate involves sampling $2pR$ uniforms and has $\mathcal{O}(p)$ cost. Coupled slice sampling with maximal couplings (Algorithm 2) or common random numbers both have expected or deterministic cost $\mathcal{O}(p)$. The remaining steps of the two-scale coupling match the one-step coupling in Algorithm 1.

We now consider the performance of the two-scale coupling on synthetic datasets. The setup is identical to that introduced in Section 2.3, where for each dimension $p$ we simulate $100$ meetings times based on independently generated synthetic datasets. We target the Horseshoe posterior ($\nu=1$) with $a_{0}=b_{0}=1$. Figure 3 shows the meeting times $\tau$ of coupled Markov chains (with a lag $L=1$), for both the one-scale coupling and the two-scale coupling with $R=1$ and threshold $d_{0}=0.5$. Figure 3 shows that the two-scale coupling can lead to orders of magnitude reductions in meeting times, compared to the one-scale coupling.

We next present some preliminary results which hint as to why the two-scale coupling leads to the drastic reduction in meeting times. We focus on the coupling of the $(\eta,\tilde{\eta})$ full conditionals, as the one-scale and two-scale algorithms differ exactly at this step. First, we show that for any monotone function $h:(0,\infty)\rightarrow\mathbb{R}$ and all components $j=1,\ldots,p$, the expected distance $\mathbb{E}[|h(\eta_{t+1,j})-h(\tilde{\eta}_{t+1,j})|]$ is the same under both common random numbers and maximal coupling.

Proposition 3.1 implies that such single-step expectations alone do not distinguish the behavior of CRN from maximal couplings. This compels us to investigate other distributional features of $|h(\eta_{t+1,j})-h(\tilde{\eta}_{t+1,j})|$ under either coupling to possibly tease apart differences in their behavior on high probability events. Focusing on the Horseshoe prior and making the choice $h(x)=\log(1+x)$, we uncover a distinction between the tail behavior of the two couplings which can substantially accumulate over $p$ independent coordinates. We offer an illustration in a stylized setting below.

Extensions of Proposition 3.2 which allow different components of $m_{t}$ and $\tilde{m}_{t}$ to take different values under some limiting assumptions are omitted here for simplicity. This stylized setting already captures why the one-scale algorithm may result in much larger meeting times in high dimensions compared to the CRN-based two-scale algorithm. We note that Proposition 3.2 remains informative whenever $|\log\mu-\log\tilde{\mu}|\ll\log\log p$, even when $\log\log p$ is not large. For example, consider $\mu=\tilde{\mu}e^{\delta}$ for some small $\delta>0$ with $\delta\ll\log\log p$. Under CRN, $\log(1+\eta_{t+1,j})$ and $\log(1+\tilde{\eta}_{t+1,j})$ remain $\delta$-close almost surely for all components $j$. In contrast, under the maximal coupling, at least one component will have larger than $\log(\alpha L)+\log\log p$ deviations with high probability. Since $\alpha$ and $L=\log(1+1/\mu)/2$ do not depend on $p$ or $\delta$, $\log(\alpha L)+\log\log p\gg\delta$ for fixed $\mu,\tilde{\mu}>0$. These larger deviations between some components of $\eta_{t+1}$ and $\tilde{\eta}_{t+1}$ under maximal couplings push the two chains much further apart when sampling $(m_{t+1},\tilde{m}_{t+1})\big{|}\eta_{t+1},\tilde{\eta}_{t+1}$, resulting in meeting probabilities that are much lower over multiple steps of the algorithm. Overall, Propositions 3.1 and 3.2 show that multi-step kernel calculations taking into account such high-probability events may be necessary to further distinguish the relative performances of one-scale and two-scale couplings. A full analytic understanding of this difference in performance remains an open problem.

Another way of understanding the benefits of two-scale coupling would be to find a distance under which the coupled chains under common random numbers exhibit a contraction. The CRN coupling part of Proposition 3.2 hints that $\bar{d}(m_{t},\tilde{m}_{t}):=\sum_{j=1}^{p}|\log m_{t,j}-\log\tilde{m}_{t,j}|$ may be a natural metric. Appendix B.1 contains discussions and related simulations comparing this metric with alternatives. Proposition 3.3 verifies that such metric bounds the total variation distance between one-step transition kernels.