Scalable Spike-and-Slab
Niloy Biswas

Harvard University

Lester Mackey<EMAIL_ADDRESS>
Microsoft Research New England

Xiao-Li Meng<EMAIL_ADDRESS>
Harvard University

Scalable Spike-and-Slab

Niloy Biswash
Lester Mackeymsr
Xiao-Li Mengh

hDepartment of Statistics, Harvard University
msrMicrosoft Research New England

Niloy<EMAIL_ADDRESS>
Machine Learning, ICML

Spike-and-slab priors are commonly used for Bayesian variable selection, due to
their interpretability and favorable statistical properties.
However, existing samplers for spike-and-slab posteriors incur prohibitive computational costs when the number of variables is large.
In this article, we propose Scalable Spike-and-Slab (), a
scalable Gibbs sampling implementation for high-dimensional
Bayesian regression with the continuous spike-and-slab
prior of George & McCulloch, 1993.
For a dataset with $n$ observations and $p$ covariates, has order $\max \{ n^2 p_t, np \}$ computational cost at iteration $t$ where $p_t$
never exceeds the number of covariates switching
spike-and-slab states between iterations $t$ and $t-1$
of the Markov chain.
This improves upon the order  $n^2p$ per-iteration cost of state-of-the-art
implementations as, typically, $p_t$ is substantially smaller than $p$.
We apply on synthetic and real-world datasets, demonstrating orders of magnitude speed-ups over existing exact samplers and significant gains in inferential quality over approximate samplers with comparable cost.

§ INTRODUCTION

§.§ Bayesian computation in high dimensions

We consider linear,  logistic, and probit regression in high dimensions,
where the number of observations $n$ is
smaller than the number of covariates $p$.
This setting is common in modern applications such as
genome-wide association studies [Guan & Stephens, 2011, Zhou et al., 2013]
and astronomy [Kelly, 2007, Sereno, 2015].
In the non-Bayesian paradigm, sparse point estimates
such as the LASSO [Tibshirani, 1996], Elastic Net [Zou & Hastie, 2005]
and SLOPE [Bogdan et al., 2015]

offer a route to variable selection. These estimates are based on
optimization based approaches, which are
computationally efficient and scale to
datasets with hundreds of thousands of covariates.

In the Bayesian paradigm, which will be our focus, one places a prior
on the unknown parameters of interest and considers the corresponding posterior distribution.
Sampling algorithms such as Markov chain Monte Carlo (MCMC) are then
used to simulate from the posterior distribution. In modern high dimensional settings,
general-purpose MCMC algorithms can have high computational cost per iteration.
This has kindled a line of work on tailored algorithms for Bayesian regression
Our manuscript participates in this wider effort to scale Bayesian inference to large data
applications. Specifically, we propose computationally efficient MCMC algorithms for
high-dimensional Bayesian linear, logistic and probit regression with spike-and-slab priors.

§.§ Variable selection with spike-and-slab priors

Consider Gaussian linear regression, logistic regression, and probit regression
with $n$ observations and $p$ covariates.  The respective likelihoods are given by
$L_{\textup{lin}}(\bm{\beta}, \sigma^2 ; \bm{X},\bm{y}) = \frac{1}{(2 \pi \sigma^2)^{n/2}} \exp \big(- \frac{\sum_{i=1}^n (y_i - \bm{x_i}^\top \bm{\beta})^2}{2\sigma^2}\big)$,
$L_{\textup{log}}(\bm{\beta} ; \bm{X},\bm{y}) = \prod_{i=1}^n \frac{\exp(y_i \bm{x_i}^\top \bm{\beta})}{1+\exp(\bm{x_i}^\top \bm{\beta})}$, and
$L_{\textup{prob}}(\bm{\beta} ; \bm{X},\bm{y}) = \prod_{i=1}^n \Phi( \bm{x_i}^\top \bm{\beta} )^{y_i} (1-\Phi( \bm{x_i}^\top \bm{\beta} ))^{1-y_i}$.
$\bm{X} \in \mathbb{R}^{n \times p}$ is the design matrix with rows $\bm{x_i}^\top$,
$\bm{y} \in \mathbb{R}^n$ (for linear regression) or $\bm{y} \in \{0,1\}^n$ (for logistic and probit regression)
is the response vector, $\bm{\beta} \in \mathbb{R}^p$ is the unknown signal,
$\sigma^2 \in (0, \infty)$ is the unknown Gaussian noise variance, and
$\Phi$ is the cumulative density function of ${\cal N}(0,1)$.

We focus on the high-dimensional
setting with $n \ll p$, where $\bm{\beta} \in \mathbb{R}^p$
is assumed to be sparse.  We use a continuous spike-and-slab
prior on $\bm{\beta}$ to capture sparsity:

σ^2 ∼(a_0/2, b_0/2 ),
(for linear regression)

σ^2 = 1,  (for logistic and probit regression)

z_j j=1,,pi.i.d.∼ (q),

β_j | z_j, σ^2 j=1,,p ind ∼ (1-z_j) 𝒩(0, σ^2 τ_0^2) + z_j 𝒩(0, σ^2 τ_1^2),

σ^2 ∼(a_0/2, b_0/2 ),
(for linear regression)

σ^2 = 1,  (for logistic and probit regression)

z_j j=1,,pi.i.d.∼ (q),

β_j | z_j,   σ^2  j=1,,pind∼  (1  -  z_j)𝒩(0,   σ^2 τ_0^2)  +  z_j𝒩(0,   σ^2 τ_1^2),

$q\in (0,1)$, $\tau_1^2 \gg \tau_0^2 >0$, and $a_0, b_0 > 0$ are
hyperparameters. Here, $\mathcal{N}(0, \sigma^2 \tau_0^2)$
and $\mathcal{N}(0, \sigma^2 \tau_1^2)$ correspond to the spike and slab parts
of the prior respectively.
In the high-dimensional setting,
a small constant $q \ll 1$ is often chosen, but the algorithms in this manuscript also
readily extend to hierarchical variants of (<ref>) with a hyperprior placed on $q$ [Scott & Berger, 2010, Castillo & van der Vaart, 2012].

Catalyzed by the works of George & McCulloch, 1993, George & McCulloch, 1997,
continuous spike-and-slab priors are now a mainstay of Bayesian variable
selection (see the recent reviews of <cit.>
and Banerjee et al., 2021).
The posterior probabilities $\mathbb{P}(z_j=1| \bm{y})$ provide a natural
interpretable approach to variable selection. The median probability model selects all covariates $j$ such that $\mathbb{P}(z_j=1| \bm{y}) > \sfrac{1}{2}$, and it is easily fitted using Monte Carlo samples and provides the optimal predictive model in the case with orthogonal design matrix, as well as extensions with certain correlated matrices [Barbieri & Berger, 2004, Barbieri et al., 2021].
Narisetty & He, 2014 have further fine-tuned the
optimal scaling of $q, \tau_0^2$, and $\tau_1^2$ with respect to the number of
covariates and sample size to establish model selection consistency for
linear regression with general design matrices in high dimensions.

One could alternatively consider point-mass spike-and-slab priors
<cit.>, where
$\beta_j | z_j \sim (1-z_j) \delta_0(\cdot) + z_j \mathcal{N}(0, \sigma^2 \tau_1^2)$
such that a degenerate Dirac distribution about zero is chosen for the spike part.
Point-mass priors have favorable statistical properties

and we hope to extend our algorithms to point-mass priors in follow-up work.

§.§ Our contributions

Our contributions are summarized below.
Throughout, we use $\mathcal{O}$ and $\Omega$ to respectively denote asymptotic upper and lower
bounds on computational complexity growth rates.

Section <ref> introduces
Scalable Spike-and-Slab (),
a computationally efficient implementation of Gibbs samplers
for linear and logistic regression with the
prior given in (<ref>).
Section <ref> investigates the
computational bottlenecks of state-of-the-art (SOTA)
implementations, which require $\Omega(n^2p)$

computational cost per iteration for datasets with $n$ observations and $p$ covariates.
Section <ref> develops ,
which overcomes existing computational bottlenecks by employing
a pre-computation based strategy and requires
$\mathcal{O}(\max\{n^2 p_t, np\})$ computational cost at iteration $t$,
where $p_t$ is no greater than the number of covariates switching
spike-and-slab states between iterations $t$ and $t-1$
of the Markov chain.
Section <ref> analyzes
the favourable computational complexity of , showing
that $p_t$ is typically much smaller than $p$ and that it can remain constant and even approach zero
under various limiting regimes as $p$ increases.

Section <ref> compares
with the SOTA exact MCMC sampler and a recently
proposed approximate MCMC sampler, which does not converge to the posterior distribution of interest.
We demonstrate that offers substantially
faster numerical runtimes compared to the SOTA exact MCMC sampler,
reporting $50\times$ speedups on synthetic datasets.
In the same experiment,
and the approximate sampler have comparable runtimes, but the asymptotically exact procedure provides more accurate variable selection.

Section <ref> demonstrates
the benefits of on a diverse suite of datasets, including two synthetic datasets and eight real-world experiments.

For example, on a genome-wide association study (GWAS) dataset with
with $n\!\approx\!2000$ and $p\!\approx\!100000$, we again observe
$50\times$ computational speedups over the SOTA exact MCMC sampler.
Finally, Section <ref> discusses directions for future work.

The open-source packages $\href{https://CRAN.R-project.org/package=ScaleSpikeSlab}{\mathrm{ScaleSpikeSlab}}$ in
$\mathrm{R}$ and $\mathrm{Python}$

implement our methods and recreate the experiments in this paper.

§ SCALABLE SPIKE-AND-SLAB

§.§ Status quo and computational bottlenecks

Gibbs samplers have long been employed to sample from the
posterior distributions corresponding to the prior in (<ref>)
<cit.>. The
computational bottleneck of existing Gibbs samplers
is linked to  sampling from the full conditional of
$\bm{\beta} \in \mathbb{R}^p$. This is given by

\begin{equation} \label{eq:full_conditional_beta}
\bm{\beta}_{t+1} | \bm{z}_t, \sigma^2_t \sim
\mathcal{N} \big(\bm{\Sigma}^{-1}_t \bm{X}^\top \bm{y}, \sigma_t^2 \bm{\Sigma}^{-1}_t \big)\quad {\rm for}\quad  \bm{\Sigma}_t = \bm{X}^\top \bm{X} + \bm{D}_t,
\end{equation}

\begin{equation} \label{eq:full_conditional_beta}
\bm{\beta}_{t+1} | \bm{z}_t, \sigma^2_t \! \sim \!
\mathcal{N} \big(\bm{\Sigma}^{-1}_t \bm{X}^\top \bm{y}, \sigma_t^2 \bm{\Sigma}^{-1}_t \big)
\end{equation}

icmlfor $\bm{\Sigma}_t \!=\! \bm{X}^\top \bm{X} + \bm{D}_t$,
where $t$ indexes the iteration of the Markov chain,

and $\bm{D}_t$ is the diagonal matrix with the vector
$\bm{z}_t \tau^{-2}_1 + (\bm{\mathrm{1}}_p-\bm{z}_t) \tau^{-2}_0$
populating its diagonal elements.

Sampling from (<ref>) using standard matrix multiplication and a generic Cholesky
decomposition that ignores the specific structure
of $\bm{\Sigma}_t$ requires $\Omega(p^3)$ computational cost,
which quickly becomes prohibitive for large $p$.
Hereafter we will refer to this generic method as the Naïve Sampler.
Ishwaran & Rao, 2005 recommend separating the components
of $\bm{\beta}$ into $B$ blocks of size $\sfrac{p}{B}$ each, and then updating each block using Gibbs
sampling, which gives a reduced $\Omega(p^3 B^{-2})$ computational cost.
However, this cost remains prohibitive for large $p$, and using a
larger number of blocks $B$ induces higher auto-correlation
between successive iterations of the Gibbs sampler.
Recently, Bhattacharya et al., 2016 developed
an algorithm based on the Woodbury matrix identity [Hager, 1989]
to sample from multivariate Gaussian distributions of the
form in (<ref>),
which requires a more favourable $\Omega(n^2p)$
computational cost and is given in Algorithm <ref>.
An $\Omega(n^2p)$ sampler of (<ref>)
[Bhattacharya et al., 2016]

Sample $\bm{r} \sim \mathcal{N}(0, \bm{I}_p)$, $\bm{\xi} \sim \mathcal{N}(0,\bm{I}_n)$.
Set $\bm{u} = \bm{D}_t^{-\frac{1}{2}}\bm{r}$ and calculate $\bm{v} = \bm{X}\bm{u} + \bm{\xi}$.
Set $\bm{v}^* = \bm{M}_t^{-1} (\frac{1}{\sigma_t} \bm{y} - \bm{v})$ for $\bm{M}_t = \bm{I}_n + \bm{X} \bm{D}_t^{-1} \bm{X}^\top$.
Return $\bm{\beta} = \sigma_t ( \bm{u} +  \bm{D}_t^{-1} \bm{X}^\top \bm{v}^* )$.

For large-scale datasets with $n$ in the thousands and $p$ in the hundreds
of thousands, as found in modern scientific applications,
the $\Omega(n^2p)$ cost per iteration is still too high.
This has spurred recent work on approximate
MCMC [Narisetty et al., 2019] for the continuous
spike-and-slab prior on logistic regression and on variational inference methods
for the point-mass spike-and-slab prior
[Titsias & Lázaro-Gredilla, 2011, Ray et al., 2020, Ray & Szabó, 2021].

Such approximate samplers can provide improved computational
speeds but do not converge to the posterior distribution of interest.

§.§ A scalable Gibbs sampler

We now develop . Our key insight is that successive pre-computation can be used to reduce the computational
cost of Algorithm <ref>. In Algorithm <ref>,
the $\Omega(n^2p)$ computational cost per iteration
arises from the calculation of the matrix $\bm{M}_t$.
Current algorithms calculate $\bm{M}_t = \bm{I}_n + \bm{X} \bm{D}_t^{-1} \bm{X}^\top$ from scratch
every iteration at $\Omega(n^2p)$ cost under standard matrix multiplication and then

solve an $n$ by $n$ linear system to obtain $v^*$
at $\Omega(n^3)$ cost.
We propose instead to use the previous state $\bm{z}_{t-1}$ and the pre-computed matrices $\bm{M}_{t-1}$ and $\bm{M}_{t-1}^{-1}$
at each step to aid the calculation of $\bm{M}_t^{-1}$. Our strategy is given below.

Specifically, denote $\bm{X} = [\bm{X}_1,\mydots,\bm{X}_p] \in \mathbb{R}^{n \times p}$ and its sub-matrices $\bm{X}_{A_t} \defeq [\bm{X}_j : j \in A_t] \in \mathbb{R}^{n \times \|\bm{z}_t\|}$
and $\bm{X}_{A^c_t} \defeq [\bm{X}_j : j \in A^c_t] \in \mathbb{R}^{n \times (p - \|\bm{z}_t\|)}$, where  $A_t \defeq \{ j: z_{j,t} = 1 \}$ and $A_t^c = \{ j: z_{j,t} = 0 \}$,  are the ordered index sets of covariates corresponding to slab states and to spike states respectively at iteration $t$, and $\|\cdot\|_1$ is the $\lp[1]$ norm on $\{0,1\}^p$. Also denote $\Delta_t \defeq \{ j: z_{j,t} \neq z_{j,t-1} \}$, the ordered index set of
covariates which switch spike-and-slab states between iterations $t$ and $t-1$;  $\delta_t \defeq \|\bm{z}_t-\bm{z}_{t-1}\|_1 = |\Delta_t|$, the number
of switches;
$\bm{D}_{\Delta_t} \defeq \Diag((\bm{D}_{t})_{j,j} : j \in \Delta_t)
\in \Diag(\mathbb{R}^{\delta_t \times \delta_t})$,
the diagonal sub-matrix of $\bm{D}_t$ composed of the diagonal entries with
ordered indices in $\Delta_t$; and $ \bm{C}_{\Delta_t} \defeq \bm{D}_{\Delta_t}^{-1} - \bm{D}_{\Delta_{t-1}}^{-1}$. Finally, let $\bm{\tilde{M}}_{\tau_0} \defeq \bm{I}_n + \tau_0^2 \bm{X} \bm{X}^T \in \mathbb{R}^{n \times n}$ and
$\bm{\tilde{M}}_{\tau_1} \defeq \bm{I}_n + \tau_1^2 \bm{X} \bm{X}^T \in \mathbb{R}^{n \times n}$, which
are fixed for all iterations.

Under this notation, there are three expressions for $\bm{M}_t$:

M_t    = M̃_τ_0 + (τ_1^2-τ_0^2) X_A_t X_A_t^⊤

= M̃_τ_1 + (τ_0^2-τ_1^2) X_A^c_t X_A^c_t^⊤

= M_t-1 + X_Δ_t  C_Δ_t X_Δ_t^⊤.

M_t    = M̃_τ_0 + (τ_1^2-τ_0^2) X_A_t X_A_t^⊤

= M̃_τ_1 + (τ_0^2-τ_1^2) X_A^c_t X_A^c_t^⊤

= M_t-1 + X_Δ_t  C_Δ_t X_Δ_t^⊤.

In (<ref>) – (<ref>),
calculating the matrix products
$\bm{X}_{A_t} \bm{X}_{A_t}^\top$, $\bm{X}_{A^c_t} \bm{X}_{A^c_t}^\top$, and $\bm{X}_{\Delta_t}  \bm{C}_{\Delta_t} \bm{X}_{\Delta_t}^\top$ requires
$\mathcal{O}(n^2 \|\bm{z}_t\|_1)$, $\mathcal{O}(n^2 (p-\|\bm{z}_t\|_1))$, and
$\mathcal{O}(n^2 \delta_t)$ cost respectively. Given
$\bm{\tilde{M}}_{\tau_0}$, $\bm{\tilde{M}}_{\tau_1}$, $\bm{M}_{t-1}$, and $\bm{z}_{t-1}$,
we evaluate whichever matrix product in
(<ref>) – (<ref>)
has minimal computational cost and thereby calculate $\bm{M}_t$ at the reduced cost of
$\mathcal{O}(n^2 p_t)$ where $p_t \defeq \min \{ \|\bm{z}_t\|_1, p-\|\bm{z}_t\|_1, \delta_t \}$.

To calculate $\bm{M}_t^{-1}$, we consider the cases $n \leq p_t$ and $p_t < n$ separately.
When $n \leq p_t$, we calculate $\bm{M}_t^{-1}$ by directly inverting the
calculated matrix $\bm{M}_t$ from (<ref>),
which requires $\mathcal{O}(n^3)$ cost.
When $p_t < n$, we apply the Woodbury matrix identity on
(<ref>). This gives

M_t^-1    = M̃_τ_0^-1 -M̃_τ_0^-1 X_A_t ( (τ_1^2-τ_0^2)^-1 I_z_t_1 +
X_A_t^⊤M̃_τ_0^-1 X_A_t )^-1 X_A_t^⊤M̃_τ_0^-1

= M̃_τ_1^-1 -M̃_τ_1^-1 X_A_t^c ( (τ_0^2-τ_1^2)^-1 I_p-z_t_1 +
X_A^c_t^⊤M̃_τ_1^-1 X_A^c_t )^-1 X_A^c_t^⊤M̃_τ_1^-1

= M_t-1^-1 - M_t-1^-1X_Δ_t (  C_Δ_t^-1 +
X_Δ_t^⊤M_t-1^-1 X_Δ_t )^-1 X_Δ_t^⊤M_t-1^-1.

M_t^1     = M̃_τ_0^1

-M̃_τ_0^1  X_A_t  (  1/τ_1^2-τ_0^2  I_z_t_1  +
X_A_t^⊤ M̃_τ_0^1  X_A_t  )^1  X_A_t^⊤ M̃_τ_0^1

= M̃_τ_1^1

-M̃_τ_1^1  X_A_t^c  (  1/τ_0^2-τ_1^2  I_p-z_t_1  +   X_A^c_t^⊤ M̃_τ_1^1  X_A^c_t  )^1  X_A^c_t^⊤ M̃_τ_1^1

= M_t1^1

-M_t1^1X_Δ_t (  C_Δ_t^1 +
X_Δ_t^⊤M_t1^1 X_Δ_t )^1 X_Δ_t^⊤M_t1^1.

Given $\bm{\tilde{M}}_{\tau_0}^{-1}$, $\bm{\tilde{M}}_{\tau_1}^{-1}$, $\bm{M}_{t-1}^{-1}$ and $\bm{z}_{t-1}$,
we evaluate whichever expression in (<ref>) – (<ref>)
has minimal computational cost to calculate $\bm{M}_t^{-1}$.
Similar to (<ref>), this requires
$\mathcal{O}(n^2 p_t)$ computational cost,
which arises from matrix inversion and multiplication.

Overall, this strategy of using the previous state $\bm{z}_{t-1}$ and the
pre-computed matrices $\bm{\tilde{M}}_{\tau_0}$, $\bm{\tilde{M}}_{\tau_0}^{-1}$,
$\bm{\tilde{M}}_{\tau_1}$, $\bm{\tilde{M}}_{\tau_1}^{-1}$, $\bm{M}_{t-1}$ and $\bm{M}_{t-1}^{-1}$,
reduces the computational cost of calculating the matrices $\bm{M}_{t}$ and $\bm{M}_{t}^{-1}$ from
$\Omega(n^2p)$ (as in all current implementations of Algorithm
<ref>) to $\mathcal{O}(n^2 p_t)$.

As we show in Sections <ref> and
<ref>, in many large-scale
applications $p_t$ is orders of magnitude smaller
than both $n$ and $p$, yielding substantial improvements in
computational efficiency.
Furthermore, we emphasize that the matrices $\bm{\tilde{M}}_{\tau_0}$,
$\bm{\tilde{M}}_{\tau_0}^{-1}$, $\bm{\tilde{M}}_{\tau_1}$, $\bm{\tilde{M}}_{\tau_1}^{-1}$
are fixed for all iterations, and the state $\bm{z}_{t-1}$ and matrices
$\bm{M}_{t-1}$ and $\bm{M}_{t-1}^{-1}$ only need to be stored temporarily to
generate samples for iteration $t$ and can be deleted after.
Therefore requires minimal additional memory compared
to current implementations.

The full Gibbs samplers for Bayesian linear, logistic,
and probit regression which make use of this pre-computation are given in Algorithms <ref> and <ref>.
The Gibbs samplers for logistic and probit regression are based on data augmentation strategies (see, e.g., O'Brien & Dunson, 2004, Narisetty et al., 2019), and the Gibbs sampler for logistic regression
requires an adjusted pre-computation strategy with $\mathcal{O} \big(\max\{ n^2 p_t, n^3, np \} \big)$ cost.
Appendix <ref> contains derivations and details of Algorithms <ref> and <ref> (the implementation of logistic regression is based on a scaled $t$-distribution
approximation to the logistic distribution, as commonly done in the literature; see [Narisetty et al., 2019]).

Bayesian linear regression with

Input: State $ \bm{C}_t \defeq (\bm{\beta}_t, \bm{z}_t, \sigma_t^2) \in
\mathbb{R}^p \times \{0,1\}^p \times (0,\infty)$,
state $\bm{z}_{t-1}$, and matrices $\bm{M}_{t-1}, \bm{M}_{t-1}^{-1}$.
[1]
Calculate $p_t$ and use $\eqref{eq:pre_compute_correction1}$ to calculate $\bm{M}_t$.

$p_t \geq n$   invert $\bm{M}_t$ to calculate $\bm{M}_t^{-1}$
use (<ref>) to calculate $\bm{M}_t^{-1}$.

Sample $\bm{\beta}_{t+1} | \bm{z}_t, \sigma_t^2$ using Algorithm <ref> from

$\mathcal{N} \big( \bm{\Sigma}^{-1}_t \bm{X}^\top \bm{y},  \sigma_t^2 \bm{\Sigma}^{-1}_t \big) \text{ for } \bm{\Sigma}_t = \bm{X}^\top \bm{X} + \bm{D}_t.$

Sample each $z_{j,t+1} | \bm{\beta}_{t+1}, \sigma_t^2$ independently from

$\Bernoulli \Big(\frac{q \mathcal{N}(\beta_{j,t+1}; 0, \sigma_t^2 \tau_1^2)}{q \mathcal{N}(\beta_{j,t+1}; 0, \sigma_t^2 \tau_1^2) + (1-q) \mathcal{N}(\beta_{j,t+1}; 0, \sigma_t^2 \tau_0^2)} \Big).$

for $j=1,\mydots,p$.
Sample $\sigma^{2}_{t+1} | \bm{\beta}_{t+1}, \bm{z}_{t+1}$ from

$\InvGamma \Big( \frac{a_0+n+p}{2}, \frac{b_0 + \| \bm{y} - \bm{X} \bm{\beta}_{t+1} \|_2^2 + \bm{\beta}_{t+1}^\top \bm{D}_{t+1} \bm{\beta}_{t+1} }{2} \Big).$

Output: $ \bm{C}_{t+1} = (\bm{\beta}_{t+1}, \bm{z}_{t+1}, \sigma^2_{t+1})$,
$\bm{z}_t$, $\bm{M}_{t}, \bm{M}_{t}^{-1}$.

Bayesian logistic & probit regression with

State $ \bm{C}_t \defeq (\bm{\beta}_t, \bm{z}_t, \bm{\tilde{y}}_t, \bm{\tilde{\sigma}}^2_t) \in \mathbb{R}^p \times \{0,1\}^p \times \mathbb{R}^n \times (0,\infty)^n$,
states $\bm{z}_{t-1}$, $\bm{\tilde{\sigma}}^2_{t-1}$, and matrices $\bm{M}_{t-1}, \bm{M}_{t-1}^{-1}$.
[1]
Logistic regression:
Use pre-computation (see Appendix <ref>) to calculate $\bm{M}_t \defeq \bm{I}_n + \bm{W}_t^{-1/2} \bm{X} \bm{D}_t^{-1} \bm{X}^\top \bm{W}_t^{-1/2}$ for
$\bm{W}_t = \Diag(\bm{\tilde{\sigma}}^2_t)$. Invert $\bm{M}_t$ to calculate $\bm{M}_t^{-1}$.

Probit regression: Calculate $p_t$ and use $\eqref{eq:pre_compute_correction1}$ to calculate $\bm{M}_t \defeq \bm{I}_n + \bm{X} \bm{D}_t^{-1} \bm{X}^\top$.
$p_t \geq n$   invert $\bm{M}_t$ to calculate $\bm{M}_t^{-1}$
use (<ref>) to calculate $\bm{M}_t^{-1}$.

Sample $\bm{\beta}_{t+1} | \bm{z}_t, \bm{\tilde{y}}_t, \bm{\tilde{\sigma}}_t^2$ using Algorithm <ref> from

$\mathcal{N} \big( \bm{\Sigma}^{-1}_t \bm{X}^\top \bm{W}^{-1}_t \bm{\tilde{y}}_t, \bm{\Sigma}^{-1}_t \big) \text{ for } \bm{\Sigma}_t\!=\!\bm{X}^\top \! \bm{W}_t^{-1} \! \bm{X}\!+\!\bm{D}_t.$

Sample each $z_{j,t+1} | \bm{\beta}_{t+1}, \bm{\tilde{y}}_t, \bm{\tilde{\sigma}}_t^2$ independently from

$\Bernoulli \Big(\frac{q \mathcal{N}(\beta_{j,t+1}; 0, \tau_1^2)}{q \mathcal{N}(\beta_{j,t+1}; 0,  \tau_1^2) + (1-q) \mathcal{N}(\beta_{j,t+1}; 0, \tau_0^2)} \Big)$

for $j=1,\mydots,p$.
Sample each $\tilde{y}_{i,t+1} | \bm{\beta}_{t+1}, \bm{z}_{t+1}, \bm{\tilde{\sigma}}_t^2$ independently from

$\mathcal{N}( \bm{x_i}^\top  \bm{\beta}_{t+1}, \tilde{\sigma}_{i,t}^2) \mathrm{I}_{[0, +\infty)} \quad \text{if } y_i=1,$ $\quad
\mathcal{N}( \bm{x_i}^\top  \bm{\beta}_{t+1}, \tilde{\sigma}_{i,t}^2) \mathrm{I}_{(-\infty, 0)} \quad \text{if } y_i=0$

for $i=1,\mydots,n$.
Logistic regression: Sample each $\tilde{\sigma}_{i,t+1}^2 | \bm{\beta}_{t+1}, \bm{z}_{t+1}, \bm{\tilde{y}}_{t+1}$ independently from

$\InvGamma \Big( \frac{v+1}{2}, \frac{w^2 \nu +
(\tilde{y}_{i,t+1} - \bm{x_i}^\top \bm{\beta}_{t+1} )^2}{2} \Big)$

for $i=1,\mydots,n$, where $\nu \defeq 7.3$ and $w^2 \defeq \sfrac{\pi^2(\nu-2)}{(3 \nu)}$ are constants.

Probit regression: Set each $\tilde{\sigma}_{i,t+1}^2=1$ for $i=1,\mydots,n$.

Output: $\bm{C}_{t+1} \!=\! (\bm{\beta}_{t+1}, \bm{z}_{t+1}, \bm{\tilde{y}}_{t+1}, \bm{\tilde{\sigma}}^2_{t+1})$, $\bm{z}_t$, $\bm{M}_{t}$, $\bm{M}_{t-1}^{-1}$

§.§ Analysis of computational complexity

We now investigate the favorable computational
complexity of Algorithms <ref> and
<ref>. Proposition <ref>,
proved in Appendix <ref>,
gives the computational cost of these Gibbs samplers for linear and
logistic regression, showing an improvement over
existing implementations which have $\Omega(n^2p)$ cost.

Algorithm <ref>
and Algorithm <ref> for
probit regression both have a computational cost of
$\mathcal{O} \big(\max\{ n^2 p_t, np \} \big)$ at iteration $t$, and
Algorithm <ref> for
logistic regression has a computational cost of
$\mathcal{O} \big(\max\{ n^2 p_t, n^3, np \} \big)$
at iteration $t$, where
$p_t = \min \{ \|\bm{z}_t\|_1, p-\|\bm{z}_t\|_1, \delta_t \}$ for

In Proposition <ref>, $\|\bm{z}_t\|_1$ and $p-\|\bm{z}_t\|_1$ are the number of
slab covariates and the number of spike covariates respectively at iteration $t$ of the Markov chain,
and $\delta_t$ is the number of covariates switching spike-and-slab states between iterations
$t$ and $t-1$ of the Markov chain. Note that $p_t \leq \min \{ \|\bm{z}_t\|_1, p-\|\bm{z}_t\|_1 \} \leq p/2$ directly. In practice, there are a variety of scenarios under which $p_t$ is significantly smaller than $p/2$.

Sparse $\bm{z}_t$. Whenever $\bm{z}_t$ is sparse relative to the full dimensionality $p$,
we have $p_t \leq \|\bm{z}_t\|_1 \ll p$. Sparsity in $\bm{z}_t$ is a common occurrence
in high-dimensional regression with a sparse signal vector $\bm{\beta}^* \in \mathbb{R}^p$, as
$\|\bm{z}_t\|_1$ often closely approximates the true sparsity $s \defeq \|\bm{\beta}^*\|_0$
with high probability. This occurs, for instance, in the settings of Narisetty & He, 2014
and Narisetty et al., 2019, where strong model selection consistency of the continuous
spike-and-slab posterior is established for linear and logistic regression respectively.

Posterior concentration.
Even when $\min \{ \|\bm{z}_t\|_1, p-\|\bm{z}_t\|_1 \}$ is comparable to $p/2$,

concentration of the spike-and-slab posterior targeted by the Gibbs sampler can lead to $\delta_t$ and hence $p_t$ remaining much smaller than $p/2$. Proposition <ref>, proved in Appendix <ref>,
calculates the expectation of $\delta_t$ explicitly in terms of
the posterior distribution that is targeted by our algorithms and
the auto-correlation of the states $(\bm{z}_t)_{t \geq 0}$.

For $\delta_t$ as given in Proposition <ref>,

𝔼[ δ_t ]     = ∑_j=1^p ℙ(z_j,t=1)ℙ(z_j,t-1=0) +

ℙ(z_j,t=0)ℙ(z_j,t-1=1) -2cov(z_j,t,z_j,t-1).

𝔼[ δ_t ] = ∑_j=1^p ℙ(z_j,t=1)ℙ(z_j,t-1=0) + ℙ(z_j,t=0)ℙ(z_j,t-1=1) -2cov(z_j,t,z_j,t-1).

Suppose the Markov chain
generated by Algorithm <ref> or <ref>
is at its stationary distribution $\pi$ at iteration $t-1$. Then,

\begin{equation} \label{eq:expected_comp_cost2}
\mathbb{E}[ \delta_t ] = 2 \sum_{j=1}^p \mathrm{var}_{\pi}(z_{j,t}) (1-\mathrm{corr}_{\pi}(z_{j,t},z_{j,t-1})).
\end{equation}

In (<ref>), note
$\mathrm{var}_{\pi}(z_{j,t}) = \mathbb{P}_{\pi}(z_{j,t}=0) \mathbb{P}_{\pi}(z_{j,t}=1) \leq \sfrac{1}{4}$
for each component $j$, with equality only when
$z_{j,t} \sim \Bernoulli(\sfrac{1}{2})$.
Therefore all components $j$ with $\mathbb{P}_{\pi}(z_{j,t}=1)$ close to $0$ or $1$

do not contribute significantly
towards $\delta_t$ in expectation.

Such posterior concentration is guaranteed whenever $\bm{z}_t$ is convergent, be it to the true model selection vector as in Narisetty & He, 2014 and Narisetty et al., 2019 or to any other value.
In such circumstances we can have $\delta_t = \| \bm{z}_t - \bm{z}_{t-1} \|_1=o(p)$ and
even $\delta_t = o(\|\bm{z}_t\|_1)$, regardless of the magnitude of $\|\bm{z}_t\|_1$.

High auto-correlation.
High auto-correlation of the Gibbs sampler can also lead
to smaller values of $\delta_t$ and hence $p_t$.
We already see from Proposition <ref> that, even for  components $j$ with bimodal
marginal posterior distributions, high auto-correlation between
successive states $z_{j,t-1}$ and $z_{j,t}$ can yield lower
$\delta_t$ in expectation.
Proposition <ref>, proved in
Appendix <ref>, provides an additional exact expression for $\delta_t$ in terms of the empirical
correlation between $\bm{z}_t$ and $\bm{z}_{t-1}$.

Let $\tau_t=\sqrt{\|\bm{z}_t\|_1(p-\|\bm{z}_t\|_1)}$ and
$\rho_t$ be the empirical correlation between $\bm{z}_t$ and $\bm{z}_{t-1}$
(that is, the correlation between $z_{J, t}$ and $z_{J, t-1}$ when $J$ is
uniform on $\{1, \ldots, p\}$). Then, for $\delta_t$ as given in Proposition <ref>,

\begin{equation} \label{eq:exactp}
\delta_t = \|\bm{z}_t\|_1+\|\bm{z}_{t-1}\|_1 - \frac{2 \|\bm{z}_t\|_1 \|\bm{z}_{t-1}\|_1 + 2 \rho_t \tau_t \tau_{t-1}}{p}.
\end{equation}

Since $|\rho_t|\le 1$, Proposition <ref> implies

δ_t    ≥(z_t_1 - z_t-1_1)^2/p + (τ_t-τ_t-1)^2/p  and

δ_t    ≤(z_t_1 - z_t-1_1)^2/p + (τ_t+τ_t-1)^2/p.

The lower bound in (<ref>) is a good approximation to $\delta_t$ when $\rho_t$ is close to one,  which is the case either when the Gibbs sampler is converging (such that $\bm{z}_t$ becomes stable) or when it gets stuck (such that $\bm{z}_t$ changes slowly with $t$). In either case, $\tau_t$ exhibits similar “stable/stuck” behavior, implying that the lower bound itself will be close to zero. This suggests $\delta_t$ and hence $p_t$ is close to zero when $\rho_t$ is close to 1, even if $\min\{\|\bm{z}_t\|_1, p-\|\bm{z}_t\|_1\}$ is not negligible.

Motivated by such discussions and theoretical analysis,
we now empirically examine how $p_t$ grows as the number of
observations $n$, the number of covariates $p$,
and the sparsity $s$ of the true signal varies.
Figure <ref> is based on synthetic
linear regression datasets. For each dataset, one Markov chain is
generated using Algorithm <ref>
to target the corresponding spike-and-slab posterior,
from which the mean and one standard error bars of $(p_t)_{t = 5000}^{10000}$ are plotted.
For Figure <ref> (Left), we consider datasets
with $n=100$, varying $p$ with $p \geq n$,
and a sparse true signal $\bm{\beta}^* \in \mathbb{R}^p$ with
components $\bm{\beta}^*_j = 2 \mathrm{I}\{j \leq s\}$

for sparsity $s=10$, and noise standard deviation $\sigma^*=2$.
For Figure <ref> (Center), we consider datasets
with $p=1000$, varying $n$ with $n \leq p$, $s=10$, and
$\sigma^*=2$. For Figure <ref> (Right), we consider datasets
with $n=10s$, $p=1000$, and $\sigma^*=2$ for varying $s \geq 1$.

Details of the synthetically generated datasets are in
Appendix <ref>.

Figure <ref> (Left) shows that
both $p_t$ is substantially smaller than both $p$ and $n$
and that it does not increase with the number of covariates $p$.
Figure <ref> (Center) shows that
$p_t$ tends to zero as $n$ increases.
Figure <ref> (Right) shows that
$p_t$ decreases

as the sparsity $s$ increases. All figures suggest that $p_t$ is
controlled by $\delta_t$ in these settings, because $\|\bm{z}_t\|_1$
takes values close to $s$, and $p-\|\bm{z}_t\|_1$ tends
to be much larger than $\|\bm{z}_t\|_1$.

Overall, Figure <ref>
highlights that not only does $p_t$ tend to be substantially
smaller than $p$, but it also tends to
be smaller than both $n$ and $s$.

By Proposition <ref>,
this showcases the substantially lower computational cost of compared to current implementations which cost
$\Omega(n^2p)$ per iteration.

The cost parameter $p_t$ for averaged over iterations $5000 < t \leq 10000$
with one standard error bars,
for synthetic linear regression datasets with varying number of
covariates $p$ (Left), varying number of observations $n$ (Center), and
varying sparsity $s$ (Right). The ground truth sparsity $s$ is also plotted for comparison.
See Section <ref> for details.

The cost parameter $p_t$ for averaged over iterations $5000 < t \leq 10000$
with one standard error bars,
for synthetic linear regression datasets with varying number of
covariates $p$ (Left), varying number of observations $n$ (Center), and
varying sparsity $s$ (Right). The ground truth sparsity $s$ is also plotted for comparison.
See Section <ref> for details.

§.§ Extensions to Scalable Spike-and-Slab

With additional memory capacity and pre-computation, we can further improve the per-iteration costs of .

For the matrices $\bm{\tilde{M}}_{\tau_0}^{-1}$ and $\bm{\tilde{M}}_{\tau_1}^{-1}$ in Section <ref>, suppose the matrices $\bm{X}^\top \bm{X}$, $\bm{X}^\top \bm{\tilde{M}}_{\tau_0}^{-1} \bm{X}$, and $\bm{X}^\top \bm{\tilde{M}}_{\tau_1}^{-1} \bm{X}$ are pre-computed. This initial step requires $\mathcal{O}(n p^2)$ computational cost and $\mathcal{O}(p^2)$ memory. Then the matrices $\bm{X}_{A_t}^\top \bm{X}_{A_t}$ and $\bm{X}_{A^c_t}^\top \bm{X}_{A^c_t}$ in (<ref>) – (<ref>) correspond to pre-computed sub-matrices of $\bm{X}^\top \bm{X}$, and
calculating $\bm{M}_t$ using (<ref>) – (<ref>)
at iteration $t$ involves matrix addition which only requires $\mathcal{O}(n^2)$ cost. Similarly, matrices $\bm{X}_{A_t}^\top \bm{\tilde{M}}_{\tau_0}^{-1} \bm{X}_{A_t}$ and $\bm{X}_{A^c_t}^\top \bm{\tilde{M}}_{\tau_1}^{-1} \bm{X}_{A^c_t}$ in (<ref>) – (<ref>) correspond to pre-computed sub-matrices of $\bm{X}^\top \bm{\tilde{M}}_{\tau_0}^{-1} \bm{X}$ and $\bm{X}^\top \bm{\tilde{M}}_{\tau_1}^{-1} \bm{X}$ respectively and do not need to be recalculated at iteration $t$. Therefore calculating $\big( (\tau_1^2-\tau_0^2)^{-1} \bm{I}_{\|\bm{z}_t\|_1} + \bm{X}_{A_t}^\top \bm{\tilde{M}}_{\tau_0}^{-1} \bm{X}_{A_t} \big)^{-1}$ or $\big( (\tau_0^2-\tau_1^2)^{-1} \bm{I}_{p-\|\bm{z}_t\|_1} +
\bm{X}_{A^c_t}^\top \bm{\tilde{M}}_{\tau_1}^{-1} \bm{X}_{A^c_t} \big)^{-1}$ in (<ref>) – (<ref>) at each iteration $t$ only requires $\mathcal{O}(\|z_t\|_1^3)$ or $\mathcal{O}((p-\|z_t\|_1)^3)$ cost respectively.

To sample from (<ref>),  consider the cases $n \leq \min \{\|z_t\|_1, p-\|z_t\|_1\}$ and $\min \{\|z_t\|_1, p-\|z_t\|_1\} < n$ separately. When $n \leq \min \{\|z_t\|_1, p-\|z_t\|_1\}$, we calculate $\bm{M}_t^{-1}$ by directly inverting the calculated matrix $\bm{M}_t$ from (<ref>), which requires $\mathcal{O}(n^3)$ cost.
When $\min \{\|z_t\|_1, p-\|z_t\|_1\} < n$, we avoid calculating $M_t^{-1}$ explicitly and instead calculate the matrix vector product $M_t^{-1}(\frac{1}{\sigma_t}\bm{y}-\bm{v})$ in Algorithm <ref> right-to-left, using whichever expression in (<ref>) – (<ref>) has minimal computational cost.

Overall, now the Gibbs samplers for linear and probit regression require only $\mathcal{O}(\max\{ \min \{\|z_t\|_1, p-\|z_t\|_1,n\}^3, np \})$ computational cost at iteration $t$. This provides lower computational cost for linear and probit regression whenever $\min \{\|z_t\|_1, p-\|z_t\|_1,n\}^3 < n^2 p_t$.
A similar extension for logistic regression requires only  $\mathcal{O}(\max\{ n^3, np \})$ computational cost at iteration $t$ and is given in Appendix <ref>.

§ COMPARISON WITH ALTERNATIVES

In this section we compare with the naïve
sampler, the SOTA exact MCMC sampler based on the sampling algorithm of
Bhattacharya et al., 2016, and the Skinny Gibbs approximate MCMC sampler of Narisetty et al., 2019 for logistic regression.
Table <ref> highlights the favorable computational cost of compared to the naïve and SOTA samplers.
The Skinny Gibbs sampler typically has lower computational cost compared to for logistic regression, and
can have lower or higher computational cost than
for probit regression depending on whether $\|\bm{z}_t\|_1^2 \le n p_t$ or not.
However, unlike , the Skinny Gibbs sampler does not converge to the correct posterior distribution.

Comparison of with alternatives
($p_t = \min \{ \|\bm{z}_t\|_1, p-\|\bm{z}_t\|_1, \delta_t \}$
for $\delta_t = \|\bm{z}_t-\bm{z}_{t-1}\|_1$).

MCMC Sampler     Cost
Converges to

Naïve     $\Omega(p^3)$     $\surd$

State-of-the-art     $\Omega(n^2p)$     $\surd$

Skinny Gibbs     $\Omega(\max \{ n \|\bm{z}_t \|^2_1, np \})$     $\times$

$\mathbf{S^3}$ (linear and probit)     ${\mathcal{O}(\max \{ n^2 p_t, np \})}$     $\mathbf{\surd}$

$\mathbf{S^3}$ (logistic)     ${\mathcal{O}(\max \{ n^2 p_t, n^3, np \})}$     $\mathbf{\surd}$

To assess the practical impact of  computational cost and asymptotic bias,
Figures <ref> and <ref>
compares the numerical runtimes and statistical performance of
with the SOTA sampler and the Skinny Gibbs sampler.
For the Skinny Gibbs sampler, we use the $\mathrm{skinnybasad}$ $\mathrm{R}$ package
of Narisetty et al., 2019, which implements only logistic regression.
We consider synthetically generated datasets
with a true signal $\bm{\beta}^* \in \mathbb{R}^p$ where
$\beta^*_j = 2 \mathrm{I}\{j \leq s\}$ for sparsity $s$.
We consider datasets with $n=10s$ observations and
$p=100s$ covariates for varying sparsity $s \geq 1$.
Details of the synthetically generated dataset are in Appendix
<ref>. For each synthetic dataset,
we run for logistic and probit regression and
the Skinny Gibbs sampler for 1000 iterations, run the
SOTA sampler for 100 iterations, and record the average
time taken per iteration. All timings were obtained using a single core of
an Apple M1 chip on a Macbook Air 2020 laptop with 16 GB RAM.

Comparison of time per iteration between , state-of-the-art (SOTA)
exact MCMC sampler and the Skinny Gibbs approximate sampler of Narisetty et al., 2019
on synthetic binary classification datasets.

See Section <ref> for details.

Comparison of time per iteration between , state-of-the-art (SOTA)
exact MCMC sampler and the Skinny Gibbs approximate sampler of Narisetty et al., 2019
on synthetic binary classification datasets.

See Section <ref> for details.

Figure <ref>
highlights that the numerical runtimes of are orders of magnitude faster than the SOTA sampler and comparable to the Skinny
Gibbs sampler. For example, for $n=4000$
observations, $p=40000$ covariates, and sparsity $s=400$,
for logistic regression requires $21500$ms per iteration on average, which is approximately $15$ times faster than the SOTA sampler
for probit regression (which requires $335000$ms per iteration on average) and $2.5$ times faster than the Skinny Gibbs sampler
(which requires $55600$ms per iteration on average), and for probit regression requires $1100$ms per iteration on average, which is approximately $50$ times faster than the SOTA sampler
for probit regression (which requires $55800$ms per iteration on average).

For larger real-life datasets with hundreds of thousands of
covariates, the numerical runtimes of are similarly
favorable compared to the SOTA sampler.
This is showcased in Section <ref>,
where for a genetics dataset, is
$50$ times faster than the SOTA sampler.

Average true positive rate (TPR) and false discovery rate (FDR) of
and the Skinny Gibbs approximate sampler [Narisetty et al., 2019] across $20$ independently generated datasets with one standard error bars.

See Section <ref> for details.

Average true positive rate (TPR) and false discovery rate (FDR) of
and the Skinny Gibbs approximate sampler [Narisetty et al., 2019] across $20$ independently generated datasets with one standard error bars.

See Section <ref> for details.

Figure <ref> plots the true positive rate (TPR)
and the false discovery rate (FDR) of variable selection based on samples
from for logistic and probit regression and Skinny Gibbs on synthetic binary classification datasets.
To assess variable selection for signals of varying magnitude,
we consider an exponentially decaying sparse true signal
$\bm{\beta}^* \in \mathbb{R}^p$ such that $\beta^*_j = 2^{\frac{9-j}{4}}$
for $j \leq s$ and $\beta^*_j=0$ for $j>s$ for sparsity $s$.
The corresponding synthetically generated datasets have $n=100$
observations, varying number of covariates $p$ with $p\geq n$, and
sparsity $s=5$. For each synthetic dataset,
we implement for logistic and probit regression
and the Skinny Gibbs sampler for $5000$ iterations with a burn-in of $1000$
iterations and calculate the TPR and FDR from the samples. We use the same prior
hyperparameters for all the algorithms, which are chosen
according to Narisetty et al., 2019.
Additional experimental details are included in Appendices <ref>.
The SOTA sampler is not shown in Figure <ref>,
as SOTA and are alternative implementations of the same Gibbs sampler
and by definition have the same statistical performance.
Figure <ref> shows that in
higher dimensions, samples from yield significantly higher TPR and lower FDR
than the Skinny Gibbs sampler. We observe similar results for other choices
of prior hyperparameters, which give to either have comparable or more favorable
statistical performance to the Skinny Gibbs sampler.

Overall, Figures <ref> and <ref> highlight
that can have comparable or even favorable computational cost to the
Skinny Gibbs sampler, whilst having the correct stationary distribution and
more favorable statistical properties in higher dimensions.
Appendix <ref> contains additional simulation results
showcasing performance for individual datasets as the chain length and the total
time elapsed varies.

§ APPLICATIONS

We now examine the benefits of on a diverse suite of regression and binary classification datasets. Table <ref>

summarizes the two synthetic and eight real-world datasets
considered, with further details in Appendix <ref>.

Synthetic and real-life datasets considered in
Section <ref>.

Dataset     $n$     $p$     Response Type

Borovecki     $31$     $22283$     Binary

Chin     $118$     $22215$     Binary

Chowdary     $104$     $22283$     Binary

Gordon     $181$     $12533$     Binary

Lymph     $148$     $4514$     Binary

Maize     $2266$     $98385$     Continuous

Malware     $373$     $503$     Binary

PCR     $60$     $22575$     Continuous

Synthetic Binary     $1000$     $50000$     Continuous

Synthetic Continuous     $1000$     $50000$     Binary

We first consider the Gordon microarray dataset [Gordon et al., 2002] with $n=181$ observations
(corresponding to a binary response vector indicating presence of lung cancer)
and $p = 12533$ covariates (corresponding to genes expression levels).
Figure <ref> shows the marginal posterior
probabilities estimated using samples from and the SOTA sampler for logistic and probit regression and
the Skinny Gibbs sampler for logistic regression, as well as the corresponding average runtimes per iteration. The marginal posterior probabilities $\pi_j\defeq \mathbb{P}_{\pi}(z_{j}=1)$ are estimated
by $\hat\pi_j=\frac{1}{I(T-S)} \sum_{i=1}^I \sum_{t=S+1}^{T} z^{(i)}_{j,t}$,
where $(\bm{z}_t^{(i)})_{t \geq 0}$ are samples from $i=1,\mydots,I$ independent Markov chains
generated using . We sample $I=5$ independent chains of length $T=5000$ iterations with a burn-in of $1000$ iterations for both and the SOTA sampler.
The average runtime per iteration with one standard error bars are calculated based on these independent chains.

Figure <ref> (Left) plots $\hat\pi_j$ against $j$ in the decreasing order of $\hat\pi_j$s.
It shows $\hat\pi_j$s based on samples from both and the SOTA sampler.

We simulate both and the SOTA sampler with the same random numbers at each iteration,
so that any differences will be due to numerical imprecision.
Figure <ref> (Left) shows that the estimates using  and the SOTA sampler are
indistinguishable. Furthermore, in this example all components of $z_t$ are identical between and the SOTA sampler
chains for all iterations $t$.
Despite producing Markov chains with indistinguishable marginal distributions and hence statistical properties,
Figure <ref> (Right) shows that has approximately $20$ and $6$ times faster runtime per iteration than SOTA for logistic and probit
regression respectively. Furthermore, for logistic regression has approximately $100$ times faster runtime per iteration than the Skinny Gibbs sampler. Overall, Figure <ref> highlights the practical value of over the SOTA sampler and the Skinny Gibbs sampler.

Comparing Bayesian logistic and probit regression with , SOTA, and Skinny Gibbs on the Gordon microarray dataset with $n=181$ observations and $p=12533$ covariates.
(Left and Middle) Marginal posterior probabilities $\mathbb{P}_{\pi}(z_{j}=1)$ estimated using samples from each chain: the recovered and SOTA probabilities are indistinguishable but differ significantly from the Skinny Gibbs probabilities.
(Right) Average runtime per sampler  iteration with one standard error bars.
See Section <ref> for details.

Figure <ref>, plotted with the y-axis on the log-scale, compares the runtimes of and the SOTA sampler for linear and probit regression on ten regression and binary classification datasets respectively.
It plots the average runtimes with one standard error bars

based on $10$ independent chains each of length $1000$ and $100$ for and the SOTA sampler respectively.
Figure <ref> shows that has lower runtimes per iteration compared to the SOTA sampler for all the datasets considered, with the most substantial speedups for larger datasets. For example, for the
Maize GWAS dataset [Romay et al., 2013, Liu et al., 2016, Zeng & Zhou, 2017]
with $n=2266$ observations (corresponding to average number of days taken for silk emergence
in different maize lines) and $p=98385$ covariates (corresponding to single nucleotide polymorphisms (SNPs) in the genome), requires $650$ms per iteration on average, which is $48$ times faster than the SOTA sampler requiring $31300$ms per iteration. For researchers, such speedups can reduce algorithm runtime from days to hours, giving substantial time and computational cost savings at no compromise to inferential quality.
Appendix <ref> contains additional results of applied to
these datasets including effective sample size (ESS) calculations, marginal posterior probabilities, and
performance under 10-fold cross-validation.

Average runtime per iteration with one standard error bars for and the SOTA sampler
for linear and probit regression applied to the ten continuous and binary response datasets.
See Section <ref> for details.

Average runtime per iteration with one standard error bars for and the SOTA sampler
for linear and probit regression applied to the ten continuous and binary response datasets.
See Section <ref> for details.

§ FURTHER WORK

The following questions arise from our work.

Extensions of to point-mass spike-and-slab
priors, as well as to non-Gaussian tails.
Whilst priors given in (<ref>) are one of the most common formulations employed in practice, a number of alternatives are available. This includes point-mass spike-and-slab priors
<cit.>, where a degenerate Dirac distribution about zero is chosen for the spike part,
and extensions which consider the heavier-tailed Laplace distribution
for the slab part instead of a Gaussian distribution
[Castillo et al., 2015, Ročková, 2018, Ray et al., 2020, Ray & Szabó, 2021].
An extension of would to be employ similar pre-computation
based strategy of Section <ref> to MCMC samplers for these alternative formulations.

Convergence complexity analysis of .
An important question that is not addressed in this article
is the number of iterations required for or other similar samplers
to converge to the their target posterior
distributions. For Gibbs samplers targeting posteriors
corresponding to continuous shrinkage priors
much theoretical progress has been made
[Pal & Khare, 2014, Qin & Hobert, 2019, Bhattacharya et al., 2022, Biswas et al., 2022].
Convergence of Gibbs samplers targeting spike-and-slab
posteriors has been less extensively studied and
requires more attention.

Diagnostics to assess the convergence of and
asymptotic variance of .
Given some time and computational budget constraints,
an immediate benefit of is the ability to run longer
Markov chains targeting spike-and-slab posteriors. This can
alleviate some concerns linked to burn-in and asymptotic
variance, but convergence and effective sample size diagnostics
[Johnson, 1998, Biswas et al., 2019, Vats & Knudson, 2021, Vehtari et al., 2021]
remain an important consideration particularly in high-dimensional
settings. We hope to investigate convergence diagnostics in
future work.

We thank Juan Shen for sharing the PCR and the Lymph Node datasets, Xiaolei Liu and Xiang Zhou for sharing the Maize GWAS dataset, and Marina Vannucci for helpful feedback. NB was supported by the NSF grant DMS-1844695, a GSAS Merit Fellowship, and a Two Sigma Fellowship Award. XLM was partially supported by the NSF grant DMS-1811308.

[Banerjee et al., 2021]
Banerjee, S., Castillo, I., and Ghosal, S.
Bayesian inference in high-dimensional models.
Springer volume on Data Science, 2021.

[Barbieri & Berger, 2004]
Barbieri, M. M. and Berger, J. O.
Optimal predictive model selection.
Annals of Statistics, 320 (3):0 870 – 897,
URL <https://doi.org/10.1214/009053604000000238>.

[Barbieri et al., 2021]
Barbieri, M. M., Berger, J. O., George, E. I., and Ročková, V.
The Median Probability Model and Correlated Variables.
Bayesian Analysis, 160 (4):0 1085 – 1112,
URL <https://doi.org/10.1214/20-BA1249>.

[Bhadra et al., 2019]
Bhadra, A., Datta, J., Polson, N. G., and Willard, B.
Lasso Meets Horseshoe: A Survey.
Statistical Science, 340 (3):0 405 – 427,
URL <https://doi.org/10.1214/19-STS700>.

[Bhattacharya et al., 2015]
Bhattacharya, A., Pati, D., Pillai, N. S., and Dunson, D. B.
Dirichlet–Laplace Priors for Optimal Shrinkage.
Journal of the American Statistical Association, 1100
(512):0 1479–1490, 2015.
URL <https://doi.org/10.1080/01621459.2014.960967>.

[Bhattacharya et al., 2016]
Bhattacharya, A., Chakraborty, A., and Mallick, B. K.
Fast sampling with Gaussian scale mixture priors in high-dimensional
Biometrika, 1030 (4):0 985–991, 2016.
ISSN 0006-3444.
URL <https://doi.org/10.1093/biomet/asw042>.

[Bhattacharya et al., 2022]
Bhattacharya, S., Khare, K., and Pal, S.
Geometric ergodicity of Gibbs samplers for the Horseshoe and its
regularized variants.
Electronic Journal of Statistics, 160 (1):0 1
– 57, 2022.
URL <https://doi.org/10.1214/21-EJS1932>.

[Biswas et al., 2019]
Biswas, N., Jacob, P. E., and Vanetti, P.
Estimating convergence of markov chains with l-lag couplings.
In Advances in Neural Information Processing Systems,
volume 32. Curran Associates, Inc., 2019.

[Biswas et al., 2022]
Biswas, N., Bhattacharya, A., Jacob, P. E., and Johndrow, J. E.
Coupling-based convergence assessment of some gibbs samplers for
high-dimensional bayesian regression with shrinkage priors.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 2022.
URL <https://doi.org/10.1111/rssb.12495>.

[Bogdan et al., 2015]
Bogdan, M., van den Berg, E., Sabatti, C., Su, W., and Candès, E. J.
SLOPE—Adaptive variable selection via convex optimization.
The Annals of Applied Statistics, 90 (3):0
1103 – 1140, 2015.
URL <https://doi.org/10.1214/15-AOAS842>.

[Carvalho et al., 2010]
Carvalho, C. M., Polson, N. G., and Scott, J. G.
The horseshoe estimator for sparse signals.
Biometrika, 970 (2):0 465–480, 04 2010.
ISSN 0006-3444.
URL <https://doi.org/10.1093/biomet/asq017>.

[Castillo & van der Vaart, 2012]
Castillo, I. and van der Vaart, A.
Needles and Straw in a Haystack: Posterior concentration for
possibly sparse sequences.
Annals of Statistics, 400 (4):0 2069 – 2101,
URL <https://doi.org/10.1214/12-AOS1029>.

[Castillo et al., 2015]
Castillo, I., Schmidt-Hieber, J., and van der Vaart, A.
Bayesian linear regression with sparse priors.
Annals of Statistics, 430 (5):0 1986–2018,
URL <https://doi.org/10.1214/15-AOS1334>.

[Dua & Graff, 2017]
Dua, D. and Graff, C.
UCI machine learning repository, 2017.
URL <http://archive.ics.uci.edu/ml>.

[Flegal et al., 2021]
Flegal, J. M., Hughes, J., Vats, D., Dai, N., Gupta, K., and Maji, U.
mcmcse: Monte Carlo Standard Errors for MCMC.
Riverside, CA, and Kanpur, India, 2021.
R package version 1.5-0.

[George & McCulloch, 1993]
George, E. I. and McCulloch, R. E.
Variable Selection via Gibbs Sampling.
Journal of the American Statistical Association, 880
(423):0 881–889, 1993.

[George & McCulloch, 1997]
George, E. I. and McCulloch, R. E.
Approaches for Bayesian Variable Selection.
Statistica Sinica, 70 (2):0 339–373, 1997.
ISSN 10170405, 19968507.
URL <http://www.jstor.org/stable/24306083>.

[Gordon et al., 2002]
Gordon, G. J. G., Jensen, R. V. R., Hsiao, L.-L. L., Gullans, S. R. S.,
Blumenstock, J. E. J., Ramaswamy, S. S., Richards, W. G. W., Sugarbaker, D.
J. D., and Bueno, R. R.
Translation of Microarray Data into Clinically Relevant Cancer
Diagnostic Tests Using Gene Expression Ratios in Lung Cancer and
Cancer Research, 620 (17):0 4963–4967,
September 2002.

[Guan & Stephens, 2011]
Guan, Y. and Stephens, M.
Bayesian variable selection regression for genome-wide association
studies and other large-scale problems.
The Annals of Applied Statistics, 50 (3):0
1780 – 1815, 2011.
URL <https://doi.org/10.1214/11-AOAS455>.

[Hager, 1989]
Hager, W. W.
Updating the inverse of a matrix.
SIAM Review, 310 (2):0 221–239, 1989.
ISSN 00361445.
URL <http://www.jstor.org/stable/2030425>.

[Hans et al., 2007]
Hans, C., Dobra, A., and West, M.
Shotgun Stochastic Search for “Large p” Regression.
Journal of the American Statistical Association, 1020
(478):0 507–516, 2007.
URL <https://doi.org/10.1198/016214507000000121>.

[Held & Holmes, 2006]
Held, L. and Holmes, C. C.
Bayesian auxiliary variable models for binary and multinomial
Bayesian Analysis, 10 (1):0 145 – 168, 2006.
URL <https://doi.org/10.1214/06-BA105>.

[Ishwaran & Rao, 2005]
Ishwaran, H. and Rao, J. S.
Spike and slab variable selection: Frequentist and Bayesian
Annals of Statistics, 330 (2):0 730 – 773,
URL <https://doi.org/10.1214/009053604000001147>.

[Johndrow et al., 2020]
Johndrow, J., Orenstein, P., and Bhattacharya, A.
Scalable Approximate MCMC Algorithms for the Horseshoe Prior.
Journal of Machine Learning Research, 210
(73):0 1–61, 2020.
URL <http://jmlr.org/papers/v21/19-536.html>.

[Johnson, 1998]
Johnson, V. E.
A coupling-regeneration scheme for diagnosing convergence in Markov
chain Monte Carlo algorithms.
Journal of the American Statistical Association, 930
(441):0 238–248, 1998.

[Johnson & Rossell, 2012]
Johnson, V. E. and Rossell, D.
Bayesian Model Selection in High-Dimensional Settings.
Journal of the American Statistical Association, 1070
(498):0 649–660, 2012.
URL <https://doi.org/10.1080/01621459.2012.682536>.

[Johnstone & Silverman, 2004]
Johnstone, I. M. and Silverman, B. W.
Needles and straw in haystacks: Empirical Bayes estimates of
possibly sparse sequences.
Annals of Statistics, 320 (4):0 1594 – 1649,
URL <https://doi.org/10.1214/009053604000000030>.

[Kelly, 2007]
Kelly, B. C.
Some Aspects of Measurement Error in Linear Regression of
Astronomical Data.
The Astrophysical Journal, 6650 (2):0
1489–1506, aug 2007.
URL <https://doi.org/10.1086/519947>.

[Liang et al., 2013]
Liang, F., Song, Q., and Yu, K.
Bayesian Subset Modeling for High-Dimensional Generalized Linear
Journal of the American Statistical Association, 1080
(502):0 589–606, 2013.
URL <https://doi.org/10.1080/01621459.2012.761942>.

[Liu et al., 2016]
Liu, X., Huang, M., Fan, B., Buckler, E. S., and Zhang, Z.
Iterative Usage of Fixed and Random Effect Models for Powerful and
Efficient Genome-Wide Association Studies.
PLOS Genetics, 120 (2):0 1–24, 2016.
URL <https://doi.org/10.1371/journal.pgen.1005767>.

[Mitchell & Beauchamp, 1988]
Mitchell, T. J. and Beauchamp, J. J.
Bayesian Variable Selection in Linear Regression.
Journal of the American Statistical Association, 830
(404):0 1023–1032, 1988.
ISSN 01621459.
URL <http://www.jstor.org/stable/2290129>.

[Narisetty & He, 2014]
Narisetty, N. N. and He, X.
Bayesian variable selection with shrinking and diffusing priors.
Annals of Statistics, 420 (2):0 789 – 817,
URL <https://doi.org/10.1214/14-AOS1207>.

[Narisetty et al., 2019]
Narisetty, N. N., Shen, J., and He, X.
Skinny Gibbs: A Consistent and Scalable Gibbs Sampler for Model
Journal of the American Statistical Association, 1140
(527):0 1205–1217, 2019.
URL <https://doi.org/10.1080/01621459.2018.1482754>.

[O'Brien & Dunson, 2004]
O'Brien, S. M. and Dunson, D. B.
Bayesian multivariate logistic regression.
Biometrics, 600 (3):0 739–746, 2004.
URL <https://doi:10.1111/j.0006-341X.2004.00224.x>.

[Pal & Khare, 2014]
Pal, S. and Khare, K.
Geometric ergodicity for Bayesian shrinkage models.
Electronic Journal of Statistics, 80 (1):0
604–645, 2014.
URL <https://doi.org/10.1214/14-EJS896>.

[Polson et al., 2013]
Polson, N. G., Scott, J. G., and Windle, J.
Bayesian Inference for Logistic Models Using Pólya–Gamma Latent
Journal of the American Statistical Association, 1080
(504):0 1339–1349, 2013.
URL <https://doi.org/10.1080/01621459.2013.829001>.

[Qin & Hobert, 2019]
Qin, Q. and Hobert, J. P.
Convergence complexity analysis of Albert and Chib’s algorithm for
Bayesian probit regression.
Annals of Statistics, 470 (4):0 2320–2347,
URL <https://doi.org/10.1214/18-AOS1749>.

[Ray & Szabó, 2021]
Ray, K. and Szabó, B.
Variational bayes for high-dimensional linear regression with sparse
Journal of the American Statistical Association, 00
(0):0 1–12, 2021.
URL <https://doi.org/10.1080/01621459.2020.1847121>.

[Ray et al., 2020]
Ray, K., Szabo, B., and Clara, G.
Spike and slab variational bayes for high dimensional logistic
In Advances in Neural Information Processing Systems,
volume 33, pp.  14423–14434. Curran Associates, Inc., 2020.

[Romay et al., 2013]
Romay, M. C., Millard, M. J., Glaubitz, J. C., Peiffer, J. A., Swarts, K. L.,
Casstevens, T. M., Elshire, R. J., Acharya, C. B., Mitchell, S. E.,
Flint-Garcia, S. A., McMullen, M. D., Holland, J. B., Buckler, E. S., and
Gardner, C. A.
Comprehensive genotyping of the USA national maize inbred seed
Genome Biology, 140 (6):0 R55, 2013.
URL <https://doi.org/10.1186/gb-2013-14-6-r55>.

[Ročková, 2018]
Ročková, V.
Bayesian estimation of sparse signals with a continuous
spike-and-slab prior.
Annals of Statistics, 460 (1):0 401 – 437,
URL <https://doi.org/10.1214/17-AOS1554>.

[Scott & Berger, 2010]
Scott, J. G. and Berger, J. O.
Bayes and empirical-Bayes multiplicity adjustment in the
variable-selection problem.
Annals of Statistics, 380 (5):0 2587 – 2619,
URL <https://doi.org/10.1214/10-AOS792>.

[Sereno, 2015]
Sereno, M.
A Bayesian approach to linear regression in astronomy.
Monthly Notices of the Royal Astronomical Society,
4550 (2):0 2149–2162, 11 2015.
ISSN 0035-8711.
URL <https://doi.org/10.1093/mnras/stv2374>.

[Tadesse & Vannucci, 2021]
Tadesse, M. G. and Vannucci, M.
Handbook of Bayesian Variable Selection.
Chapman and Hall/CRC, 2021.
URL <https://doi.org/10.1201/9781003089018>.

[Tibshirani, 1996]
Tibshirani, R.
Regression Shrinkage and Selection Via the Lasso.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 580 (1):0 267–288, 1996.
URL <https://doi.org/10.1111/j.2517-6161.1996.tb02080.x>.

[Titsias & Lázaro-Gredilla, 2011]
Titsias, M. and Lázaro-Gredilla, M.
Spike and slab variational inference for multi-task and multiple
kernel learning.
In Advances in Neural Information Processing Systems,
volume 24. Curran Associates, Inc., 2011.

[Vats & Knudson, 2021]
Vats, D. and Knudson, C.
Revisiting the Gelman–Rubin Diagnostic.
Statistical Science, 360 (4):0 518 – 529,
URL <https://doi.org/10.1214/20-STS812>.

[Vats et al., 2019]
Vats, D., Flegal, J. M., and Jones, G. L.
Multivariate output analysis for Markov chain Monte Carlo.
Biometrika, 1060 (2):0 321–337, 2019.
ISSN 0006-3444.
URL <https://doi.org/10.1093/biomet/asz002>.

[Vehtari et al., 2021]
Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P.-C.
Rank-Normalization, Folding, and Localization: An Improved
$\widehat{R}$ for Assessing Convergence of MCMC (with Discussion).
Bayesian Analysis, 160 (2):0 667 – 718, 2021.
URL <https://doi.org/10.1214/20-BA1221>.

[Yang et al., 2016]
Yang, Y., Wainwright, M. J., and Jordan, M. I.
On the computational complexity of high-dimensional Bayesian
variable selection.
Annals of Statistics, 440 (6):0 2497–2532,
URL <https://doi.org/10.1214/15-AOS1417>.

[Zeng & Zhou, 2017]
Zeng, P. and Zhou, X.
Non-parametric genetic prediction of complex traits with latent
Dirichlet process regression models.
Nature Communications, 80 (1):0 456, 2017.
URL <https://doi.org/10.1038/s41467-017-00470-2>.

[Zhou et al., 2013]
Zhou, X., Carbonetto, P., and Stephens, M.
Polygenic Modeling with Bayesian Sparse Linear Mixed Models.
PLOS Genetics, 90 (2):0 1–14, 02 2013.
URL <https://doi.org/10.1371/journal.pgen.1003264>.

[Zou & Hastie, 2005]
Zou, H. and Hastie, T.
Regularization and variable selection via the elastic net.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 670 (2):0 301–320, 2005.
URL <https://doi.org/10.1111/j.1467-9868.2005.00503.x>.

§ PROOFS

Consider Step $1$ of Algorithm <ref> and
Algorithm <ref> for probit regression. Given pre-computed matrices $\bm{\tilde{M}}_{\tau_0}$,
$\bm{\tilde{M}}_{\tau_0}$, $\bm{M}_{t-1}$, $\bm{M}_{t-1}^{-1}$
and state $\bm{z}_{t-1}$, calculating $\bm{M}_{t}$, $\bm{M}_{t}^{-1}$ requires
$\mathcal{O}(n^2 p_t)$ cost

by (<ref>) and (<ref>),
where $p_t = \min \{ \|\bm{z}_t\|_1, p-\|\bm{z}_t\|_1, \delta_t \}$ for $\delta_t = \|\bm{z}_t - \bm{z}_{t-1}\|_1$.

Consider Step $1$ of Algorithm <ref> for logistic regression. This requires $\mathcal{O}( \max\{n^2 p_t, n^3\} )$ cost, where the $\mathcal{O}(n^2p_t)$ cost arises from the calculation of $\bm{M}_t$ using (<ref>) and the $\mathcal{O}(n^3)$ cost arises from inverting $\bm{M}_t$ to calculate $\bm{M}_t^{-1}$.

Given $\bm{M}_{t}^{-1}$, Step $2$ of Algorithms <ref> and
<ref> then requires $\mathcal{O}(n p)$ cost, which arises from the matrix vector product $\bm{X}\bm{D}_t^{-\frac{1}{2}}r$ for $r \sim \mathcal{N}(0, \bm{I}_p)$ in Algorithm <ref>.
By component-wise independence, Step $3$ of Algorithms <ref> and
<ref> costs $\mathcal{O}(p)$ and
Step $4$ of Algorithm <ref> cost $\mathcal{O}(n)$.
Step $4$ of Algorithm <ref> and
Step $5$ of Algorithm <ref> for probit regression both cost $\mathcal{O}(1)$.
Step $5$ of Algorithm <ref> for logistic regression both costs $\mathcal{O}(n)$.

This gives an overall computational cost of $\mathcal{O}(\max \{n^2 p_t, np \})$ for
Algorithm <ref> and Algorithm <ref> for probit regression at
iteration $t$, and a cost of $\mathcal{O}(\max \{n^2 p_t, n^3, np \})$ for
<ref> for logistic regression.

By linearity, $\mathbb{E}[\delta_t] = \sum_{j=1}^p \mathbb{P}(z_{j,t} \neq z_{j,t-1})$
where for each component $j$ the random variables $z_{j,t}$ and $z_{j,t-1}$ are on $\{0,1\}$. For each $j$, we obtain

ℙ(z_j,t ≠z_j,t-1)    = ℙ(z_j,t=1, z_j,t-1=0) + ℙ(z_j,t=0, z_j,t-1=1)

= ( ℙ(z_j,t=1) - ℙ(z_j,t=1, z_j,t-1=1) ) +
( ℙ(z_j,t-1=1) - ℙ(z_j,t=1, z_j,t-1=1) )

= ( ℙ(z_j,t=1) - cov(z_j,t, z_j,t-1) - ℙ(z_j,t=1)ℙ(z_j,t-1=1) ) +

( ℙ(z_j,t-1=1) - cov(z_j,t, z_j,t-1) - ℙ(z_j,t=1)ℙ(z_j,t-1=1) )

= ℙ(z_j,t=1)ℙ(z_j,t-1=0) + ℙ(z_j,t=0)ℙ(z_j,t-1=1)
- 2 cov(z_j,t, z_j,t-1).

When $z_{j, t-1}$ follows the stationary $\pi$, $z_{j,t} \sim z_{j, t-1}$ and $\mathrm{var}(z_{j,t}) = \mathbb{P}(z_{j,t}=1)\mathbb{P}(z_{j,t-1}=0)$. Consequently,

ℙ(z_j,t ≠z_j,t-1) = 2 var_π(z_j,t) - 2 cov(z_j,t, z_j,t-1) = 2 var_π(z_j,t)(1 - corr_π(z_j,t, z_j,t-1)).

Note that $a^2=a$ if $a$ only takes the value $0$ and $1$. This gives

\begin{equation}\label{eq:identity}
{\rm I}\{z_{j, t}\not= z_{j, t-1}\}=(z_{j, t}-z_{j, t-1})^2=z_{j, t}+z_{j, t-1}-2z_{j, t}z_{j, t-1}.
\end{equation}

Let $J$ be the uniform random variable on the integers $\{1, \ldots, p\}$. Then,

\begin{equation}\label{eq:ptexact}
\delta_t = p \big( \EJ(z_{J, t})+\EJ(z_{J, t-1}) - 2 \EJ(z_{J, t}z_{J, t-1}) \big),
\end{equation}

where the expectation is taken with respect to the random index $J$. Note that
$\EJ(z_{J, t})=\|\bm{z}_t\|_1/p$ and $\VJ(z_{J, t})=(\|\bm{z}_t\|_1/p)(1-\|\bm{z}_t\|_1/p)=\tau_t^2/p^2$.
We obtain

\begin{align}
\EJ(z_{J, t}z_{J, t-1}) = \CJ(z_{J, t}, z_{J, t-1})+ \EJ(z_{J, t})\EJ(z_{J, t-1})
= \frac{\rho_t \tau_t \tau_{t-1} + \|\bm{z}_t\|_1 \|\bm{z}_{t-1}\|_1}{p^2},\label{eq:etexact}
\end{align}

where $\rho_t=\cJ(z_{J, t}, z_{J, t-1})$. Combining (<ref>)-(<ref>) yields (<ref>).

§ ALGORITHM DERIVATIONS

§.§ Linear regression with spike-and-slab priors

For linear regression with the continuous spike-and-slab priors in (<ref>),
the posterior density of $(\bm{\beta}, \bm{z}, \sigma^2) \in \mathbb{R}^p \times \{0,1\}^p \times (0,\infty)$ is given by

π(β, z, σ^2 | y) ∝     𝒩 (y; Xβ, σ^2) (σ^2 ; a_0/2, b_0/2 )

∏_j=1^p ( q 𝒩(β_j; 0, σ^2 τ_1^2) )^z_j ((1-q) 𝒩(β_j; 0, σ^2 τ_0^2))^1-z_j.

From (<ref>), we can calculate the conditional distributions. We have

π(β | z, σ^2, y)    ∝𝒩 (y; Xβ, σ^2 I_n) ∏_j=1^p 𝒩(β_j; 0, σ^2 τ_1^2)^z_j  𝒩(β_j; 0, σ^2 τ_0^2)^1-z_j

∝𝒩 (y; Xβ, σ^2 I_n) 𝒩(β; 0, σ^2 D^-1)  for  D (z τ_1^-2 + (1_p-z) τ_0^-2)

∝exp{ -1/2σ^2 I_n ( β^⊤X^⊤X β - 2β^⊤X^⊤y + β^⊤D β ) }

∝𝒩(β; Σ^-1 X^⊤y, σ^2 Σ^-1)  for  Σ = X^⊤X + D,

π(z | β, σ^2, y)    ∝∏_j=1^p ( q 𝒩(β_j; 0, σ^2 τ_1^2) )^z_j ((1-q) 𝒩(β_j; 0, σ^2 τ_0^2))^1-z_j

∝∏_j=1^p ( z_i; q 𝒩(β_j; 0, σ^2 τ_1^2)/q 𝒩(β_j; 0, σ^2 τ_1^2) + (1-q) 𝒩(β_j; 0, σ^2 τ_0^2) ),  and

π(σ^2 | β, z, y)    ∝𝒩 (y; Xβ, σ^2) (σ^2 ; a_0/2, b_0/2 ) 𝒩 (β ; 0, σ^2 D^-1 )

∝(1/σ^2)^n/2 exp{ -1/2 σ^2  y - X β _2^2 } ( 1/σ^2)^a_0/2 exp{ -1/2 σ^2 b_0 } (1/σ^2)^p/2 exp{-1/2 σ^2 β^⊤D β }

∝( σ^2 ; a_0+n+p/2, b_0 +  y - X β _2^2 + β^⊤D β/2 )

as given in Algorithm <ref>.

§.§ Probit regression with spike-and-slab priors

Consider the probit regression likelihood, where for each observation $i=1,\mydots,n$,
$\mathbb{P}( y_i = 1 | \bm{x_i}, \bm{\beta}) = 1 - \mathbb{P}( y_i = 0 | \bm{x_i}, \bm{\beta}) =
\Phi( \bm{x_i}^\top \bm{\beta} )$
for $\bm{x_i}^\top$ the $i$-th row of the design matrix $\bm{X}$.
and $\Phi$ that cumulative density function of a univariate Normal distribution.
We obtain $y_i = \mathrm{I} \{ \tilde{y}_i>0 \}$ for
$\tilde{y}_i | \bm{\beta} \sim \mathcal{N}(\bm{x_i}^\top \bm{\beta}, 1)$.
The Bayesian probit regression model is then given by

z_j    i.i.d.∼ (q)    for all  j=1,,p

β_j | z_j      ind  ∼ (1-z_j) 𝒩(0, τ_0^2) + z_j 𝒩(0, τ_1^2)
for all  j=1,,p

ỹ_i | β      ind  ∼ 𝒩(x_i^⊤β, 1)
for all  i=1,,n

y_i           = I{ ỹ_i > 0 }    for all  i=1,,n.
For the prior and likelihood in (<ref>), the posterior density
of $(\bm{\beta}, \bm{z}, \bm{\tilde{y}}) \in \mathbb{R}^p \times \{0,1\}^p \times \mathbb{R}^p$ is given by

π(β, z, ỹ | y) ∝    ( ∏_i=1^n I { I {ỹ_i>0} = y_i } 𝒩 (ỹ_i; x_i^⊤β, 1) )

∏_j=1^p ( q 𝒩(β_j; 0, τ_1^2) )^z_j ((1-q) 𝒩(β_j; 0, τ_0^2))^1-z_j.

From (<ref>), we can calculate the conditional distributions. We obtain

π(β | z, ỹ, y)     ∝𝒩(ỹ; X β,  I_n) 𝒩(β; 0, D^-1)  for  D (z τ_1^-2 + (1_p-z) τ_0^-2)

∝𝒩(β; Σ^-1 X^⊤ỹ, Σ^-1)  for  Σ = X^⊤X + D,

π(z | β, ỹ, y)    ∝∏_j=1^p ( z_i; q 𝒩(β_j; 0, τ_1^2)/q 𝒩(β_j; 0, τ_1^2) + (1-q) 𝒩(β_j; 0,  τ_0^2) ),   and

π(ỹ | β, z, y)    ∝∏_i=1^n 𝒩 (ỹ_i; x_i^⊤β, 1) I { I {ỹ_i>0} = y_i }.

as required for probit regression in Algorithm <ref>.

§.§ Logistic regression with spike-and-slab priors

We first describe the Bayesian logistic regression model considered.
Consider the logistic regression likelihood, where for each observation $i=1,\mydots,n$,
$\mathbb{P}( y_i = 1 | \bm{x_i}, \bm{\beta}) = 1 - \mathbb{P}( y_i = 0 | \bm{x_i}, \bm{\beta}) = \frac{\exp(\bm{x_i}^\top \bm{\beta})}{1+\exp(\bm{x_i}^\top \bm{\beta})}$
for $\bm{x_i}^\top$ the $i$-th row of the design matrix $\bm{X}$.
We obtain $y_i = \mathrm{I} \{ \tilde{y}_i>0 \}$ where $\tilde{y}_i \overset{ind}{\sim} \mathrm{Logistic}(\bm{x_i}^\top \bm{\beta}, 1)$, corresponding to the logistic distribution centered about $\bm{x_i}^\top \bm{\beta}$ and scale parameter $1$.

§.§.§ Student's $t$-distribution based approximation of the logistic regression likelihood.

Following O'Brien & Dunson, 2004 and Narisetty et al., 2019, we can approximate $\mathrm{Logistic}(\bm{x_i}^\top \bm{\beta}, 1)$ with $\bm{x_i}^\top \bm{\beta} + w t_{\nu}$, where $t_{\nu}$ denotes a $t$-distribution with $\nu$ degrees of freedom and $w$ is a multiplicative factor. The constants $\nu \defeq 7.3$ and $w^2 \defeq \frac{\pi^2(\nu-2)}{3 \nu}$ are chosen following O'Brien & Dunson, 2004,
in order to match the variance of the logistic distribution and to minimize the integrated squared distance between the respective densities. The Gaussian scale representation of this $t$-distribution is

\begin{equation}
\tilde{y}_i | \bm{x_i}, \bm{\beta}, \tilde{\sigma}_i  \sim \mathcal{N}(\bm{x_i}^\top \bm{\beta}, \tilde{\sigma}_i^2),  \quad \tilde{\sigma}_i^2 \sim \InvGamma \Big(\frac{v}{2}, \frac{w^2 v}{2} \Big), \label{eq:logistic_approx}
\end{equation}

where each $\tilde{\sigma}_i^2$ is an augmented variable. The Bayesian logistic regression model is then given by

z_j    i.i.d.∼ (q)    for all  j=1,,p

β_j | z_j      ind  ∼ (1-z_j) 𝒩(0, τ_0^2) + z_j 𝒩(0, τ_1^2)
for all  j=1,,p

σ̃^2_i    i.i.d.∼ (ν/2, w^2ν/2)

ỹ_i | β, σ̃^2_i      ind  ∼ 𝒩(x_i^⊤β, σ̃_i^2)
for all  i=1,,n

y_i           = I{ ỹ_i > 0 }    for all  i=1,,n.
Let $\bm{\tilde{\sigma}}^2$ denote the vector with entries $\tilde{\sigma}_i^2$ for $i=1,\mydots,n$. For the  prior and likelihood in (<ref>), the posterior density of $(\bm{\beta}, \bm{z}, \bm{\tilde{y}}, \bm{\tilde{\sigma}}^2)$ on $\mathbb{R}^p \times \{0,1\}^p \times \mathbb{R}^n \times (0,\infty)^n$ is given by

π(β, z, ỹ, σ̃^2 | y)
∝    ∏_j=1^p ( q 𝒩(β_j; 0, τ_1^2) )^z_j ((1-q) 𝒩(β_j; 0, τ_0^2))^1-z_j

∏_i=1^n I { I {ỹ_i>0} = y_i } 𝒩 (ỹ_i; x_i^⊤β, σ̃_i^2) (σ̃_i^2 ; ν/2, w^2 ν/2 ).

From (<ref>), we can calculate the conditional distributions. Let $\bm{W}=\Diag(\bm{\tilde{\sigma}}^2)$. We obtain

π(β | z, ỹ, σ̃^2, y)     ∝𝒩(ỹ; X β, W) 𝒩(β; 0, D^-1)  for  D (z τ_1^-2 + (1_p-z) τ_0^-2)

∝𝒩(β; Σ^-1 X^⊤W^-1 ỹ, Σ^-1)  for  Σ = X^⊤W^-1 X + D,

π(z | β, ỹ, σ̃^2, y)    ∝∏_j=1^p ( z_i; q 𝒩(β_j; 0, τ_1^2)/q 𝒩(β_j; 0,  τ_1^2) + (1-q) 𝒩(β_j; 0, τ_0^2) ),

π(ỹ | β, z, σ̃^2, y)    ∝∏_i=1^n 𝒩 (ỹ_i; x_i^⊤β, σ̃_i^2) I { I {ỹ_i>0} = y_i }  ,  and

π(σ̃^2 | β, z, ỹ, y)    ∝∏_i=1^n 𝒩 (ỹ_i; x_i^⊤β, σ̃_i^2)
(σ̃_i^2 ; ν/2, w^2 ν/2 )

∝∏_i=1^n ( σ̃_i^2 ; ν+1/2, w^2 ν+ (ỹ_i - x_i^⊤β )^2/2 ).

as required for logistic regression in Algorithm <ref>.

A scalable Gibbs sampler for logistic regression.
The computational bottleneck of existing Gibbs samplers for logistic regression is linked to sampling from the full conditional of $\bm{\beta} \in \mathbb{R}^p$. This is given by

\begin{equation} \label{eq:full_conditional_beta_logistic}
\bm{\beta}_{t+1} | \bm{z}_t, \bm{\tilde{\sigma}}^2_t \sim
\mathcal{N} \big(\bm{\Sigma}^{-1}_t \bm{X}^\top \bm{W}_t^{-1} \bm{\tilde{y}}, \bm{\Sigma}^{-1}_t \big)\quad {\rm for}\quad  \bm{\Sigma}_t = \bm{X}^\top \bm{W}_t^{-1} \bm{X} + \bm{D}_t,
\end{equation}

where $t$ indexes the iteration of the Markov chain,
$\bm{W}_t$ is the diagonal matrix with the vector
$\bm{\tilde{\sigma}}^2_t$ populating its diagonal elements, and $\bm{D}_t$ is the diagonal matrix with the vector $\bm{z}_t \tau^{-2}_1 + (\bm{\mathrm{1}}_p-\bm{z}_t) \tau^{-2}_0$ populating its diagonal elements. To sample from (<ref>), we can use
the $\Omega(n^2p)$ sampler of Bhattacharya et al., 2016, which is given in Algorithm <ref>.

An $\Omega(n^2p)$ sampler of (<ref>)
[Bhattacharya et al., 2016]

Sample $\bm{r} \sim \mathcal{N}(0, \bm{I}_p)$, $\bm{\xi} \sim \mathcal{N}(0,\bm{I}_n)$.
Set $\bm{u} = \bm{D}_t^{-\frac{1}{2}}\bm{r}$ and calculate $\bm{v} = \bm{W}_t^{-1/2}\bm{X}\bm{u} + \bm{\xi}$.
Set $\bm{v}^* = \bm{M}_t^{-1} ( \bm{W}_t^{-1/2} \bm{\tilde{y}} - \bm{v})$ for $\bm{M}_t = \bm{I}_n + \bm{W}_t^{-1/2} \bm{X} \bm{D}_t^{-1} \bm{X}^\top \bm{W}_t^{-1/2}$.
Return $\bm{\beta} = \bm{u} +  \bm{D}_t^{-1} \bm{X}^\top \bm{W}_t^{-1/2} \bm{v}^*$.

Following the strategy in Section <ref>, for logistic regression uses pre-computation to reduce the computational cost of Algorithm <ref>. Using the notation from Section <ref> with $\bm{M}_t \defeq \bm{I}_n + \bm{W}_t^{-1/2} \bm{X} \bm{D}_t^{-1} \bm{X}^\top \bm{W}_t^{-1/2}$, we note

M_t    = I_n +
W_t^-1/2 (

M̃_τ_0 - I_n
+ (τ_1^2-τ_0^2)X_A^c_t X_A^c_t^T
) W_t^-1/2

= I_n +
W_t^-1/2 (

M̃_τ_1 - I_n
+ (τ_0^2-τ_1^2)X_A^c_t X_A^c_t^T
) W_t^-1/2

= I_n +
W_t^-1/2 (
W_t-1^1/2(M_t-1 - I_n)W_t-1^1/2 + X_Δ_t C_Δ_t X_Δ_t^T
) W_t^-1/2 .

In (<ref>) – (<ref>),
calculating the matrix products
$\bm{X}_{A_t} \bm{X}_{A_t}^\top$, $\bm{X}_{A^c_t} \bm{X}_{A^c_t}^\top$, and $\bm{X}_{\Delta_t}  \bm{C}_{\Delta_t} \bm{X}_{\Delta_t}^\top$ requires
$\mathcal{O}(n^2 \|\bm{z}_t\|_1)$, $\mathcal{O}(n^2 (p-\|\bm{z}_t\|_1))$, and
$\mathcal{O}(n^2 \delta_t)$ cost respectively. Given
$\bm{\tilde{M}}_{\tau_0}$, $\bm{\tilde{M}}_{\tau_1}$, $\bm{M}_{t-1}$, and $\bm{z}_{t-1}$,
we evaluate whichever matrix product in
(<ref>) – (<ref>) has minimal computational cost and thereby calculate $\bm{M}_t$ at the reduced cost of
$\mathcal{O}(n^2 p_t)$ where $p_t \defeq \min \{ \|\bm{z}_t\|_1, p-\|\bm{z}_t\|_1, \delta_t \}$.

To calculate $\bm{M}_t^{-1}$, we calculate $\bm{M}_t^{-1}$ by directly inverting the calculated matrix $\bm{M}_t$, which requires $\mathcal{O}(n^3)$ cost. Overall, this strategy

reduces the computational cost of calculating the matrices $\bm{M}_{t}$ and $\bm{M}_{t}^{-1}$ from
$\Omega(n^2p)$ to $\mathcal{O}(\max\{n^2 p_t,n^3\})$.

Extensions to Scalable Spike-and-Slab for logistic regression.
Suppose the matrices $\bm{X}^\top \bm{X}$ is pre-computed. This initial step requires $\mathcal{O}(n p^2)$ computational cost and $\mathcal{O}(p^2)$ memory. Then the matrices $\bm{X}_{A_t}^\top \bm{X}_{A_t}$ and $\bm{X}_{A^c_t}^\top \bm{X}_{A^c_t}$ in (<ref>) – (<ref>) correspond to pre-computed sub-matrices of $\bm{X}^\top \bm{X}$, and
calculating $\bm{M}_t$ using (<ref>) – (<ref>)
each iteration $t$ involves matrix addition and diagonal matrix multiplication which only requires $\mathcal{O}(n^2)$ cost. To sample from (<ref>), we calculate $\bm{M}_t^{-1}$ by directly inverting the calculated matrix $\bm{M}_t$ from (<ref>), which requires $\mathcal{O}(n^3)$ cost.

Overall, now the Gibbs samplers for logistic regression requires $\mathcal{O}(\max\{ n^3, np \})$ computational cost at iteration $t$,
which is an improvement compared to .

§ EXPERIMENT DETAILS

Figure <ref> of Section <ref>.
In Figure <ref>, we use the same prior
hyperparameters for all the algorithms. Following Narisetty et al., 2019,
we choose $\tau_0^2 = \frac{1}{n}$, $\tau_1^2 = \max \{ \frac{p^{2.1}}{100n},1 \}$ and
$q = \mathbb{P}(z_j=1)$ such that
$\mathbb{P}(\sum_{j=1}^p \mathrm{I}\{z_j=1\} > K)=0.1$ for
$K=\max \{ 10, \log n \}$.

The true positive rate (TPR) and the false discovery rate (FDR) correspond to the proportion of non-zero and zero components of $\beta^*_j$ that are correctly selected respectively. They are calculated as
$\frac{1}{s} \sum_{j=1}^s \mathrm{I} \{ \mathbb{P}_{\pi}(z_{j}=1) > 0.5 \}$
and $ \frac{1}{p-s} \sum_{j=s+1}^p \mathrm{I} \{ \mathbb{P}_{\pi}(z_{j}=1) > 0.5 \}$ respectively, where the marginal posterior probabilities $\pi_j \defeq \mathbb{P}_{\pi}(z_{j}=1)$ are estimated by $\hat{\pi}_j \defeq \frac{1}{4000} \sum_{t=\mathrm{1001}}^{\mathrm{5000}} z_{j,t}$
for sample points $(\bm{z}_{t})_{t \geq 0}$ generated using or Skinny Gibbs.
The lines in Figure <ref> correspond to the average TPR and FDR across $20$ independently generated datasets, and the grey bands correspond to one standard error of the averages.

§ DATASET DETAILS

Synthetic continuous response dataset in Section <ref>.
In Figure <ref>, synthetic linear regression datasets
are considered. For number of observations $n$ and number of covariates $p$,
we generate a design matrix $\bm{X} \in \mathbb{R}^{n \times p}$ such that
each $[\bm{X}]_{i,j} \overset{i.i.d.}{\sim} \mathcal{N}(0, 1)$ for all
$1 \leq i \leq n$ and $1 \leq j \leq p$, which is then scaled to
ensure each column has a mean of $0$ and a standard error of $1$.
We choose the true signal $\bm{\beta}^* \in \mathbb{R}^p$ such that
$\bm{\beta}^*_j = 2 \mathrm{I}\{ j \leq s \}$, where $s$ is the
sparsity parameter corresponding to the number of non-zero components.
Given $\bm{X}$ and $\bm{\beta}^*$, we generate $\bm{y} = \bm{X}\bm{\beta}^* + \sigma^* \bm{\epsilon}$ for
$\bm{\epsilon} \sim \mathcal{N}(0,\bm{I}_n)$, where
$\sigma^*=2$ is the Gaussian noise standard deviation.

Synthetic binary response dataset in
Section <ref>.
In Figures <ref> and <ref>,
synthetic binary classification datasets are considered.
For number of observations $n$ and number of covariates $p$,
we generate a design matrix $\bm{X} \in \mathbb{R}^{n \times p}$ such that
each $[\bm{X}]_{i,j} \overset{i.i.d.}{\sim} \mathcal{N}(0, 1)$ for all
$1 \leq i \leq n$ and $1 \leq j \leq p$, which is then scaled to
ensure each column has a mean of $0$ and a standard error of $1$.
We choose the true signal $\bm{\beta}^* \in \mathbb{R}^p$ such that
$\bm{\beta}^*_j = 2^{\frac{9-j}{4}} \mathrm{I}\{ j \leq s \}$,

where $s$ is the sparsity parameter corresponding to the number of non-zero components.
Given $\bm{X}$ and $\bm{\beta}^*$, we generate $y_i = \mathrm{I}\{ \tilde{y}_i > 0\}$ for
$\tilde{y}_i \sim \mathrm{Logistic}(\bm{x_i}^\top \bm{\beta}^*, 1)$ for $i=1\mydots,n$,
where $\bm{x_i}^\top$ is the $i$-th row of $\bm{X}$ and
$\mathrm{Logistic}(\bm{x_i}^\top \bm{\beta}^*, 1)$ is the Logistic distribution
with mean $\bm{x_i}^\top \bm{\beta}^*$ and scale parameter $1$.

Datasets in Section <ref>.

The Malware detection dataset from the UCI machine learning repository
[Dua & Graff, 2017] has $n=373$ observations with binary responses and
$p=503$ covariates, and is publicly available on

The Borovecki, Chowdary, Chin and Gordon datasets are all high-dimensional
microarray datasets. They are publicly available
on the $\mathrm{datamicroarray}$ package in $\mathrm{R}$.
The Borovecki dataset has $n=31$ observations with binary responses and $p=22283$
The Chowdary dataset has $n=104$ observations with binary responses
and $p=22283$ covariates.
The Chin dataset has $n=118$ observations with binary
responses and $p=22215$ covariates.
The Gordon dataset has $n=181$ observations with binary
responses and $p=12533$ covariates.

The PCR GWAS dataset has $n=60$ observations with continuous responses
and $p=22575$ covariates, and is publicly available on
<www.ncbi.nlm.nih.gov/geo> (accession number $GSE3330$).
The Lymph Node GWAS dataset has $n=148$ observations with binary responses
and $p=4514$ covariates, and has been previously considered
[Hans et al., 2007, Liang et al., 2013, Narisetty et al., 2019].
The Maize GWAS dataset has $n=2266$ observations with continuous responses
and $p=98385$ covariates, and has been previously considered
[Romay et al., 2013, Liu et al., 2016, Zeng & Zhou, 2017].
The Lymph Node GWAS and the Maize GWAS datasets are not publicly available.

The synthetic continuous dataset has $n=1000$ observations
and $p=50000$ covariates. The design matrix $\bm{X}$ is generates such that
each $[\bm{X}]_{i,j} \overset{i.i.d.}{\sim} \mathcal{N}(0, 1)$ for all
$1 \leq i \leq n$ and $1 \leq j \leq p$, which is then scaled to
ensure each column has a mean of $0$ and a standard error of $1$.
The true signal $\bm{\beta}^* \in \mathbb{R}^p$ is chosen such that
$\beta^*_j = 2^{\frac{9-j}{4}}\mathrm{I}\{ j \leq s \}$,

where $s$ is the sparsity parameter corresponding to the number of non-zero components.
Given $\bm{X}$ and $\bm{\beta}^*$, we generate $\bm{y} = \bm{X}\bm{\beta}^* + \sigma^* \bm{\epsilon}$ for
$\bm{\epsilon} \sim \mathcal{N}(0,\bm{I}_n)$, where
$\sigma^*=2$ is the Gaussian noise standard deviation.
The synthetic binary classification dataset is generated as in
Section <ref>, with $n=1000$ observations
and $p=50000$ covariates.

§ ADDITIONAL EXPERIMENTS

Variable selection performance as a function of time or number of iterations.

Avrage rue positive rate (TPR) and false discovery rate (FDR) plotted against
the number of iterations and the total time elapsed in seconds. We consider
, SOTA, and the Skinny Gibbs approximate sampler [Narisetty et al., 2019]
applied to a synthetic binary classification dataset with $n=200$ observations and
$p=1000$ covariates. The TPR and FDR are averaged over
$10$ independent chains, and one standard error bars are shown on the left and center-left plots and omitted on the right and center-right plots for visibility.
The SOTA sampler is omitted from the Left and Center-Left plots as its output has the same marginal
distribution and statistical performance as .

See Section <ref> for details.

Figure <ref> plots the average true positive rate (TPR)
and the false discovery rate (FDR) of variable selection based on samples
from and Skinny Gibbs as the
length of the chains are varied. The TPR and FDR are averaged over
$10$ independent chains, and one standard error bars are shown.
We consider a synthetic binary classification dataset
generated using a logistic regression model as in Section <ref>,
with $n=200$ observations, $p=1000$ covariates, sparsity $s=10$, and an
exponentially decaying sparse true signal $\bm{\beta}^* \in \mathbb{R}^p$ such that
$\beta^*_j = 2^{\frac{9-j}{4}}$ for $j \leq s$ and $\beta^*_j=0$ for $j>s$.
The TPR and FDR are calculated as in Section <ref> with
the marginal posterior probabilities $\pi_j \defeq \mathbb{P}_{\pi}(z_{j}=1)$
now estimated by $\hat{\pi}_j \defeq \frac{1}{T-999} \sum_{t=\mathrm{1000}}^{\mathrm{T}} z_{j,t}$
for a burn-in of $1000$, a varying chain length $T \geq 1000$, and sample points $(\bm{z}_{t})_{t \geq 0}$ generated using or Skinny Gibbs. We use the same prior
hyperparameters for all the algorithms, which are chosen
according to Narisetty et al., 2019.

Figure <ref> Left and Center-Left plot the
TPR and FDR against the chain length $T$. It shows that for both logistic
and probit regression have higher TPR and lower FDR than Skinny Gibbs for all chain
lengths. Furthermore, for logistic regression has higher TPR and lower FDR than for probit regression,
which is expected as the synthetic dataset for this example is generated using a
logistic regression model. The SOTA sampler is omitted from the Left and Center-Left plots as its output has the same marginal
distribution and statistical performance as .

Figure <ref> Center-Right and Right
plot the TPR and FDR against total time elapsed in seconds to generate samples
using , SOTA, or Skinny Gibbs chains with a burn-in of $1000$ iterations. The standard error bars are now omitted for better visibility.
For each time budget, we observe better variable selection performance from when compared with the slower SOTA implementation or with Skinny Gibbs.

Effective Sample Size of for the datasets in Section <ref>.
Figure <ref> shows the Effective Sample Size per iteration
and per unit of time (in seconds) of and the SOTA sampler for the datasets in Section <ref>.
The ESS is calculated using the $\mathrm{mcmcse}$ package [Flegal et al., 2021, Vats et al., 2019]
for one chain of length $10000$ iterations with a burn-in of $1000$ iterations for each dataset.
The average ESS of the $\beta$ components are then plotted.
Figure <ref> Right shows that has significantly higher ESS per second compared
to the corresponding SOTA sampler for all the datasets considered.

Effective sample size (ESS) per iteration and per second of and the SOTA sampler for some of the datasets in Section <ref>. The
ESS is calculated using one chain of length $10000$ iterations with a burn-in of $1000$ iterations.
The ESS per iteration of the SOTA sampler is omitted from the Left plot as it implements the same Gibbs sampler as .

Performance metrics for the datasets in Section <ref>.
Figures <ref> – <ref> show various performance metrics of
for some of the datasets considered in Section <ref>.
Figures <ref> – <ref> (Left) plot the marginal posterior
probability estimates $\hat\pi_j$ against $j$ in the decreasing order of $\hat\pi_j$s,
following the setup in Figure <ref>.
For datasets with continuous valued responses, $\hat\pi_j$s are based on
samples from for linear regression.
For datasets with binary valued responses, $\hat\pi_j$s are based on
samples from for logistic and probit regression, and the Skinny Gibbs
sampler from logistic regression.
We use samples from $5$ independent
chains of length $10000$ iterations with a burn-in of $1000$ iterations.

Estimates based on samples from the SOTA sampler are not shown,
as they implement the same Gibbs sampler as (other than possible numerical discrepancies, as discussed in Section

Figures <ref> – <ref> (Center) show the
average time taken per iteration with one standard error bars for , the SOTA sampler,
and the Skinny Gibbs sampler based on $5$ independent chains of length $10000$ iterations.

Figures <ref> – <ref> (Right) show the 10-fold
cross-validation average root-mean-square error (RMSE) against the total time elapsed to run
one and one SOTA chain.
To compute this evaluation, we partition the observed dataset into $10$ folds uniformly at random and, for each fold $k$, run a chain conditioned on all data outside of fold $k$ and evaluate its performance on the held-out data in the $k$-th fold.
The average RMSE is calculated as
$ \frac{1}{10} \sum_{k=1}^{10} r_k$, where $r_k$ is the RMSE for the
$k^{th}$ fold. For datasets with continuous valued responses, the
quantities $r_k$ for linear regression are calculated as
$(\frac{1}{|D_k|} \sum_{i \in D_k} (y_i - \hat{y}_i)^2)^{1/2}$ where
$D_k$ is the $k^{th}$ fold, $\hat{y}_i \defeq \frac{1}{T-1000} \sum_{t=1001}^{T} x_i^T \beta_{t}$
are the predicted responses, and $(\beta_{t})_{t \geq 0}$ are samples from and SOTA
targeting the posterior distribution of the $k^{th}$ training set.
For datasets with binary valued responses, the quantities $r_k$
are calculated as $(\frac{1}{|D_k|} \sum_{i \in D_k} (y_i - \hat{p}_i)^2)^{1/2}$,
where $D_k$ is the $k^{th}$ fold,
$\hat{p}_i \defeq \frac{1}{T-1000} \sum_{t=1001}^{T} \mathrm{Logistic}(x_i^T \beta_{t})$
and $\hat{p}_i \defeq \frac{1}{T-1000} \sum_{t=1001}^{T} \Phi(x_i^T \beta_{t})$ are the
predicted probabilities for logistic and probit regression respectively, and $(\beta_{t})_{t \geq 0}$ are samples
from and SOTA targeting the posterior distribution of the $k^{th}$ training set.
Figures <ref> – <ref> (Right) plot the average RMSE
against total time elapsed in seconds to generate samples using or SOTA chains
with a burn-in of $1000$ iterations.
The RMSE of the Skinny Gibbs sampler is not available,
as the $\mathrm{skinnybasad}$ package does not output the full chain trajectories required
for RMSE calculations.

Borovecki dataset with $n=31$ observations, $p=22283$ covariates and binary valued responses.

Chin dataset with $n=118$ observations, $p=22215$ covariates and binary valued responses.

Chowdary dataset with $n=104$ observations, $p=22283$ covariates and binary valued responses.

Gordon dataset with $n=181$ observations, $p=12533$ covariates and binary valued responses.

Lymph dataset with $n=148$ observations, $p=4514$ covariates and binary valued responses.

Malware dataset with $n=373$ observations, $p=503$ covariates and binary valued responses.

PCR dataset with $n=60$ observations, $p=22575$ covariates and continuous valued responses.