$C_{t+1} \defeq (\beta_{t+1}, \eta_{t+1}, \sigma^2_{t+1}, \xi_{t+1})$ and

$\tilde{C}_{t+1} \defeq (\tilde{\beta}_{t+1}, \tilde{\eta}_{t+1}, \tilde{\sigma}_{t+1}^2, \tilde{\xi}_{t+1})$.

§.§ Approximate MCMC for high-dimensional logistic regression

Section <ref> considers Bayesian logistic regression
with spike and slab priors applied to a
a malware detection dataset and a lymph node GWAS dataset.
The Malware detection dataset from the UCI machine learning repository [Dua and Graff, 2017] has $n=373$ observations (corresponding to a binary response vector indicating whether a file is malicious or non-malicious) and $d=503$ covariates (corresponding to features of the files), and
is publicly available on <kaggle.com/piyushrumao/malware-executable-detection>.
The lymph node GWAS dataset [Hans et al., 2007, Liang et al., 2013, Narisetty et al., 2019] has $n=148$ observations (corresponding to a binary response
vector indicating high or low risk status of the lymph node that is related to
breast cancer) and $d=4514$ covariates (corresponding to SNPs in the genome)
is not publicly available.

The logistic regression likelihood is given by
$L(\beta; y, X) = \prod_{i=1}^n (1+\exp(-y_i x_i^\top \beta))^{-1}$ where
$y \in \{-1,1\}^n$ is the response vector,
$X \in \mathbb{R}^{n \times d}$ is the scaled design matrix with rows $x_i^\top$,
and $\beta \in \mathbb{R}^d$ is an unknown signal vector.
The spike and slab prior is given by

\begin{equation} \label{eq:spike_slab}
Z_j \overset{i.i.d.}{\sim} \Bernoulli(q), \qquad \beta_j | Z_j=0 \sim \mathcal{N}(0, \tau^2_{0}),
\qquad \beta_j | Z_j=1 \sim \mathcal{N}(0, \tau^2_{1})
\end{equation}

for $j=1, \mydots, d$ where $q \in (0,1), \tau_0>0,$ and $\tau_1>0$ are hyper-parameters
with $\tau_0 \ll \tau_1$ such that $Z_i=0$ and $Z_i=1$ correspond to null and
non-null components of $\beta_j$ respectively.
By considering the posterior distribution of each variable $Z_j$ on
$\{0,1\}$, spike and slab priors provide an interpretable method for
Bayesian variable selection <cit.>.

The state-of-the-art exact MCMC algorithms to sample from posteriors
corresponding to the prior in (<ref>)
are Gibbs samplers which cost $\mathcal{O}(n^2 d)$ per iteration
[Bhattacharya et al., 2016].
Narisetty et al., 2019 have recently developed approximate MCMC methods for
this setting. Their approximate MCMC algorithm, called Skinny Gibbs,
is based on matrix approximations of the form

\begin{equation}
\begin{pmatrix}
X_A^\top X_A + \tau_{1}^{-2}I & X_A^\top X_{A^c} \\
X_{A^c}^\top X_A & X_{A^c}^\top X_{A^c} + \tau_{0}^{-2}I \\
\end{pmatrix} \approx \begin{pmatrix}
X_A^\top X_A + \tau_{1}^{-2}I & 0 \\
0 & ((n-1) + \tau_{0}^{-2})I \\
\end{pmatrix}
\end{equation}

where $A= \{j: Z_{j}=1\}$, $X_A$ is an $n \times |A|$ matrix corresponding to the active columns
$j \in A$ of the design matrix, and $X_{A^c}$ is an $n \times (d-|A|)$ matrix

corresponding to the inactive columns $j \notin A$.

This gives an overall computation cost of $\mathcal{O}(n \min \{d, |A|^2 \})$ per iteration.

In Section <ref>, we use couplings to assess the
quality of the Skinny Gibbs algorithm.
The upper bounds in Figure <ref> are given by our estimator $\cub[2]$
(<ref>) with $S=1000,  T=3000,$ and $I=100$ for both the malware and lymph node GWAS datasets, where these values were chosen based on initial runs.
The lower bounds in Figure <ref> are estimated using (<ref>)
based on the same samples from the coupled chains used to
calculate the upper bound estimate.
We consider a CRN coupling between one marginal chain corresponding to the exact MCMC kernel
and another corresponding to the Skinny Gibbs kernel.
The CRN coupled kernel is given in Algorithm <ref>.

exact chain current state $C_t \defeq (\beta_t, z_t, e_t, w_t) \in \mathbb{R}^d \times \{0,1\}^d \times \mathbb{R}^n \times \mathbb{R}^n$ and

approximate chain current state $\tilde{C}_t \defeq (\tilde{\beta}_t, \tilde{z}_t, \tilde{e}_t, \tilde{w}_t) \in \mathbb{R}^d \times \{0,1\}^d \times \mathbb{R}^n \times \mathbb{R}^n$.

* Sample $(\beta_{t+1}, \tilde{\beta}_{t+1}) | z_t, e_t, w_t, \tilde{z}_t, \tilde{e}_t, \tilde{w}_t$
with common random numbers and the fast sampling algorithms of Bhattacharya et al., 2016, marginally targeting

(β_A,t+1, β_A^c,t+1) | z_t,  e_t, w_t ∼𝒩 ( Σ^-1 X^⊤W y,  Σ^-1 )     for

X_A^⊤W X_A + τ_1^-2I     X_A^⊤W X_A^c

X_A^c^⊤W X_A     X_A^c^⊤W X_A^c + τ_0^-2I

(β̃_Ã, t+1, β̃_Ĩ, t+1) | z̃_t,  ẽ_t, w̃_t ∼𝒩 ( Σ̃^-1 X^⊤W̃ y,  Σ̃^-1 )
X_Ã^⊤W̃ X_Ã + τ_1^-2 I     0

0     ((n-1) + τ_0^-2) I

where $W = \Diag( w_t )$ and $\tilde{W} = \Diag( \tilde{w}_t )$,  $A= \{j: z_{j,t}=1\}$ and $\tilde{A}= \{j: \tilde{z}_{j,t}=1\}$ are the index sets of active components,
$X_A$ and $X_{\tilde{A}}$ are matrices corresponding to the active (or inactive) columns of $X$
with columns $j \in A$ and $j \in \tilde{A}$ respectively,
$\beta_{A, t+1}$ and $\tilde{\beta}_{\tilde{A}, t+1}$ are vectors of active components of $\beta_{t+1}$ and $\tilde{\beta}_{t+1}$ respectively.

* Sample $(z_{t+1}, \tilde{z}_{t+1})$ given $\beta_{t+1}, \tilde{\beta}_{t+1}, e_t, \tilde{e}_t, w_t,  \tilde{w}_t$ with common random numbers sequentially in order for $j=1, \mydots, p$ such that each $z_{j,t+1}$ and $\tilde{z}_{j,t+1}$ are Bernoulli random variables with odds

q 𝒩(β_j,t+1, 0, τ_1^2)/(1-q) 𝒩(β_j,t+1, 0, τ_0^2)  and

q 𝒩(β̃_j,t+1, 0, τ_1^2)/(1-q) 𝒩(β̃_j,t+1, 0, τ_0^2)

exp( β̃_j,t+1 X_j^⊤W̃ ( Y - X_C_j β_C_j, t+1 ) + 1/2X_j^⊤(I-W̃) X_j β_j,t+1^2 )

respectively where $\mathcal{N}(\cdot; \mu, \Sigma)$ is the probability density of the normal distribution with mean $\mu$ and variance $\Sigma$,  $C_j \defeq \{k:  \tilde{z}_{k,t+1}=1 \text{ for } k <j \text{ or } \tilde{z}_{k,t}=1 \text{ for } k > j \}$ is the index set of active components in $\{1, \mydots, p\} \backslash \{j\}$,
$X_{C_j}$ is a matrix of the columns of $X$ which correspond to indices in $C_j$,  and $\tilde{\beta}_{C_j, t+1}$ is a vector of the components of $\tilde{\beta}_{t+1}$
which correspond to indices in $C_j$.

Common random numbers coupling of an exact and an approximate
Markov chain for Bayesian logistic regression with spike and slab priors.

* Sample $(e_{t+1}, \tilde{e}_{t+1}) | \beta_{t+1}, \tilde{\beta}_{t+1}, z_{t+1}, \tilde{z}_{t+1}, w_t, \tilde{w}_t$ with common random numbers component-wise independently such that for each $i=1,\mydots, n$

e_i,t+1    ∼{
𝒩( x_i^⊤β_t+1,  w_i,t^-1) I_[0, ∞)     if  y_i=1

𝒩( x_i^⊤β_t+1,  w_i,t^-1) I_(-∞, 0)     if  y_i=0

ẽ_i,t+1    ∼{
𝒩( x_Ã,i^⊤β̃_Ã,t+1,  w̃_i,t^-1) I_[0, ∞)     if  y_i=1

𝒩( x_Ã,i^⊤β̃_Ã,t+1,  w̃_i,t^-1) I_(-∞, 0)     if  y_i=0

where $x_i^\top $ and $x_{\tilde{A}, i}^\top $ are the $i^{th}$ row of the $X$ and $X_{\tilde{A}}$ respectively.

* Sample $(w_{t+1}, \tilde{w}_{t+1}) | \beta_{t+1}, \tilde{\beta}_{t+1}, z_{t+1}, \tilde{z}_{t+1}, e_{t+1}, \tilde{e}_{t+1}$.  We take this variable to be fixed,  and set $w_{i,t} = \tilde{w}_{i,t} = 3/\pi^2$ for all $i=1, \mydots, n$ and $t \geq 0$, where $3/\pi^2$ is the precision of the logistic distribution.  In the case this variable can vary,  they can be sampled using common random numbers such that for each $i=1, \mydots, n$,

w_i,t+1 ∼Γ( ν+1/2,  K (y_i- x_i^⊤β_t+1)^2/2 )  and
w̃_i,t+1 ∼Γ( ν+1/2,  K(y_i- x_Ã,i^⊤β̃_Ã,t+1)^2/2 )

where $\nu=7.3$, $K \defeq (\pi^2 (\nu-2)/3)$ are fixed constants as given in Narisetty et al., 2019.

$C_t \defeq (\beta_t, z_t, e_t, w_t)$ and $C_t \defeq (\tilde{\beta}_t, \tilde{z}_t, \tilde{e}_t, \tilde{w}_t)$.

§ ADDITIONAL ALGORITHMS

$(X_t, Y_t)$, unnormalized densities $p$ and $q$
of $P$ and $Q$ respectively, step sizes $\sigma_P$ and $\sigma_Q$
Sample $\epsilon_{CRN} \sim \mathcal{N}(0,I_d)$. Calculate proposals

$$X^* \defeq X_t + \frac{1}{2} \sigma_P^2 \nabla \log p(X_{t}) + \sigma_P \epsilon_{CRN} \text{ and } Y^* \defeq Y_t + \frac{1}{2} \sigma_Q^2 \nabla \log q(Y_{t}) + \sigma_Q \epsilon_{CRN}$$

Sample $U_{CRN} \sim \Uniform([0,1])$
if $U_{CRN} \leq \frac{p(X^*) \mathcal{N}( X^* ; X_t + \frac{1}{2} \sigma_P^2 \nabla \log p(X_{t}), \sigma_P^2I_d )}{p(X_t) \mathcal{N}( X_t ; X^* + \frac{1}{2} \sigma_P^2 \nabla \log p(X^*), \sigma_P^2I_d )}  $, then set \( X_{t+1} = X^* \) ; else set \( X_{t+1} = X_{t} \)
if $U_{CRN} \leq \frac{q(Y^*) \mathcal{N}( Y^* ; Y_t + \frac{1}{2} \sigma_Q^2 \nabla \log q(Y_{t}), \sigma_Q^2I_d )}{q(Y_t) \mathcal{N}( Y_t ; Y^* + \frac{1}{2} \sigma_Q^2 \nabla \log q(Y^*), \sigma_Q^2I_d )}  $, then set \( Y_{t+1} = Y^* \) ; else set \( Y_{t+1} = Y_{t} \)
$(X_{t+1}, Y_{t+1})$
Common random numbers coupling of two MALA kernels marginally
targetting distributions $P$ and $Q$ respectively

$(X_t, Y_t)$, unnormalized densities $p$ and $q$
of $P$ and $Q$ respectively, step sizes $\sigma_P$ and $\sigma_Q$
Sample $\epsilon_{CRN} \sim \mathcal{N}(0,I_d)$. Calculate proposals

$$X^* \defeq X_t + \frac{1}{2} \sigma_P^2 \nabla \log p(X_{t}) + \sigma_P \epsilon_{CRN} \text{ and } Y^* \defeq Y_t + \frac{1}{2} \sigma_Q^2 \nabla \log q(Y_{t}) + \sigma_Q \epsilon_{CRN}. $$

Sample $U \sim \Uniform([0,1])$
if $U \leq \frac{p(X^*) \mathcal{N}( X^* ; X_t + \frac{1}{2} \sigma_P^2 \nabla \log p(X_{t}), \sigma_P^2I_d )}{p(X_t) \mathcal{N}( X_t ; X^* + \frac{1}{2} \sigma_P^2 \nabla \log p(X^*), \sigma_P^2I_d )}  $, then set \( X_{t+1} = X^* \) ; else set \( X_{t+1} = X_{t} \)
Set \( Y_{t+1} = Y^* \)
$(X_{t+1}, Y_{t+1})$
Common random numbers coupling of a MALA kernel and an ULA kernel
marginally targeting distributions $P$ and $Q$ respectively