We consider Markov chain Monte Carlo (MCMC) algorithms for Bayesian high-dimensional regression with
continuous shrinkage priors. A common challenge with these
algorithms is the choice of the number of iterations to perform. This is
critical when each iteration is expensive, as is the
case when dealing with modern data sets, such as genome-wide association
studies with thousands of rows and up to hundred of thousands of columns.  We
develop coupling techniques tailored to the setting of high-dimensional
regression with shrinkage priors, which enable practical, non-asymptotic
diagnostics of convergence without relying on traceplots or long-run
asymptotics. By establishing geometric drift and minorization
conditions for the algorithm under consideration, we prove that
the proposed couplings have finite expected meeting time.
Focusing on a class of shrinkage priors which includes the
“Horseshoe”, we empirically demonstrate the scalability of the proposed
couplings. A highlight of our findings is that less
than 1000 iterations can be enough for a Gibbs sampler to reach stationarity
in a regression on $100,000$ covariates.  The numerical results also
illustrate the
impact of the prior on the computational efficiency of the coupling,
and suggest the use of priors where the local precisions are
Half-t distributed with degree of freedom larger than one.

§ INTRODUCTION

§.§ Iterative computation in high dimensions

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
<cit.> or

Elastic Net <cit.>, 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 <cit.>, 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
Comparatively less attention has been put on the design of reliable stopping criteria.
Stopping criteria for MCMC, such as
univariate effective sample size  <cit.>, or the $\hat{R}$

convergence diagnostic <cit.>,
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 <cit.> and use
coupled lagged chains to monitor convergence.

§.§ Bayesian shrinkage regression with Half-$t(\nu)$ priors

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

\begin{equation}
\label{eq:linearregression}
L(\beta, \sigma^2 ; X, y) = \frac{1}{(2 \pi \sigma^2)^{n/2}} \exp \Big( - \frac{1}{2 \sigma^2} \| y - X \beta \|^2 \Big),
\end{equation}

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

β_j |σ^2, ξ, ηj=1,…,pind∼ 𝒩 (0, σ^2/ξη_j ),  ξ∼π_ξ(·),  η_j j=1,…,pi.i.d.∼ π_η(·),  σ^2 ∼(a_0/2, b_0/2 ),

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
<cit.>, 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 \ge 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

\begin{equation}\label{eq:half_t_eta}
\pi_{\eta}(\eta_j) \propto \frac{1}{\eta_j^{\frac{2-\nu}{2}}(1+\nu \eta_j)^{\frac{\nu + 1}{2}}} \, \mathbbm{1}_{(0, \infty)}(\eta_j).
\end{equation}

The case $\nu = 1$ corresponds to the Horseshoe prior
<cit.> and has received
overwhelming attention in the literature among this prior class; see
<cit.> for a recent overview. The use of Half-$t(\nu)$
priors on scale parameters in hierarchical Gaussian models was popularized by
<cit.>. However, the subsequent
literature has largely gravitated towards the default choice
<cit.> 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
<cit.> 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
<cit.> 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
<cit.> 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
<cit.> based on couplings
of lagged chains, applied to the question of convergence in
<cit.>. 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.

§.§ Our contribution

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
<ref>) for Bayesian linear regression
focusing on Half-$t(\nu)$ local shrinkage priors, extending
the algorithm of <cit.> for the Horseshoe.
Our algorithm has an overall computation cost of $\mathcal{O}(n^2p)$ 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
<ref>) and show that it
results in meeting times that have exponential tails, which in turn
validates their use in convergence diagnostics <cit.>.
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
Our proofs utilize a uniform bound on
generalized ridge regression estimates, which enables us to avoid a technical
assumption of <cit.>, 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, <cit.> 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 <ref>).
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 <ref>) that the proposed
coupled MCMC algorithm is applicable to big data settings.
Figure <ref> 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.

Upper bounds on the total variation between the marginal distributions of the chain
and the target distribution for Bayesian regression on a GWAS dataset with
$(n,p) \approx (2000, 100000)$.

Scripts in R <cit.> are available from
<https://github.com/niloyb/CoupledHalfT> to reproduce the figures of the article.

§ COUPLED MCMC FOR REGRESSION WITH HALF-$T(\NU)$ PRIORS

We develop MCMC techniques for Bayesian shrinkage regression with
Half-$t(\nu)$ priors. The model is given by (<ref>)–(<ref>)
and $\eta^{-1/2}\sim t(\nu)$ for some $\nu \geq 1$ as in (<ref>).
Posterior distributions resulting from such heavy-tailed shrinkage priors
have desirable statistical properties, but their features pose challenges to
generic-purpose MCMC algorithms. Characterization of the marginal prior and
posterior densities of $\beta$ on $\mathbb{R}^p$ given $\xi$ and $\sigma^2$, and
further comments on these challenges
are presented in Appendix <ref>.
Specifically, the resulting posterior distributions
present 1) multimodality, 2) heavy tails and 3) infinite density at zero.
This hints at a trade-off between statistical accuracy and computational feasibility, since the very
features that present computational difficulties are crucial for optimal
statistical performance across sparsity levels and signal strengths.

§.§ A blocked Gibbs sampler for Half-$t(\nu)$ priors

Blocked Gibbs samplers are popularly used for Bayesian regression with
global-local shrinkage priors, due to the analytical tractability of full
conditional distributions
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$), <cit.> have
developed exact and approximate MCMC algorithms for
high-dimensional settings,
building on the algorithms of <cit.> and
<cit.>. We extend that sampler
to general Half-$t(\nu)$ priors,
and summarize it in (<ref>).
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.

π(η_t+1 | β_t, σ_t^2, ξ_t)

∝∏_j=1^p e^-m_t,j η_t+1,j/η_t+1,j^1-ν/2 (1+ νη_t+1,j)^ν+1/2   for  m_t,j = ξ_t β_t,j^2/2 σ_t^2,

π(ξ_t+1 | η_t+1) ∝L(y | ξ_t+1, η_t+1 ) π_ξ(ξ_t+1)

σ_t+1^2 | η_t+1, ξ_t+1 ∼( a_0+n/2, y^T M^-1 y + b_0/2 )

β_t+1 | η_t+1, ξ_t+1, σ_t+1^2 ∼𝒩(Σ^-1 X^Ty, σ_t+1^2 Σ^-1)

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^TX + \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}}$  <cit.>.

Full details of our sampler are in Algorithms <ref> and <ref>
of Appendix <ref>, and derivations are in Appendix

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^2p)$ 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 <cit.> 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 (<ref>)
has $\mathcal{O}(n^2p)$ computation cost.

§.§ Coupling the proposed Gibbs sampler

We develop couplings for the blocked Gibbs sampler presented in Section
<ref>. 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 <cit.>. 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 <cit.>
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

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 <cit.>. Such tail condition also
ensures that the expected computation time of unbiased estimators as in
<cit.>, 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 <ref>, 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 <ref> is geometrically
We then illustrate difficulties encountered by such scheme in high
dimensions, which will motivate alternate strategies
described in Section <ref>. For simplicity, we mostly omit the lag $L$ from the notation.

§.§ One-scale coupling

For the blocked Gibbs sampler in Section <ref>, 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
<cit.> and
<cit.>). It is included in Algorithm
<ref> of Appendix <ref>, and has an expected
computation cost of two units <cit.>.

Our initial coupled MCMC kernel is given in Algorithm <ref>,
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
<ref> 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 <ref>, 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.

$ C_t := (\beta_t, \eta_t, \sigma^2_t, \xi_t) $ and $\tilde{C}_t := (\tilde{\beta}_t, \tilde{\eta}_t, \tilde{\sigma}^2_t, \tilde{\xi}_t)$.

$C_t = \tilde{C}_t$
Sample $C_{t+1} | C_{t}$ as in (<ref>) and set $\tilde{C}_{t+1} = C_{t+1}$.

* Sample $(\eta_{t+1}, \tilde{\eta}_{t+1}) \big| C_t, \tilde{C}_t$ component-wise independently using Algorithm <ref>.

* Sample $( (\xi_{t+1}, \sigma_{t+1}^2, \beta_{t+1} ), (\tilde{\xi}_{t+1}, \tilde{\sigma}_{t+1}, \tilde{\beta}_{t+1} ))$ given $\eta_{t+1}, \tilde{\eta}_{t+1}$ as follows:

* Sample $(\xi_{t+1}, \tilde{\xi}_{t+1}) \big| \eta_{t+1}, \tilde{\eta}_{t+1}, \xi_{t}, \tilde{\xi}_{t}$ using coupled Metropolis–Rosenbluth–Teller–Hastings (Algorithm <ref> in Appendix <ref>).

* Sample $( \sigma_{t+1}^2, \tilde{\sigma}_{t+1}^2 ) | \xi_{t+1}, \eta_{t+1}, \tilde{\xi}_{t+1}, \tilde{\eta}_{t+1}$ from a maximal coupling of two Inverse Gamma distributions (Algorithm <ref> in Appendix <ref>).

* Sample $( \beta_{t+1}, \tilde{\beta}_{t+1} ) | \sigma_{t+1}^{2}, \xi_{t+1}, \eta_{t+1}, \tilde{\sigma}_{t+1}^{2}, \tilde{\xi}_{t+1}, \tilde{\eta}_{t+1}$ using common random numbers for Gaussian scale mixture priors (Algorithm <ref> in Appendix <ref>).

$(C_{t+1} , \tilde{C}_{t+1})$.
A one-scale coupled MCMC kernel for Half-$t(\nu)$ priors.

A coupling of slice samplers marginally targeting
$p( \cdot | m_j)$ and $p( \cdot | \tilde{m}_j)$, where
$p(\eta_j | m) \propto (\eta_j^{\frac{1-\nu}{2}}(1+\nu \eta_j)^{\frac{\nu + 1}{2}})^{-1} e^{-m \eta_j}$ on $(0,\infty)$.

$\eta_{t,j}, \tilde{\eta}_{t,j}, m_{t,j} := \xi_{t} \beta_{t,j}^2 / ( 2 \sigma_{t}^2 ), \tilde{m}_{t,j} := \tilde{\xi}_{t} \tilde{\beta}_{t,j}^2 / (2 \tilde{\sigma}_{t}^2) >0 $.

* Sample
$U^{crn}_{j,*} \sim \Uniform(0,1)$, and set $U_{j,*} := U^{crn}_{j,*} (1+\nu \eta_{t,j})^{-\frac{\nu + 1}{2}} $ and
$\tilde{U}_{j,*} := U^{crn}_{j,*} (1+\nu \tilde{\eta}_{t,j})^{-\frac{\nu + 1}{2}} $.

* Sample
$(\eta_{t+1,j}, \tilde{\eta}_{t+1,j})  | U_{j,*}, \tilde{U}_{j,*}$ from a maximal coupling of distributions
$P_j$ and $\tilde{P}_j$ using Algorithm <ref> in Appendix <ref>,

where $P_{j}$ and $\tilde{P}_{j}$ have unnormalized densities $\eta \mapsto \eta^{s-1} e^{ - m_j \eta}$
and $\eta \mapsto \eta^{s-1} e^{ - \tilde{m}_j \eta}$ on $(0, T_{j})$
and $(0, \tilde{T}_{j})$ respectively, where $T_{j}=(U_{j,*}^{-2/(1+\nu)}-1)/\nu$,
$\tilde{T}_{j}=(\tilde{U}_{j,*}^{-2/(1+\nu)}-1)/\nu$ and $s=(1+\nu)/2$.

$(\eta_{t+1,j}, \tilde{\eta}_{t+1,j})$.
Coupled Slice Sampler for Half-$t(\nu)$ priors using maximal couplings.

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 <ref>
and an application of <cit.>.
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 <cit.>
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 (<ref>)
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 <ref>
and <ref> of Appendix
<ref> for completeness; see also <cit.> 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
computational cost for the full blocked Gibbs sampler,

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 <ref> gives a geometric drift condition for this variant of the blocked Gibbs sampler.

(Geometric drift).  Consider the Markov chain $(\beta_t, \xi_t, \sigma^2_t)_{t\geq 0}$ generated from the blocked Gibbs sampler in (<ref>), where
each component $\eta_{t+1,j} | \beta_{t,j}, \xi_t, \sigma_t^2$ and $\xi_{t+1} | \eta_{t+1}$
are sampled perfectly such that $(\beta_t, \xi_t, \sigma^2_t) \mapsto \eta_{t+1} \mapsto  (\beta_{t+1}, \xi_{t+1}, \sigma^2_{t+1})$ with $\eta_{t+1}$ intermediary for each $t\geq0$. Assume the global shrinkage prior $\pi_\xi(\cdot)$ has a compact support.
For $m_j = \xi \beta_j^2/(2 \sigma^2)$, $c\in (0,1/2)$ and $d \in (0,1)$, define

\begin{equation} \label{eq:lyapunov_function}
V(\beta, \xi, \sigma^2) = \sum_{j=1}^p m_j^{-c} + m_j^{d}.
\end{equation}

Then for each $\nu \geq 1$, there exist some
$c\in (0,1/2), d \in (0,1)$ such that $V$ is a drift function, i.e.,
there exists $\gamma \in (0,1)$ and $K \in (0, \infty)$ such that
for all $\beta_{t} \in \mathbb{R}^p$, $\xi_t >0$ and $\sigma_t^2 >0$,

𝔼[V(β_t+1, ξ_t+1, σ_t+1^2) | β_t, ξ_t, σ_t^2] ≤γV(β_t, ξ_t, σ_t^2) + K.

The drift (or “Lyapunov”) function in (<ref>),
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
helps verify that the meeting times under the one-scale coupling
have exponential tails and hence finite expectation.

Consider the blocked Gibbs sampler in (<ref>).
In the setup of Proposition <ref>, assume that the global shrinkage
prior $\pi_\xi(\cdot)$ has a compact support. Write $Z_t = (\beta_t, \xi_t, \sigma^2_t)$ and
$\tilde{Z}_t = ( \tilde{\beta}_t, \tilde{\xi}_t, \tilde{\sigma}^2_t)$.
Consider the one-scale coupling given by
$(Z_t, \tilde{Z}_t) \mapsto (\eta_{t+1}, \tilde{\eta}_{t+1}) \mapsto (Z_{t+1}, \tilde{Z}_{t+1})$
where $(\eta_{t+1}, \tilde{\eta}_{t+1})$ are maximally coupled component-wise, and
$(Z_{t+1}, \tilde{Z}_{t+1})$ are coupled using common random numbers. Denote
the meeting time by $\tau := \inf \{ t \geq 0: Z_t = \tilde{Z}_t \}$. Then

$\mathbb{P}( \tau > t ) \leq A_0 \kappa_0^t$

for some constants $A_0 \in (0, \infty)$ and $\kappa_0 \in (0,1)$, and for all $t \geq 0$.

Exponentially-tailed meeting times immediately imply that the blocked Gibbs sampler in Section
<ref> is geometrically ergodic. This follows
from noting that <cit.> is closely linked to
a minorization condition <cit.>, and details
are included in Proposition <ref> of Appendix <ref>.
For a class of Gibbs samplers targeting shrinkage
priors including the Bayesian LASSO, the Normal-Gamma prior <cit.>, and
Dirichlet-Laplace prior <cit.>, <cit.> and
<cit.> have proven geometric ergodicity based on drift and
minorization arguments <cit.>. For the Horseshoe prior,
<cit.> has recently established geometric ergodicity. In
the high-dimensional setting with $p > n$ for the Horseshoe, the
proof of <cit.> 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,
<cit.> have also established geometric ergodicity
for the Horseshoe prior without requiring such truncation.

Our proofs of exponentially-tailed meeting times and geometric ergodicity
generalize to a larger class of priors satisfying some moment conditions on the full
conditionals of $\eta_j$.
Consider a compactly supported prior $\pi_\xi$ on $\xi$, and a prior
$\pi_{\eta}$ on each $\eta_j$. Then
the unnormalized density of the full conditionals of each $\eta_j$ is given by
$\pi \big(\eta_{j} \mid  m_j = \xi \beta_j^2/(2 \sigma^2) \big) \propto \eta_j^{1/2} e^{- m_j \eta_j} \pi_{\eta}(\eta_j)$.

Consider the following assumptions on $\pi_\eta$.

* For some $0<c < 1/2$, there exists $0 < \epsilon < \Gamma(\frac{1}{2}-c)^{-1}\sqrt{\pi}$ and $K^{(1)}_{c, \epsilon} < \infty$ such that
$\mathbb{E}[ \eta_{j}^c | m_j] \leq \epsilon m_j^{-c} + K^{(1)}_{c, \epsilon}$ for all $m_j > 0$.

* For some $0 < d < 1$, there exists $0 < \epsilon < 2^d $ and $K^{(2)}_{d, \epsilon} < \infty$ such that
$\mathbb{E}[ \eta_{j}^{-d} | m_j] \leq \epsilon m_j^{d} + K^{(2)}_{d, \epsilon}$ for all $m_j > 0$.

Then the blocked Gibbs sampler in (<ref>) is geometrically ergodic.

We now consider the performance of Algorithm
<ref> on synthetic datasets. In particular, we empirically illustrate that the rate $\kappa_0$ in Proposition <ref> 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 <ref> 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
<ref>. 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
Figure <ref>, with a vertical axis on a logarithmic scale, shows
that the meeting times increase exponentially and have heavier tails with increasing dimension.

Meeting times against dimension for posteriors associated with a Horseshoe prior under the one-scale coupling of Algorithm <ref>,
with $n=100, s=10, \sigma_*=0.5$.

§ COUPLING STRATEGIES FOR HIGH-DIMENSIONAL SETTINGS

To address the scaling issues observed with the one-scale coupling in
Section <ref>,

we develop a two-scale coupling in Section
<ref> that delivers better empirical performance.
We also offer some formal results supporting the empirical improvements.
In Section <ref> we develop other related coupling strategies
which have exponentially-tailed meeting times and have similar empirical
performance as the two-scale coupling.
In Section <ref>, we investigate the
impact of the degree of freedom $\nu$ for Half-$t(\nu)$ priors on our two-scale
coupling. We find that Half-$t(\nu)$ priors with higher degrees of freedom
$\nu > 1$ can give similar statistical performance to the Horseshoe ($\nu=1$),
whilst enjoying orders of magnitude computational improvements through shorter
meeting times.

§.§ A two-scale coupling

We develop a two-scale coupling algorithm,
which slightly differs from the terminology used in stochastic dynamics

and in analysis of pre-conditioned HMC
It is also closely linked to the delayed one-shot coupling
of <cit.>, 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
<ref>, 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
<ref> and try to obtain exact meetings.
When $d(C_t, \tilde{C}_t) >d_0$, we instead
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

\begin{equation} \label{eq:distance_metric}
d(C_t, \tilde{C}_t) = \mathbb{P}\big( \eta_{t+1} \neq \tilde{\eta}_{t+1} | C_t, \tilde{C}_t \big),
\end{equation}

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 <ref>.  As the coupling in Algorithm
<ref> is independent component-wise, we can obtain a simplified expression involving only univariate probabilities as

d(C_t, C̃_t)    =  1 - ∏_j=1^p ℙ ( η_t+1,j = η̃_t+1,j |  η_t,j, η̃_t,j, m_t,j, m̃_t,j ) ,

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 <ref> since the maximal coupling slice sampler (Algorithm <ref>) 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

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

\begin{equation} \label{eq:distance_metric_estimate1}
\widehat{d_R}^{(1)}(C_t, \tilde{C}_t) = \frac{1}{R} \sum_{r=1}^R \mathbbm{1} \{  \eta^{(r)}_{t+1} \neq \tilde{\eta}^{(r)}_{t+1} \},
\end{equation}

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

\begin{equation} \label{eq:distance_metric_estimate3}
\widehat{d_R}^{(2)}(C_t, \tilde{C}_t) = 1 - \prod_{j=1}^p \Big(\frac{1}{R} \sum_{r=1}^R \mathbb{P}\Big( \eta_{t+1,j} = \tilde{\eta}_{t+1,j} \Big| U^{(r)}_{j}, \tilde{U}^{(r)}_{j}, m_{t,j}, \tilde{m}_{t,j} \Big) \Big),
\end{equation}

where $\big(U^{(r)}_{j}, \tilde{U}^{(r)}_{j} \big)$ are sampled independently for
$r=1,\ldots,R$ using Algorithm <ref>
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 <ref>.

Further metric calculation details are in Appendix <ref>.
For a number of samples $R$, estimates (<ref>)
and (<ref>) both have computation cost
of order $pR$. Compared to the estimate in (<ref>),
the estimate in (<ref>)
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 (<ref>) and (<ref>)
having low variance. Henceforth, we will use the estimate
in (<ref>)
in our experiments with two-scale couplings, unless specified otherwise.

Overall the two-scale coupling kernel
has twice the $\mathcal{O}(n^2p)$ cost of the single chain
kernel in (<ref>).
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
<ref>) 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

We now consider the performance of the two-scale coupling
on synthetic datasets. The setup is identical to that introduced in Section
<ref>, 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 <ref> 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 <ref> shows that the two-scale coupling
can lead to orders of magnitude reductions in meeting times, compared to the one-scale

Meeting times against dimension for posteriors associated with a Horseshoe prior under
one-scale and two-scale couplings,
with $n=100,s=10, \sigma_*=0.5$.

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.

Consider the setup of Proposition <ref>,
such that the $\eta$ full conditionals are sampled perfectly. Then
for any monotone function $h:(0, \infty) \rightarrow \mathbb{R}$ and
all components $j=1, \ldots, p$,

\begin{equation}
\mathbb{E}_{ \text{max} }
\big[ | h(\eta_{t+1,j}) - h(\tilde{\eta}_{t+1,j}) | \big| m_{t}, \tilde{m}_{t} \big] =
\mathbb{E}_{ \text{crn} }
\big[ | h(\eta_{t+1,j}) - h(\tilde{\eta}_{t+1,j}) | \big| m_{t}, \tilde{m}_{t} \big],
\end{equation}

where $\mathbb{E}_\text{max}$ and $\mathbb{E}_\text{crn}$ correspond to
expectations under the maximal coupling and CRN coupling respectively of
$(\eta_{t+1}, \tilde{\eta}_{t+1})$ given $m_t$ and $\tilde{m}_t$.

Proposition <ref> 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.
For the Horseshoe prior $(\nu=1)$,
consider when $m_t = \mu \mathbf{1}_p \in (0,\infty)^p$ and
$\tilde{m}_t= \tilde{\mu} \mathbf{1}_p \in (0,\infty)^p$ for some positive scalars
$\mu$ and $\tilde{\mu}$ with $\tilde{\mu} \neq \mu$. Then under a CRN coupling,

\begin{equation} \label{eq:crn}
\mathbb{P} \big( \max_{j=1}^p | \log(1+\eta_{t+1,j}) - \log(1+\tilde{\eta}_{t+1,j}) | > | \log \mu - \log \tilde{\mu} |
\ \big| \ m_t, \tilde{m}_t \big) = 0
\end{equation}

for all $p \geq 1$. Define a positive constant  $L=\log(1+1/\max \{\mu , \tilde{\mu} \})/2$. Under maximal coupling, for
any $\alpha \in (0,1)$ and any $\alpha' \in (\alpha,1)$, there
exists a constant $D>0$ which does not depend on $p$ such that for all $p>1$,

ℙ (max_j=1^p | log(1+η_t+1,j) - log(1+η̃_t+1,j) |
> log( αL ) + loglogp   |  m_t, m̃_t ) >
1 - e^-D p^1-α'.

Extensions of Proposition <ref> 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 <ref> 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 <ref> and
<ref> 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 <ref> 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 <ref> contains
discussions and related simulations comparing this metric with
alternatives. Proposition <ref>
verifies that such metric bounds the total variation distance between one-step transition kernels.

Let $\mathcal{P}$ denote the Markov transition kernel associated with the update from $(\beta_t,\xi_t, \sigma_t^2)$ to $(\beta_{t+1},\xi_{t+1}, \sigma_{t+1}^2)$ for the Gibbs sampler in (<ref>).
Let $Z = (\beta, \xi, \sigma^2)$, $\tilde{Z} = (\tilde{\beta}, \tilde{\xi}, \tilde{\sigma}^2)$,
$m_{j} = \xi \beta_{j}^2/(2 \sigma^2)$, and $\tilde{m}_{j} = \tilde{\xi} \tilde{\beta}_{j}^2/(2 \tilde{\sigma}^2)$
for $j=1, \ldots, p$. Then,

\begin{equation} \label{eq:tv_metric_upper_bound}
\text{TV} \big( \mathcal{P}( Z, \cdot), \mathcal{P}( \tilde{Z}, \cdot) \big)^2
\leq
% \frac{(1+\nu)(e-1)}{4}  \sum_{j=1}^p | \log(m_j) - \log(\tilde{m}_j) | =
\frac{(1+\nu)(e-1)}{4} \bar{d}(Z, \tilde{Z}).
\end{equation}

By Proposition <ref>, it suffices to show that
$\big(\bar{d}(Z_t, \tilde{Z}_t)\big)_{t \geq 0}$ contracts under CRN  to a
sufficiently small value such that the upper bound in
(<ref>) is strictly less that $1$. This would ensure
that CRN coupling will bring the marginal chains close enough that the probability of one-step max coupling resulting in a meeting is high. Note that (<ref>) corresponds to the total
variation distance of the full chain, so under this setup we would attempt to
maximally couple the full vectors $\eta_{t+1}$ and $\tilde{\eta}_{t+1}$ (rather
than maximally coupling component-wise) when close. Analytical calculations to
establish such a contraction, and to understand how the contraction rate scales with
dimension and other features of the problem such as the choice of
prior, sparsity, and signal to noise ratio
requires further work.

§.§ Alternative coupling strategies

The rapid meeting of the two-scale coupling in simulations recommends its use. However, we were unable to establish finite expected meeting times with this coupling, which prevents its use in constructing unbiased estimates. We now describe several small modifications of the two-scale coupling that have similar empirical performance but for which we can prove that the meeting time distribution has exponential tails.

Our basic approach is to

maximally couple each component $(\eta_{t,j}, \tilde{\eta}_{t,j})$ until a failed meeting
$\{ \eta_{t,j} \neq \tilde{\eta}_{t,j} \}$ occurs, after which we
switch to CRN for the remaining components.
The ordering of the components $j=1,\ldots,p$ can
be deterministic or randomized at each iteration. Under this construction,
only one component of $(\eta_{t}, \tilde{\eta}_{t})$, corresponding to
the failed meeting, will evolve independently; the other components will either meet
or evolve under CRN. Therefore, we expect to obtain benefits similar to the two-scale
coupling in high dimensions. An apparent advantage is that we can bypass the choice
of metric and associated threshold. While the switch-to-CRN coupling has the same $\mathcal{O}(n^2p)$
cost per iteration as the two-scale coupling, it can give faster run-times in practice when the
dimension $p$ is large compared to the number of observations $n$ as it avoids repeated calculation of a metric.

We report the performance of the switch-to-CRN coupling
with the ordering of the components
randomized at each iteration.

The setup is identical to that in Section
<ref>, where for each dimension $p$
we simulate $100$ meetings times based
on independently generated synthetic datasets.
Figure <ref> shows that the switch-to-CRN coupling
leads to similar meeting times as the two-scale coupling.

Meeting times for posteriors from the Horseshoe prior under the two-scale and the switch-to-CRN coupling algorithms.
$n=100,s=10, \sigma_*=0.5$.

We can establish that the meeting times under
this switch-to-CRN coupling have exponential tails and hence
finite expectation.

Consider the blocked Gibbs sampler in (<ref>).
We follow the setup in Proposition <ref>, and assume that the global shrinkage
prior $\pi_\xi(\cdot)$ has a compact support. Denote $Z_t = (\beta_t, \xi_t, \sigma^2_t)$ and
$\tilde{Z}_t = ( \tilde{\beta}_t, \tilde{\xi}_t, \tilde{\sigma}^2_t)$.
Consider the switch-to-CRN coupling given by
$(Z_t, \tilde{Z}_t) \mapsto (\eta_{t+1}, \tilde{\eta}_{t+1}) \mapsto (Z_{t+1}, \tilde{Z}_{t+1})$
where $(\eta_{t+1}, \tilde{\eta}_{t+1})$ are maximally coupled component-wise
(under any fixed or random ordering of the components) until the
first failed meeting, after which common random numbers are employed, and
$(Z_{t+1}, \tilde{Z}_{t+1})$ are coupled using common random numbers. Denote
the meeting time by
$\tau := \inf \{ t \geq 0: Z_t = \tilde{Z}_t \}$. Then

$\mathbb{P}( \tau > t ) \leq A_1 \kappa_1^t$

for some constants $A_1 \in (0, \infty)$ and $\kappa_1 \in (0,1)$, and for all $t \geq 0$.

As in Proposition <ref>, our proof of Proposition
<ref> implies a rate $\kappa_1$ which tends to $1$
exponentially as dimension $p$ increases. To obtain more favorable rates with
respect to the dimension as hinted by Figure
<ref> would require a better
understanding of the CRN coupling, for example establishing a
contraction in some parts of the state space as discussed above.

Overall Figures <ref> and <ref>
show significant empirical improvements in
using the two-scale coupling and related strategies, compared to the one-scale
coupling of Section <ref>. Therefore, we recommend the use
of two-scales couplings and variants such as the switch-to-CRN coupling
for practical high-dimensional applications. Furthermore, in settings where estimating
the metric $d$ is computationally expensive or analytically intractable,
variants such as the switch-to-CRN coupling may offer faster numerical run-times per iteration.
Henceforth, we will use the two-scale coupling of Section
<ref> in our experiments
and GWAS dataset applications, unless specified otherwise.

§.§ Computational and statistical impact of degree of freedom $\nu$

We empirically investigate the impact of the degree of freedom $\nu \geq 1$ for
Half-$t(\nu)$ priors.  Higher degrees of freedom $\nu$ corresponds to priors on
$\beta$ which have stronger shrinkage towards zero (Proposition
<ref>).  We consider meeting times resulting from
the two-scale coupling of Section <ref>, as well as
the statistical performance of the posteriors.  On the computation side, recent
work on convergence analysis <cit.>
has highlighted the impact of shrinkage on convergence of simpler MCMC
algorithms. Under a common random numbers coupling, <cit.>
have studied the convergence of the Gibbs sampler
<cit.> for Bayesian probit regression, proving dimension-free convergence rates for priors with
sufficiently strong shrinkage towards zero. On the statistical estimation side,
Half-$t(\nu)$ priors have been long proposed <cit.> for Bayesian hierarchical
models. <cit.> have established
(near) minimax optimal estimation for the Horseshoe ($\nu=1$) in the Normal
means model. <cit.> have extended the results of
<cit.> to show minimax optimality of a wider class of
priors including Half-$t(\nu)$ priors. Recently, <cit.>
has established that the degree of freedom $\nu$ can impact the multiplicative
constant in the posterior contraction rate for the Normal means model.
Posterior contraction rates of continuous shrinkage priors in the
regression setting beyond the Normal means model would deserve further

Meeting times

Mean Squared Error

Meeting times under the two-scale coupling algorithm and
mean squared error with Half-$t(\nu)$ priors and LASSO estimates.
In Figure <ref> $n=100,s=20$ and $\sigma_*=2$;
in Figure <ref> $n=100,p=500,s=20, \sigma_*=2$.

We consider the meeting times obtained with Half-$t(\nu)$ priors under the
two-scale coupling (with $R=1$ and threshold $d_0=0.5$) on synthetic datasets.
The synthetic datasets are generated as per Figure <ref>
of Section <ref>.
Figures <ref> is then based on $20$
synthetic datasets generated independently for different degrees of freedom
$\nu \in [1,2]$.  Figure <ref> plots
meeting times $\tau$ of $1$-lag coupled Markov chains against dimension.  It
indicates that even a small increase in the degree of freedom $\nu > 1$
compared to the Horseshoe ($\nu=1$) can lead to orders of magnitude
computational improvements through shorter meeting times.

We also consider the statistical performance of Bayesian shrinkage regression with
Half-$t(\nu)$ priors on synthetic datasets in Figure <ref>.
The synthetic datasets are generated as per Figure <ref> of Section
<ref>. We use the blocked Gibbs sampler of Algorithm
<ref> to draw samples from the posteriors corresponding to
different degrees of freedom $\nu$, by simulating chains of length $1000$, with
a burn-in period chosen using the coupling-based total variation distance upper bounds
of <cit.>.
Specifically we employ a burn-in of $300$ steps for $\nu \geq 1.2$ and of $600$ steps for $\nu = 1$.

Figure <ref> shows similar mean squared error (MSE) for different
values of $\nu \geq 1$. For a fixed synthetic dataset, the MSE in Figure
<ref> is calculated as $\| \beta_* - \hat{\beta} \|_{2}^2/p$
$\hat{\beta} := \sum_{t=1}^{1000} \beta_t / 1000$,
is the MCMC estimator obtained from the Markov chain after burn-in.
The error bars in Figure <ref> are generated by simulating synthetic datasets
independently $100$ times and each time calculating the corresponding MSE from a Markov chain
for different values of $\nu \in [1,2]$.
The MSE of the LASSO is also included,
with regularization parameter chosen with cross-validation using the
$\texttt{glmnet}$ package <cit.>.  For $\nu \in [1,2]$,
we observe a lower MSE using Half-$t(\nu)$ priors compared to the LASSO.
The grey plots in Figures <ref> show that much larger
$\nu \in \{4,6,8\}$ can lead to higher MSE, as the corresponding posteriors
are strongly concentrated about zero and less able to identify non-null signals.
Overall, Figure <ref> suggests
that Half-$t(\nu)$ priors with degrees of freedom $\nu \in [1,2]$
can result in comparable statistical performance
while much larger values of $\nu$ should be discouraged.

§ RESULTS ON GWAS DATASETS

Section <ref> suggests
that Half-$t(\nu)$ priors with higher degrees of freedom $\nu > 1$
can give similar statistical performance to the Horseshoe ($\nu = 1$),
whilst allowing orders of magnitude computational improvements.

Motivated by this observation, we apply our algorithms to genome-wide association study
(GWAS) datasets using $\nu = 2$. We consider a bacteria dataset <cit.>
with $n=71$ observations (corresponding to production of the vitamin riboflavin),
and $p=4,088$ covariates (corresponding to single nucleotide polymorphisms (SNPs) in the genome),
and a maize dataset <cit.>
with $n=2,266$ observations (corresponding to average number of days taken for silk emergence
in different maize lines), and $p=98,385$ covariates (corresponding to
SNPs). Bayesian methods are well-suited to such GWAS datasets, as they
provide interpretable notions of uncertainty through marginal posterior
probabilities which allow for multimodality,
enable the use of prior information, and account for the uncertainty
in parameter estimates <cit.>.
Furthermore, heavy-tailed continuous shrinkage priors can be particularly effective in the
GWAS setting, as the associations between SNPs and the phenotype are expected to be sparse,
and such priors have shown competitive empirical performance in polygenic prediction

For both datasets, we target the Half-$t(2)$ posterior with $a_0=b_0=1$.
We use the two-scale coupling with $R=1$ and threshold $d_0= 0.5$,
and initialize chains independently from the prior.
For the $\xi$ updates, we use a Metropolis–Rosenbluth–Teller–Hastings
step-size of $\sigma_{\text{MRTH}}=0.8$.
Based on 100 independent coupled chains,
Figure <ref> shows
upper bounds on the total variation distance to stationarity, $\text{TV}(\pi_t, \pi)$
as a function of $t$,
using $L$-lag couplings
with $L=750$ and $L=200$ for the maize and riboflavin datasets respectively.
The results indicate that these Markov chains converge to the corresponding
posterior distributions in less than $1000$ and $500$ iterations respectively.
We note that the lags $L=750$ and $L=200$ were selected based on
some preliminary runs of the coupled chains to obtain informative total variation bounds,
which incurs an additional preliminary cost. Any choice of lag results in valid upper bounds,
but only large enough lags lead to upper bounds that are in the interval $[0,1]$ even for small iterations $t$.

In these GWAS examples, the run-time per iteration may be significant.
For the maize dataset, on a 2015 high-end laptop,
each iteration of the coupled chain takes approximately $60$ seconds, and running
one coupled chain until meeting can take more than one day. This run-time is dominated
by the calculation of the weighted matrix cross-product $X \text{Diag}(\eta)^{-1} X^T$,
with an $\mathcal{O}(n^2p)$ cost. In this setting,
a scientist wanting to run chains for $10^4$ or $10^5$ iterations may have
to wait for weeks or months.
Using the proposed couplings, a scientist with access to parallel computers
can be confident that samples obtained after $10^3$ iterations
would be indistinguishable from perfect samples from the target;
for a similar cost they can also obtain unbiased estimators
as described in <cit.>.
For the riboflavin dataset, each iteration of the coupled chain takes
approximately $0.1$ seconds, and running one coupled chain until meeting takes
only $10$ to $20$ seconds.  The run-time there is dominated by the component-wise
$\eta$ updates which have $\mathcal{O}(p)$ cost.  In the riboflavin example,
the coupling-based diagnostics enable the scientist to perform reliable Bayesian computation
directly from personal computers. Overall, the GWAS examples emphasize that
couplings of MCMC algorithms can aid practitioners in high-dimensional

Riboflavin dataset with $n=71$ and $p=4,088$

Maize dataset with $n=2,266$ and $p=98,385$

$L$-Lag coupling-based upper bounds on $\text{TV}(\pi_t, \pi)$,
the distance between the chain at time $t$ and its limiting distribution,
in Bayesian shrinkage regression with Half-$t(2)$ prior on GWAS datasets.
Maize dataset: $n=2,266$ and $p=98,385$; Riboflavin dataset: $n=71$ and $p=4,088$.

§ DISCUSSION

We have introduced couplings of Gibbs samplers
for Bayesian shrinkage regression with Half-$t(\nu)$ priors.
The proposed two-scale
coupling is operational in realistic
settings, including a
GWAS setting with $n \approx 2,000, p \approx 100,000$.

Firstly, our work participates in a wider effort to apply MCMC algorithms
and coupling techniques to ever more challenging settings
<cit.>. In particular, the short meeting times
obtained here vindicate the use of coupling techniques in
multimodal, high-dimensional sampling problems.
This relates to questions raised in some comments of <cit.>
and <cit.> in the discussion of
<cit.>, and shows that the dimension
of the state space may not necessarily be the most important driver of
the performance of MCMC algorithms and couplings thereof.

Secondly, we have applied $L$-lag couplings
<cit.> to obtain upper bounds on the total variation
distance between the Markov chain at a finite time and its stationary distribution.
This allows empirical investigations of the mixing time
and how it varies
with the inputs of the problem, including number of observations,
number of covariates, signal to noise ratio and sparsity.
We find that coupling techniques constitute
a convenient  non-asymptotic tool to monitor the performance of MCMC algorithms.

Thirdly, we observe that Half-$t(\nu)$ priors with degrees of freedom $\nu$
higher than one give similar statistical estimation performance than the Horseshoe
whilst providing significant computational advantages.
This contributes towards the discussion on the impact of the prior
on the trade-off between statistical estimation and computational
feasibility, in the setting of high-dimensional Bayesian regression.

The following questions arise from our work.

* Convergence of the blocked Gibbs sampler.
Short meeting times suggest that the blocked Gibbs sampler converges quickly,
even in high dimensions. This motivates a more formal
study of the convergence
of that Markov chain, to understand how the convergence rate
varies with features of
the data generating process and of the prior. The Lyapunov function
in Proposition <ref> and the two-scale coupling may prove useful for such analysis.
Our initial work in this area suggests that finding a convenient metric that gives
sharp bounds on the metric $d$ used here, while simultaneously being amenable to
theoretical analysis, will be a key step.

* Alternative coupling algorithms.

Section <ref> indicates that many coupling strategies
are available in the present setting, and more generally for various
Gibbs samplers.
For example, we could combine the switch-to-CRN coupling with the
two-scale coupling to form other coupling strategies. Under the
two-scale coupling setup, we could apply the switch-to-CRN coupling only when
the chains are close with respect to a chosen metric, and
employ CRN couplings when far away. Alternatively, we could always apply the
switch-to-CRN coupling, and estimate component-wise meeting probabilities to
select the order of the updates, for example such that components more
likely to meet are updated first. Some couplings may result in shorter meeting times
for the Horseshoe prior. This may allow the Horseshoe prior to
remain competitive with Half-$t(\nu)$ priors
with higher degrees of freedom. In our experiments
we did not find ways to obtain shorter meetings
for the Horseshoe, but we certainly cannot rule out that
possibility. If one could obtain short meetings for the Horseshoe
in high dimensions, the apparent trade-off between statistical performance and computing cost
identified in the present work would disappear.

* Interplay between posterior concentration and MCMC convergence.
For Bayesian regression with spike-and-slab priors, <cit.> and
<cit.> have shown that
posterior contraction can aid the convergence of MCMC in
high dimensions. The performance of the proposed couplings
motivates similar investigations for continuous shrinkage priors.

The authors thank Xiaolei Liu and Xiang Zhou for sharing the Maize GWAS data,
and Yves F. Atchadé, Qian Qin and Neil Shephard for helpful comments.
Niloy Biswas gratefully acknowledges support from a GSAS Merit Fellowship and
a Two Sigma Fellowship Award.
Anirban Bhattacharya gratefully acknowledges support
by the National Science Foundation through an NSF CAREER award (DMS-1653404).
Pierre E. Jacob gratefully acknowledges support by the National Science Foundation
through grants DMS-1712872 and DMS-1844695.
The maize dataset GWAS simulations in this article were run on the FASRC Cannon cluster supported
by the FAS Division of Science Research Computing Group at Harvard University.

In Section <ref>, we describe some features of the posterior distributions arising in high-dimensional regression with shrinkage priors.
In Section <ref> we
describe metrics to monitor the distance between pairs of chains and
provide details on the calculation of metrics to use in two-scale coupling strategies.
In Section <ref> we provide the proofs for all the theoretical results in this document.

Section <ref> provides pseudo-codes for various algorithms mentioned in this document. Section <ref> provides some derivations of certain steps in the Gibbs sampler.

§ FEATURES OF THE TARGET COMPLICATE USE OF GENERIC MCMC

While posterior distributions resulting from heavy-tailed shrinkage priors

have desirable statistical properties, their
features pose challenges to generic MCMC

For Half-$t(\nu)$ priors, we show that the resulting posterior distributions
present 1) multimodality, 2) heavy tails and 3) poles at zero. This hints at a trade-off
between statistical accuracy and computational difficulty, since the very
features that present computational challenges are crucial for optimal
statistical performance across sparsity levels and signal strengths.

Proposition <ref> gives the marginal prior and
posterior densities of $(\beta, \sigma^2, \xi)$ up to normalizing constants,
i.e. integrating over $\eta$. For any real-valued functions $f$ and $g$, we
write $f(x) \asymp g(x)$ for $x \rightarrow x_0$ when $m < \liminf_{x
\rightarrow x_0} \big| f(x)/g(x) \big| \leq \limsup_{x \rightarrow x_0}
\big| f(x)/g(x) \big| < M $ for some $0<m\leq M < \infty$.

The marginal prior of component $\beta_j$ on $\mathbb{R}$ given $\xi$ and $\sigma^2$ has density

\begin{equation} \label{eq:marginal_beta_prior}
\pi(\beta_j | \xi, \sigma^2) \propto U \Big( \frac{1+\nu}{2}, 1, \frac{\xi \beta_j^2}{2 \sigma^2 \nu} \Big) \asymp \begin{cases}
- \log|\beta_j| & \text{ for } \; |\beta_j| \rightarrow 0 \\
|\beta_j|^{-(1+\nu)} & \text{ for } \; |\beta_j| \rightarrow +\infty
\end{cases},
\end{equation}

where $\Gamma(\cdot)$ is the Gamma function, and $U(a,b,z) := \Gamma(a)^{-1} \int_0^\infty x^{a-1}(1+x)^{b-a-1}e^{-zx}dx$ is the confluent hypergeometric function of the second kind <cit.> for any $a, b, z>0$.
The marginal posterior density of $\beta$ on $\mathbb{R}^p$ given $\xi$ and $\sigma^2$ is

π( β| σ^2, ξ,  y )    ∝𝒩( y; X β, σ^2 I_n ) ∏_j=1^p U ( 1+ν/2, 1, ξβ_j^2/2 σ^2 ν ),

thus $ \pi(\beta | \sigma^2, \xi,  y ) \overset{\|\beta \| \rightarrow 0}{\asymp} - \prod_{j=1}^p \log ( |\beta_j| ) $ and $- \frac{\partial }{\partial \beta_j} \log \pi(\beta | \sigma^2, \xi,  y ) \overset{\beta_j \rightarrow 0}{\asymp} - (\beta_j \log |\beta_j|)^{-1}  $.

Figure <ref> gives an example of such marginal posterior
in a stylized setting. Here

\[n=2, \quad p=3, \quad X = \begin{pmatrix} 1 & 1
& 0\\ 1 & 0 & 1 \end{pmatrix},\quad y = X  (1\; 0\; 0)^T = (1\; 1)^T,\]

with $\nu=2$ and $\xi = \sigma^2=0.25$.

Component $\beta_1$ in Figure <ref> illustrates
the potential multimodality in the posterior. The posterior distribution
also has polynomial tails in high dimensions, as stated in Corollary

Posterior densities $\pi( \beta_1, \beta_2 | \sigma^2, \xi,  y ), \pi( \beta_1, \beta_3 | \sigma^2, \xi,  y )$ and $\pi( \beta_2, \beta_3 | \sigma^2, \xi,  y )$.

Let $e^{\perp} \in \mathbb{R}^p$ be a unit vector such that $X e^{\perp} = 0$
and let $\lambda>0$. Assume entries $e^{\perp}_j \in \mathbb{R} $ are non-zero
for all $j=1,\ldots,p$, which is necessary to ensure that
$\pi( \lambda e^{\perp} | \sigma^2, \xi,  y )$ is
finite for all $\lambda>0$. Then

\begin{equation}
\pi( \lambda e^{\perp} | \sigma^2, \xi,  y ) \asymp \lambda^{-p(1+\nu)}
\end{equation}

as $\lambda \rightarrow + \infty$.
That is, the marginal posterior of $\beta$ on $\mathbb{R}^p$
has polynomial tails along all directions in the null space $\ker(X)$ of $X$.

In our stylized example, $e^{\perp} = 3^{-1/2} (1\, -1 \, -1)^T$
satisfies $X e^{\perp} = 0$. It follows from Corollary
<ref> that $ \pi( \lambda e^{\perp}
| \sigma^2, \xi,  y ) \asymp \lambda^{-3(1+\nu)} \asymp \lambda^{-9}) $ for
large $\lambda$. More generally, whenever $p > n$ and the marginal prior on
$\beta$ has polynomial tails, the subspace $\ker(X)$ of $\mathbb{R}^p$ has
vector space dimension at least $(p-n)$, and therefore the posterior has
polynomial tails over a non-trivial subspace. For such heavy-tailed targets,
the most popular “generic” MCMC algorithms encounter significant theoretical
and practical challenges. Here, we use the term generic for algorithms that
are not explicitly designed to target this class of distributions, but
instead rely on a general strategy to explore state spaces.

Specifically, for heavy-tailed targets, lack of geometric ergodicity
has been established for Random Walk Metropolis–Hastings (RWMH)
<cit.>, Metropolis-Adjusted Langevin (MALA)
<cit.> and Hamiltonian Monte Carlo (HMC)
and <cit.> show lack of geometric ergodicity for
targets which have heavy-tails in every direction.  Their arguments can be
modified to hold for targets which are heavy-tailed along any one direction as
in Corollary <ref>, by considering a
projection of the full chain along that direction.]. Better performance
can sometimes be obtained using heavy-tailed proposals
<cit.> or
variable transformations when available <cit.>.

The posterior density in (<ref>)
also approaches infinity when any component $\beta_j$ approaches zero, as
$\lim_{\beta_j \rightarrow 0} U( (1+\nu)/2,1, (\xi \beta_j^2)/(2 \sigma^2 \nu)) = \infty$
for any $\nu, \xi, \sigma^2> 0$. This is highlighted by the modes at the origin in Figure <ref>.
The pole at zero corresponds to a discontinuity in the gradients of the negative log-density, as given in Proposition
Such diverging gradients can lead to numerically unstable
leapfrog integration for HMC samplers

§ DETAILS ON THE CHOICE AND CALCULATION OF METRIC

§.§ Comparison with other metrics

Figure <ref> plots
$d(C_t, \tilde{C}_t)$ alongside the $L_1$ metric and the log-transformed
$L_1$ metric on parameters $m_{t,j}$ and $\tilde{m}_{t,j}$ as in
(<ref>). It is based on one trajectory of our
two-scale coupled chain under common random numbers (threshold $d_0=0$),
for a synthetic dataset generated as per Section <ref> targeting the Horseshoe
posterior ($\nu=1$) with $a_0=b_0=1$. Figure <ref> shows
that $d(C_t, \tilde{C}_t)$ is bounded above by both the $L_1$ metric and
the $L_1$ metric between the logarithms of $m_j, \tilde{m}_j$. Therefore, we
could consider the capped $L_1$ metric ($\min \{ \sum_j | m_{t,j} -
\tilde{m}_{t,j} | , 1 \}$) or the capped $L_1$ metric on the logs ($\min \{
\sum_j | \log( |m_{t,j}/ \tilde{m}_{t,j}|) | , 1 \}$) to obtain upper
bounds on $d$. Finding metrics which are easily calculable and accurate
could be investigated further and  may prove valuable
in the theoretical analysis of the proposed algorithms.

Metric trajectories for the posterior associated with a Horseshoe prior
under common random numbers coupling
with $n=100$, $p=1000$, $s=10$ and $\sigma_*=0.5$.

§.§ Choice of metric threshold for the two-scale coupling

For different values of $d_0$, Figure
<ref> shows averaged trajectories $(d(C_t,
\tilde{C}_t)) _{t \geq 0}$. When $d_0=1$ (one-scale coupling), metric
$d(C_t, \tilde{C}_t)$ remains near $1$, and when $d_0=0$ (common
random numbers coupling) metric $d(C_t, \tilde{C}_t)$ gets very close to
but does not exactly equal $0$. For all values of the threshold shown,
the chains exactly meet and $d(C_t, \tilde{C}_t)$ reaches zero at similar
times. Figure <ref> considers higher dimensional settings. It
plots the histograms of $(d(C_t, \tilde{C}_t))_{t=1}^{2000}$ when
$d_0=0$. The histograms show that $d(C_t, \tilde{C}_t)$ takes
values close to either $1$ or $0$ for different dimensions $p$. This

may explain why different threshold values sufficiently away from $0$ and $1$ give
similar meeting times. Figure <ref> shows that different thresholds
give similar meeting times for the Horseshoe.

Two-scale coupling algorithm performance for the Horseshoe
prior under different values of $d_0$.
Unless specified otherwise, $n=100,p=500,s=20$ and $\sigma_*=2$.

§.§ Calculating metric $d$

We present details of estimating the metric $d$ defined in (<ref>).
As the coupling in Algorithm <ref> is independent component-wise, we have

d(C_t, C̃_t)    = 1 - ℙ( η_t+1 = η̃_t+1 | C_t, C̃_t ) = 1 - ∏_j=1^p ℙ ( η_t+1,j = η̃_t+1,j |  η_t,j, η̃_t,j, m_t,j, m̃_t,j )

for $m_j = \xi \beta_j^2 / (2 \sigma^2)$
and $\tilde{m}_j = \tilde{\xi} \tilde{\beta}_j^2 / (2 \tilde{\sigma}^2)$.

§.§.§ An integration based approximation

We consider the approximation

\begin{equation}
\mathbb{P} \big( \eta_{t+1,j} \neq \tilde{\eta}_{t+1,j} |  \eta_{t,j}, \tilde{\eta}_{t,j}, m_{t,j}, \tilde{m}_{t,j} \big) \approx \text{TV}(M_{t,j}, \tilde{M}_{t,j}),
\end{equation}

where $M_{t,j}$ and $\tilde{M}_{t,j}$ correspond to target distributions of the
marginal slice samplers in Algorithm <ref>.
This corresponds to approximating the probability of not meeting under
Algorithm <ref>
with the total variation distance (i.e. probability of not meeting under a maximal coupling)
between the marginal target distributions $M_{t,j}$ and $\tilde{M}_{t,j}$ of the slice samplers.

We can evaluate $\text{TV}(M_{t,j}, \tilde{M}_{t,j})$ analytically using Proposition <ref>.
This motivates the approximation
d̃(C_t, C̃_t) := 1 - ∏_j=1^p ( 1 - TV(M_t,j, M̃_t,j) )

for $d(C_t, \tilde{C}_t)$.

Consider distributions $M$ and $\tilde{M}$ on $[0,\infty)$ with densities

\begin{equation} \label{eq:slice_sampler_target_density}
p(\eta_j | m) \propto \frac{1}{\eta_j^{\frac{1-\nu}{2}} (1+ \nu \eta_j)^{\frac{\nu+1}{2}}} e^{-m \eta_j} \text{ and }
p(\eta_j | \tilde{m}) \propto \frac{1}{\eta_j^{\frac{1-\nu}{2}} (1+ \nu \eta_j)^{\frac{\nu+1}{2}}} e^{- \tilde{m} \eta_j},
\end{equation}

respectively on $[0,\infty)$, such that $M$ and $\tilde{M}$ are parameterized by the fixed constant $m > 0$ and $\tilde{m} > 0$.

Then, $\text{TV}(M, \tilde{M}) = 0$ when $m=\tilde{m}$. When $m \neq \tilde{m}$,

TV(M, M̃) = | U_1+ν/2, 1, m/ν(νK)/U(1+ν/2, 1, m/ν) - U_1+ν/2, 1, m̃/ν(νK)/U(1+ν/2, 1, m̃/ν) |,

where $U_{a,b,z}(t) := \frac{1}{\Gamma(a)} \int_0^t x^{a-1}(1+x)^{b-a-1} e^{-zx} dx$ is defined as the lower incomplete confluent hypergeometric function of the second kind, such that $U_{a,b,z}(\infty) = U(a,b,z)$ gives the confluent hypergeometric function of the second kind, and

\begin{equation}
K := \frac{\log U(\frac{1+\nu}{2}, 1, \frac{m}{\nu})  - \log U(\frac{1+\nu}{2}, 1, \frac{\tilde{m}}{\nu}) }{\tilde{m}-m}.
\end{equation}

We have

\begin{equation}
p(\eta_j | m) = \frac{1}{\eta_j^{\frac{1-\nu}{2}} (1+ \nu \eta_j)^{\frac{\nu+1}{2}}} e^{-m \eta_j} \frac{1}{Z}
\ \text{ and } \
p(\eta_j | \tilde{m}) = \frac{1}{\eta_j^{\frac{1-\nu}{2}} (1+ \nu \eta_j)^{\frac{\nu+1}{2}}} e^{-\tilde{m} \eta_j} \frac{1}{\tilde{Z}},
\end{equation}

Z    = ∫_0^∞1/η_j^1-ν/2 (1+ νη_j)^ν+1/2 e^-m η_j dη_j = ν^-ν+1/2 Γ( ν+1/2 ) U (ν+1/2, 1, m/ν ),

Z̃    = ∫_0^∞1/η_j^1-ν/2 (1+ νη_j)^ν+1/2 e^-m̃ η_j dη_j = ν^-ν+1/2 Γ( ν+1/2 ) U (ν+1/2, 1, m̃/ν ).

Case $m = \tilde{m}$. $\text{TV}(M, \tilde{M})=0$ is immediate.

Case $m < \tilde{m}$. Note that

p(η_j | m) ≤p(η_j | m̃)    ⇔exp( (m̃-m) η) ≤Z/Z̃

⇔η≤log( U(1+ν/2, 1, m/ν) ) - log( U(1+ν/2, 1, m̃/ν) ) /m̃-m =: K.

1 - TV(M, M̃)    = ∫_0^∞min(p( η_j | m), p( η_j | m̃)) d η_j

= ∫_0^K p( η_j | m) d η_j + ∫_K^∞p( η_j | m̃) d η_j

= 1/Z ∫_0^K 1/η_j^1-ν/2 (1+ νη_j)^ν+1/2 e^-m η_j d η_j

+ 1/Z̃ ∫_K^∞1/η_j^1-ν/2 (1+ νη_j)^ν+1/2 e^-m̃ η_j d η_j

= U_1+ν/2, 1, m/ν(νK)/U(1+ν/2, 1, m/ν) + 1 - U_1+ν/2, 1, m̃/ν(νK)/U(1+ν/2, 1, m̃/ν)

= 1 - | U_1+ν/2, 1, m̃/ν(νK)/U(1+ν/2, 1, m̃/ν) - U_1+ν/2, 1, m/ν(νK)/U(1+ν/2, 1, m/ν) |.

Equation (<ref>) directly follows.

Case $m > \tilde{m}$. Follows from the $m < \tilde{m}$ case by symmetry.

§.§.§ A Rao–Blackwellized estimator

We have

ℙ  ( η_t+1,j = η̃_t+1,j |  η_t,j, η̃_t,j, m_t,j, m̃_t,j )

= 𝔼 [ℙ ( η_t+1,j = η̃_t+1,j | U_j,*, Ũ_j,*, m_t,j, m̃_t,j ) | η_t,j, η̃_t,j, m_t,j, m̃_t,j ],

where $m_{t,j} = \xi_{t} \beta_{t,j}^2 / (2 \sigma_{t}^2)$ and
$\tilde{m}_{t,j} = \tilde{\xi}_{t} \tilde{\beta}_{t,j}^2 / (2 \tilde{\sigma}_{t}^2)$,
$(U_{j,*}, \tilde{U}_{j,*}) | m_{t,j}, \tilde{m}_{t,j}$ are sampled using common random numbers as in Algorithm <ref>,
and $(\eta_{t+1,j}, \tilde{\eta}_{t+1,j}) | U_{j,*}, \tilde{U}_{j,*}, m_{t,j}, \tilde{m}_{t,j}$ are sampled as in Algorithm <ref>.
We can analytically evaluate probability $\mathbb{P} \big( \eta_{t+1,j} = \tilde{\eta}_{t+1,j} \Big| U_{j,*}, \tilde{U}_{j,*}, m_{t,j}, \tilde{m}_{t,j} \big)$, which motivates the Rao-Blackwellized estimator from (<ref>),

\begin{equation}
\widehat{d_R}^{(2)}(C_t, \tilde{C_t}) := 1 - \prod_{j=1}^p \Big(\frac{1}{R} \sum_{r=1}^R \mathbb{P} \Big( \eta_{t+1,j} = \tilde{\eta}_{t+1,j} \Big| U_{j,r}, \tilde{U}_{j,r}, m_{t,j}, \tilde{m}_{t,j} \Big) \Big),
\end{equation}

where $R$ is the number of samples. For each $r=1,\ldots,R$, $(U_{j,r}, \tilde{U}_{j,r})$ is sampled independently as in Algorithm <ref> and maximal coupling probability $\mathbb{P} \big( \eta_{t+1,j} = \tilde{\eta}_{t+1,j} \big| U_{j,r}, \tilde{U}_{j,r}, m_{t,j}, \tilde{m}_{t,j} \big)$ is calculated analytically using Proposition <ref>.

Suppose $\eta_j \big| U_{j,*}, m_{j} \sim P_j$ and $\tilde{\eta}_j \big| \tilde{U}_{j,*}, \tilde{m}_{j} \sim \tilde{P}_j$
where marginal distributions ${P}_j$ and $\tilde{P}_j$ correspond to the distribution in the second step of the slice sampler
in Algorithm <ref>. Then,

ℙ     ( η_j = η̃_j | U_j,*, Ũ_j,*, m_j, m̃_j )
= {
γ_s(m_j K̃_j)/γ_s(m_j T_j,*) + γ_s(m̃_j (T_j,* ∧T̃_j,*))-γ_s(m̃_j K̃_j)/γ_s(m̃_j T̃_j,*)     for  m_j < m̃_j

γ_s(m_j (T_j,* ∧T̃_j,*))/γ_s(m_j (T_j,* ∨T̃_j,*))     for  m_j = m̃_j

γ_s(m̃_j K̃_j)/γ_s(m̃_j T̃_j,*) + γ_s(m_j (T_j,* ∧T̃_j,*))-γ_s(m_j K̃_j)/γ_s(m_j T_j,*)     for  m_j > m̃_j

for $m_{t,j} = \xi_{t} \beta_{t,j}^2 / (2 \sigma_{t}^2)$ and
$\tilde{m}_{t,j} = \tilde{\xi}_{t} \tilde{\beta}_{t,j}^2 / (2 \tilde{\sigma}_{t}^2)$,
$T_{j,*} :=(U_{j,*}^{-\frac{2}{1+\nu}}-1)/\nu$ and
$\tilde{T}_{j,*} :=(\tilde{U}_{j,*}^{-\frac{2}{1+\nu}}-1)/\nu$,
$s = \frac{1 + \nu}{2}$,
$\gamma_{s}(x) = \frac{1}{\Gamma(s)} \int_0^x t^{s-1} e^{-t} dt \in [0,1]$ the regularized incomplete lower Gamma function, $a \wedge b := \min \{ a, b \} $, $a \vee b := \max \{ a, b \} $, and
K̃_j := ( 0 ∨(  log((m̃_j/m_j)^s γ_s(m_j T_j,*)/γ_s(m̃_j T̃_j,*)))/m̃_j-m_j ) ) ∧( T_j,* ∧T̃_j,* ).

We work component-wise and drop the subscripts for ease of notation. Then distributions
$P$ and $\tilde{P}$ have densities

p(η; m, T)    = η^s-1 e^-m η m^s/Γ(s) γ_s(m T)    on  (0, T]

p(η; m̃, T̃)    = η^s-1 e^- m̃ η  m̃^s/Γ(s) γ_s( m̃ T̃)
on  (0, T̃]

respectively where $T :=(U_{*}^{-\frac{2}{1+\nu}}-1)/\nu$ and $\tilde{T} :=(\tilde{U}_{*}^{-\frac{2}{1+\nu}}-1)/\nu$,
$s = \frac{1+\nu}{2}$ and \(\gamma_{s}(x) = \frac{1}{\Gamma(s)} \int_0^x t^{s-1} e^{-t} dt \in [0,1] \) is the regularized incomplete lower Gamma function. Note that

ℙ_max( η= η̃ | T, T̃, m, m̃ ) = ∫_0^T ∧T̃ p(η; m, T) ∧p(η; m̃, T̃) d η.

Case $m = \tilde{m}$.

∫_0^T ∧T̃ p(η; m, T) ∧p(η; m̃, T̃) d η   = ∫_0^T ∧T̃ η^s-1 (m)^s/Γ(s) ( 1/γ_s(m T) ∧1/γ_s(m T̃) ) d η

= ∫_0^T ∧T̃ η^s-1 (m)^s/Γ(s) 1/γ_s(m(T ∨T̃) )

= γ_s(m(T ∧T̃) )/γ_s(m(T ∨T̃) )

as required.

Case $m < \tilde{m}$.  For $0 < \eta \leq T \wedge \tilde{T}$,

p(η; m, T) ≤p(η; m̃, T̃)    ⇔exp((m̃-m)η) ≤(m̃/m)^s γ_s(m T)/γ_s(m̃ T̃)
⇔η≤ log((m̃/m)^s γ_s(m T)/γ_s(m̃ T̃)))/m̃-m =: K.

Denote $\tilde{K} := \big( 0 \vee K \big) \wedge \big( T \wedge \tilde{T} \big) $, such that

K̃ := {
0     when  K ≤0

K      when  0 < K ≤T ∧T̃

T ∧T̃     when   T ∧T̃ < K.

This gives

∫_0^T ∧T̃ p(η; m, T) ∧p(η; m̃, T̃) d η   = ∫_0^K̃ p(η; m, T) d η+ ∫_K̃^T ∧T̃ p(η; m̃, T̃) d η

= ( Γ(s) γ_s(m K̃)/(m)^s ) ( (m)^s/Γ(s) γ_s(m T))

+ ( Γ(s) (γ_s(m̃ (T ∧T̃))-γ_s(m̃ K))/(m̃)^s ) ( (m̃)^s/Γ(s) γ_s(m̃ T̃)))

= γ_s(m K̃)/γ_s(m T) + γ_s(m̃ (T ∧T̃))-γ_s(m̃ K̃)/γ_s(m̃ T̃).

as required.

Case $m > \tilde{m}$.  Follows from the $m < \tilde{m}$ case by symmetry.

§ PROOFS

§.§ Marginal $\beta$ prior and posterior densities

Marginal Prior on $\beta$.
Let $\pi_{\eta}$ denote the component-wise prior density of each $\eta_j$ based on the Half-$t(\nu)$ distribution, and
$\mathcal{N}( \cdot; \mu, \Sigma )$ the density of a multivariate Normal distribution with mean $\mu$ and $\Sigma$.
It suffices to calculate the marginal prior for each $\beta_j$ given $\xi, \sigma^2$, by marginalizing over each $\eta_j$.
We obtain,

π(β_j | ξ, σ^2)    = ∫_0^∞𝒩 ( β_j; 0, σ^2/ξη_j ) π_η(η_j) dη_j

∝∫_0^∞ √(ξη_j) /√(2 πσ^2) exp( - ξβ_j^2/2σ^2 η_j ) 1/η_j^2-ν/2(1+νη_j)^ν+ 1/2 dη_j

∝∫_0^∞1/η_j^1-ν/2(1+νη_j)^ν+ 1/2 exp( - ξβ_j^2/2σ^2 η_j )  dη_j

= ∫_0^∞1/(t/ν)^1-ν/2(1+t)^ν+ 1/2 exp( -ξβ_j^2/2σ^2 ν t) dt 1/ν

= Γ(ν+ 1/2)/ν^ν+ 1/2 U ( ν+1/2, 1, ξβ_j^2/2σ^2 ν )

∝U ( ν+1/2, 1, ξβ_j^2/2σ^2 ν ) ,

where (<ref>) follows from the definition of the confluent hypergeometric function of the second kind: $U(a,b,z) := \frac{1}{\Gamma(a)} \int_0^\infty x^{a-1}(1+x)^{b-a-1}e^{-zx}dx$
for any $a, b, z>0$. By the asymptotic expansion of the confluent hypergeometric function of the second kind <cit.>, we obtain

\begin{align} \label{eq:prop1proofeqmarginal_beta_prior}
\pi(\beta_j | \xi, \sigma^2) \propto U \Big( \frac{\nu+1}{2}, 1, \frac{\xi \beta_j^2}{2\sigma^2 \nu} \Big) \asymp \begin{cases}
-\frac{1}{\Gamma(\frac{1+\nu}{2})} \log \Big( \frac{\xi \beta_j^2}{2 \sigma^2 \nu} \Big) & \asymp - \log (|\beta_j|) \text{ for } |\beta_j| \rightarrow 0 \\
\Big( \frac{\xi \beta_j^2}{2 \sigma^2 \nu} \Big)^{-\frac{1+\nu}{2}} & \asymp  |\beta_j|^{-(1+\nu)} \text{ for } |\beta_j| \rightarrow +\infty
\end{cases}
\end{align}

as required for (<ref>).

Marginal Posterior of $\beta$.
Let $\pi_\xi$ denote the prior density on $\xi$ and $\pi_{\sigma^2}$ the $\InvGamma(\frac{a_0}{2}, \frac{b_0}{2})$ prior density on $\sigma^2$. Then,

π(β| ξ, σ^2, y) ∝𝒩( y; Xβ, σ^2 I_n ) ( ∏_j=1^p π( β_j | ξ, σ^2 ) ) ∝𝒩( y; Xβ, σ^2 I_n ) ∏_j=1^p U ( 1+ν/2, 1, ξβ^2_j/2 σ^2 ν )

as required for (<ref>).

Equations (<ref>)
and (<ref>) directly give

\begin{equation}
\pi (\beta | \xi, \sigma^2, y) \asymp \prod_{j=1}^p U \Big( \frac{\nu+1}{2}, 1, \frac{\xi \beta_j^2}{2\sigma^2 \nu} \Big)  \asymp - \prod_{j=1}^p \log | \beta_j | \text{ for } \ \|\beta \| \rightarrow 0.
\end{equation}

Gradient of negative log-density of posterior of $\beta$.
From (<ref>), we obtain

-∂/∂β_j ( logπ(β| ξ, σ^2, y) )    =
-∂/∂β_j ( - 1/2 σ^2 y-X β_2^2 + ∑_j=1^p logU ( 1+ν/2, 1, ξβ^2_j/2 σ^2 ν )  )

= -[1/σ^2 X^T (y-X β)]_j -   d/d β_j U ( 1+ν/2, 1, ξβ^2_j/2 σ^2 ν )/U ( 1+ν/2, 1, ξβ^2_j/2 σ^2 ν ).

Recall  <cit.> for $a,b > 0$,

\begin{align}
\frac{d}{dz} U(a,b,z) \overset{ \forall z > 0}{=} -a U(a+1,b+1,z), \
U(a,2,z) \overset{z \rightarrow 0^+}{\asymp} \frac{1}{\Gamma(a)} \frac{1}{z}, \
U(a,1,z) \overset{z \rightarrow 0^+}{\asymp} -\frac{1}{\Gamma(a)} \log(z).
\end{align}

This gives

d/d β_j U ( 1+ν/2, 1, ξβ^2_j/2 σ^2 ν )/U ( 1+ν/2, 1, ξβ^2_j/2 σ^2 ν )    =
- 1+ν/2 U ( 1+ν/2+1, 1+1, ξβ^2_j/2 σ^2 ν )/U ( 1+ν/2, 1, ξβ^2_j/2 σ^2 ν ) ξβ_j/ σ^2 ν β_j →0≍ 1/ β_j log| β_j |.

Therefore by (<ref>), $-\frac{\partial}{\partial \beta_j} \Big( \log \pi(\beta | \xi, \sigma^2, y) \Big) \overset{\beta_j \rightarrow 0}{\asymp} - \frac{1}{ \beta_j \log | \beta_j |}$ as required.

By Proposition <ref>,

π(β| ξ, σ^2, y) = 1/Z_ξ, σ^2, ν 𝒩( y; Xβ, σ^2 I_n ) ∏_j=1^p U ( 1+ν/2, 1, ξβ^2_j/2 σ^2 ν ),

$Z_{\xi, \sigma^2, \nu}:= \int_{\mathbb{R}^p} \mathcal{N}( y; X\beta, \sigma^2 I_n ) \prod_{j=1}^p U \Big( \frac{1+\nu}{2}, 1, \frac{\xi \beta^2_j}{2 \sigma^2 \nu} \Big) d \beta$
is the normalizing constant. As $p >n$, there exists some $e^{\perp} \in \mathbb{R}^p$ such that $Xe^{\perp}=0$.
When $e^{\perp}_j=0$ for some $j=1,\ldots,p$, the posterior density is infinite
by Proposition <ref>.
Consider $e^{\perp} \in \mathbb{R}^p$
such that each $e^{\perp}_j \in \mathbb{R} $ is non-zero for $j=1,\ldots,p$.
For any such fixed such $e^{\perp}$ and large $\lambda>0$, we obtain

π(λe^⊥ | ξ, σ^2, y)    = 1/Z_ξ, σ^2, ν 𝒩( y; 0, σ^2 I_n ) ∏_j=1^p U ( 1+ν/2, 1, ξλ^2 (e^⊥_j)^2/2 σ^2 ν )

∝∏_j=1^p U ( 1+ν/2, 1, ξλ^2 (e^⊥_j)^2/2 σ^2 ν )

≍∏_j=1^p ( ξλ^2 (e^⊥_j)^2/2 σ^2 ν )^-1+ν/2

where (<ref>) follows from the asymptotic expansion of the
confluent hypergeometric function of the second kind <cit.>.

§.§ Geometric Ergodicity

Initial technical results.
We first record some technical results.
Consider any matrix $A$ such that $A=\begin{bmatrix} B \\ C \end{bmatrix}$
for submatrices $B$ and $C$. Then $\|A\|_2 \leq \|B\|_2 + \|C\|_2$, where $\| \cdot \|_2$ is the spectral norm of a matrix.

Take any vector $v$ of length equal to the number of columns of $A$ with $\| v \|_2 =1$. Then

Av _2^2 =   B
C  v _2^2 =  Bv _2^2 +  Cv _2^2 ≤B_2^2 + C_2^2.

Taking the supremum over all such $v$ with $\| v \|_2 =1$ gives
$\|A\|_2^2 \leq \|B\|_2^2 + \|C\|_2^2 \leq (\|B\|_2 + \|C\|_2)^{2}$.

(Sherman–Morrison–Woodbury formula <cit.>)
For any $p\times p$ invertible matrix $A$, $p\times n$ matrix $U$
and $n\times p$ matrix $V$,

\begin{equation} \label{eq:woodbury}
(A+UV)^{-1}=A^{-1} - A^{-1}U(I_n+V A^{-1} U )^{-1} V A^{-1}.
\end{equation}

When $U=V^T=x \in \mathbb{R}^p$, we obtain
$(A+xx^T)^{-1}=A^{-1} - A^{-1} x  x^T A^{-1}/(1+x^T A^{-1} x)$.

(A uniform bound on the full conditional mean of $\beta$).
Let $X$ be an $n\times p$ matrix with columns $x_i \in \mathbb{R}^n$ for $i=1, \ldots, p$,
and let $R=\text{Diag}(r_1,\ldots,r_p)$ be a $p\times p$ diagonal matrix with positive entries. Then

\begin{equation}
\sup_{r_1, \ldots, r_p>0} \|(X^TX + R)^{-1} X^T\|_2 \leq  C_p,
\end{equation}

for some $C_p < \infty$, where $\| \cdot \|_2$ is the spectral norm of a matrix.
That is, the spectral norm of $(X^TX + R)^{-1} X^T$ is uniformly bounded over all positive diagonal matrices $R$.

In an earlier version of the present manuscript, we had explicitly
characterized the uniform upper bound as $\| X^\dagger y \|_2$ for
$X^{\dagger}$ the Moore-Penrose pseudoinverse of $X$, which is not true in
general. Existence of an uniform upper bound has been also established by
<cit.> using block matrix inversions and induction.

Fix any $n \geq 1$. We first consider the case $p=1$, such that $X=x_1 \in \mathbb{R}^n$, $R=r_1$ for some scalar $r_1>0$. We obtain

\begin{equation} \label{eq:one_dim_bound}
\| (X^TX + R)^{-1} X^T \|_2 = \frac{1}{r_1+x_1^T x_1} \| x_1^T \|_2 = \frac{\| x_1 \|_2}{r_1+ \|x_1\|_2^2} \leq
\frac{\mathrm{I}_{x_1 \neq 0}}{\|x\|_2},
\end{equation}

which is an uniform bound for all $r_1>0$.

We now consider $p > 1$. We have $X=[x_1 \ldots x_p]$ with $x_i \in \mathbb{R}^n$ for $i=1, \ldots, p$ and
$R=\text{Diag}(r_1,\ldots,r_p)$ with $r_i >0$ for $i=1, \ldots, p$. For any
$p\times p$ permutation matrix $P$, the matrices $XP^T$ and $PRP^T$ correspond to a reordering of the columns of $X$ and the diagonal entries of $R$. We

\begin{equation}
\| \big( (XP^T)^T(XP^T) + PRP^T \big)^{-1} (XP^T)^T\|_2 = \| P(X^TX + R)^{-1}X^T \|_2 = \| (X^TX + R)^{-1}X^T \|_2.
\end{equation}

This implies that $\| (X^TX + R)^{-1}X^T \|_2$ is invariant under any
reordering of the $p$ columns of $X$
and the $p$ diagonal entries of $R$. By considering such reordering, we can assume the diagonal entries of
$R$ are non-decreasing without loss of generality; that is $R=\text{Diag}(r_1,\ldots,r_p)$ with $0<r_1\leq \ldots\leq r_p$.
Denote $A_p :=(X^TX + R)^{-1}X^T$. Then,

A_p    = ( R^-1 - R^-1 X^T ( I_n + X R^-1 X^T )^-1 X R^-1 ) X^T  by Lemma <ref>

= R^-1 X^T - R^-1 X^T ( I_n - ( I_n + X R^-1 X^T )^-1 )

= R^-1 X^T (I_n + X R^-1 X^T)^-1.

It therefore suffices to show that

\begin{equation} \label{eq:Cp_bound}
C_{p} := \sup_{0<r_1 \leq \ldots \leq r_{p}} \| A_p \|_2 = \sup_{0<r_1 \leq \ldots \leq r_{p}} \|{R}^{-1} {X}^T (I_n+ {X} {R}^{-1} {X}^T)^{-1}\|_2
\end{equation}

is finite for any $p > 1$. For $k=1, \ldots, p-1$, define

\begin{equation} \label{eq:Ck_bound}
C_{k} := \sup_{0<r_1 \leq \ldots \leq r_{k}} \|A_{k}\|_2 = \sup_{0<r_1 \leq \ldots \leq r_{k}} \|{R}_k^{-1} {X}_k^T (I_n+ {X}_k {R}_k^{-1} {X}_k^T)^{-1}\|_2,
\end{equation}

where $X_k=[x_1 \ldots x_k]$ is the $n\times k$ matrix corresponding to the first $k$ columns of $X$
and ${R}_k = \text{Diag}(r_1, \ldots, r_k)$ with $0<r_1\leq \ldots \leq r_k$
is the $k\times k$ diagonal matrix corresponding to the top $k$ diagonal entries of $R$,
and $A_{k} := R_{k}^{-1} X_{k}^T (I_n + X_{k} R_{k}^{-1} X_{k}^T)^{-1}$.
Then $X_{p-1} R_{p-1}^{-1} X_{p-1}^T = \sum_{i=1}^{p-1} x_i x_i^T/r_i$.
By Lemma <ref> with $U=V^T=x_{p} / r_p^{1/2}$,

(I_n+ X R^-1 X^T)^-1    = ( (I_n+ X_p-1 R_p-1^-1 X_p-1^T) + x_p x_p^T/r_p )^-1

=(I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1

- (I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1 x_p x_p^T (I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1 /r_p + x_p^T (I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1 x_p.

Substituting this decomposition into (<ref>), we obtain

A_p =
R_p-1^-1 X_p-1^T

r_p^-1 x_p^T
(I_n+ X R^-1 X^T)^-1 =

where, after some simplifications, $v_p \in \mathbb{R}^n$ and $V_{p-1}$ is a $p-1$ by $n$ matrix such that

v_p^T    := r_p^-1 x_p^T(I_n+ X R^-1 X^T)^-1 =
x_p^T (I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1/ r_p + x_p^T (I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1 x_p ,

V_p-1    := R_p-1^-1 X_p-1^T (I_n+ X R^-1 X^T)^-1 = A_p-1 ( I_n- x_p v_p^T ).

We have $\|V_{p-1}\|_2 \leq C_{p-1} (1 + \|x_p\|_2 \|v_p^{T}\|_2 )$ by (<ref>), subadditivity
and submultiplicativity. Then by Lemma <ref>,

A_p _2    ≤V_p-1_2 + v_p^T_2 ≤C_p-1 (1 + x_p_2 v_p^T_2 ) + v_p^T_2.

Therefore to show $C_p$ in (<ref>) is finite, it suffices to bound $\|v_p^{T}\|_2$ uniformly over all $R=\text{Diag}(r_1,\ldots,r_p)$ with
$0<r_1\leq \ldots\leq r_p$. Note that $x_{p} = X_{p-1} a_{p} + b_{p}$ for some $a_{p} \in \mathbb{R}^{p-1}$ and some
$b_{p} \in \mathbb{R}^n$ such that $X_{p-1}^T b_{p} = 0$. Here
$X_{p-1} a_{p}$ corresponds to the projection of $x_{p}$ to the column space of $X_{p-1}$, and
$b_{p}$ corresponds to the component of $x_{p}$ that is perpendicular to the column space of $X_{p-1}$.
Note that $b_{p} = (I_n+ X_{p-1} R_{p-1}^{-1} X_{p-1}^T)b_{p}$. This implies

(I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1 b_p    = b_p,

b_p^T (I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1 X_p-1 a_p     = b_p^T X_p-1 a_p = 0 .

Therefore by (<ref>),

v_p^T _2    =  a_p^T X_p-1^T (I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1 + b_p^T/ r_p + a_p^T X_p-1^T (I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1 X_p-1 a_p + b_p^T b_p  _2

≤ a_p^T X_p-1^T (I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1/ r_p + a_p^T X_p-1^T (I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1 X_p-1 a_p _2 +
b_p^T/ r_p + b_p^T b_p _2   by subadditivity

≤ a_p^T X_p-1^T (I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1/ r_p  _2 + I_b_p ≠0/b_p_2   as  (I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1  is p.s.d.

=  (R_p-1 a_p)^T R_p-1^-1 X_p-1^T(I_n+ X_p-1 R_p-1^-1 X_p-1^T)^-1/ r_p  _2 + I_b_p ≠0/b_p_2

≤ R_p-1/r_p _2  a_p _2 C_p-1
+ I_b_p ≠0/b_p_2   by submultiplicativity and (<ref>)

≤ a_p _2 C_p-1 + I_b_p ≠0/b_p_2
as $r_{p} \geq r_j$ for all $j = 1,\ldots,p-1$.

Note that $x_p, a_p$ and $b_p$  do not depend on entries of the diagonal matrix $R$.
Taking the supremum over all $R=\text{Diag}(r_1,\ldots,r_p)$ with
$0<r_1\leq \ldots\leq r_p$ in (<ref>), by (<ref>) we obtain the recurrence relation

C_p ≤C_p-1 (1 + x_p_2 (  a_p _2 C_p-1 + I_b_p ≠0/b_p_2 ) ) + (  a_p _2 C_p-1 + I_b_p ≠0/b_p_2 )

for all $p \geq 1$. As $C_1$ is finite by (<ref>),
(<ref>) implies that $C_{p}$ is finite for all $p$.

(Moments of the full conditionals of each $\eta_j$).
For fixed constants $m>0$ and $\nu>0$, define a distribution on $(0, \infty)$
with probability density function

\begin{equation} \label{eq:eta_full_conditional_pdf}
p(\eta | m) \propto \frac{1}{\eta^{\frac{1-\nu}{2}} (1+ \nu \eta)^{\frac{\nu+1}{2}}} e^{-m \eta}.
\end{equation}

Then for any $c > \max \{-\frac{\nu+1}{2}, -1 \}$ and $\eta \sim p(\cdot | m) $, we have

𝔼[η^c | m] =1/ν^c  Γ(1+ν/2+c)/ Γ(1+ν/2) U(1+ν/2+c, 1+c, m/ν)/U(1+ν/2, 1, m/ν).

Recall the definition of the confluent hypergeometric function of the second kind $U(a,b,z) := \frac{1}{\Gamma(a)} \int_0^\infty x^{a-1}(1+x)^{b-a-1} e^{-zx} dx $. We obtain

∫_0^∞η^c/η^1-ν/2 (1+ νη)^ν+1/2 e^-m η d η   = ν^-ν+1/2-c ∫_0^∞x^ν+1/2+c-1(1+x)^-ν+1/2 e^- m/ν x dx

= ν^-ν+1/2-c Γ(ν+1/2+c ) U (ν+1/2+c, 1+c, m/ν ).

𝔼[η^c | m]    =  ν^-ν+1/2-c Γ(ν+1/2+c ) U (ν+1/2+c, 1+c, m/ν )/ν^-ν+1/2 Γ(ν+1/2 ) U (ν+1/2, 1, m/ν )

= 1/ν^c Γ(ν+1/2+c ) U (ν+1/2+c, 1+c, m/ν )/Γ(ν+1/2 ) U (ν+1/2, 1, m/ν ).

Ratios of confluent hypergeometric function evaluations, such as the ones
arising in Lemma <ref>,
can be upper-bounded using the following result.

Let $U$ denote the confluent hypergeometric function of the second kind, and fix $r>0$.

* Let $c>0$. Then

\begin{equation} \label{eq:confluent_ratio_bound1}
\frac{U \big(r+c, 1+c, x \big)}{ U \big(r, 1, x \big)} < x^{-c} \quad \text{for all } x>0.
\end{equation}

Furthermore, for every $\epsilon > 0$, there exists a $K^{(1)}_{\epsilon, c, r} < \infty$ such that

\begin{equation} \label{eq:confluent_ratio_bound2}
\frac{U \big(r+c, 1+c, x \big)}{ U \big(r, 1, x \big)} < \epsilon x^{-c} + K^{(1)}_{\epsilon, c, r} \quad \text{for all } x>0.
\end{equation}

* Let $\max \{ -1, -r \} <c<0$. Then there exists some $K^{(2)}_{c, r} < \infty$ such that

\begin{equation} \label{eq:confluent_ratio_bound3}
\frac{ U \big(r+c, 1+c, x \big)}{U \big(r, 1, x \big)} < x^{-c} + K^{(2)}_{c, r} \quad \text{for all } x>0.
\end{equation}

Case: $c>0$.
By definition, we have

x^c  Γ(r+c) U (r+c, 1+c, x )/ Γ(r) U (r, 1, x )     = x^c ∫_0^∞t^r+c-1 (1+t)^-r e^-xt dt/∫_0^∞t^r-1 (1+t)^-r e^-xt dt

∫_0^∞s^r+c-1 (x+s)^-r e^-s ds/∫_0^∞s^r-1 (x+s)^-r e^-s ds  for  s=xt

= Γ(r+c)/Γ(r) 𝔼[ (x+S)^-r ]/𝔼[ (x+S̃)^-r ]  for  S ∼Gamma(r+c, 1)  and  S̃ ∼Gamma(r, 1)

< Γ(r+c)/Γ(r)

where the final equality follows as the distribution $\Gamma(r+c, 1)$ stochastically dominates the distribution $\Gamma(r, 1)$, and hence
$\mathbb{E}[ (x+S)^{-r} ] < \mathbb{E}[ (x+\tilde{S})^{-r} ]$ for all $x>0$ as $r>0$ and $c>0$. This gives $U \big(r+c, 1+c, x \big) / U \big(r, 1, x \big) < x^{c}$ as required for (<ref>).

We now prove (<ref>). By definition, we have

Γ(r+c) U (r+c, 1+c, x )/ Γ(r) U (r, 1, x )  = ∫_0^∞t^r+c-1 (1+t)^-r e^-xt dt/∫_0^∞t^r-1 (1+t)^-r e^-xt dt = 𝔼[T_x^c ] = x^-c 𝔼[Y_x^c ]

where $T_x$ is a random variable with density $p_{T_x}(t) \propto t^{r-1}(1+t)^{-r} e^{-xt}$ on $[0, \infty)$, and $Y_x:= xT_x$ is a random variable with density $p_{Y_x}(t) \propto y^{r-1}(x+y)^{-r} e^{-y}$ on $[0, \infty)$.

* Consider $\mathbb{E}[Y_x^c ]$ as $x \rightarrow 0^+$. We obtain

lim_x →0^+ 𝔼[Y_x^c ] = lim_x →0^+  ∫_0^∞y^c + r-1(x+y)^-r e^-y dy  / ∫_0^∞y^r-1(x+y)^-r e^-y dy  = Γ(c) /  lim_x →0^+ ∫_0^∞( y/x+y )^r  y^-1 e^-y dy  = 0.

* Consider $\mathbb{E}[T_x^c ]$ as a function of $x$ for $x >0$. From the density of $T_x$, we note that for all $0 < x_1 < x_2$, $T_{x_1}$ stochastically dominates $T_{x_2}$, and as $c >0$, $T^c_{x_1}$ also stochastically dominates $T^c_{x_2}$. Therefore, $\mathbb{E}[T_x^c ]$ is a decreasing function of $x$ for $x > 0$.

As $\lim_{x \rightarrow 0^+} \mathbb{E}[Y_x^c ]=0$, for any $\epsilon>0$, there exists $\delta>0$ such that $\mathbb{E}[Y_x^c ] < \frac{\Gamma(r+c) }{\Gamma(r)} \epsilon$ for all $x \leq \delta$. This gives $x^{-c}\mathbb{E}[Y_x^c ] < \frac{\Gamma(r+c) }{\Gamma(r)} \epsilon x^{-c}$ for all $x \leq \delta$. Also for all $x > \delta$, we have $x^{-c}\mathbb{E}[Y_x^c ] = \mathbb{E}[T_x^c ] \leq \mathbb{E}[T_\delta^c ] < \infty $ as $\mathbb{E}[T_x^c ]$ is a decreasing function of $x$ for $x > 0$. Overall, we obtain

Γ(r+c) U (r+c, 1+c, x )/ Γ(r) U (r, 1, x )     = x^-c 𝔼[Y_x^c ] < Γ(r+c) /Γ(r) ϵx^-c + 𝔼[T_δ^c ]

⇒U (r+c, 1+c, x )/ U (r, 1, x )    < ϵx^-c + K^(1)_ϵ, c, r

for $K^{(1)}_{\epsilon, c, r} = \frac{\Gamma(r)}{\Gamma(r+c)} \mathbb{E}[T_\delta^c ] = \frac{ U \big(r+c, 1+c, \delta \big)}{U \big(r, 1, \delta \big)} < \infty $ as required.

Case: $\max \{ -1, -r \} <c<0$. Note that $U \big(r+c, 1+c, x \big)$ and $U \big(r, 1, x \big)$ are well-defined as $r, r+c, c+1>0$.

* Consider $\frac{ U \big(r+c, 1+c, x \big)}{U \big(r, 1, x \big)} - x^{-c}$ as $x \rightarrow \infty$. We first recall the asymptotic expansion of the confluent hypergeometric function of the second kind <cit.> at infinity,

\begin{equation} \label{eq:asympt_U_expansion}
U(a,b,z) = z^{-a} - a(a-b+1)z^{-(a+1)} + \mathcal{O}( z^{-(a+2)} )
\end{equation}

as $z \rightarrow \infty$ for any $a,b > 0$.
For large $x > 0$, this gives

x^c+1 (  U (r+c, 1+c, x )/U (r, 1, x ) - x^-c )    = x^c+1 (  x^r U (r+c, 1+c, x ) - x^r-c U (r, 1, x ) /x^r U (r, 1, x ) )

=x (  x^r+c  U (r+c, 1+c, x ) - x^r U (r, 1, x ) /x^r U (r, 1, x ) )

=x ( (1 - (r+c)r/x) - (1 - r^2/x ) /1 + 𝒪(x^-2) )

= -rc + 𝒪(x^-1)

where (<ref>) follows from the asymptotic expansion in (<ref>). As $c+1>0$, we obtain

lim_x →∞ (  U (r+c, 1+c, x )/U (r, 1, x ) - x^-c ) = 0 .

* Consider $\frac{ U \big(r+c, 1+c, x \big)}{U \big(r, 1, x \big)} - x^{-c}$ as $x \rightarrow 0^+$. By definition of confluent hypergeometric function of the second kind, we have

lim_x →0^+  Γ(r+c) U (r+c, 1+c, x )/ Γ(r) U (r, 1, x )  = lim_x →0^+ ∫_0^∞t^r+c-1(1+t)^-r e^-xt dt/∫_0^∞t^r-1(1+t)^-r e^-xt dt.

Note that as $r+c >0$ and $c-1 < -1$,

\begin{equation}
\int_0^\infty t^{r+c-1}(1+t)^{-r} dt = \int_0^\infty \big(\frac{t}{1+t}\big)^r t^{c-1} dt \leq \int_0^1 t^{r+c-1} dt + \int_1^\infty t^{c-1} dt = \frac{1}{r+c} - \frac{1}{c} < \infty .
\end{equation}

Therefore, $\lim_{x \rightarrow 0^+} \int_0^\infty t^{r+c-1}(1+t)^{-r} e^{-xt}dt = \frac{1}{r+c} - \frac{1}{c}  < \infty$. However, we have

\begin{equation}
\lim_{x \rightarrow 0^+} \int_0^\infty t^{r-1}(1+t)^{-r} e^{-xt} dt  = \lim_{x \rightarrow 0^+} \int_0^\infty \big(\frac{t}{1+t}\big)^r t^{-1} e^{-xt} dt  = \infty .
\end{equation}

Therefore $\lim_{x \rightarrow 0^+} \frac{ \Gamma(r+c) U \big(r+c, 1+c, x \big)}{ \Gamma(r) U \big(r, 1, x \big)}  = 0$. As $c < 0$, we obtain

\begin{equation} \label{eq:hyper_geom_asymp2}
\lim_{x \rightarrow 0^+} \frac{ U \big(r+c, 1+c, x \big)}{U \big(r, 1, x \big)} - x^{-c} = 0.
\end{equation}

From (<ref>) and (<ref>), we have $ \lim_{x \rightarrow \infty} \frac{ U \big(r+c, 1+c, x \big)}{U \big(r, 1, x \big)} - x^{-c} = 0$ and $ \lim_{x \rightarrow 0^+} \frac{ U \big(r+c, 1+c, x \big)}{U \big(r, 1, x \big)} - x^{-c} = 0$. As $x \mapsto \frac{ U \big(r+c, 1+c, x \big)}{U \big(r, 1, x \big)} - x^{-c}$ is a continuous function on $(0, \infty)$, this implies that $\frac{ U \big(r+c, 1+c, x \big)}{U \big(r, 1, x \big)} - x^{-c}$ is bounded above for $x \in (0, \infty)$. In particular, there exists some $K^{(2)}_{c, r}  < \infty$ such that

\begin{equation}
\frac{ U \big(r+c, 1+c, x \big)}{U \big(r, 1, x \big)} - x^{-c} < K^{(2)}_{c, r}
\end{equation}

for all $x > 0$ as required for (<ref>).

(Bound on moments of the full conditionals of each $\eta_j$).
Consider a random variable $\eta$ on $(0, \infty)$ with density as in (<ref>)
in Lemma <ref> for some $\nu >0$.

Consider when $c>0$. There for every $\epsilon > 0$,
there exists some $K^{(1)}_{\epsilon, c, \nu} < \infty$ such that

𝔼[η^c | m] ≤ϵm^-c + K^(1)_ϵ, c, ν   for all  m>0.

Consider when $\max \{-1, - \frac{1+\nu}{2} \} < c < 0$. Then, there exists some
$K^{(2)}_{c, \nu} < \infty$ such that

𝔼[η^c | m] ≤Γ(1+ν/2+c)/ Γ(1+ν/2) m^-c + K^(2)_c, ν
for all  m>0.

By Lemma <ref>, for all $\max \{-1, - \frac{1+\nu}{2} \} < c$,

𝔼[η^c | m] =1/ν^c  Γ(1+ν/2+c)/ Γ(1+ν/2) U(1+ν/2+c, 1+c, m/ν)/U(1+ν/2, 1, m/ν).

Consider when $c >0$. We apply Lemma <ref> with $r = \frac{1+\nu}{2}$, $x=\frac{m}{\nu}$, and $\tilde{\epsilon}=\epsilon \Big( \frac{\Gamma(\frac{1+\nu}{2}+c)}{ \Gamma(\frac{1+\nu}{2})} \Big)^{-1}$. Then there exists some $ \tilde{K}^{(1)}_{\epsilon, c, \nu} < \infty$ such that

𝔼[η^c | m] ≤1/ν^c Γ(1+ν/2+c)/ Γ(1+ν/2) (ϵ̃ (m/ν )^-c + K̃^(1)_ϵ, c, ν ) = ϵm^-c + K^(1)_ϵ, c, ν

for $K^{(1)}_{{\epsilon}, c, \nu} = \frac{1}{\nu^{c}}  \frac{\Gamma(\frac{1+\nu}{2}+c)}{ \Gamma(\frac{1+\nu}{2})} \tilde{K}^{(1)}_{\epsilon, c, \nu} < \infty $.

Consider when $\max \{-1, - \frac{1+\nu}{2} \} < c < 0$. By Lemma <ref> with $r = \frac{1+\nu}{2}$ and $x=\frac{m}{\nu}$, there exists some $\tilde{K}^{(2)}_{c, \nu} < \infty$ such that

𝔼[η^c | m] ≤1/ν^c  Γ(1+ν/2+c)/ Γ(1+ν/2) ( ( m/ν )^-c + K̃^(2)_c, ν) = Γ(1+ν/2+c)/ Γ(1+ν/2) m^-c + K^(2)_c, ν

for $K^{(2)}_{c, \nu} = \frac{1}{\nu^{c}}  \frac{\Gamma(\frac{1+\nu}{2}+c)}{ \Gamma(\frac{1+\nu}{2})} \tilde{K}^{(2)}_{c, \nu} < \infty $.

Let $X \sim \mathcal{N}(\mu, \sigma^2)$ on $\mathbb{R}$. Then, for $c \in (0,1)$,

\begin{equation}
\mathbb{E}[ |X|^{-c} ] \leq \frac{1}{\sigma^{c} 2^{c/2}} \frac{\Gamma(\frac{1-c}{2})}{\sqrt{\pi}}
\end{equation}

We follow the proof in <cit.>, which is
included here for completeness. Recall the inverse moments of the univariate Normal distribution,

\begin{equation}
\mathbb{E}[ |X|^{-c} ]  = \frac{1}{\sigma^{c} 2^{c/2}} \frac{\Gamma(\frac{1-c}{2})}{\sqrt{\pi}} \exp( - \frac{\mu^2}{2 \sigma^2}) M \Big( \frac{1-c}{2}, \frac{1}{2}, \frac{\mu^2}{2 \sigma^2} \Big)
\end{equation}

where $M(a, b , z) := {}_{1}{F}_{1}(a ; b ; z) $ is the confluent hypergeometric function of the first kind.

For $a<b$, we know $M(a, b , z) = \frac{\Gamma(b)}{\Gamma(a) \Gamma(b-a)} \int_{0}^{1} e^{zt} t^{a-1}(1-t)^{b-a-1}dt$.
For $a=\frac{1-c}{2} < b = \frac{1}{2}$ and $z=\frac{\mu^2}{2 \sigma^2}$, this gives

exp(-z) M(a, b, z)    = Γ(b)/Γ(a) Γ(b-a) ∫_0^1 e^-z(1-t) t^a-1(1-t)^b-a-1dt

= Γ(b)/Γ(a) Γ(b-a) ∫_0^1 e^-zt (1-t)^a-1 t^b-a-1 dt.

Therefore, $\exp(-z) M(a, b , z)$ is decreasing for $z\geq 0$ and $\exp(-z) M(a, b , z) \leq \exp(-0) M(a, b , 0)=1$ for all $z\geq 0$. This gives

\begin{equation}
\mathbb{E}[ |X|^{-c} ] = \frac{1}{\sigma^{c} 2^{c/2}} \frac{\Gamma(\frac{1-c}{2})}{\sqrt{\pi}} \exp( - \frac{\mu^2}{2 \sigma^2}) M \Big( \frac{1-c}{2}, \frac{1}{2}, \frac{\mu^2}{2 \sigma^2} \Big) \leq \frac{1}{\sigma^{c} 2^{c/2}} \frac{\Gamma(\frac{1-c}{2})}{\sqrt{\pi}}.
\end{equation}

(A uniform bound on the full conditional mean of $1/\sigma^2$).
Consider the full conditional distribution of $\sigma^2$ given $\eta, \xi$ in Algorithm <ref>. Given $\eta \in \mathbb{R}_{>0}^p, \xi>0$, we have

\begin{equation}
\sigma^2 | \xi, \eta \sim \InvGamma \bigg( \frac{a_0+n}{2}, \frac{y^T M_{\xi, \eta}^{-1} y + b_0}{2} \bigg),
\end{equation}

where $M_{\xi, \eta} = I_n + \xi^{-1} X \text{Diag}(\eta^{-1}) X^T$. Then

\begin{equation} \label{eq:uniform_bound_mean_sigma2}
\mathbb{E} \Big[ \frac{1}{\sigma^2} \Big| \xi, \eta \Big] = \frac{a_0+n}{y^T M_{\xi, \eta}^{-1} y + b_0} < \frac{a_0+n}{b_0} < \infty.
\end{equation}

Recall that the mean of a $\mathrm{Gamma}(a,b)$ distribution (under shape and rate parameterization) is $a/b$. Equation (<ref>) now directly follows as $M_{\xi, \eta}^{-1}$ is positive semi-definite.

(Maximal coupling probability of each component $\eta_j$).
For fixed constants $m>0$ and $\nu>0$, define a distribution on $(0, \infty)$
with probability density function

\begin{equation} \label{eq:eta_full_conditional_pdf_1}
p(\eta | m) \propto \frac{1}{\eta^{\frac{1-\nu}{2}} (1+ \nu \eta)^{\frac{\nu+1}{2}}} e^{-m \eta}.
\end{equation}

Consider random variables $\eta \sim p(\cdot | m) $, and $\tilde{\eta} \sim p(\cdot | \tilde{m})$
for some fixed $\nu>0$. Fix some $0<a<b<\infty$, and take $m, \tilde{m} \in (a,b)$.
Then under a maximal coupling,

\begin{equation} \label{eq:eta_coupling_bound}
\mathbb{P}( \eta = \tilde{\eta} | m, \tilde{m} ) =
\int_0^\infty \min \{ p(\eta^*|m), p(\eta^*|\tilde{m}) \} d\eta^* \geq
U \Big( \frac{1+\nu}{2}, 1, \frac{b}{\nu}\Big) \Big/ U \Big( \frac{1+\nu}{2}, 1, \frac{a}{\nu}\Big).
\end{equation}

Fix some $0<a<b<\infty$, and take $m, \tilde{m} \in (a,b)$. Then

ℙ( η= η̃ | m, m̃ )     =
∫_0^∞min{ p(η^*|m), p(η^*|m̃) } dη^* , as  (η, η̃) | m, m̃  are maximally coupled

= ∫1/(η^*)^1-ν/2 (1+ νη^*)^ν+1/2  min{ e^-m η^*/U ( 1+ν/2, 1, m/ν) , e^-m̃ η^*/U (1+ν/2, 1, m̃/ν)  } d η^*.

Recall $m, \tilde{m} \in (a, b) $. Note also that $x \mapsto U \big(\frac{1+\nu}{2}, 1, x\big)$ is a decreasing function for $x>0$. This gives

\begin{equation}
\min \Big\{ \frac{e^{-m \eta^*}}{U \big( \frac{1+\nu}{2}, 1, \frac{m}{\nu}\big)} , \frac{e^{-\tilde{m} \eta^*}}{U \big(\frac{1+\nu}{2}, 1, \frac{\tilde{m}}{\nu}\big)}  \Big\}
> \frac{e^{-b \eta^*}}{U \big( \frac{1+\nu}{2}, 1, \frac{a}{\nu}\big)}.
\end{equation}

Therefore for all $m, \tilde{m} \in (a,b)$,

ℙ( η= η̃ | m, m̃ )
> ∫1/(η^*)^1-ν/2 (1+ νη^*)^ν+1/2  e^-b η^*/U ( 1+ν/2, 1, a/ν) d η^*
= U ( 1+ν/2, 1, b/ν)/U ( 1+ν/2, 1, a/ν).

Equation (<ref>) follows from taking an infimum
over all $m, \tilde{m} \in (a,b)$.

Drift and Minorization conditions.

For ease of notation, denote $\Sigma_{t+1}:= X^TX + \xi_{t+1} \text{Diag}(\eta_{t+1})$ and let $e_j$ be the $j^{th}$ unit vector of the standard basis on $\mathbb{R}^p$. We first consider $\sum_{j=1}^p m_j^{d}$. Note that

𝔼 [ ∑_j=1^p m_t+1,j^d | η_t+1, ξ_t+1, σ_t+1^2 ]    = ( ξ_t+1/2 σ_t+1^2 )^d ∑_j=1^p  𝔼 [ | β_t+1,j |^2d | η_t+1, ξ_t+1, σ_t+1^2 ]

≤( ξ_t+1/2 σ_t+1^2 )^d ∑_j=1^p  𝔼 [ β^2_t+1,j | η_t+1 , ξ_t+1, σ_t+1^2 ]^d

= ( ξ_t+1/2 σ_t+1^2 )^d ∑_j=1^p  ( (e_j^T Σ_t+1^-1 X^Ty)^2 + σ_t+1^2 (e_j^T Σ_t+1^-1 e_j) )^d

≤( ξ_t+1/2 σ_t+1^2 )^d ∑_j=1^p  ( |e_j^T Σ_t+1^-1 X^Ty|^2d + σ_t+1^2d (e_j^T Σ_t+1^-1 e_j)^d )

≤( ξ_t+1/2 σ_t+1^2 )^d ( ∑_j=1^p (|e_j^T Σ_t+1^-1 X^Ty|^2d)^1/d )^d p^1-d

+ σ_t+1^2d ( ξ_t+1/2 σ_t+1^2 )^d ∑_j=1^p (e_j^T ( ξ_t+1
Diag(η_t+1) )^-1 e_j)^d

= ( ξ_t+1/2 σ_t+1^2 )^d ( Σ_t+1^-1 X^Ty ^2_2 )^d p^1-d + 1/2^d ∑_j=1^p 1/η_t+1,j^d

≤( ξ_t+1/2 σ_t+1^2 )^d ( C_p^2  y_2^2 )^d p^1-d + 1/2^d ∑_j=1^p 1/η_t+1,j^d,

where (<ref>) follows from Jensen's inequality, (<ref>) follows as $(x+y)^d \leq x^d + y^d$ for $x, x \geq 0$ and $d\in (0,1)$, (<ref>) follows from Hölder's inequality and
from $e_j^T \Sigma_{t+1}^{-1}e_j \leq e_j^T\big( \xi_{t+1} \text{Diag}(\eta_{t+1}) \big)^{-1} e_j$ as $X^TX$ is positive-definite matrix, and (<ref>) follows from Lemma <ref> for some finite constant $C_p $ which does not depend on $\eta_{t+1}$ or $\xi_{t+1}$.

As the prior $\pi_{\xi}(\cdot)$ on $\xi$ has compact support, we have $b \leq \xi_{t+1} \leq B$ for all $t \geq 0$ for some $0<b \leq B < \infty$. This upper bound on $\xi_{t+1}$, Lemma <ref> and Jensen's inequality gives

𝔼 [( ξ_t+1/2 σ_t+1^2 )^d | η_t+1 ] ≤B^d/2^d 𝔼 [ ( 1/σ_t+1^2 ) | η_t+1 ]^d ≤B^d/2^d (a_0+n/b_0)^d = ( B (a_0+n)/2 b_0 )^d .

Therefore, by (<ref>) and (<ref>), we obtain

𝔼 [ ∑_j=1^p m_t+1,j^d | η_t+1 ] = 𝔼 [ 𝔼 [ ∑_j=1^p m_t+1,j^d | η_t+1, ξ_t+1, σ_t+1^2 ] | η_t+1 ] ≤1/2^d ∑_j=1^p 1/η_t+1,j^d + K_1

for $K_1 := \Big( \frac{B (a_0+n)}{2 b_0} \Big)^d \Big( C_p^{2} \| y\|_2^2 \Big)^d p^{1-d} < \infty$. By (<ref>) and Corollary <ref>, we obtain

𝔼 [ ∑_j=1^p m_t+1,j^d | β_t, ξ_t, σ_t^2 ]    = 𝔼 [ 𝔼 [ ∑_j=1^p m_t+1,j^d | η_t+1 ] | β_t, ξ_t, σ_t^2 ]

≤1/2^d ∑_j=1^p 𝔼 [ η_t+1,j^-d | β_t, ξ_t, σ_t^2 ] + K_1

≤1/2^d ( ∑_j=1^p Γ(1+ν/2-d)/ Γ(1+ν/2)  m_t,j^d + K^(2)_d, ν ) + K_1

≤1/2^d Γ(1+ν/2-d)/ Γ(1+ν/2) ∑_j=1^p   m_t,j^d + K_2

for $K_2 := \frac{p}{2^{d}} K^{(2)}_{d, \nu} + K_1 < \infty$. Note that for every $\nu \geq 1$, there exists $d \in (0,1)$
such that $\frac{1}{2^{d}} \frac{\Gamma(\frac{\nu+1}{2}-d)}{\Gamma(\frac{\nu+1}{2})} < 1$. We can verify this by considering the
function $f(d) := \frac{1}{2^{d}} \frac{\Gamma(\frac{\nu+1}{2}-d)}{\Gamma(\frac{\nu+1}{2})}$. It suffices to note that
$f(0) = 1$ and $f'(0)= -\psi(\frac{\nu+1}{2}) - \log(2) < 0$ for $\nu \geq 1$, where $\psi(x):= \frac{\Gamma'(x)}{\Gamma(x)}$
is the Digamma function.

We now consider $\sum_{j=1}^p m_j^{-c}$ for $c\in (0,1/2)$. We obtain

𝔼 [ ∑_j=1^p     m_t+1,j^-c | σ^2_t+1, ξ_t+1, η_t+1 ] = ( ξ_t+1/2 σ_t+1^2 )^-c ∑_j=1^p  𝔼 [ | β_t+1,j |^-2c | σ^2_t+1, ξ_t+1, η_t+1  ]

≤( ξ_t+1/2 σ_t+1^2 )^-c 1/2^c  Γ(1-2c/2)/√(π)  ∑_j=1^p ( σ_t+1^2 e_j^T Σ_t+1^-1 e_j )^-c

≤( ξ_t+1/2 σ_t+1^2 )^-c 1/2^c  Γ(1-2c/2)/√(π)  ∑_j=1^p ( σ_t+1^2 e_j^T (X_2^2 I_p + ξ_t+1 Diag(η_t+1))^-1 e_j )^-c

= ( ξ_t+1/2 σ_t+1^2 )^-c 1/2^c  Γ(1-2c/2)/√(π)  ∑_j=1^p (X_2^2 + ξ_t+1 η_t+1,j)^c/σ_t+1^2c

≤Γ(1/2-c)/√(π)  ∑_j=1^p ( (X_2^2/ξ_t+1 )^c + η_t+1,j^c )

≤Γ(1/2-c)/√(π)  ∑_j=1^p η_t+1,j^c + p Γ(1/2-c)/√(π) (X_2^2/b )^c

for some $b>0$. Here (<ref>) follows from Lemma <ref>, (<ref>) follows as $e_j^T ( \|X\|_2^2 I_p + \xi_{t} \text{Diag}(\eta_t) )^{-1} e_j \leq e_j^T \Sigma_{t+1}^{-1} e_j$ for $\| \cdot \|_2$ the $L_2$ operator norm of a matrix, (<ref>) follows as $(x+y)^c \leq x^c + y^c$ for all $x,y>0$ and $0< c <1/2$, and (<ref>) follows as the prior on $\xi$ has compact support, so $\xi_{t+1} \geq b$ for some $b>0$. Therefore, by (<ref>),

𝔼 [ ∑_j=1^p m_t+1,j^-c | η_t+1 ] = 𝔼 [ 𝔼 [ ∑_j=1^p m_t+1,j^-c | σ^2_t+1, ξ_t+1, η_t+1 ] | η_t+1 ] ≤Γ(1/2-c)/√(π)  ∑_j=1^p η_t+1,j^c + K_3.

for $K_3 := p \frac{\Gamma(\frac{1}{2}-c)}{\sqrt{\pi}} \Big(\frac{\|X\|_2^{2}}{b} \Big)^c < \infty$. By (<ref>) and Corollary <ref> with $r=\frac{1+\nu}{2}$ and any $\epsilon < \frac{\sqrt{\pi}}{\Gamma(\frac{1}{2}-c)}$, there exists $K^{(1)}_{\epsilon, c, \nu} < \infty$ such that

𝔼 [ ∑_j=1^p m_t+1,j^-c | σ^2_t, ξ_t, β_t  ]    = 𝔼 [  𝔼 [ ∑_j=1^p m_t+1,j^-c | η_t+1 ] | σ^2_t, ξ_t, β_t  ]

≤Γ(1/2-c)/√(π) ∑_j=1^p 𝔼 [ η_t+1,j^c | σ^2_t, ξ_t, β_t  ] + K_3

≤Γ(1/2-c)/√(π) ∑_j=1^p ( ϵm_t,j^-c +K^(1)_ϵ, c, ν  ) + K_3

= ( Γ(1/2-c)/√(π) ϵ) ∑_j=1^p m_t,j^-c  + K_4,

where $\frac{\Gamma(\frac{1}{2}-c)}{\sqrt{\pi}}  \epsilon < 1$ and
$K_4 = p \frac{\Gamma(\frac{1}{2}-c)}{\sqrt{\pi}}K^{(1)}_{\epsilon, c, \nu} + K_3$.

We can now combine (<ref>) and (<ref>) together. Then, for some fixed $d \in (0,1)$ such that $\frac{1}{2^{d}} \frac{\Gamma(\frac{\nu+1}{2}-d)}{\Gamma(\frac{\nu+1}{2})} < 1$ and any $c \in (0,1/2)$, we obtain

𝔼 [ ∑_j=1^p m_t+1,j^-c + m_t+1,j^d | m_t ] ≤γ_drift ( ∑_j=1^p m_t,j^-c + m_t,j^d ) + K_drift

for some

0 < γ_drift    := max {Γ(1/2-c)/√(π)  ϵ, 1/2^d Γ(ν+1/2-d)/Γ(ν+1/2) } <1

0 < K_drift    := K_2 + K_4 < ∞.

Therefore, for any such $d \in (0,1)$ and $c\in (0,1/2)$,

\begin{equation}
V(\beta) = \sum_{j=1}^p m_j^{-c} + m_j^{d}
\end{equation}

for $m_j = \frac{\xi \beta_j^2}{2 \sigma^2}$ is a Lyapunov function.

(Minorization condition).
For $R>0$, let
$S(R) = \{ (\beta, \xi, \sigma^2) \in \mathbb{R}^p \times \mathbb{R}_{>0} \times \mathbb{R}_{>0} : V(\beta, \xi, \sigma) < R \}$
denote the sub-level sets of the Lyapunov function in Proposition
<ref>. Let $\mathcal{P}$ denote the Markov transition kernel
associated with the update from $(\beta_t,\xi_t, \sigma_t^2)$ to $(\beta_{t+1},\xi_{t+1}, \sigma_{t+1}^2)$
implied by the update rule in (<ref>).
Let $Z = (\beta, \xi, \sigma^2)$ and $\tilde{Z} = (\tilde{\beta}, \tilde{\xi}, \tilde{\sigma}^2)$.
Then for all $R > 0$,

there exists $\epsilon \in (0,1)$ such that

\begin{equation} \label{eq:minorization}
\text{TV} \Big( \mathcal{P} ( Z, \cdot ), \mathcal{P} (\tilde{Z}, \cdot) \Big) < 1 - \epsilon \quad \text{ for all } Z, \tilde{Z} \in S(R),
\end{equation}

where $\text{TV}$ stands for the total variation distance between two probability distributions.
In particular for $R > 1$, $\epsilon = ( U \big( \frac{1+\nu}{2}, 1, \frac{R^{1/d}}{\nu}\big) / U \big( \frac{1+\nu}{2}, 1, \frac{R^{-1/c}}{\nu}\big))^p$

Denote $Z = (\beta, \xi, \sigma^2)$ and $\tilde{Z} =
(\tilde{\beta}, \tilde{\xi}, \tilde{\sigma}^2)$, and let $Z, \tilde{Z} \in S(R)$
for some $R>1$. We obtain

TV     ( 𝒫(Z, ·), 𝒫(Z̃, ·) )
≤ℙ ( Z_1 ≠Z̃_1 | Z, Z̃ ) < 1 - ϵ,

where $(Z_{1}, \tilde{Z}_{1})$ corresponds to one-scale coupling of
$\mathcal{P}(Z, \cdot)$ and $\mathcal{P}(\tilde{Z}, \cdot)$ as in
Proposition <ref>, the first inequality follows from the
coupling inequality, and the second inequality follows from the proof of Proposition
<ref> for
$\epsilon = \Big( U \big( \frac{1+\nu}{2}, 1, \frac{R^{1/d}}{\nu}\big) \Big/ U \big( \frac{1+\nu}{2}, 1, \frac{R^{-1/c}}{\nu}\big) \Big)^p \in (0,1)$
as $R>1$. Note that $\epsilon$ does not depend on $Z$ and $\tilde{Z}$
as required.

Geometric ergodicity follows from Proposition
<ref> and Proposition <ref> with any $R > \max \{1,2K/(1-\gamma)\}$

§.§ Coupled Kernel Proofs

For $R>0$, let
$S(R) := \{ (\beta, \xi, \sigma^2) \in \mathbb{R}^p \times \mathbb{R}_{>0} \times \mathbb{R}_{>0} : V(\beta, \xi, \sigma) < R \}$
denote the sub-level sets of the Lyapunov function in Proposition <ref>, with
corresponding contraction constant $\gamma \in (0,1)$ and additive constant $K>0$.
Let $ Z_t := (\beta_t, \xi_t, \sigma_t^2)$ and $\tilde{Z}_t = (\tilde{\beta}_t, \tilde{\xi}_t, \tilde{\sigma}_t^2)$
for all $t \geq 0$. By <cit.>, it suffices to prove that
for any $R>0$ sufficiently large such that $\gamma + K/(1+R) < 1$, there exists
$\epsilon \in (0,1)$ such that

\begin{equation} \label{eq:one_scale_condition}
\inf_{Z_t, \tilde{Z}_t \in S(R)}
\mathbb{P} \big( Z_{t+1} = \tilde{Z}_{t+1} \big| Z_{t}, \tilde{Z}_{t} \big) > \epsilon,
\end{equation}

under the one-scale coupling. We now prove (<ref>) directly for the
one-scale coupling in Proposition <ref>.
Fix some $R>1$ such that $\gamma + K/(1+R) < 1$, and take any $Z_t, \tilde{Z}_t \in S(R)$.

ℙ     ( Z_t+1 = Z̃_t+1 | Z_t, Z̃_t)

= 𝔼 [ ℙ ( Z_t+1 = Z̃_t+1 | η_t+1, η̃_t+1, Z_t, Z̃_t )
| Z_t, Z̃_t ]

= 𝔼 [ ℙ ( Z_t+1 = Z̃_t+1 | η_t+1, η̃_t+1 )
| Z_t, Z̃_t ]

= 𝔼 [ 1_ { η_t+1 = η̃_t+1 }  | Z_t, Z̃_t ]
, as  (Z_t+1,Z̃_t+1) | (η_t+1, η̃_t+1)  follows CRN coupling

= ℙ ( η_t+1 = η̃_t+1 | Z_t, Z̃_t )

= ∏_j=1^p ℙ ( η_t+1,j = η̃_t+1,j | Z_t, Z̃_t )
, as  (η_t+1, η̃_t+1) | Z_t, Z̃_t  is maximally coupled componentwise

= ∏_j=1^p   ℙ ( η_t+1,j = η̃_t+1,j | m_j, m̃_j )
, for  m_j = ξ_t β_t,j^2/2σ_t^2  and
m̃_j =  ξ̃_t β̃_t,j^2/2 σ̃_t^2.

As $Z, \tilde{Z} \in S(R)$ for some $R>1$, $m_j, \tilde{m}_j \in (R^{-1/c}, R^{1/d}) $ for all $j=1,\ldots, p$.
By Lemma <ref>,

ℙ ( η_t+1,j = η̃_t+1,j | m_j, m̃_j ) ≥U ( 1+ν/2, 1, R^1/d/ν) / U ( 1+ν/2, 1, R^-1/c/ν)

for all $j=1,\ldots, p$. Therefore,

ℙ     ( Z_t+1 = Z̃_t+1 | Z_t, Z̃_t)
≥∏_j=1^p U ( 1+ν/2, 1, R^1/d/ν)/U ( 1+ν/2, 1, R^-1/c/ν)
= ϵ
for $\epsilon := ( U ( \frac{1+\nu}{2}, 1, \frac{R^{1/d}}{\nu}) / U ( \frac{1+\nu}{2}, 1, \frac{R^{-1/c}}{\nu}))^p \in (0,1)$ as $R>1$.

As $\eta_{t+1}$ given $m_t$ and $\tilde{\eta}_{t+1}$ given $\tilde{m}_t$ are
conditionally independent component-wise, it suffices to consider
the univariate case where $p=1$. Henceforth we assume $p=1$
and we drop the component subscripts for ease of notation.
Without loss of generality, take $\tilde{m}_t \leq m_t$
and $h:(0, \infty)$ to be monotone increasing.
Then the densities of $\eta_{t+1}|m_t$ and $\tilde{\eta}_{t+1}|\tilde{m}_t$
on $(0, \infty)$ are given by

p(η| m_t) ∝e^-m_t η/η^1-ν/2 (1+νη)^ν+1/2
p(η̃ | m̃_t) ∝e^-m̃_t η/η^1-ν/2 (1+νη)^ν+1/2

respectively. As $\tilde{m}_t \leq m_t$, the distribution of
$\tilde{\eta}_{t+1}$ stochastically dominates the distribution of $\eta_{t+1}$.
Therefore under a CRN coupling of ${\eta}_{t+1}$ and $\tilde{\eta}_{t+1}$
(specifically a CRN coupling of Algorithm <ref>), $\tilde{\eta}_{t+1} \geq \eta_{t+1}$
almost surely, and $h(\tilde{\eta}_{t+1}) \geq h(\eta_{t+1})$ almost surely.
We obtain

\begin{equation}
\mathbb{E}_{ \text{crn} } \big[ | h(\eta_{t+1}) - h(\tilde{\eta}_{t+1}) | \big| m_{t}, \tilde{m}_{t} \big] =
\mathbb{E} \big[h(\tilde{\eta}_{t+1}) - h(\eta_{t+1}) \big| m_{t}, \tilde{m}_{t} \big].
\end{equation}

Now consider a maximal coupling of $\eta_{t+1}$ and $\tilde{\eta}_{t+1}$ with independent residuals
<cit.>. Let $q=\text{TV}(p(\cdot | m_t) , p(\cdot | \tilde{m}_t))$.
Then $\eta_{t+1}=\tilde{\eta}_{t+1}$ with probability
$1-q$. Otherwise with probability $q$,
$(\eta_{t+1}, \tilde{\eta}_{t+1}) = (\eta^*_{t+1}, \tilde{\eta}^*_{t+1})$
where $\eta^*_{t+1}$ and $\tilde{\eta}^*_{t+1}$ are independent random variables
with densities

p(η^*_t+1 | m_t+1, m̃_t+1)    =
max{ p(η^*_t+1 | m_t) - p(η^*_t+1 | m̃_t), 0 }/q  and

p(η̃^*_t+1 | m_t+1, m̃_t+1)    =
max{ p(η̃^*_t+1 | m̃_t) - p(η̃^*_t+1 | m_t), 0 }/q

respectively on $(0, \infty)$.
This implies that the supports of $\eta^*_{t+1}$ and $\tilde{\eta}^*_{t+1}$ are disjoint.

By (<ref>), we further have $\eta^*_{t+1} \leq \tilde{\eta}^*_{t+1}$ almost surely.
Overall, we obtain $\tilde{\eta}_{t+1} \geq \eta_{t+1}$ almost surely under a maximal
coupling, and hence

\begin{equation}
\mathbb{E}_{ \text{max} } \big[ | h(\eta_{t+1}) - h(\tilde{\eta}_{t+1}) | \big| m_{t}, \tilde{m}_{t} \big] =
\mathbb{E} \big[ h(\tilde{\eta}_{t+1}) - h(\eta_{t+1}) \big| m_{t}, \tilde{m}_{t} \big] =
\mathbb{E}_{ \text{crn} } \big[ | h(\eta_{t+1}) - h(\tilde{\eta}_{t+1}) | \big| m_{t}, \tilde{m}_{t} \big]
\end{equation}

as required.

CRN coupling result. To prove (<ref>), it suffices to show that under CRN,

\begin{equation} \label{eq:crn_as_bound}
| \log (1+\eta_{t+1,j}) - \log (1+\tilde{\eta}_{t+1,j}) | \leq | \log m_{t, j} - \log \tilde{m}_{t, j} |,
\end{equation}

almost surely for all $j=1, \ldots, p$. Henceforth, to prove (<ref>)
we work component-wise and drop the subscripts for ease of notation.
For the Horseshoe prior, we then have $p(\eta | m) \propto \exp(-m\eta)/(1+\eta)$
on $(0, \infty)$. The cumulative density function is then given by
$F_m(t) = 1 - E_1(m(t+1))/E_1(m)$ where $E_1$ is the exponential integral function. Therefore,
under CRN coupling, we obtain

1 + η= 1/m E_1^-1(U^* E_1(m))  and
1 + η̃ = 1/m̃ E_1^-1(U^* E_1(m̃)),

for $U^* \sim \Uniform(0,1)$. Then

d/d logm ( log(1 + η) )    = -1 + d/d logm ( log( E_1^-1(U^* E_1(m)) ) )

= -1 + m/ E_1^-1(U^* E_1(m))  ( d/d m ( E_1^-1(U^* E_1(m)) ) )

= -1 + m/ E_1^-1(U^* E_1(m))  U^* E_1'(m) /  E_1'( E_1^-1(U^* E_1(m)) )
, where  E'_1(z) = -e^-z/z

= -1 + m/ m(1+η)  U^* E_1'(m) /  E_1'( m(1+η) )
, as  m(1+η) = E_1^-1(U^* E_1(m))

=  -1 + U^* e^m η , as  E'_1(z) = -e^-z/z

= -1 +E_1(m(1+η)) e^m (1+η)/E_1(m) e^m , as  U^* = E_1(m(1+η))/E_1(m)

∈(-1,0)  as  x ↦e^x E_1(x)  is positive and decreasing on  (0, ∞).

Therefore, as $|\frac{d \log(1 + \eta)}{d \log m}| \leq 1$ and by the mean value theorem we obtain

\begin{equation}
| \log(1 + \eta) - \log(1 + \tilde{\eta}) | \leq | \log m - \log \tilde{m} |,
\end{equation}

almost surely as required for (<ref>). Finally, (<ref>)
implies (<ref>) with $m_{t, j} = \mu$ and $\tilde{m}_{t, j} = \tilde{\mu}$ for all $j = 1, \ldots, p$.

Maximal coupling result. To prove (<ref>), take $\tilde{\mu} < \mu$
without loss of generality. For $c>0$, denote the event

$$E_c := \big\{ \max_{j=1}^p | \log(1+\eta_{t+1,j}) - \log(1+\tilde{\eta}_{t+1,j}) | > c \big\},$$

where $(\eta_{t+1,j}, \tilde{\eta}_{t+1,j})$ are i.i.d. random variables jointly distributed
according to a maximal coupling.

By component-wise independence,

ℙ(E_c | μ, μ̃) =
1 - ∏_j=1^p ℙ ( | log(1+η_t+1,j) - log(1+η̃_t+1,j) | ≤c | μ, μ̃ )
= 1 - x_c^p

for $x_c=\mathbb{P} \big( | \log(1+\eta_{t+1,1}) - \log(1+\tilde{\eta}_{t+1,1}) | \leq c \big| \mu, \tilde{\mu} \big)$.
Henceforth, we focus on this quantity $x_c$ and drop the subscripts for ease of notation.
For the Horseshoe prior, the marginal densities densities of $\eta|\mu$ and $\tilde{\eta}|\tilde{\mu}$
are given by

p(η| μ) = e^-μη/1+η 1/e^μE_1(μ)  and
p(η| μ̃) = e^-μ̃ η/1+η 1/e^μ̃ E_1(μ̃)
on  (0, ∞).

Consider a maximal coupling of $\eta$ give $\mu$ and $\tilde{\eta}$ given $\tilde{\mu}$
with independent residuals <cit.>. Let $q=\text{TV}(p(\cdot | \mu) , p(\cdot | \tilde{\mu}))$.
Then $\eta=\tilde{\eta}$ with probability $1-q$. Otherwise with probability $q$,
$(\eta, \tilde{\eta}) = (\eta^*, \tilde{\eta}^*)$ where $\eta^*$ and $\tilde{\eta}^*$
are independent random variables on $(0, \infty)$ with densities

p(η^* | μ, μ̃)    = max{ p(η^* | μ) - p(η^* | μ̃), 0 }/q  and

p(η̃^* | μ, μ̃)    = max{ p(η̃^* | μ̃) - p(η̃^* | μ), 0 }/q,

respectively. Note that $p(\eta^* | \mu) \geq p(\eta^* | \tilde{\mu})$ if and only if  $1 + \eta^* \leq K$ for constant
$K:=(\log E_1(\mu) - \log E_1(\tilde{\mu}))/(\tilde{\mu}-\mu)$. Therefore, $1+\eta^*$ has support $[1, K)$
and $1+\tilde{\eta}^*$ has support $[K, \infty)$.
We obtain,

x_c    = (1-q) + q ℙ ( | log(1+η^*) - log(1+η̃^*) | ≤c )

= 1 - q ℙ ( | log(1+η^*) - log(1+η̃^*) | > c )

= 1 - q ℙ ( log(1+η̃^*) - log(1+η^*) > c )  , as  η^* ≤η̃^*

= 1 - q ℙ ( 1+η̃^* > (1+η^*) e^c )

≤1 - q ℙ ( 1+η̃^* > t )  for  t=K e^c , as  1+η̃^* ≤K  almost surely

= 1 - ∫_t^∞1/1+η̃^*  ( e^-μ̃η̃^*/e^μ̃ E_1(μ̃) -
e^-μη̃^*/ e^μ E_1(μ) ) dy , by  (<ref>)

= 1 - ∫_t^∞1/y ( e^-μ̃y/E_1(μ̃) -
e^-μy/E_1(μ) ) dy  with  y=1+η̃^*

= 1  - ( E_1(μ̃ t)/E_1(μ̃) - E_1(μt)/E_1(μ)).

We now consider this upper bound on $x_p$ as $t \rightarrow \infty$.
Note that for $\mu, \tilde{\mu}$ fixed with $\tilde{\mu}<\mu$,

lim_t →∞ ( E_1(μ̃) μ̃t e^μ̃t ( E_1(μ̃t)/E_1(μ̃) - E_1(μt)/E_1(μ) ) )
= lim_t →∞ ( μ̃t e^μ̃t E_1(μ̃t) ) lim_t →∞ ( 1 - E_1(μt) E_1(μ̃)/E_1(μ̃ t) E_1(μ) )

= lim_t →∞ ( μ̃t e^μ̃t E_1(μ̃t) )  as  lim_t →∞ E_1(μt)/E_1(μ̃ t) = 0

= 1  as  lim_x →∞ ( x e^x E_1(x) ) = 1.

Therefore, there exists some constant $D_1>0$ such that

\begin{equation} \label{eq:prob_limit_bound}
\Big( \frac{E_1(\tilde{\mu }t)}{E_1(\tilde{\mu })} - \frac{E_1(\mu t)}{E_1(\mu )} \Big) > \frac{D_1}{ \tilde{\mu }t e^{\tilde{\mu }t}}
\end{equation}

for all $t>0$. By (<ref>) and (<ref>), we obtain

x_c^p    < (1 - D_1/t e^μ̃t)^p

= ( (1 - 1/D_1^-1 μ̃t e^μ̃t)^D_1^-1 μ̃t e^μ̃t )^p/(D_1^-1 μ̃t e^μ̃t)

< exp(-p/(D_1^-1 μ̃t e^μ̃t) )  as   (1-1/x)^x < 1/e  for  x ≥1.

For $p > 1$, now take $c= \log \big( \frac{\alpha \log p}{\tilde{\mu } K} \big)$ for any fixed constant $\alpha \in (0,1)$.
Recall $t=K e^c$, such that $\tilde{\mu }t = \alpha \log p$. Then

x_c^p    < exp(-p/(D_1^-1 μ̃t e^μ̃t) ) = exp( -D_1 p^1-α/(αlogp) )
< exp( -D p^1-α' )

for any $\alpha' \in (\alpha, 1)$ and some constant $D$ which does not depend on $p$.

ℙ (E_c |  m_t, m̃_t )    > 1 - exp( -D p^1-α' )

for $c = \log (\alpha/(\tilde{\mu } K)) + \log \log p$.
Recall that $K:=(\log E_1(\mu ) - \log E_1(\tilde{\mu }))/(\tilde{\mu }-\mu )$.
We obtain

1/μ̃ K    = μ-μ̃/μ̃ (logE_1(μ̃) - logE_1(μ))

= μ/μ̃-1/logE_1(μ̃) - logE_1(μ)

≥logμ- logμ̃ /logE_1(μ̃) - logE_1(μ)  as  x-1≥log(x)  for  x>0

= e^μ̂ E_1(μ̂)  for some  μ̂ ∈[μ̃, μ]  as  d/d (logμ) (-logE_1(μ)) = 1/e^μ E_1(μ)

≥e^μ E_1(μ)  as  x ↦e^x E_1(x)  is decreasing for all  x>0

≥1/2 log(1 + 1/μ)  as  e^x E_1(x) ≥1/2 log(1+1/x)  for all  x>0

= 1/2 log(1 + 1/max{μ, μ̃ } )

=: L.

Therefore, $c \geq \log (\alpha L) + \log \log p$.
By (<ref>),

ℙ (max_j=1^p | log(1+η_t+1,j) - log(1+η̃_t+1,j) |
> log(αL) + loglogp   |  m_t, m̃_t )

> 1 - exp( -D p^1-α').

The case $\tilde{\mu } > \mu $ follows by symmetry to give (<ref>).

Fix some $\nu \geq 1$. For any $m > 0$, let $p(\cdot | m)$ denote the component-wise
full conditional density of $\eta$ given $m$. This is given by

\begin{equation} %\label{eq:eta_full_conditional_pdf_2}
p(\eta | m) = \frac{\exp(-m \eta)}{\eta^{\frac{1-\nu}{2}} (1+ \nu \eta)^{\frac{\nu+1}{2}}} \frac{1}{Z_m}
\end{equation}

on $(0, \infty)$, where $Z_m$ is the normalization constant. We will first show

\begin{equation} \label{eq:tv_metric_upper_bound_componentwise}
\text{TV} \big( p(\cdot | m), p(\cdot | \tilde{m}) \big)^2 \leq \frac{(1+\nu)(e-1)}{4} | \log(m) - \log(\tilde{m}) |.
\end{equation}

By a generalization of Pinkser's inequality <cit.>, we obtain

\begin{equation}
\text{TV} \big( p(\cdot | m), p(\cdot | \tilde{m}) \big)^2 \leq \frac{1}{2 \alpha} D_{\alpha}( p(\cdot | m) \| p(\cdot | \tilde{m}) )
\end{equation}

for all $\alpha \in (0,1]$, where

$$D_{\alpha}(p(\cdot | m) \| p(\cdot | \tilde{m})) := \frac{1}{\alpha -1}
\log \Big( \int p(\eta | m)^\alpha p(\eta | \tilde{m})^{1-\alpha} d \eta \Big)$$

is the Rényi divergence between probability densities $p$ and $q$ with respect to the
Lebesgue measure. For $\alpha=1/2$, we obtain

D_α(p(·| m)  p(·| m̃))    = -2 log(
∫_0^∞exp(-m + m̃/2 η)/η^1-ν/2 (1+ νη)^ν+1/2 1/√(Z_m Z_m̃) dη)
= log( Z_m Z_m̃/Z_(m+m̃)/2^2  )

Let $L(m)= \log Z_m$. Then for $\eta | m \sim p(\cdot | m)$, note that

d L(m)/dm = d Z_m/dm/Z_m = - 𝔼[η| m] ≤0,  and
d^2 L(m)/dm^2 = d^2 Z_m/dm^2 Z_m - (d Z_m/dm)^2/Z_m^2 =
𝔼[η^2 | m] - 𝔼[η| m]^2 ≥0,

such that $L(m)$ is a decreasing convex function. For $\tilde{m} \geq m$, this gives

L(m) + L(m̃) - 2 L ( m+m̃/2 )
≤( L'(m̃)-L'(m) ) (m̃-m/2 )

= ( 𝔼[η| m]-𝔼[η| m̃] ) (m̃-m/2).

Furthermore, note that by Lemmas <ref> and <ref>,

𝔼[η| m] =1/ν  Γ(1+ν/2+1)/ Γ(1+ν/2) U(1+ν/2+1, 1+1, m/ν)/U(1+ν/2, 1, m/ν) ≤1/ν  Γ(1+ν/2+1)/ Γ(1+ν/2) ( m/ν )^-1 = 1+ν/2 m^-1

where $U(a,b,z)$ corresponds to the confluent hypergeometric function of
the second kind.
By (<ref>), we therefore obtain

D_α(p(·| m)  p(·| m̃))
≤(𝔼[η| m]-𝔼[η| m̃] ) (m̃-m/2)

≤1+ν/2 (m̃-m/2 m )

= 1+ν/4  ( exp( logm̃ - logm ) - 1) ),

where $ \exp( |\log \tilde{m} - \log m |) - 1 \leq (e-1) |\log m - \log \tilde{m}|$
for $|\log m - \log \tilde{m}| \leq 1$. This gives

\begin{equation} \label{eq:renyi_div}
\text{TV} \big( p(\cdot | m), p(\cdot | \tilde{m}) \big)^2 \leq
\frac{(1+\nu)(e-1)}{4} \big( \log(\tilde{m}) - \log(m) \big),
\end{equation}

for $\tilde{m} \geq m$. The case $\tilde{m} \leq m$ follows by symmetry,
giving (<ref>).

We now extend (<ref>) to the
multivariate setting. Let $\pi(\eta | m):= \prod_{j=1}^p p(\eta_j| m_j)$ now denote the
full conditional density of $\eta \in (0, \infty)^p$ given $m \in (0, \infty)^p$. Then we obtain

TV ( π(·| m), π(·| m̃) )^2
≤D_1/2( π(·| m)  π(·| m̃) )

= ∑_j=1^p D_1/2( p(·| m_j)  p(·| m̃_j) )  by conditional independence

= (1+ν)(e-1)/4  ∑_j=1^p | log(m_j) - log(m̃_j) |  by (<ref>).

Finally, let $Z = (\beta, \xi, \sigma^2)$ and $\tilde{Z} = (\tilde{\beta}, \tilde{\xi}, \tilde{\sigma}^2)$, and
$m_{j} = \xi \beta_{j}^2/(2 \sigma^2)$ and $\tilde{m}_{j} = \tilde{\xi} \tilde{\beta}_{j}^2/(2 \tilde{\sigma}^2)$
for $j=1, \ldots, p$. For any such $Z, \tilde{Z}$, we obtain

TV ( 𝒫(Z, ·), 𝒫(Z̃ | ·) )
= 1/2 ∫| p( η^* | Z) - p(η^*|Z̃) | p( Z^* | η^* ) dη^* dZ^*

= 1/2 ∫| p( η^* | Z) - p(η^* | Z̃) | dη^*

= 1/2 ∫| π( η^* | m) - π(η^* | m̃) | dη^*

= TV ( π(·| m), π(·| m̃) ).

Therefore, by (<ref>), we obtain

TV ( 𝒫(Z, ·), 𝒫(Z̃ | ·) )^2
≤(1+ν)(e-1)/4  ∑_j=1^p | log(m_j) - log(m̃_j)|.

This proposition follows from a modification of
the proof of
Proposition <ref>. Following the setup and notation
in the proof of Proposition <ref>, we have

ℙ     ( Z_t+1 = Z̃_t+1 | Z_t, Z̃_t)

= 𝔼 [ ℙ ( Z_t+1 = Z̃_t+1 | η_t+1, η̃_t+1, Z_t, Z̃_t )
| Z_t, Z̃_t ]

= 𝔼 [ ℙ ( Z_t+1 = Z̃_t+1 | η_t+1, η̃_t+1 )
| Z_t, Z̃_t ]

= 𝔼 [ 1_ { η_t+1 = η̃_t+1 }  | Z_t, Z̃_t ]
, as  (Z_t+1,Z̃_t+1) | (η_t+1, η̃_t+1)  follows common random numbers coupling

= ℙ ( η_t+1 = η̃_t+1 | Z_t, Z̃_t )

= ℙ ( η_1,t+1 = η̃_1,t+1 | Z_t, Z̃_t )
∏_j=2^p ℙ ( η_t+1,j = η̃_t+1,j |
1 ≤i ≤j-1 ⋂ { η_i,t+1 = η̃_i,t+1 }, Z_t, Z̃_t )

under the switch-to-CRN coupling

≥∏_j=1^p U ( 1+ν/2, 1, R^1/d/ν)/U ( 1+ν/2, 1, R^-1/c/ν)
by Lemma <ref>, as in the proof of Proposition <ref> for  Z_t, Z̃_t ∈S(R),

= ( U ( 1+ν/2, 1, R^1/d/ν)/U ( 1+ν/2, 1, R^-1/c/ν) )^p =: ϵ,

where $\epsilon \in (0,1)$ since $R > 1$.

§ ADDITIONAL ALGORITHMS

$ C_t := (\beta_t, \eta_t, \sigma_t^2, \xi_t) \in \mathbb{R}^p \times \mathbb{R}^p_{> 0} \times \mathbb{R}_{> 0} \times \mathbb{R}_{> 0}$.

* Sample $\eta_{t+1} | \beta_{t}, \sigma^2_{t}, \xi_{t}$ component-wise
independently using Algorithm <ref>, targeting

\begin{equation} \label{eq:eta_full_conditionals}
\pi(\eta_{t+1} | \beta_{t} , \sigma^2_{t}, \xi_{t}) \propto \prod_{j=1}^p
\pi\big( \eta_{t+1,j} | m_{t,j} := \frac{\xi_{t} \beta_{t,j}^2}{2  \sigma^2_{t}} \big) \propto \prod_{j=1}^p \frac{e^{-m_{t,j} \eta_{t+1,j}}}{\eta_{t+1,j}^{\frac{1-\nu}{2}}(1+\nu \eta_{t+1,j})^{\frac{\nu + 1}{2}}} .
\end{equation}

* Sample $\xi_{t+1},\sigma^2_{t+1}, \beta_{t+1}$ given $\eta_{t+1}$ as follows:

* Sample $\xi_{t+1}$ using Metropolis–Rosenbluth–Teller–Hastings with step-size $\sigma_{\text{MRTH}}$:
Given $\xi_{t}$, propose $\log(\xi^*) \sim \mathcal{N}(\log(\xi_{t}), \sigma^2_{\text{MRTH}})$.
Calculate acceptance probability

\[q =  \frac{L(y | \xi_*, \eta_{t+1} ) \pi_{\xi}(\xi_*)  }{L(y | \xi_{t}, \eta_{t+1} ) \pi_{\xi}(\xi_{t})} \frac{\xi^*}{\xi_{t}},\]

where $\pi_{\xi}(\cdot) $ is the prior for $\xi$,
$ M_{\xi,\eta} := I_n + \xi^{-1} X\, \text{Diag}(\eta^{-1})\, X^T $ and

\begin{equation} \log(L(y | \xi, \eta )) =  - \frac{1}{2} \log (|M_{\xi,\eta}|) - \frac{a_0 + n}{2} \log( b_0 + y^T M_{\xi,\eta}^{-1} y ). \label{eq:Half_t_Exact_algo_ll}
\end{equation}

Set $\xi_{t+1} := \xi^*$ with probability $\min(1,q)$, otherwise set $\xi_{t+1} := \xi_{t}$.

* Sample $\sigma^2_{t+1}$ given $\xi_{t+1},\eta_{t+1}$:

\begin{equation}
\sigma^2_{t+1} | \xi_{t+1}, \eta_{t+1} \sim \InvGamma\bigg(\frac{a_0+n}{2}, \frac{y^T M_{\xi_{t+1},\eta_{t+1}}^{-1} y + b_0}{2}\bigg).
\end{equation}

* Sample $\beta_{t+1}$ given $\sigma^2_{t+1}, \xi_{t+1}, \eta_{t+1}$ using Algorithm <ref> in Appendix <ref>, targeting

\begin{equation}
\beta_{t+1} | \sigma^2_{t+1}, \xi_{t+1}, \eta_{t+1} \sim \mathcal{N} \big( \Sigma^{-1} X^T y,  \sigma^2_{t+1} \Sigma^{-1} \big)
\text{ for } \Sigma = X^TX + \xi_{t+1} \text{Diag}(\eta_{t+1} ).
\end{equation}

$ C_{t+1} := (\beta_{t+1}, \eta_{t+1}, \sigma^2_{t+1}, \xi_{t+1})$.
MCMC for Bayesian shrinkage regression with Half-$t(\nu)$ priors.

Slice sampler targeting $p(\eta_j | m_j) \propto (\eta_j^{\frac{1-\nu}{2}}(1+\nu \eta_j)^{\frac{\nu + 1}{2}})^{-1} e^{-m_j \eta_j}$ on $(0, \infty)$.
$m_j > 0$, state $\eta_{t,j} > 0$.

* Sample $U_{t,j} | \eta_{t,j} \sim \Uniform(0, (1+\nu \eta_{t,j})^{-\frac{\nu + 1}{2}})$.

* Sample $\eta_{t+1,j} | U_{t,j} \sim P_j$, where $P_j$
has unnormalized density $\eta \mapsto \eta^{s-1} e^{ - m_j \eta}$ on
$(0, T_{t,j})$, with $T_{t,j}=(U_{t,j}^{-2/(1+\nu)}-1)/\nu$ and $s=(1+\nu)/2$.