##  1 Introduction

There is a large literature on Bayesian sparse models and the inspired
algorithms for estimating the posterior. Originally motivated for solving the
variable selection problem in linear regression, spike-and-slab priors (or,
“discrete spike-and-slab”, [tadesse2021handbook]) have the marginal form of a
two-component mixture for each parameter element: one component (spike) from a
point mass at zero, and the other (slab) from a continuous distribution
(mitchell1988bayesian); the continuous elements can be independent or
dependent a priori. For linear regression with Gaussian errors, very efficient
Markov chain Monte Carlo (MCMC) algorithms have been developed. When the slab
prior distribution follows a Gaussian, the Stochastic Search Variable
Selection (SSVS) algorithm (george1995stochastic) exploits the posterior
conjugacy and samples from the marginal posterior of the binary inclusion
variables. SSVS is a Gibbs sampler that tries to flip each binary inclusion
variable one at a time from a Bernoulli full conditional distribution; at the
end of each iteration, it samples the regression coefficients given the vector
of binary inclusion variables. As an alternative to using the marginal
posterior, the Orthogonal Data Augmentation (ODA) algorithm (ghosh2011rao)
introduces an augmented design matrix (along with augmented responses) to
append the observed design matrix, such that the Gram matrix becomes diagonal
hence easily invertible, enabling the use of two-block Gibbs sampler. The ODA
algorithm can be extended to some generalized linear models if the latent
Gaussian has a fixed variance, such as probit regression (albert1993bayesian).
Further, for design matrix that is high-dimensional or contains highly
correlated predictors, various algorithms have been proposed such as shotgun
algorithm (hans2007shotgun), parallel tempering (bottolo2010evolutionary),
correlation-based search (kwon2011efficient), two-parameter flipping
Metropolis-Hastings algorithm under $g$-prior slab (yang2016computational). In
the meantime, there is a comparably vast literature on continuous shrinkage
priors with excellent performance for the task of variable selection
(george1995stochastic; rovckova2018spike; polson2010shrink;
carvalho2010horseshoe; piironen2017sparsity; armagan2013generalized;
bhattacharya2015dirichlet; bai2018beta). Since our focus is on the posterior
with exact zeros, for brevity, we will skip the detail of continuous
shrinkage.

With the rich literature, there is a recent interest in structured sparsity
(hoff2017lasso; griffin2023structured) that has inspired new extensions of
sparsity priors. Specifically, the sparsity is “structured” in the sense that:
(i) the occurrences of zeros could be dependent, according to some temporal,
spatial, or group structure; (ii) the non-zeros could have some correlation
structure, such as smoothness over some spatial domain. We now provide a few
examples. In the task of change-point detection, one may model a time series
as having the mean increments to be sparse over a continuous period of time so
that the mean function would become a step function (tibshirani2005sparsity;
betancourt2017bayesian). In the scalar-on-image regression, one may model the
regression coefficients to be spatially smooth for those non-zeros, while
being zero over several continuous regions (kang2018scalar). For these models,
a critical computational challenge arises that the above existing algorithms
for spike-and-slab priors cannot be applied directly here, due to the lack of
conjugacy or fixed variance for latent Gaussian. Because of the correlations
among the elements of the parameter, updating the full conditional of one
element at a time suffers from slow mixing of Markov chains. To facilitate the
model development and posterior estimation under structured sparsity,
xu2020bayesian propose the L1-ball priors: starting with a continuous random
variable $\beta\in\mathbb{R}^{p}$ from a “precursor” distribution (that could
be correlated or smooth over some input space), one projects $\beta$ to the
boundary of the L1 ball. This projection produces $\theta$ containing exact
zeros with positive probability (in both prior and posterior). To explain the
gist of the idea, we focus on the equivalent generative process using soft-
thresholding transform:

$\beta\sim\pi^{\beta}_{0},\qquad\theta=\text{sign}(\beta)\circ(|\beta|-\kappa)_{+},$

where the operation is conducted element-wise in the second term, and
$(\cdot)_{+}$ is the truncation function to $[0,\infty)$. In the projection to
L1-ball with L1-norm $\theta\in\mathbb{R}^{p}:\|\theta\|_{1}\leq r$ for some
$r>0$, $\kappa>0$ is a scalar output from a deterministic transformation of
$(\beta,r)$, as shown in xu2020bayesian. On the other hand, one could
generalize $\kappa$ to be another free non-negative parameter or a vector of
non-negative values. Therefore, we refer to the generalized class of
l1ball_type_prior as the L1-ball-type priors or models. As one motivating
example (that we will revisit later), the soft-thresholded Gaussian process
prior (kang2018scalar) is built with $\beta$ from a Gaussian process and
$\kappa$ as a scalar. This prior enjoys a nice property that $\theta$ is
continuous over the input domain with high sparsity at the same time, hence is
particularly useful for selecting sub-regions of an image for regression or
smoothing.

Focusing on the computational aspect, the soft-thresholding transform is
differentiable almost everywhere with respect to $\pi_{0}^{\beta}$. This means
we can use off-the-shelf gradient-based MCMC algorithms (duane1987hybrid;
girolami2011riemann; hoffman2014no; livingstone2022barker) for its posterior
estimation. The strength of these algorithms, besides low implementation cost
due to the good accessibility of software (carpenter2017stan;
bingham2019pyro), is in the rapid convergence to the region near the posterior
mode and a high acceptance rate for changing the zero/non-zero status of
multiple elements at the same time. On the other hand, notice that if the
likelihood is parameterized via $\theta$ and $\kappa$ but not dependent on
$\beta$ a priori, then at the state with some $\theta_{j}=0$, the partial
derivative of the log-posterior density with respect to $\beta_{j}$ is zero.
As a consequence, the algorithm relying on one-step diffusion [such as
Metropolis-adjusted Langevin algorithm (MALA) (rossky1978brownian)] would be
not efficient to explore changing $\theta_{j}$ to a non-zero state. That is
why multiple-step diffusion algorithms such as Hamiltonian Monte Carlo
(neal2011mcmc) and No-U-Turn sampler (hoffman2014no) are preferable.

Since multiple evaluations of the gradient can be expensive, there is room for
improvement that can be made in terms of gaining computational efficiency.
This motivates us to do exploration in the class of Gibbs samplers, which
tends to have a very low computing cost per iteration and requires almost no
tuning. To avoid one-element-at-a-time Gibbs sampling, we introduce an
augmented latent variable that allows us to cancel out the interaction term
among the parameter elements. This enables us to update the elements in a
block in each iteration. We show that this new algorithm achieves much higher
efficiency in terms of effective sample size per time unit. Further, we
establish the geometric ergodicity guarantee in linear regression.
Interestingly, our proposed latent variable is broadly applicable beyond the
L1-ball-type model setting and can be applicable in the challenging cases of
estimating latent Gaussian models, such as those involving multivariate
truncated Gaussian or latent Gaussian process.

##  2 Blocked Gibbs Sampling with Anti-correlation Latent Gaussian

###  2.1 Motivating Problems

We now introduce the motivating sampling problem and provide the necessary
notations. Let $\theta\in\mathbb{R}^{p}$ be the parameter of interest,
associated with precursor random variable $\beta\in\mathbb{R}^{p}$ via
l1ball_type_prior. We use $\mathcal{Y}$ to denote the data, where
$\mathcal{Y}_{i}=(x_{i},y_{i})$ with $y_{i}$ as outcome and $x_{i}$ as
predictor. With some loss of generality, we focus on the following form on the
conditional posterior:

$&\Pi(\theta,\beta\mid
M,\phi,H,\psi,\kappa,\mathcal{Y})\propto\exp\bigg{[}-\frac{1}{2}(\theta^{\prime}M\theta-2\phi^{\prime}\theta)\bigg{]}\exp\bigg{[}-\frac{1}{2}(\beta^{\prime}H\beta-2\psi^{\prime}\beta)\bigg{]},\\\
\theta=\text{sign}(\beta)\circ(|\beta|-\kappa)_{+},$

where $M$ and $H$ are both positive semi-definite matrix of size $p\times p$,
$\phi\in\mathbb{R}^{p},\psi\in\mathbb{R}^{p}$. And we focus on vector
$\kappa\in\mathbb{R}^{p}$ from now on. We treat $\kappa$ as given for now, but
will discuss the case when $\kappa$ is assigned with another prior.

Although the first line in eq:posterior_form resembles a product of two latent
Gaussian densities, the equality constraint in the second line makes it a
degenerate density. We now provide three concrete examples to illustrate the
generality of the above type of problems.

In a linear regression model for $y=X\theta+\epsilon$, with a predictor matrix
$X\in\mathbb{R}^{n\times p}$, $\epsilon\sim\text{N}(0,\Omega^{-1})$ for some
positive definite $\Omega$, and independent Gaussian precursor
$\beta_{j}\stackrel{{\scriptstyle indep}}{{\sim}}\text{N}(0,\tau_{j})$ for
some given hyper-parameter $\tau_{j}>0$. We have eq:posterior_form with
$M=X^{\prime}\Omega X$, $\phi=X^{\prime}\Omega y$, $H=\text{diag}(1/\tau_{j})$
and $\psi=0$.

In a logistic regression with binary
$y_{i}\sim\text{Bernoulli}[1/(1+\exp(-x_{i}^{\prime}\theta)]$, we can use the
Polya-Gamma data augmentation to create a Gaussian augmented likelihood (see
polson2013bayesian for detail). Under independent Gaussian precursor
$\beta_{j}\stackrel{{\scriptstyle indep}}{{\sim}}\text{N}(0,\tau_{j})$, we
have eq:posterior_form with $M=X^{\prime}\Omega X$,
$\phi=X^{\prime}\Omega(y-1/2)$, $H=\text{diag}(1/\tau_{j})$, $\psi=0$, and
$\Omega=\text{diag}(\omega_{j})$ with $\omega_{j}$ marginally from the Polya-
Gamma distribution $\text{PG}(1,0).$

In a sparse smoothing model $y=\theta+\epsilon$, with
$\epsilon\sim\text{N}(0,I{\sigma^{2}})$, $y_{i}$ associated with some spatial
coordinate $s_{i}$, for which the goal is to obtain a spatially correlated
$\theta_{i}$, while making some $\theta_{i}$’s zero in some contiguous regions
of $s_{i}$’s. We can use a spatially correlated Gaussian precursor
$\beta\sim\text{N}[0,K(s,s)]$ with $K$ some covariance kernel such as the
squared exponential kernel. We have eq:posterior_form with $M=I\sigma^{-2}$,
$\phi=y\sigma^{-2}$, $H=[K(s,s)]^{-1}$ and $\psi=0$.

For any non-diagonal $M$ or $H$, the quadratic term $\theta^{\prime}M\theta$
or $\beta^{\prime}H\beta$ in the exponent of eq:posterior_form makes it
difficult to explore a large change in the parameter. This motivates us to use
some latent variables to cancel out those terms.

###  2.2 Anti-correlation Gaussian

If $M$ is non-diagonal and we want to cancel out $\theta^{\prime}M\theta$,
consider a Gaussian latent variable $r\in\mathbb{R}^{p}$:

$&(r\mid\theta,M)\sim\No\big{[}(dI-M)\theta,\;\;(dI-M)\big{]},$

where $d>0$ is a chosen constant to make $dI-M$ positive definite. Similarly,
if $H$ is non-diagonal, we use
$(t\mid\beta,H)\sim\No[(eI-H)\beta,\;\;(eI-H)]$, with $e>0$ making $eI-H$
positive definite. In practice, we find using smaller $d$ and $e$ lead to
better mixing performance, hence we recommend choosing them to be as small as
possible.

We refer to eq:anti-quad-r as an “anti-correlation Gaussian”, because it
cancels out the quadratic correlation terms. For brevity, from now on, we
focus on the case when we use both $r$ and $t$, and we can see that the
conditional posterior of $(\beta_{j},\theta_{j})$ becomes independent over
$j$:

$&\Pi(\beta,\theta\mid
M,\phi,H,r,t,\kappa,\mathcal{Y})\propto\prod_{j=1}^{p}\exp\bigg{\\{}-\frac{1}{2}[d\theta^{2}_{j}-2(\phi_{j}+r_{j})\theta_{j}+e\beta^{2}_{j}-2(\psi_{j}+t_{j})\beta_{j}]\bigg{\\}},\\\
\theta_{j}=\text{sign}(\beta_{j})(|\beta_{j}|-\kappa_{j})_{+}.$

The conditional independence allows us to draw $\beta_{j}$’s in a block. It is
not hard to see that each $\beta_{j}$ follows a three-component mixture,
corresponding to the cases when $\theta_{j}>0$, $\theta_{j}=0$ or
$\theta_{j}<0$, for which we use a discrete variable $b_{j}$ taking value
$1,0$ or $-1$, respectively. We have $b_{j}$ following:

$&\Pi(b_{j}=0\mid\cdot)=c_{j}^{-1}m^{-1}_{0}\exp\bigg{[}\frac{1}{2}\frac{(\psi_{j}+t_{j})^{2}}{e}\bigg{]},\\\
\Pi(b_{j}=1\mid\cdot)=c_{j}^{-1}m^{-1}_{1}\exp\bigg{[}\frac{1}{2}\frac{(\phi_{j}+r_{j}+\psi_{j}+t_{j}+d\kappa_{j})^{2}}{d+e}-\frac{1}{2}d\kappa^{2}_{j}-(\phi_{j}+r_{j})\kappa_{j}\bigg{]}\\\
\Pi(b_{j}=-1\mid\cdot)=c_{j}^{-1}m^{-1}_{-1}\exp\bigg{[}\frac{1}{2}\frac{(\phi_{j}+r_{j}+\psi_{j}+t_{j}-d\kappa_{j})^{2}}{d+e}-\frac{1}{2}d\kappa^{2}_{j}+(\phi_{j}+r_{j})\kappa_{j}\bigg{]},$

where we use $\Pi(b\mid\cdot)$ as a shorthand notation for the full
conditional posterior of $b$ from now on. In the above, $c_{j}^{-1}$ is a
normalizing constant to make the three probabilities add to one; $m_{0}$,
$m_{1}$ and $m_{-1}$ are the constants in the following truncated Gaussian
densities eq:theta_b that multiplied to the exponential terms:

$&(\beta_{j}\mid
b_{j}=0,\cdot)\sim\No_{(-\kappa_{j},\kappa_{j})}(\frac{\psi_{j}+t_{j}}{e},\frac{1}{e}),\\\
(\beta_{j}\mid
b_{j}=1,\cdot)\sim\No_{(\kappa_{j},\infty)}(\frac{\phi_{j}+r_{j}+\psi_{j}+t_{j}+d\kappa_{j}}{d+e},\frac{1}{d+e}),\\\
(\beta_{j}\mid
b_{j}=-1,\cdot)\sim\No_{(-\infty,-\kappa_{j})}(\frac{\phi_{j}+r_{j}+\psi_{j}+t_{j}-d\kappa_{j}}{d+e},\frac{1}{d+e}),$

where the subscript $(\tilde{\alpha}_{1},\tilde{\alpha}_{2})$ denotes the
support of the truncated Gaussian. Therefore, we can draw each $b_{j}$ first,
and then $\beta_{j}$ from the truncated Gaussian. Note that $\theta$ is
completely determined once $\beta$ and $\kappa$ are given. As a result,
besides the other parameters, we have a two-block update scheme based on
sampling $(r,t\mid\beta,\cdot)$, and sampling $(\beta\mid r,t,\cdot)$.

######  Remark 2.1.

In the above, we present the simple case when each $\kappa_{j}$ is given. How
to estimate $\kappa_{j}$ in MCMC depends on the structure and prior one
imposes on the vector $\kappa$. In the examples of this article, we use a
common $\kappa_{j}=\kappa_{0}$ for all $j$, and prior $\pi_{0}(\kappa_{0})$.
We use slice sampling step to update $\kappa_{0}$, with
$\Pi(\kappa_{0}\mid\cdot)\propto\pi_{0}(\kappa_{0})\Pi(\theta,\beta\mid
M,\phi,H,\psi,1\kappa_{0},\mathcal{Y})$ where the second term from
eq:posterior.

###  2.3 Efficient Sampling of the Anti-Correlation Gaussian in Regression

We now focus on making the sampling of the anti-correlation Gaussian more
efficient for large $p$. Motivated by regression settings and with some loss
of generality, we focus on the scenario where $M$ is updated as a three-matrix
product $M=X^{\prime}\Omega X$, where $X$ is fixed but $\Omega$ can change in
each iteration. For simplicity, we focus on sampling of
$r\sim\No[(dI-M)\theta,\;(dI-M)]$, and the method is trivially extensible to
the sampling of $t$, if needed.

Canonically, sampling this multivariate Gaussian involves decomposing the
covariance matrix $(dI-X^{\prime}\Omega X)$, which is an $O(p^{3})$ operation.
In the case of simple regression with homoscedasticity $\Omega=I\sigma^{-2}$
(such as linear regression with homoscedasticity, or the latent Gaussian in
probit regression with $\sigma^{2}=1$ (albert1993bayesian)), we could use
$d=(\lambda_{p}(X^{\prime}X)+\tilde{\epsilon})\sigma^{-2}$ with a small
$\tilde{\epsilon}>0$ and $\lambda_{p}(\cdot)$ denoting the spectral norm, then
we only need to decompose
$(d\sigma^{2}I-X^{\prime}X)=\mathcal{L}\mathcal{L}^{\prime}$ for one time
before the start of the Markov chain. In each Markov chain iteration, we can
sample $r$ efficiently via $r=\sigma^{-1}\mathcal{L}\gamma_{1}+(dI-M)\theta$
with $\gamma_{1}\sim\text{N}(0,I_{p})$.

On the other hand, when $\Omega$ is more complicated than $I\sigma^{-2}$, and
is updated in every iteration (such as in logistic regression where $\Omega$
is updated from a Polya-Gamma distribution), then we would need to compute the
decomposition each time. Similar issues have been noted by
bhattacharya2016fast and nishimura2022prior during the sampling of a Gaussian
with covariance $(X^{\prime}\Omega X+D)^{-1}$ with $D$ positive definite, the
first group of authors exploits the Woodbury matrix identity and involves an
inversion of an $n\times n$ matrix, and the second group use a preconditioned
conjugate gradient method to iteratively solve a linear system. Here our
problem is slightly different, in the form of sampling of a Gaussian with
covariance $(dI-X^{\prime}\Omega X)$ (or equivalently, a transform
$r=(dI-X^{\prime}\Omega X)\gamma_{2}$, with
$\gamma_{2}\sim\text{N}[\theta,(dI-X^{\prime}\Omega X)^{-1}]$). The
subtraction in the covariance/precision matrix creates a new challenge, for
which we develop a new non-iterative algorithm.

For ease of presentation, we now focus on the $p\geq n$ case, and one can
modify it to accommodate the $p\leq n$ case easily. Before the start of the
Markov chain, we pre-compute the singular value decomposition (SVD) of
$X=U_{X}\Lambda_{X}V^{\prime}_{X}$, with $U_{X}$ an $n\times n$ matrix,
$V_{X}$ an $p\times n$ matrix, $\Lambda_{X}$ an $n\times n$ diagonal matrix
with diagonal entries $\Lambda_{X,(i,i)}\geq 0$ for $i=1,\ldots,n$. We denote
$\Lambda^{*}_{X}:=\Lambda_{X,(1,1)}$ as the largest singular value, and
$b_{\Omega}:=[\lambda_{n}(\Omega)]^{-1}$ with $\lambda_{n}(\Omega)$ the
largest eigenvalue of $\Omega$. We assume
$d>\lambda_{n}(\Omega)(\Lambda^{*}_{X})^{2}$, which is sufficient to ensure
positive definiteness of $dI-X^{\prime}\Omega X$. Further, we obtain a matrix
$V_{X}^{\dagger}\in\mathbb{R}^{p\times(p-n)}$ such that
$[V_{X}\;V^{\dagger}_{X}]$ forms an orthonormal matrix, using the full SVD of
$X$. With the above, we use the following sampling scheme.

  1. 1.

Sample $\gamma_{1}\sim\No(0,dI_{n})$, $\gamma_{2}\sim\No(0,dI_{p-n})$,
$\quad\gamma_{3}\sim\No[\Lambda_{X}\gamma_{1}/d,\;b_{\Omega}I_{n}-(\Lambda_{X})^{2}/d]$;

  2. 2.

Sample $\eta\sim\No(0,\Omega^{-1}-b_{\Omega}I_{n})$;

  3. 3.

Set
$r=V_{X}\gamma_{1}+V^{\dagger}_{X}\gamma_{2}-X^{\prime}\Omega(U_{X}\gamma_{3}+\eta)+(dI-X^{\prime}\Omega
X)\theta$.

In the above, the Gaussian could be degenerate with zero variance for some
elements.

######  Theorem 2.2.

The above algorithm produces a random sample from the anti-correlation
Gaussian, $r\sim\No\big{[}(dI-X^{\prime}\Omega
X)\theta,\;\;(dI-X^{\prime}\Omega X)\big{]}.$

######  Proof 2.3.

It is not hard to see that $\bmatrix\gamma_{1}\\\ \gamma_{2}\\\
\gamma_{3}\sim\text{N}(\bmatrix 0\\\ 0\\\ 0,\bmatrix dI_{n}&O\Lambda_{X}\\\
OdI_{p-n}O\\\ \Lambda_{X}Ob_{\Omega}I_{n}),\;\text{and}\;\bmatrix
V_{X}\gamma_{1}+V^{\dagger}_{X}\gamma_{2}\\\
U_{X}\gamma_{3}+\eta\sim\text{N}(\bmatrix 0\\\ 0,\bmatrix dI_{p}&X^{\prime}\\\
X\Omega^{-1}),$ where $O$ in the above denotes a conformable matrix filled
with zeros. The first covariance is positive definite because
$db_{\Omega}>\Lambda^{2}_{X,i,i}$ for any $i=1,\ldots,n$. Then we see that
$\text{VAR}[V_{X}\gamma_{1}+V^{\dagger}_{X}\gamma_{2}-X^{\prime}\Omega(U_{X}\gamma_{3}+\eta)]=dI+X^{\prime}\Omega\Omega^{-1}\Omega
X-2X^{\prime}\Omega X=dI-X^{\prime}\Omega X$. Adding $(dI-X^{\prime}\Omega
X)\theta$ yields the result.

We now discuss the complexity of the above sampling algorithm. The sampling of
$\gamma_{1}$, $\gamma_{2}$, and $\gamma_{3}$ can be carried out very
efficiently due to the diagonal covariance. Often $\Omega$ is a diagonal
matrix, hence sampling of $\eta$ can be quite efficient as well. For the other
cases with non-diagonal $\Omega$ at a large $n$, one often adopts a structured
(such as close-to-low-rank) $\Omega$ or $\Omega^{-1}$ under which $\eta$ can
be sampled efficiently (ghosh2009default; chandra2021bayesian).

######  Remark 2.4.

With the above sampler for $r$ and diagonal $\Omega$, no matrix decomposition
or inversion is needed during the Markov chain Monte Carlo. With vectorized
computation on matrix-vector products, each iteration of posterior sampling
has only a linear cost in $\max(n,p)$.

##  3 Geometric Ergodicity Guarantee

Now we establish the convergence guarantee of our proposed Gibbs sampling
algorithm for linear models. For a clear presentation, we focus on the theory
under models with fixed $(M,\phi,H,\psi,\kappa)$, which is associated with the
two-block updating scheme as presented above; and one could potentially use
our result and further extend to more complicated models.

Our target distribution is associated with the joint posterior density for
$(\theta,\beta,r,t)$: $\Pi(r,t\mid M,\theta,H,\beta)\Pi(\theta,\beta\mid
M,\phi,H,\psi,\kappa,\mathcal{Y}),$ and we denote the associated measure by
$\mu_{\theta,\beta,r,t}(\cdot)$. In MCMC, we start with initial values for
those random variables drawn from a certain initial distribution, and we
denote them by $(\theta^{0},\beta^{0},r^{0},t^{0})$, then we update their
values iterative via the following Markov transition kernel:
$&\mathcal{K}(\theta^{m+1},\beta^{m+1},r^{m+1},t^{m+1}\mid r^{m},t^{m})\\\
=\Pi_{\mathcal{K}}(r^{m+1}\mid\theta^{m+1})\Pi_{\mathcal{K}}(t^{m+1}\mid\beta^{m+1})\Pi_{\mathcal{K}}(\beta^{m+1},\theta^{m+1}\mid
r^{m},t^{m})$ for $m=0,1,2,\ldots$, where the first two densities are based on
eq:anti-quad-r and that of $t$, and the last one is based on eq:b and
eq:theta_b, with all the fixed parameters omitted. We denote the associated
measure for $(r^{m},t^{m},\theta^{m},\beta^{m})$ by
$P^{m}[(r^{0},t^{0},\theta^{0},\beta^{0}),\cdot]$.

In order to show that
$P_{\theta,\beta,r,t}^{m}[(r^{0},t^{0},\theta^{0},\beta^{0}),\cdot]$ converges
to $\mu_{\theta,\beta,r,t}(\cdot)$, we can simply show that
$P^{m}_{r,t}[(r^{0},t^{0}),\cdot]$ converges to $\mu_{r,t}(\cdot)$, with
$\mu_{r,t}(\cdot)$ the marginal posterior measure integrated over
$(\theta,\beta)$, and $P^{m}_{r,t}[(r^{0},t^{0},\theta^{0},\beta^{0}),\cdot]$
the marginal measure of $(r^{m},t^{m})$ after $m$ MCMC iterations after
integrating out $(\theta^{m},\beta^{m})$. Note that the information of
$(\theta^{0},\beta^{0})$ has never been used, as long as we take the first
step in drawing from $\Pi_{\mathcal{K}}(\beta^{1},\theta^{1}\mid
r^{0},t^{0})$. Therefore, we will use $P^{m}_{r,t}[(r^{0},t^{0}),\cdot]$ as a
simplified notation. A similar technique was used by hobert2012tvgs for a two-
part update scheme. With those ingredients in place, we are ready to state the
geometric ergodicity results.

######  Theorem 3.1.

There exists a real-valued function $C_{1}(r^{0},t^{0})$ and $0<\gamma<1$ such
that for all $(r,t)$,

$\|P^{m}_{r,t}[(r^{0},t^{0}),\cdot]-\mu_{r,t}(\cdot)\|_{\text{TV}}\leq
C_{1}(r^{0},t^{0})\gamma^{m}.$

where $\|\cdot\|_{TV}$ denotes the total variation norm. We defer the proof to
Appendix A. Further, we have the following result.

######  Theorem 3.2.

The Markov chain generated by
$\mathcal{K}(\theta^{m+1},\beta^{m+1},r^{m+1},t^{m+1}\mid r^{m},t^{m})$ is
geometrically ergodic. Specifically, there exists a real-valued function
$C_{2}(\theta^{0},\beta^{0},r^{0},t^{0})$ and $0<\gamma<1$ such that for all
$(\theta^{0},\beta^{0},r^{0},t^{0})$,

$\|P^{m}_{(\theta,\beta,r,t)}[(\theta^{0},\beta^{0},r^{0},t^{0}),\cdot]-\mu_{(\theta,\beta,r,t)}(\cdot)\|_{\text{TV}}\leq
C_{2}(\theta^{0},\beta^{0},r^{0},t^{0})\gamma^{m},$

where $\gamma$ is the same rate as the one in eq:marginal_ge.

The proof is a direct extension of Theorem 3.1, since
$\Pi_{\mathcal{K}}(\theta,\beta\mid r,t)$ in the Markov transition kernel
after $m$ iterations and $\Pi(\theta,\beta\mid
r,t,M,\phi,H,\psi,\kappa,\mathcal{Y})$ in the posterior distribution are the
same thing. Therefore, the convergence rate remains the same. See
DiaconisKK2008GibbsSE; juns1994covariance; robert1995convergencecontrol;
rosenthal2001deinitializing for detail.

##  4 Extensions

The data augmentation using anti-correlation Gaussian can be employed beyond
the scope of L1-ball-type models. We now discuss two useful classes of
extensions.

The first extension is on the sampling of latent Gaussian models
(robert_monte_2004; fahrmeir_multivariate_2010; gelman_bayesian_2014), whole
likelihood generally takes the form:

$L(y;\theta)\propto\int|{M}({\theta})|^{1/2}\exp\left[-\frac{1}{2}z^{T}M({\theta})z\right]\prod_{i=1}^{n}g\left(y_{i}\mid
z_{i},\theta\right)dz,$

with $g$ some density/mass function for the data $y$, and $M(\theta)$ is a
positive definite matrix that depends on $\theta$ as the parameter of
interest. This framework covers a wide range of applications such as spatial,
regression or dynamic modeling. Although it seems intuitive to simply sample
$z$ as a latent variable during MCMC, the presence of non-Gaussian $g$ often
disrupts posterior conjugacy, creating difficulty in efficiently updating $z$
especially when $z$ is in high dimension. This is one major motivation for
posterior approximation methods such as integrated nested Laplace
approximation (INLA) (rue2009approximate) as alternative to MCMC.

Using the anti-correlation Gaussian $(r\mid z,\theta)\sim
N\\{[dI-M(\theta)]z,dI-M(\theta)\\}$, we can obtain conditionally independent
$z_{i}$ that can be updated efficiently in a block.

$\Pi(z\mid r,\theta,y)\propto\prod_{i=1}^{n}g\left(y_{i}\mid
z_{i},\theta\right)\exp(-\frac{1}{2}dz_{i}^{2}+z_{i}r_{i})$

Due to the algorithmic similarity to the sampling of the soft-thresholded
Gaussian process model to be presented in the numerical experiments, we omit
further detail on the sampling for latent Gaussian models.

The second extension is on the sampling of truncated multivariate Gaussian:

$\Pi(\theta\mid\mu,\Sigma,R)\propto\exp\left[-\frac{1}{2}(\theta-\mu)^{\rm
T}\Sigma^{-1}(\theta-\mu)\right]1(\theta\in R)$

where $\mu\in\mathbb{R}^{p}$, $\Sigma$ is positive definite, and $R$ some
constrained set of dimension $p$. Despite being commonality in statistical
application, the truncated multivariate Gaussian can pose a challenge for
sampling, especially in high dimension (Wilhelm2010tmvtnormAP;
chow2014preconditioned), due to having correlation and truncation at the same
time. Using the anti-correlation Gaussian $(r\mid\theta,\mu,\Sigma)\sim
N[(dI-\Sigma^{-1})(\theta-\mu),dI-\Sigma^{-1}]$, we have

$\Pi(\theta\mid\mu,r)\propto 1(\theta\in
R)\prod_{j=1}^{p}\exp[-\frac{1}{2}d(\theta_{j}-\mu_{j})^{2}+(\theta_{j}-\mu_{j})r_{j}].$

When the constraints in $R$ are separable over each sub-dimension, then
$\theta_{j}$ is conditionally independent over $j$. For example, with simple
box constraint $R=\cap_{j=1}^{p}\\{x:l_{j}<x_{j}<r_{j}\\}$, then each
$\theta_{j}$ follows a univariate truncated Gaussian and $\theta$ can be
sampled in a block.

##  5 Numerical Examples

We now use several numerical examples to show the computational efficiency of
the anti-correlation blocked Gibbs sampler.

###  5.1 Linear Regression

We first consider the setting where $M$ is non-diagonal in eq:posterior_form.
We use linear regression with $y_{i}\sim\text{N}(x_{i}^{\rm
T}\theta,\sigma^{2}),$ $x_{i}\in\mathbb{R}^{p}$, $i=1,\ldots,n$, $x_{i}$’s
simulated from a multivariate Gaussian with mean zero, and autoregressive
correlation $0.5^{|j-j^{\prime}|}$ between $x_{i,j}$ and $x_{i,j^{\prime}}$.
For the ground truth, we simulate the first $10$ elements of $\theta_{j}$ from
$\No(5,0.5^{2})$ , and set the rest $\theta_{j^{\prime}}$ as zero, and
$\sigma^{2}=10^{-4}$. For the prior specification, we use
$\beta_{j}\stackrel{{\scriptstyle indep}}{{\sim}}\text{N}(0,\tau_{j})$, with
$\tau_{j}\sim\text{Inverse-Gamma}(a_{j},b_{j})$, and
$\kappa_{0}\sim\text{Exp}(\lambda)$; for the noise variance, we use
$\sigma^{2}\sim\text{Inverse-Gamma}(a_{\sigma},b_{\sigma})$. In these
experiments, we use $p=500$, $n=400$, $a_{j}=5$, $b_{j}=1$, $a_{\sigma}=1$,
$b_{\sigma}=1$ and $\lambda=0.5$. In this article, we use rate
parameterization for all inverse-gamma and exponential distributions.

For benchmarking, we compare our algorithm against three other samplers: (i)
the No-U-Turn HMC sampler (hoffman2014no); (ii) the component-wise slice
sampler (neal2003slice) that updates one element at a time; (iii) the
empirical principal component-based slice sampler (thompson2011slice) with
randomized direction. For (iii), we use an empirical estimate of the
covariance between all parameters (based on running 100 iterations of the
component-wise slice sampler during an adaptation period, and the adaptation
time cost is excluded from the calculation of computing efficiency), take the
first $10$ principal components, and run the slice sampler along each
direction under slight perturbation (adding Gaussian noise $\text{N}(0,1/4p)$
to each element of direction, then rescaling the perturbed direction vector to
have unit 2-norm). For (ii) and (iii), we use the stepping out and shrinkage
procedures (neal2003slice) to draw an appropriate subset of slices.

We repeat the experiments 10 times, each under a different random seed number.
For each algorithm, we run MCMC for 10,000 iterations, with the first 2,000 as
warm-ups. We record the total wall clock time (in seconds), and calculate the
effective sample size. Since the popular NUTS sampler, Stan
(carpenter2017stan) is coded in C++, while all the other three algorithms are
done in R, to have a fair comparison on the algorithmic efficiency, we
estimate the running time of NUTS under a pure R implementation from the
‘hmclearn’ package (thomas2021learning). Specifically, we first obtain the
total number of leap-frog steps from the diagonistics of Stan program, then
run the same number of leap-frog steps in hmclearn and record time.

Algorithm | Anti-corr Gaussian | NUTS | Comp-wise slice | EPC slice
---|---|---|---|---
Time spent on 1000 iters | 7.52 | 6659.63 | 938.29 | 20.56
Table 1:  Running time for four algorithms. The time unit is in second based
on pure R implementation for each algorithm.

[t].24 [width=]Images/linear_regression/tp_theta1_AntiCorr.png Iteration

Figure 1: Trace plot of $\theta_{1}$ for Anti-corr Gaussian.

[t].24 [width=]Images/linear_regression/tp_theta1_NUTS.png Iteration

Figure 2: Trace plot of $\theta_{1}$ for NUTS.

[t].24 [width=]Images/linear_regression/tp_theta1_EPCSlice.png Iteration

Figure 3: Trace plot of $\theta_{1}$ for EPC slice.

[t].24 [width=]Images/linear_regression/tp_theta1_CompWiseSlice.png Iteration

Figure 4: Trace plot of $\theta_{1}$ for Comp-wise slice.

[t].49 [width=]Images/linear_regression/new_ESS_s_first10.png Algorithm

ESS/s

Figure 5: Average ESS/s for the first 10 entries of $\theta$.

[t].49 [width=keepaspectratio]Images/linear_regression/new_ESS_s_rest.png
Algorithm

ESS/s

Figure 6: Average ESS/s for the last 490 entries of $\theta$. Figure 7:
Traceplots and effective sample size per computing time (ESS/s) for the four
algorithms: anti-correlation Gaussian (Anti-corr Gaussian), No-U-Turn Sampler
(NUTS), empirical principal component-based slice sampler (EPC slice), and
component-wise slice sampler (Comp-wise slice).

In the result, all algorithms have the Markov chain sample concentrated near
the ground truth. As shown in Table 1, there is a clear difference in the
running time. The anti-correlation Gaussian is the fastest, followed by EPC
slice, component-wise slice, and NUTS. The NUTS algorithm has the fastest
mixing; however, this is at the cost of expensive computation per iteration.
As shown in Figure 7, the anti-correlation Gaussian algorithm has the highest
effective sample size per computing time.

###  5.2 Application on Soft-Thresholded Gaussian Process

We next consider the setting where $H$ is non-diagonal in eq:posterior_form.
We use soft-thresholded Gaussian process (kang2018scalar) as an example, which
takes the form of:
$\theta_{s}=\text{sign}(\beta_{s})\circ(|\beta_{s}|-\kappa_{0})_{+},\qquad\beta\sim\text{GP}[0,K(\cdot,\cdot)],$
where $s$ is the location in some input location space $\mathbb{S}$,
$\kappa_{0}>0$ is a scalar, and $\text{GP}$ denotes a Gaussian process from
which any finite-dimensional realization follows a multivariate (potentially
degenerate) Gaussian with zero mean, and covariance parameterized by
$\text{Cov}(\beta_{s},\beta_{s^{\prime}})=K(s,s^{\prime})$. With an
appropriate choice of $K$, we can obtain useful properties on a GP
realization, such as continuity or smoothness over $s$. After applying soft-
thresholding, we preserve the continuity over $\mathbb{S}$, and smoothness in
each open set where $\theta_{s}\neq 0$. Therefore, the above is useful for
obtaining a “smooth and sparse” parameterization for $\theta$. In this
article, we consider a simple but useful application in image smoothing, in
the form of:
$y_{s}=\theta_{s}+\epsilon_{s},\qquad\epsilon_{s}\stackrel{{\scriptstyle
iid}}{{\sim}}\text{N}(0,\sigma^{2}),$ where $s$ is the pixel location,
$s=(i,j)$ with $i=1,\ldots,n_{1}$ and $j=1,\ldots,n_{2}$. For the covariance
function, we use
$K(s,s^{\prime})=\tau\exp[-{\|s-s^{\prime}\|^{2}_{2}}/(2\xi^{2})]$.

In application, we use the functional magnetic resonance imaging scan of one
human subject who was performing a motor task. We take a sectional view along
the anterior-posterior axis of the brain, corresponding to $91\times 91$
pixels per frame, and over a time period, we collect $280$ frames. Using
superscript $f$ to index each frame, we assume
$y^{f}=\\{y_{s}^{f}\\}_{\text{all}s}$ to be independent over $f$ with frame-
varying mean $\theta^{f}=\\{\theta_{s}^{f}\\}_{\text{all}s}$, but the other
parameters $\sigma^{2}$, $\tau$, $\xi$ and $\kappa_{0}$ are invariant over
$f=1,\ldots,280$. To finish the prior specification, we use
$\tau\sim\text{Inverse-Gamma}(0.1,0.1)$, $\sigma^{2}\sim\text{Inverse-
Gamma}(0.1,0.1)$ and $\kappa_{0}\sim\text{Exp}(0.5)$. tikhonov2020joint has
previously shown that when the goal of using GP is to obtain smoothing of
data, estimating the bandwidth $\xi$ under a reasonable degree of precision is
empirically adequate. Specifically, they impose a simple discrete prior on
$\xi$ supported on a finite set of values, so that all possible precision
matrix $H$ and its determinant could be pre-computed, hence further
facilitating the computation. We follow their solution here, with
$\xi=0.5\tilde{t}$ and $\tilde{t}$ a discrete uniform prior supported on
$\\{1,\ldots,20\\}$.

When using our proposed Gibbs sampler for posterior estimation, we sample
$\beta$ using the anti-correlation Gaussian data augmentation technique,
sample $\sigma^{2}$ and $\tau$ from their corresponding inverse-gamma full
conditional distributions, and sample each of $\xi$ and $\kappa_{0}$ using the
slice sampling algorithm. We run $10,000$ iterations with the first $2,000$ as
burn-ins, and it takes about 270 minutes on a MacBook Pro with a 10-core CPU.
We report the effective sample size per hour for four randomly selected
locations of $\theta^{100}$ at $f=100$: $315.14$ at $s=3300$, $315.44$ at
$s=3320$, $268.95$ at $s=3720$, and $227.16$ at $s=6700$. To further assess
the mixing, we plot both the traceplots at two locations and an
autocorrelation function (ACF) plot in Figure 11. The traceplots show fast
mixing and the ACF drops quickly.

[t]0.27 [width=]Images/STGP/tp_theta3320.png Iteration

Figure 8: Traceplot of $\theta^{100}_{3320}$.

[t]0.27 [width=]Images/STGP/tp_theta6700.png Iteration

Figure 9: Traceplot of $\theta^{100}_{6700}$.

[t]0.41 [width=]Images/STGP/acf_theta.png Lag

ACF

Figure 10: ACF plot of $\theta^{100}$. Figure 11: Traceplots and
autocorrelation function (ACF) plot. Each box on the ACF plot incorporates the
ACF for all pixels of an image.

In Figure 22, (a)-(c) show the posterior mean estimates for three frames, and
(d)-(f) do the point-wise standard deviations. We can clearly see both
smoothness and sparsity in the mean estimates, as well as how they change over
time. To show the benefits of using a soft-thresholded GP, we also plot the
results of simple thresholding on the raw images from (g)-(j), which lead to
non-smooth estimates and scattered noisy non-zero points on the estimates.

To show how our proposed algorithm improves the computational efficiency of
the soft-thresholded GP method, we also run a No-U-Turn sampler on the same
model using Stan. We try NUTS with max tree depth equal to 7, it takes 68
hours to finish $10,000$ iterations, and the effective sample size per hour
for estimating $\theta^{100}$ is $8.1$ on average, $9.2$ at $s=3300$, $9.3$ at
$s=3320$, $7.4$ at $s=3720$, and $6.7$ at $s=6700$. Therefore, our algorithm
improves the computing speed substantially.

[t].32 [width=]Images/STGP/post_est_frame50.png [t].32
[width=]Images/STGP/post_est_frame100.png [t].32
[width=]Images/STGP/post_est_frame200.png [t].32
[width=]Images/STGP/sd_frame50.png [t].32 [width=]Images/STGP/sd_frame100.png
[t].32 [width=]Images/STGP/sd_frame200.png

Figure 12: Posterior estimate $\hat{\theta}^{50}$. Figure 13: Posterior
estimate $\hat{\theta}^{100}$. Figure 14: Posterior estimate
$\hat{\theta}^{200}$. Figure 15: Point-wise standard deviations for
$\hat{\theta}^{50}$. Figure 16: Point-wise standard deviations for
$\hat{\theta}^{100}$. Figure 17: Point-wise standard deviations for
$\hat{\theta}^{200}$.

[t].24 [width=]Images/STGP/rawimage_0.png [t].24
[width=]Images/STGP/rawimage_1.png [t].24 [width=]Images/STGP/rawimage_3.png
[t].24 [width=]Images/STGP/rawimage_5.png

Figure 18: Raw image $y^{100}$. Figure 19: Raw image $y^{100}$ after the soft-
thresholding transform with threshold=1. Figure 20: Raw image $y^{100}$ after
the soft-thresholding transform with threshold=3. Figure 21: Raw image
$y^{100}$ after the soft-thresholding transform with threshold=5. Figure 22:
Posterior estimates, point-wise standard deviations, and direct soft-
thresholding transform on a raw image. The posterior estimates by STGP are
smooth and sparse with point-wise uncertainty estimated whereas the direct
soft-thresholding transform only leads to noisy estimates.

##  6 Discussion

The latent Gaussian form is very common in Bayesian posterior estimation. In
this article, we show the advantages of using anti-correlation Gaussian to
cancel out the quadratic term in the exponent. We expect our technique to be
broadly applicable beyond our presented scenarios. For example, one could
potentially extend this technique to address sampling problems of Bingham
distribution in Stiefel manifold (hoff2009simulation), the costly evaluation
of determinant in determinantal point process (kulesza2012determinantal), and
so forth.

When it comes to linear regression variable selection using discrete spike-
and-slab prior, there is some close connection of our algorithm to the ODA
algorithm (ghosh2011rao), which uses an augmented design matrix to make the
complete gram matrix diagonal. Since one could reparameterize a discrete
spike-and-slab prior as a special case of the L1-ball prior by setting
$\kappa$ according to a quantile of $\pi_{0}(\beta)$, we expect our geometric
ergodicity result could be extended to analyzing the ODA algorithm. In
comparison, our algorithm has the advantages of not having to create an
augmented design matrix and easy application in non-linear models. We expect
some nice properties that are further exploited in the ODA algorithm, such as
the collapsed sampling step based on marginalizing out $\theta_{j}$ given if
$b_{j}=1$ or $0$, could be obtained in some specific form of latent Gaussian
model after we augment the anti-correlation Gaussian.

##  Appendix A Proof of Geometric Ergodicity in Linear Models

new_ergodicity.tex

##  Appendix B Additional Results

###  B.1 Alternative Algorithm for Sampling Anti-correlation Gaussian

We now describe another efficient algorithm to sample from the anti-
correlation Gaussian when $M=X^{\prime}\Omega X$. Since $dI-X^{\prime}\Omega
X$ is positive definite, we have the following converging series:
$(I-d^{-1}X^{\prime}\Omega X)^{-1}=\sum_{k=0}^{\infty}(d^{-1}X^{\prime}\Omega
X)^{k}.$

When truncating the right-hand side at $K$, we have:

$\sum_{k=0}^{K}(d^{-1}X^{\prime}\Omega X)^{k}=[I-(d^{-1}X^{\prime}\Omega
X)^{K+1}](I-d^{-1}X^{\prime}\Omega X)^{-1},$

where the relative error term $\|(d^{-1}X^{\prime}\Omega
X)^{K+1}\|_{2}\leq\rho^{K}$, with $\rho=\lambda_{p}(X^{\prime}\Omega X)/d$.
Therefore, with a selected $D$ and a numeric tolerance $\epsilon$, we can
control $\rho$ and find out a suitable truncation
$\hat{K}=\log_{\rho}\epsilon$. In our code, we use $\epsilon=10^{-8}$, and
$\rho=2/3$ leading to $\hat{K}=45$.

To generate a Gaussian with mean $0$ and covariance as in eq:truncated_cov, we
can run the following algorithm, starting from a $p$-element $\gamma=0$.

  1. 1.

Simulate $\eta_{1}\sim\text{N}(0,I_{p})$, $\eta_{2}\sim\text{N}(0,I_{n}),$

  2. 2.

Update
$\gamma\leftarrow\eta_{1}+d^{-1}X^{\prime}\Omega^{1/2}\eta_{2}+(d^{-1}X^{\prime}\Omega
X)\gamma$.

And it is not hard to see that each iteration increases $K$ by $2$, therefore,
we only need to run the above algorithm for $23$ iterations. Afterward, we set
$r=(dI-X^{\prime}\Omega X)d^{-1/2}\gamma+(dI-X^{\prime}\Omega X)\theta,$ which
yields a Gaussian with mean $(dI-X^{\prime}\Omega X)\theta$, variance
numerically very close to $(dI-X^{\prime}\Omega X)$.

###  B.2 Benchmark: Sampling Truncated Multivariate Gaussian using Anti-
correlation Gaussian

A comparison regarding the algorithm efficiency is conducted between the anti-
correlation Gaussian and the ‘rtmvnorm()’ function in the R package ‘tmvtnorm’
(See Wilhelm2010tmvtnormAP for details). The argument “algorithm ” in the
function is set to be ‘gibbsR’ for a fair comparison in R language; otherwise,
it uses Fortran to accelerate the function. We sample from $p$-dimensional
truncated multivariate Gaussian distribution with zero mean, identity
covariance matrix, and truncation area $R=\prod_{i=1}^{p}(-4,-3]$. For such
specifications, rejection sampling breaks down even when $p$ is as small as
10. We conducted $20$ simulations with $p=10$ for each algorithm. For each
simulation, we do 12,000 iterations with 2,000 burn-ins and record the total
clock wall time (in seconds). Figure 26 presents the trace plots of one
randomly selected element of the parameter and compares the effective sample
size per computing time (ESS/s) for anti-correlation Gaussian and rtmvnorm.
Both algorithms enjoy fast mixing; however, anti-correlation Gaussian has an
average ESS/s more than twice the average ESS/s for rtmvnorm across the 20
replications.

[t]0.32 [width=]Images/TruncMVN/tp_2_AntiCorr.png Iteration

Figure 23: Traceplot of the second element of the parameter using Anti-Corr
Gaussian.

[t]0.32 [width=]Images/TruncMVN/tp_2_builtin.png Iteration

Figure 24: Traceplot of the second element of the parameter using rtmvnorm.

[t]0.31 [width=]Images/TruncMVN/Avg_ESS_s.png Algorithm

ESS/s ($\times$1,000)

Figure 25: Average ESS/s. Figure 26: Traceplots and effective sample size per
computing time (ESS/s) for anti-correlation Gaussian Gibbs sampler and
rtmvnorm.