# Bayesian $L_{\frac{1}{2}}$ regression

Xiongwen Ke and Yanan Fan
UNSW Data Science Hub,
School of Mathematics and Statistics,
UNSW Sydney, 2052 Australia  Corresponding author<EMAIL_ADDRESS>
###### Abstract

It is well known that bridge regression enjoys superior theoretical properties
than traditional LASSO. However, the current latent variable representation of
its Bayesian counterpart, based on the exponential power prior, is
computationally expensive in higher dimensions. In this paper, we show that
the exponential power prior has a closed-form scale mixture of normal
decomposition for $\alpha=(\frac{1}{2})^{\gamma},\gamma\in\mathbb{N}^{+}$. We
develop a partially collapsed Gibbs sampling scheme, which outperforms
existing Markov chain Monte Carlo strategies, we also study theoretical
properties under this prior when $p>n$. In addition, we introduce a non-
separable bridge penalty function inspired by the fully Bayesian formulation
and a novel, efficient, coordinate-descent algorithm. We prove the algorithm’s
convergence and show that the local minimizer from our optimization algorithm
has an oracle property. Finally, simulation studies were carried out to
illustrate the performance of the new algorithms.

Keywords: Bridge shrinkage; High dimensional regression; MCMC; Sparse
optimization.

##  1 Introduction

Consider the linear regression problem

$\boldsymbol{Y}=\boldsymbol{X}\boldsymbol{\beta}+\sigma\boldsymbol{\epsilon}$
(1)

where $\boldsymbol{Y}=(y_{1},\ldots,y_{n})$ is an n-dimensional response
vector, assumed to have been centered to 0 to avoid the need for an intercept.
$\boldsymbol{X}$ is an $n\times p$ design matrix consisting of $p$
standardised covariate measurements.
$\boldsymbol{\beta}=(\beta_{1},\ldots,\beta_{p})^{T}$ is a $p\times 1$ vector
of unknown coefficients,
$\boldsymbol{\epsilon}=(\epsilon_{1},\ldots,\epsilon_{n})$ is a vector of
i.i.d noise with mean $0$ and variance 1, and $\sigma>0$ is the error standard
deviation. The bridge estimator (Frank and Friedman,, 1993) is the solution to
the objective function of the form

$\arg\min_{\boldsymbol{\beta}}\frac{1}{2}\|\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}\|^{2}+\lambda\sum_{j=1}^{p}\left|\beta_{j}\right|^{\alpha}$

with $0<\alpha<1,\lambda>0$. When $\alpha=1$, the bridge estimator is the same
as the LASSO (whose Bayesian counterpart is also known as laplace or double
exponential prior) (Park and Casella,, 2008). The bridge estimator produces a
sparser solution than the LASSO penalty ($\alpha=1$), while the case of
$\alpha=0$ is the NP-hard problem of subset selection.

Several papers have studied statistical properties of bridge regression
estimators and their optimization strategies from a frequentist point of view
(Knight et al.,, 2000; Huang et al.,, 2008; Zou and Li,, 2008). From a
Bayesian perspective, the bridge penalty induces the exponential power prior
distribution for $\boldsymbol{\beta}$ of the form

$\pi(\boldsymbol{\beta}|\lambda,\alpha)\propto\prod_{j=1}^{p}\frac{\lambda^{1/\alpha}}{2\Gamma(1+1/\alpha)}\exp\left(-\lambda\left|\beta_{j}\right|^{\alpha}\right),\quad
0<\alpha<1.$

This prior can be expressed as a scale mixture of normals (Polson et al.,,
2014; Armagan,, 2009; West,, 1987)

$\pi\left(\boldsymbol{\beta}|\lambda,\alpha\right)\propto\prod_{j=1}^{p}\int_{0}^{\infty}\mathcal{N}\left(0,\tau_{j}^{2}\lambda^{-2/\alpha}\right)\pi\left(\tau_{j}^{2}\right)d\tau_{j}^{2}$

where $\pi(\tau_{j}^{2})$ is the density for scale parameter $\tau_{j}^{2}$.
However, there is no closed form expression for $\pi(\tau_{j}^{2})$. Polson et
al., (2014) proposed to work with the conditional distribution of
$\boldsymbol{\tau}^{2}|\boldsymbol{\beta}$ in their MCMC strategy:

$\pi(\tau^{2}|\boldsymbol{\beta})\sim\prod_{j=1}^{p}\frac{\exp\left(-\lambda^{2/\alpha}\left|\beta_{j}\right|^{2}\tau_{j}^{2}\right)\pi\left(\tau_{j}^{2}\right)}{E_{\pi\left(\tau_{j}^{2}\right)}\left\\{\exp\left(-\lambda^{2/\alpha}\left|\beta_{j}\right|^{2}\tau_{j}^{2}\right)\right\\}}$

which is an exponentially-tilted stable distribution. They suggested the use
of the double rejection sampling algorithm of Devroye, (2009) to sample this
conditional distribution. However, the approach is very complicated and not
easy to scale to high dimensions.

A second approach proposed by Polson et al., (2014), is based on a scale
mixture of triangular (SMT) representation. More recently, an alternative
uniform-gamma representation was introduced by Mallick and Yi, (2018). When
the design matrix $\boldsymbol{X}$ is orthogonal, both the scale mixture of
triangles (Polson et al.,, 2014) and the scale mixture of uniforms (Mallick
and Yi,, 2018) work well, but both of them suffer from poor mixing when the
design matrix is strongly collinear. This is because both sampling schemes
need to generate from truncated multivariate normal distributions, which are
still difficult to sample from efficiently in higher dimensions.

In this paper, we propose a new data-augmentation strategy, which provides us
with a laplace-gamma mixture representation of the bridge penalty at
$\alpha=\frac{1}{2}$, or $L_{(\frac{1}{2})^{\gamma}}$ prior for $\gamma=1$,
which is further decomposed into a normal-exponential-gamma mixture. This
procedure can be generalised to $\alpha=(\frac{1}{2})^{\gamma}$ for any
$\gamma\in\mathbb{N}^{+}$. This closed form representation introduces extra
latent variables, allowing us to circumvent sampling difficult conditional
distribution $\boldsymbol{\tau}^{2}|\boldsymbol{\beta}$ directly. We further
leverage the conjugacy structure in our model by removing certain conditional
components in the Gibbs sampler without disturbing the stationary
distribution. This is done by using the partially collapsed Gibbs sampling
strategy (Park and Van Dyk,, 2009; Van Dyk and Park,, 2008; Van Dyk and Jiao,,
2015). Thus improving the speed of convergence of the Markov chain Monte Carlo
(MCMC) algorithm.

While continuous shrinkage priors have obvious computational advantages over
spike-and-slab priors (which generally involve a combinatorial search over all
$2^{p}$ possible models), their resulting posterior inference from MCMC does
not automatically provide a sparse solution of the coefficients, requiring
extra and ad hoc steps for variable selection. On the other hand, it is well
known that the penalized likelihood estimators have a Bayesian interpretation
as posterior modes under the corresponding prior. Thus, it is possible to
construct very flexible sparse penalty functions from a Bayesian perspective,
whose solution is the generalized thresholding rule (Ročková and George,,
2016). We consider a full Bayesian formulation of the exponential power prior,
which introduces a non separable bridge (NSB) penalty

${\bf
pen}(\boldsymbol{\beta})=\log\int\prod_{j}\pi(\beta_{j}|\lambda)\pi(\lambda)d\lambda$

with the hyper-parameter $\lambda$ being integrated out with respect to some
suitable choice of hyper-prior. Polson and Scott, (2012); Song and Liang,
(2017) showed that if the global-local shrinkage prior has a heavy and flat
tail, and allocates a sufficiently large probability mass in a very small
neighborhood of zero, then its posterior properties are as good as those of
the spike-and-slab prior. However these properties often lead to unbounded
derivatives for the penalty function at zero. In particular, our non-separable
bridge penalty is a non-convex, non-separable and non-lipschitz function,
posing a very challenging problem for optimization. For the second main
innovation of this paper, we propose a simple CD algorithm which is capable of
producing sparse solutions. We further show that the convergence of the
algorithm to the local minima is guaranteed theoretically.

The paper is structured as follows: In Section 2, we present the normal
exponential gamma mixture representation for $L_{\frac{1}{2}}$ prior for
$\gamma=1$ and its generalization to $\gamma>1$. We also discuss the full
Bayesian approach for hyper-parameter selection. In Section 3, we produce
theoretical conditions for posterior consistency and contraction rates. In
Section 4, we propose a efficient posterior sampling scheme based on the
partially collapsed Gibbs sampler.In Section 5, we develop a coordinate
dsecent optimization algorithm for variable selection and study its oracle
property. In Section 6, we provide some numerical results from simulation
studies. Finally in Section 7, we conclude with some discussions.

##  2 The $L_{\frac{1}{2}}$ prior

In this section, we begin with the decomposition of the exponential power
prior with $L_{(\frac{1}{2})^{\gamma}}$. We first present results for the case
$\gamma=1$, and then generalise the decompositions to $\gamma>1$. Throughout
this work we refer to the $L_{(\frac{1}{2})^{\gamma}}$ prior as the
$L_{\frac{1}{2}}$ prior for simplicity.

###  2.1 The Laplace mixture representation

The scale mixture of representation was first found by West, (1987), who
showed that a function $f(x)$ is completely monotone if and only if it can be
represented as a Laplace transform of some function $g(\cdot)$:

$f(x)=\int_{0}^{\infty}\exp(-sx)g(s)\mathrm{d}(s).$

To represent the exponential power prior with $0<\alpha<1$ as a gaussian
mixture, we let $x=\frac{t^{2}}{2}$,

$\exp\left(-|t|^{\alpha}\right)=\int_{0}^{\infty}\exp\left(-st^{2}/2\right)g(s)\mathrm{d}s$

where $\exp\left(-|t|^{\alpha}\right)$ is the the Laplace transform of $g(s)$
evaluated at $\frac{t^{2}}{2}$ (Polson and Scott,, 2012; Polson et al.,,
2014). Unfortunately, there is no general closed form expression for $g(s)$.
However, for the special case $\alpha=\frac{1}{2^{\gamma}}$ where $\gamma$ is
any positive integer, we can construct a data augmentation scheme to represent
$g(s)$. We begin with special case $\alpha=\frac{1}{2}$ in Lemma 2.1.

######  Lemma 2.1.

The exponential power distribution with $\alpha=\frac{1}{2}$, of the form
$\pi_{\lambda}(\beta)=\frac{\lambda^{2}}{4}\exp(-\lambda|\beta|^{\frac{1}{2}}),$
can be decomposed as $\beta|s\sim\mathrm{DE}(0,s)$,
$s\sim\mathrm{Gamma}(\frac{3}{2},\frac{\lambda^{2}}{4})$ mixture or
equivalently, $\beta|v\sim\mathrm{DE}(0,\frac{v}{\lambda^{2}})$,
$v\sim\mathrm{Gamma}(\frac{3}{2},\frac{1}{4})$ mixture.

Here $\mathrm{DE}(b)$ denotes a laplace (also known as double exponential
distribution with mean 0 and variance $2b^{2}$. Since the Laplace distribution
can be represented by the normal-exponential mixture (Andrews and Mallows,,
1974), we therefore obtain the normal-exponential-gamma mixture representation
for the exponential power prior with $\alpha=\frac{1}{2}$. Lemma 2.2 provides
a recursive relation for the more general $\alpha=\frac{1}{2^{\gamma}}$ case.

######  Lemma 2.2.

The exponential power distribution with $\alpha=(\frac{1}{2})^{\gamma}$,
($\gamma\geq 1$, $\gamma$ is positive integer) of the form
$\pi_{\lambda}(\beta)\propto\frac{\lambda^{2^{\gamma}}}{2(2^{\gamma}!)}\exp\left(-\lambda\left|\beta\right|^{\frac{1}{2^{\gamma}}}\right)$
can be represented as the mixture of the exponential power distribution with
$\alpha=(\frac{1}{2})^{\gamma-1}$ and Gamma distribution. More specifically,
we have

$\pi(\beta|s_{\gamma})\propto\frac{1}{2(2^{\gamma-1}!)s_{\gamma}^{2^{\gamma-1}}}\exp\left(-\frac{|\beta|^{\frac{1}{2^{\gamma-1}}}}{s_{\gamma}}\right),\quad
s_{\gamma}\sim{\mathrm{Gamma}\left(\frac{2^{\gamma}+1}{2},\frac{\lambda^{2}}{4}\right)}$

or equivalently

$\pi(\beta|v_{\gamma},\lambda)\propto\frac{\lambda^{2^{\gamma}}}{2(2^{\gamma-1}!)v_{\gamma}^{2^{\gamma-1}}}\exp\left(-\frac{\lambda^{2}|\beta|^{\frac{1}{2^{\gamma-1}}}}{v_{\gamma}}\right),\quad
v_{\gamma}\sim\mathrm{Gamma}\left(\frac{2^{\gamma}+1}{2},\frac{1}{4}\right).$

Applying Lemmas 2.1 and 2.2, the main Theorem 2.3 provides the analytic
expressions under the scale mixture of normal representation for the
exponential power prior.

######  Theorem 2.3.

The exponential power prior with $\alpha=(\frac{1}{2})^{\gamma}$, ($\gamma\geq
1$ and $\gamma$ is positive integer) can be represented as laplace mixture
with $i=1,\ldots,\gamma-1$,

$\beta|v_{1},\lambda\sim\mathrm{DE}\left(\frac{v_{1}}{\lambda^{2^{\gamma}}}\right),\quad
v_{i}|v_{i+1}\sim\mathrm{Gamma}\left(\frac{2^{i}+1}{2},\frac{1}{4v_{i+1}^{2}}\right),\quad
v_{\gamma}\sim\mathrm{Gamma}\left(\frac{2^{\gamma}+1}{2},\frac{1}{4}\right)\\\
$

which can further be represented as global-local shrinkage prior:

$\displaystyle\beta|\tau^{2},\lambda\sim
N\left(\mathbf{0},\frac{\tau^{2}}{\lambda^{2^{\gamma+1}}}\right),\quad$
$\displaystyle\tau^{2}|v_{1}\sim\mathrm{Exp}\left(\frac{1}{2v_{1}^{2}}\right),$
(2) $\displaystyle
v_{i}|v_{i+1}\sim\mathrm{Gamma}\left(\frac{2^{i}+1}{2},\frac{1}{4v_{i+1}^{2}}\right),\quad$
$\displaystyle
v_{\gamma}\sim\mathrm{Gamma}\left(\frac{2^{\gamma}+1}{2},\frac{1}{4}\right),$

where $i=1,\ldots,\gamma-1$, $\tau$ is the local shrinkage parameter and
$\frac{1}{\lambda^{2^{\gamma}}}$ is the global shrinkage parameter. The proof
can be found in section 1 of supplementary materials.

As we can see, the marginal distribution for the local shrinkage parameter
$\tau$ has no closed form expression, leading to computational difficulties in
existing MCMC schemes (Polson et al.,, 2014; Mallick and Yi,, 2018). However,
we show that for the cases $\alpha=(\frac{1}{2})^{\gamma}$, introducing latent
variables $v_{i}$ leads to a computationally efficient decomposition. The
augmented representation now makes it easier to create efficient computational
algorithms such as MCMC or the EM algorithm for finding posterior modes.

###  2.2 Choice of hyper-prior

The hyper-parameter $\lambda$ controls global shrinkage, and is critical to
the success of the $L_{\frac{1}{2}}$ prior. It is possible to fix $\lambda$ by
selecting a value for it via an empirical Bayes approach, e.g., marginal
maximum likelihood (Casella,, 2001). However, in extremely sparse cases, the
empirical Bayes estimate of the global shrinkage parameter might collapse to 0
( Scott and Berger, (2010); Datta et al., (2013)). Assigning a hyper-prior to
$\lambda$ allows the model to achieve a level of self-adaptivity and can boost
performance(Scott and Berger,, 2010; Ročková and George,, 2018). Here we
propose to assign a half Cauchy prior to $\frac{1}{\sqrt{\lambda}}$. That is

$\frac{1}{\sqrt{\lambda}}\sim\mathrm{Cauchy_{+}}(0,1),$

which has the following scale mixture representation:

$\lambda\sim\mathrm{Gamma}(\frac{1}{2},\frac{1}{b}),\quad
b\sim\mathrm{InvGamma}(\frac{1}{2},1),$ (3)

This allows us to keep the conjugate structure of our Gibbs sampling scheme
later.

The main motivation for using this prior is that the density of the global
shrinkage parameter $t=\frac{1}{\lambda^{2^{\gamma}}}$ approaches infinity at
the origin. This prior yields a proper posterior, allowing for strong
shrinkage, while its thick tail can accommodate a wide range of values. This
can be seen by the simple calculation:

$\pi(t)=\frac{1}{2^{\gamma+1}(1+t^{\frac{1}{2^{\gamma+1}}})t^{1-\frac{1}{2^{\gamma+2}}}};\quad\pi(z)=\frac{1}{2^{\gamma+1}(1+z^{\frac{1}{2^{\gamma+1}}})z^{1-\frac{1}{2^{\gamma+2}}}}$

where $z=\frac{1}{t}=\lambda^{2^{\gamma}}$. Thus
$\lim_{t\rightarrow{0}}\pi(t)=\infty$ and
$\lim_{z\rightarrow{0}}\pi(z)=\infty$.

##  3 Posterior consistency and contraction rate

We focus on the high-dimensional case where the number of predictors $p$ is
much larger than the number of observations $n$ ($p>n$) and most of the
coefficients in the parameter vector
$\boldsymbol{\beta}=(\beta_{1},\ldots,\beta_{p})$ are zero. We consider the
exponential power prior with $\alpha=\frac{1}{2^{\gamma}}$ on
$\boldsymbol{\beta}$ under the following set up,

$\boldsymbol{Y}|\boldsymbol{\beta},\sigma^{2}\sim\mathcal{N}_{n}\left(\boldsymbol{X}\beta,\sigma^{2}\mathbf{I}_{\mathbf{n}}\right),\quad\boldsymbol{\beta}|\lambda\propto\prod_{j=1}^{P}\frac{\lambda^{2^{\gamma}}}{2(2^{\gamma}!)}\exp\left(-\lambda\left|\beta_{j}\right|^{\frac{1}{2^{\gamma}}}\right),\quad\pi\left(\sigma^{2}\right)\propto\sigma^{-2}.$

Mallick and Yi, (2018) showed strong consistency of the posterior under
exponential power prior in the case $p<n$. When $p>n$, the non-invertibility
of the design matrix complicates analysis. In the following, we show that
under the local invertibility assumption of the Gram matrix
$\frac{\boldsymbol{X}^{T}\boldsymbol{X}}{n}$, the $L_{2}$ contraction rates
for the posterior of $\boldsymbol{\beta}$ are nearly optimal and no worse than
the convergence rates achieved by spike-and-slab prior (Castillo et al.,,
2015). Therefore,there is no estimation performance loss due to switching from
spike-and-slab prior to the $L_{\frac{1}{2}}$ prior.

In what follows, we rewrite the dimension $p$ of the model by $p_{n}$ to
indicate that the number of predictors can increase with the sample size $n$,
and similarly we rewrite $\lambda$ by $\lambda_{n}$. We use $\beta_{0}$ and
$\sigma_{0}$ to indicate true regression coefficients and the true standard
deviation. Let $S_{0}$ be the set containing the indices of the true nonzero
coefficients, where $S_{0}\subseteq\left\\{1,..,p_{n}\right\\}$, then
$s_{0}=|S_{0}|$ denote the size of the true model. We use $\Pi_{n}$ to
emphasis that the prior distribution is sample size dependent. We now state
our main theorem on the posterior contraction rates for the exponential power
prior based on the following assumptions:
A1. The dimensional is high $p_{n}\succ n$ and $\log(p_{n})=o(n)$.
A2. The true number of nonzero $\beta_{0}$ satisfies $s_{0}=o(n/\log p_{n})$.
A3. All the covariates are uniformly bounded. In other words, there exists a
constant $k>0$ such that
$\lambda_{\max}(\boldsymbol{X}^{T}\boldsymbol{X})\leq{kn^{\alpha}}$ for some
$\alpha\in{[1,+\infty)}$.
A4. There exist constants $v_{2}>v_{1}>0$ and an integer $\widetilde{p}$
satisfying $s_{0}=o(\widetilde{p})$ and $\widetilde{p}=o(s_{0}\log n)$, so
that
$nv_{1}\leq{\lambda_{min}(\boldsymbol{X}_{s}^{T}\boldsymbol{X}_{s})}\leq{\lambda_{max}(\boldsymbol{X}_{s}^{T}\boldsymbol{X}_{s})}\leq{nv_{2}}$
for any model of size $|s|\leq{\widetilde{p}}$.
A5. $||\beta_{0}||_{\infty}={E}$ and E is some positive number independent of
$n$

######  Theorem 3.1.

(Posterior contraction rates). Let $\epsilon_{n}=\sqrt{s_{0}\log p_{n}/n}$ and
suppose that assumptions A1-A5 hold. Under the linear regression model with
unknown Gaussian noise, we endow $\boldsymbol{\beta}$ with the exponential
power prior
$\Pi_{n}(\boldsymbol{\beta})=\frac{\lambda_{n}^{\frac{1}{\alpha}}}{2\Gamma(1+\frac{1}{\alpha})}\exp{\left(-\lambda_{n}|\boldsymbol{\beta}|^{\alpha}\right)}$,
where $0<\alpha\leq\frac{1}{2}$ and $\lambda_{n}\propto(\log p_{n})$ and
$\sigma^{2}\sim\mathrm{InvGamma}(\delta,\delta)$. Then

$\displaystyle\Pi_{n}\left(\boldsymbol{\beta}:\left\|\boldsymbol{\beta}-\boldsymbol{\beta}_{0}\right\|_{2}\gtrsim\epsilon_{n}\mid\boldsymbol{Y}\right)\rightarrow
0$
$\displaystyle\Pi_{n}\left(\boldsymbol{\beta}:\left\|\boldsymbol{X}\boldsymbol{\beta}-\boldsymbol{X}\boldsymbol{\beta}_{0}\right\|_{2}\gtrsim\sqrt{n}\epsilon_{n}\mid\boldsymbol{Y}\right)\rightarrow
0$
$\displaystyle\Pi_{n}\left(\sigma^{2}:\left|\sigma-\sigma_{0}\right|\gtrsim\sigma_{0}\epsilon_{n}\mid
Y\right)\rightarrow 0$

where $\Pi_{n}(\cdot\mid\boldsymbol{Y})$ denote the posterior distribution
under the prior $\Pi_{n}$.

######  Theorem 3.2.

(Dimensionality). We define the generalized dimension as

$\gamma_{j}(\beta_{j})=I(|\beta_{j}|>a_{n})\,\,\,\text{and}\,\,\,|\gamma(\boldsymbol{\beta})|=\sum_{j=1}^{p_{n}}\gamma_{j}(\beta_{j})$

where $a_{n}=\frac{\sigma_{0}}{\sqrt{2k}}\sqrt{s_{0}\log p_{n}/n}/p_{n}$ and
$\lim_{n\rightarrow{\infty}}a_{n}=0$.Then under the assumptions A1-A5, for
sufficient large $M_{3}>0$, we have

$\sup_{\theta_{0}}P_{\theta_{0}}^{n}\Pi\left(\boldsymbol{\beta}:|\mathcal{\gamma}(\boldsymbol{\beta})|>M_{3}s_{0}\mid\boldsymbol{Y}\right)\rightarrow
0$

.

Details of proofs for Theorems 3.1 and 3.2 can be found in section 6 of the
supplementary materials.

##  4 MCMC for Bayesian $L_{\frac{1}{2}}$ regression

For the linear model with gaussian error of the form (1), together with the
global-local $L_{\frac{1}{2}}$ shrinkage prior of the form (2), and the hyper-
prior (3), the full conditional posterior distributions are given by

$\displaystyle\boldsymbol{\beta}|\sigma^{2},\boldsymbol{\tau}^{2},\boldsymbol{v},\lambda,b,\boldsymbol{Y}\sim\mathrm{N_{p}}\left((\boldsymbol{X}^{T}\boldsymbol{X}+\sigma^{2}\lambda^{4}\boldsymbol{D}_{\boldsymbol{\tau}^{2}}^{-1})^{-1}\boldsymbol{X}^{T}\boldsymbol{Y},\sigma^{2}(\boldsymbol{X}^{T}\boldsymbol{X}+\sigma^{2}\lambda^{4}\boldsymbol{D}_{\boldsymbol{\tau}^{2}}^{-1})^{-1}\right)$
$\displaystyle
r_{j}=\frac{1}{\tau_{j}^{2}}|\boldsymbol{\beta},\boldsymbol{v},\lambda,b,\boldsymbol{Y}\sim\mathrm{InvGaussian}\left(\frac{1}{\lambda^{2}v_{j}|\beta_{j}|},\frac{1}{v_{j}^{2}}\right),\quad
j=1,...,p$
$\displaystyle\sigma^{2}|\boldsymbol{\beta},\boldsymbol{\tau}^{2},\boldsymbol{v},\lambda,b,\boldsymbol{Y}\sim\mathrm{InvGamma}\left(\frac{n}{2},\frac{1}{2}(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta})^{T}(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta})\right)$
$\displaystyle\lambda^{4}|\boldsymbol{\beta},\boldsymbol{\tau}^{2},\boldsymbol{v},\sigma^{2},b,\boldsymbol{Y}\propto(\lambda^{4})^{\frac{2p-3.5}{4}}\exp\left(-\lambda^{4}\sum_{j=1}^{p}\frac{\beta_{j}^{2}}{2\tau_{j}^{2}}-\frac{\lambda}{b}\right)$
$\displaystyle
b|\boldsymbol{\beta},\boldsymbol{\tau}^{2},\boldsymbol{v},\lambda,\sigma^{2},\boldsymbol{Y}\sim\mathrm{InvGamma}(1,1+\lambda)$
$\displaystyle
v_{j}|\boldsymbol{\tau}^{2},\boldsymbol{\beta},\sigma^{2},\lambda,\boldsymbol{Y}\propto(v_{j})^{-\frac{3}{2}}\exp\left(-\frac{\tau_{j}^{2}}{2v_{j}^{2}}-\frac{1}{4}v_{j}\right),\quad
j=1,...,p$

where
$\boldsymbol{D}_{\boldsymbol{\tau}^{2}}=\mathrm{diag}(\tau_{1}^{2},...,\tau_{p}^{2})$.
Here the full conditional distributions for $\boldsymbol{v}$ and $\lambda$ do
not follow any standard distributions, so would require additional sampling
strategies, such as the Metropolised-Gibbs sampler or adaptive rejection
sampling (Gilks and Wild,, 1992). The partially collapsed Gibbs sampler (PCG)
(Park and Van Dyk,, 2009) speeds up convergence by reducing the conditioning
in some or all of the component draws of its parent Gibbs sampler, which also
provides a closed-form conditional posterior in our case. In this section, we
develop a partially collapsed Gibbs sampler based on the approach of (Park and
Van Dyk,, 2009), first for the case with $\gamma=1$, then separately for the
case with $\gamma>1$.

###  4.1 Case 1: $L_{\frac{1}{2}}$ prior with $\gamma=1$

Van Dyk and Park, (2008) showed that, by applying three basic steps:
marginalize, permute and trim, one can construct the PCG sampler which
maintain the target joint posterior distribution as the stationary
distribution.

Following the prescription in Van Dyk and Park, (2008), the first step of the
PCG sampler is marginalisation, where we sample the $\tau_{j}^{2}$ and
$v_{j}$’s multiple times in the Gibbs cycle without perturbing the stationary
distribution. In the second step we permute the ordering of the update which
helps us identify trimming strategies. Based on the final trimming step, we
are able to construct a sequential sampling scheme to sample the posterior
where each update is obtained from a closed-form distribution with some
conditional components being removed. We provide the exact sampling scheme
below in Algorithm 1:

Initialise the parameters
$\boldsymbol{\beta},\lambda,\boldsymbol{v},\boldsymbol{\tau}^{2},\sigma^{2}$
and set $t=1$, sequentially update the parameters until at $t=T$, via
$\displaystyle{\bf S1.}$ $\displaystyle\mbox{Sample
}\boldsymbol{\beta}^{\prime}\sim\mathrm{N_{p}}\left(({\boldsymbol{X}}^{T}{\boldsymbol{X}}+\sigma^{2}\lambda^{4}{\boldsymbol{D}}_{\boldsymbol{\tau}^{2}}^{-1})^{-1}{\boldsymbol{X}}^{T}\boldsymbol{Y},\sigma^{2}(\boldsymbol{X}^{T}\boldsymbol{X}+\sigma^{2}\lambda^{4}\boldsymbol{D}_{\boldsymbol{\tau}^{2}}^{-1})^{-1}\right)$
$\displaystyle{\bf S2.}$ $\displaystyle\mbox{Sample
}\lambda^{\prime}\sim\mathrm{Gamma}\left(2p+0.5,\sum_{j=1}^{p}|\beta^{\prime}_{j}|^{\frac{1}{2}}+\frac{1}{b}\right)$
$\displaystyle{\bf S3.}$ $\displaystyle\mbox{Sample
}\frac{1}{v^{\prime}_{j}}\sim\mathrm{InvGaussian}\left(\sqrt{\frac{1}{4\lambda^{\prime
2}|\beta^{\prime}_{j}|}},\frac{1}{2}\right),\quad j=1,...,p$
$\displaystyle{\bf S4.}$ $\displaystyle\mbox{Sample
}\frac{1}{{\tau^{\prime}}_{j}^{2}}\sim\mathrm{InvGaussian}\left(\frac{1}{{\lambda^{\prime}}^{2}{v^{\prime}}_{j}|\beta^{\prime}_{j}|},\frac{1}{{v^{\prime}}_{j}^{2}}\right),\quad
j=1,...,p$ $\displaystyle{\bf S5.}$ $\displaystyle\mbox{Sample }\sigma^{\prime
2}\sim\mathrm{InvGamma}\left(\frac{n}{2},\frac{1}{2}(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}^{\prime})^{T}(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}^{\prime})\right)$
$\displaystyle{\bf S6.}$ $\displaystyle\mbox{Sample
}b^{\prime}\sim\mathrm{InvGamma}\left(1,1+\lambda^{\prime}\right)$
$\displaystyle{\bf S7.}$ $\displaystyle\mbox{Set
}\boldsymbol{\beta},\lambda,\boldsymbol{v},\boldsymbol{\tau}^{2},\sigma^{2}\mbox{
to
}\boldsymbol{\beta}^{\prime},\lambda^{\prime},\boldsymbol{v}^{\prime},\boldsymbol{\tau}^{\prime
2},\sigma^{\prime 2}\mbox{ and increment }t=t+1.$

Algorithm 1 PCG $L_{\frac{1}{2}},\gamma=1$ update algorithm

The details of the three steps and derivations for
$\pi(\lambda|\boldsymbol{\beta},\sigma^{2},\boldsymbol{Y})$ and
$\pi(h_{j}=\frac{1}{v_{j}}|\beta_{j},\sigma^{2},\lambda,\boldsymbol{Y})$ are
given in Section 2 of Supplementary Materials.

###  4.2 Case 2: $L_{\frac{1}{2}}$ prior with $\gamma\in\mathbb{N}^{+}$

We can similarly design a PCG sampler for the more general
$\alpha=(\frac{1}{2})^{\gamma},\gamma\in\mathbb{N}^{+}$ case. The construction
here is slightly more tedious but within the same spirit as $\gamma=1$ case.
Under $\gamma>1$, we have an expanded parameter space with
$\boldsymbol{v}=\\{v_{i,j}\\},i=1,\ldots,\gamma,j=1,\ldots,p$, We omit the
details and refer the reader to section 2 of supplementary materials for the
intuition under $\gamma=1$. We describe the algorithm for $\gamma>1$ below

Initialise the parameters
$\boldsymbol{\beta},\lambda,\boldsymbol{v},\boldsymbol{\tau}^{2},\sigma^{2}$
and set $t=1$, sequentially update the parameters until at $t=T$, via
$\displaystyle{\bf S1.}$ $\displaystyle\mbox{Sample
}\boldsymbol{\beta}^{\prime}\sim\mathrm{N}_{p}\left((\boldsymbol{X}^{T}\boldsymbol{X}+\sigma^{2}\lambda^{2^{\gamma+1}}\boldsymbol{D}_{\boldsymbol{\tau}^{2}}^{-1})^{-1}\boldsymbol{X}^{T}\boldsymbol{Y},\sigma^{2}(\boldsymbol{X}^{T}\boldsymbol{X}+\sigma^{2}\lambda^{2^{\gamma+1}}\boldsymbol{D}_{\tau^{2}}^{-1})^{-1}\right)$
$\displaystyle{\bf S2.}$ $\displaystyle\mbox{Sample
}\lambda^{\prime}\sim\mathrm{Gamma}\left(2^{\gamma}p+0.5,\sum_{j=1}^{p}\left|\beta^{\prime}_{j}\right|^{\frac{1}{2^{\gamma}}}+\frac{1}{b}\right)$
$\displaystyle{\bf S3.}$ $\displaystyle\mbox{Sample
}\frac{1}{{v^{\prime}}_{\gamma,j}}\sim\mathrm{InvGaussian}\left(\frac{1}{2\lambda^{\prime}|\beta^{\prime}_{j}|^{\frac{1}{2^{\gamma}}}},\frac{1}{2}\right),\quad
j=1,...,p$ $\displaystyle\mbox{For }i=1,\ldots,\gamma-1,$
$\displaystyle\mbox{Sample
}\frac{1}{{v^{\prime}}_{i,j}}\sim\mathrm{InvGaussian}\left(\frac{1}{2v^{\prime}_{i+1,j}\lambda|\beta^{\prime}_{j}|^{\frac{1}{2^{i}}}},\frac{1}{2{v^{\prime}}_{i+1,j}^{2}}\right),\quad
j=1,...,p$ $\displaystyle{\bf S4.}$ $\displaystyle\mbox{Sample
}\frac{1}{{\tau^{\prime}}_{j}^{2}}\sim\mathrm{InvGaussian}\left(\frac{1}{{\lambda^{\prime}}^{2^{\gamma}}{v^{\prime}}_{1,j}|\beta^{\prime}_{j}|},\frac{1}{{v^{\prime}}_{1,j}^{2}}\right),\quad
j=1,...,p$ $\displaystyle{\bf S5.}$ $\displaystyle\mbox{Sample
}{\sigma^{\prime}}^{2}\sim\mathrm{InvGamma}\left(\frac{n}{2},\frac{1}{2}(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}^{\prime})^{T}(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}^{\prime})\right)$
$\displaystyle{\bf S6.}$ $\displaystyle\mbox{Sample
}b^{\prime}\sim\mathrm{InvGamma}\left(1,1+\lambda’\right)$ $\displaystyle{\bf
S7.}$ $\displaystyle\mbox{Set
}\boldsymbol{\beta},\lambda,\boldsymbol{v},\boldsymbol{\tau}^{2},\sigma^{2}\mbox{
to
}\boldsymbol{\beta}^{\prime},\lambda^{\prime},\boldsymbol{v}^{\prime},\boldsymbol{\tau}^{\prime
2},\sigma^{\prime 2}\mbox{ and increment }t=t+1.$

Algorithm 2 PCG $L_{\frac{1}{2}},\gamma>1$ update algorithm

Note that this scheme also works for $\gamma=0$ (Bayesian LASSO) case. In this
case, step 3 is removed and step 4 is replaced by

${\bf S4}.\quad\mbox{Sample
}\frac{1}{{\tau^{\prime}}_{j}^{2}}\sim\mathrm{InvGaussian}\left(\frac{1}{\lambda|\beta^{\prime}_{j}|},1\right),\quad
j=1,...,p.$

##  5 Variable selection: posterior mode search via coordinate desecent
optimization

Whilst the PCG sampler presented in section 4 produces full posterior
distributions for all the parameters, it does not produce exact zero’s for the
regression coefficients $\boldsymbol{\beta}$, making variable selection
difficult. It may be possible to set a threshold value to separate the
important and unimportant variables, however it is not clear how the threshold
can be chosen optimally. If the goal of inference is to identify important
variables very quickly, we propose below a fast optimisation strategy which
searches for the posterior modes and is capable of producing sparse solutions
suitable for variable selection.

Coordinate descent (CD) type of algorithms have been developed for bridge
penalties previously, see for example, Marjanovic and Solo, (2012, 2013,
2014). However, the bridge penalty is limited by its lack of ability to adapt
to the sparsity pattern across the coordinates. Without knowing the true
sparsity, specifying the global shrinkage parameter $\lambda$ can be
challenging. In the spirit of the full Bayesian treatment, by assigning a
suitable hyper-prior for $\lambda$ to the bridge penalty, we can work with the
marginalized the hyper-parameter $\lambda$, which can achieve a level of
adaptiveness. A similar approach was used by Ročková and George, (2016, 2018)
for the spike and slab LASSO penalty.The marginalised log posterior
distribution can be written as:

$\displaystyle\log\pi\left(\boldsymbol{\beta},b,\sigma^{2}\right|\boldsymbol{Y})=$
$\displaystyle-\frac{1}{2\sigma^{2}}\left\|Y-X\boldsymbol{\beta}\right\|_{2}^{2}-(n+2)\log\sigma+\log\int\pi(\boldsymbol{\beta}|\lambda)\pi(\lambda|b)d\lambda+\log\pi(b)$
(4) $\displaystyle=$
$\displaystyle-\frac{1}{2\sigma^{2}}\left\|Y-X\boldsymbol{\beta}\right\|_{2}^{2}-(n+2)\log\sigma-\mathbf{pen_{b}}(\boldsymbol{\beta})-\log
C_{b}+\log\pi(b)$

where $b$ is the hyper-parameter in (3) for $\frac{1}{\sqrt{\lambda}}$ and
$\lambda$ is marginalised out. Finally, Jeffreys’ prior is used for
$\sigma^{2}$ in (4). The prior $\pi(\boldsymbol{\beta}|b)$ is now written in
terms of the penalty function
$\mathbf{pen_{b}}(\boldsymbol{\beta})=-\log\int\pi(\boldsymbol{\beta}|\lambda)\pi(\lambda|b)d\lambda-\log
C_{b}$, where $C_{b}$ is the normalizing constant which depends on $b$. The
posterior mode can then be obtained by maximising (4) with respect to
$\boldsymbol{\beta}$, $b$ and $\sigma^{2}$

$(\boldsymbol{\beta}_{MAP},b_{MAP},\sigma_{MAP}^{2})=\arg\max_{\boldsymbol{\beta},b,\sigma^{2}}\log\pi(\boldsymbol{\beta},b,\sigma^{2}|\boldsymbol{Y}),$

which can be achieved by iteratively updating the
$\sigma^{2},\boldsymbol{\beta},b$ in turn. In fact,the parameter $b$ is not of
interest to be estimated. It is also very dangerous to do that in high
dimension and very sparse problem(see section 5.1), however it is not possible
to analytically integrate it out as for $\lambda$. In section 5.4, We will
discuss a strategy to avoid iteratively updating $b$.

When $\sigma^{2}$ is unknown, the non-convexity of the gaussian negative log-
likelihood is a well-known difficulty for optimisation (Bühlmann and Van De
Geer,, 2011). For our model, the iterative optimisation algorithm appear to
work well in low dimensional cases $n>p$, but often shrinks everything to zero
when $n<p$. However, it appears to works very well when $\sigma^{2}$ is known.
It is not clear why variance estimation can sometimes fail in high dimensions.
Here we suggest the following simple strategy to tackle this issue.

If we treat $\sigma_{0}^{2}$ as an unknown constant in the penalized
regression with separable penalty, the objective function can be re-written as

$\arg\min_{\boldsymbol{\beta}}\frac{1}{2}\|\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}\|^{2}+\sigma_{0}^{2}\lambda\sum_{j=1}^{p}\mathbf{pen}(\beta_{j})$

and define $\lambda^{{}^{\prime}}=\sigma_{0}^{2}\lambda$ as a new penalty
parameter. We assign the same prior to $\lambda^{{}^{\prime}}$ instead of
$\lambda$ in (3). Then, if we marginalize over $\lambda^{{}^{\prime}}$, we
have a simpler objective function to (4) with the same penalty function.

$\log\pi\left(\boldsymbol{\beta},b\right|\boldsymbol{Y})=-\frac{1}{2}\left\|Y-X\boldsymbol{\beta}\right\|_{2}^{2}-\mathbf{pen_{b}}(\boldsymbol{\beta})-\log
C_{b}+\log\pi(b)$ (5)

where
$\mathbf{pen_{b}}(\boldsymbol{\beta})=-\log\int\pi(\boldsymbol{\beta}|\lambda^{\prime})\pi(\lambda^{\prime}|b)d\lambda^{\prime}-\log
C_{b}=-\log\int\pi(\boldsymbol{\beta}|\lambda)\pi(\lambda|b)d\lambda-\log
C_{b}$ and there is no parameter $\sigma^{2}$ in the objective function. We
then follow a two stage approach.In stage one, we find
$\hat{\boldsymbol{\beta}}$ and $\hat{s}$ (the number of nonzero elements in
$\hat{\boldsymbol{\beta}}$) by maximising (5). In stage two, after finishing
the optimization, we consider the variance estimator:

$\hat{\sigma}^{2}=\frac{\|Y-X\hat{\beta}\|_{2}^{2}}{n-\hat{s}}.$ (6)

Fan and Lv, (2011) studied oracle properties of the non-concave penalized
likelihood estimator in the high dimensional setting. Fan et al., (2012)
showed that the variance can be consistently and efficiently estimated by
$\hat{\sigma}_{\lambda}^{2(SCAD)}$ of the form in (6). The main difference in
our approach here is that rather than using cross-validation to determine
$\lambda^{\prime}$, we marginalised over $\lambda^{\prime}$.

###  5.1 The KKT condition and coordinate-descent algorithm

For the optimisation of $\boldsymbol{\beta}$, Ročková and George, (2016) gave
the following necessary Karush-Kuhn-Tucker (KKT) condition for the global
mode:

$\beta_{j}=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\left(|z_{j}|-\frac{\partial\operatorname{\bf
pen}(\boldsymbol{\beta})}{\partial|\beta_{j}|}\right)_{+}\mathrm{sign}(z_{j})$
(7)

where
$z_{j}=\boldsymbol{X}_{j}^{T}(\boldsymbol{Y}-\boldsymbol{X}_{-j}\beta_{-j})$
and $(z)_{+}=max(0,z)$. The derivatives of $\operatorname{\bf
pen}(\boldsymbol{\beta})$ plays a crucial role here. For example, for
separable priors where $\lambda$ is treated as fixed, then under the Laplace
priors, the solution is the traditional LASSO estimator
$\beta_{j}=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\left(|z_{j}|-\sigma^{2}\lambda\right)_{+}\mathrm{sign}(z_{j})$.
For the non-separable prior, we integrate out the hyper-parameter $\lambda$
with respect to the $L_{(\frac{1}{2})^{\gamma}}$ prior, which gives us

$\displaystyle\pi(\beta_{1},...,\beta_{p}|b)$
$\displaystyle=\frac{(2^{\gamma}p-0.5)!}{\sqrt{\pi
b}2^{p}[(2^{\gamma})!]^{p}}\frac{1}{(\sum_{j=1}^{p}|\beta_{j}|^{\frac{1}{2^{\gamma}}}+1/b)^{2^{\gamma}p+0.5}}.$

Then

$\frac{\partial\operatorname{\bf
pen}_{b}(\boldsymbol{\beta})}{\partial|\beta_{j}|}=-\frac{\partial\log\pi(\boldsymbol{\beta}|b)}{\partial|\beta_{j}|}=\frac{(2^{\gamma}p+0.5)|\beta_{j}|^{\frac{1}{2^{\gamma}}-1}}{2^{\gamma}(\sum_{j=1}^{p}|\beta_{j}|^{\frac{1}{2^{\gamma}}}+\frac{1}{b})}=\left(p+\frac{1}{2^{\gamma+1}}\right)\frac{1}{|\beta_{j}|^{1-\frac{1}{2^{\gamma}}}}\frac{1}{\sum_{j=1}^{p}|\beta_{j}|^{\frac{1}{2^{\gamma}}}+\frac{1}{b}},$

and thus we have

$\beta_{j}=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\left(|z_{j}|-\frac{C_{1}}{|\beta_{j}|+C_{2}|\beta_{j}|^{1-\frac{1}{2^{\gamma}}}}\right)_{+}\mathrm{sign}(z_{j})$
(8)

where we set

$C_{1}=\left(p+\frac{1}{2^{\gamma+1}}\right),\quad C_{2}=\sum_{i\neq
j}|\beta_{i}|^{\frac{1}{2^{\gamma}}}+\frac{1}{b}.$ (9)

For $\gamma>0$, the shrinkage term $\frac{\partial\operatorname{\bf
pen}_{b}(\boldsymbol{\beta})}{\partial|\beta_{j}|}$ is global-local adaptive,
it depends on the information both from itself and from the other coordinates.
This is consistent with the framework of global-local shrinkage priors (Polson
and Scott,, 2012). While for $\gamma=0$, the shrinkage term
$\frac{\partial\operatorname{\bf
pen}_{b}(\boldsymbol{\beta})}{\partial|\beta_{j}|}=\frac{1}{\sum_{j=1}^{p}|\beta_{j}|+1/b}$,
does not have the local adaptive properties.

However, the KKT condition given by Equation (7) does not apply to penalties
with unbounded derivatives. In fact, the set of $\hat{\beta}$ satisfying
Equation (7) is only a superset of the solution set for our non-separable
penalty. We can see that for non-separable bridge penalty with $\gamma>0$, the
solution $\beta_{j}=0$ always satisfies equation (8). Here, we derive the true
KKT condition for our penalty, which serves as a basis for the CD algorithm.

######  Theorem 5.1.

Let the constants $\tilde{\beta}_{j}$ and $h(\tilde{\beta}_{j})$, be such that
they satisfy the following conditions

$\displaystyle 1$
$\displaystyle=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}C_{1}\frac{1+(1-\frac{1}{2^{\gamma}})C_{2}\tilde{\beta}_{j}^{-\frac{1}{2^{\gamma}}}}{\tilde{\beta}_{j}^{2}\left(1+C_{2}\tilde{\beta}_{j}^{-\frac{1}{2^{\gamma}}}\right)^{2}}$
(10) $\displaystyle h(\tilde{\beta}_{j})$
$\displaystyle=\tilde{\beta}_{j}+(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\tilde{\beta}_{j}+C_{2}\tilde{\beta}_{j}^{(1-\frac{1}{2^{\gamma}})}}.$

Then the solutions satisfy

$\beta_{j}=T(\boldsymbol{\beta}_{-j},\beta_{j})=\left\\{\begin{array}[]{ll}0&\text{
if
}(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|<h(\tilde{\beta}_{j})\\\
0&\text{ if }(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|\geq
h(\tilde{\beta}_{j})\quad\text{and}\quad\delta_{-j}(\beta_{j}^{{}^{\prime\prime}})>0\\\
\operatorname{sign}(z_{j})\beta_{j}^{{}^{\prime\prime}}I(\beta_{j}\neq
0)&\text{ if }(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|\geq
h(\tilde{\beta}_{j})\quad\text{and}\quad\delta_{-j}(\beta_{j}^{{}^{\prime\prime}})=0\\\
\operatorname{sign}(z_{j})\beta_{j}^{{}^{\prime\prime}}&\text{ if
}(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|\geq
h(\tilde{\beta}_{j})\quad\text{and}\quad\delta_{-j}(\beta_{j}^{{}^{\prime\prime}})<0\end{array}\right.$
(11)

where
$\delta_{-j}(\beta_{j})=L_{-j}(\beta_{j})-L_{-j}(0)$,$\beta_{j}^{{}^{\prime\prime}}\in[\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]$,
$C_{1}$ and $C_{2}$ are as given in (9) and $L_{-j}$ denotes the loss function
Equation (5) with all except the $j$th $\boldsymbol{\beta}$ fixed.

The proof of the theorem is provided in the section 3 of the supplementary
materials, where in the step one and step two of the proof, we also show the
following corollary:

######  Corollary 5.1.1.

$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|\geq
h(\tilde{\beta}_{j})$ if and only if $\beta_{j}^{{}^{\prime\prime}}$ exists
and can be computed by fixed point iteration:
$\beta_{j}^{(k+1)}=\rho(\beta_{j}^{(k)})$, where

$\rho(\beta_{j})=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|-(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\beta_{j}+C_{2}\beta_{j}^{1-\frac{1}{2^{\gamma}}}}$
(12)

with the initial condition
$\beta_{j}^{(0)}\in[\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]$.
When
$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|=h(\tilde{\beta}_{j})$,
$\beta_{j}^{{}^{\prime\prime}}=\tilde{\beta}_{j}$.

Remark 1: From Theorem 5.1, we see that the global mode is a blend of hard
thresholding and nonlinear shrinkage. $h(\tilde{\beta}_{j})$ is the selection
threshold where $\tilde{\beta}_{j}$ needs to be solved numerically according
to Equation (10). Rather than evaluating the selection threshold by solving
this difficult non-linear Equation directly, we suggest selecting the variable
by checking the convergence of the fixed point iterations.By corollary 5.1.1,
we can run the fixed point iteration with an initial guess
$\beta_{j}^{(0)}=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|$. If
it fails to converge after a long iteration or $\rho(\beta_{j}^{(k)})<0$, we
conclude that $\beta_{j}=0$.

However, this is still a computationally intensive. A single block of update
requires running the non-parallel fixed point iteration algorithms $p$ times
sequentially. The next lemma allows us to run the fixed point iteration only
on a small subset of variables. It is also a useful auxiliary result for us to
show the oracle property of the estimator later.

######  Lemma 5.2.

If the regressors haves been centered and standardized with $\|X_{j}\|^{2}=n$,
for $1\leq j\leq p$, then by Equation (10), for $\gamma\in\mathbb{N}^{+}$, we
always have $2\tilde{\beta}_{j}\leq h(\tilde{\beta}_{j})\leq
3\tilde{\beta}_{j}$ and

$\tilde{\beta}_{j}^{(2-\frac{1}{2^{\gamma}})}=\frac{C_{1}M^{\prime}(\gamma)}{\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j}}\frac{1}{\tilde{\beta}_{j}^{\frac{1}{2^{\gamma}}}+C_{2}}=M^{\prime}(\gamma)\frac{p+\frac{1}{2^{\gamma+1}}}{n\tilde{\beta}_{j}^{\frac{1}{2^{\gamma}}}+nC_{2}}$

for $M^{\prime}(\gamma)\in(\frac{1}{2},1]$. If the estimator satisfies
$\sum_{i=1}^{p_{n}}|\hat{\beta}_{j}|^{\frac{1}{2^{\gamma}}}=O(s_{0})$ and
$\|\hat{\beta}\|_{\infty}<\infty$, then
$\tilde{\beta}_{j}^{(2-\frac{1}{2^{\gamma}})}=O(\frac{p_{n}}{s_{0}n+nb^{-1}})$.
In addition, there exists another lower bound for the selection threshold
$u(\beta_{-j})<h(\tilde{\beta}_{j})$, which is not a function of
$\tilde{\beta}_{j}$.

$u(\beta_{-j})=2\left\\{\frac{C_{1}(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}}{2C_{2}+2[(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]^{\frac{1}{2^{\gamma}}}}\right\\}^{\frac{1}{2-\frac{1}{2^{\gamma}}}}.$
(13)

The proof of Lemma 5.2 is in section 3 of supplementary materials.
Remark 2: Lemma 5.2 provides an explicit lower bound for the selection
threshold $h(\tilde{\beta}_{j})$, we see that $\hat{\beta}_{j}=0$ if
$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|\leq u(\beta_{-j})$.
Therefore, we only need to run the fixed point iterations in
$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|>u(\beta_{-j})$ cases.
To further speed up, we suggest using the null model $\beta_{Initial}=0$ as
initialization. Then in the early stages of the updates, the lower bound (13)
will be large, which means we only need to run fixed point iteration in a very
small set of variables.

We iteratively updating $\hat{b}$ by taking
$\frac{\partial\log\pi\left(\hat{\boldsymbol{\beta}},b\right|\boldsymbol{Y})}{\partial
b}=0$ , which leads to solving the nonlinear equation:

$\frac{1}{b}=\frac{1}{2b^{2}}+\frac{2^{\gamma}p+0.5}{2b^{2}\sum_{j=1}^{p}|\hat{\beta_{j}}|^{\frac{1}{2^{\gamma}}+}+2b}.$

Since the equation above implies the inequality

$\frac{2^{\gamma}p-1.5}{2\sum_{j=1}^{p}|\hat{\beta_{j}}|^{\frac{1}{2^{\gamma}}}}\leq
b\leq\frac{1}{2}+\frac{2^{\gamma}p+0.5}{2\sum_{j=1}^{p}|\hat{\beta_{j}}|^{\frac{1}{2^{\gamma}}}}$

we can approximate the solution by

$\hat{b}^{(k)}\approx\frac{2^{\gamma-1}p}{\sum_{j=1}^{p}|\hat{\beta_{j}}^{(k)}|^{\frac{1}{2^{\gamma}}}}.$
(14)

Remark 3: This approach can be viewed as an empirical Bayes approach as the
hyper-parameter $b$ is learned from the data. Of most concern is that the
empirical Bayes estimator has the risk of a degenerate solution (e.g., setting
all $\beta_{j}=0$) in the high dimensional and very sparse setting, resulting
in an inappropriate statement of the sparsity of the regression model (Scott
and Berger,, 2010; Polson and Scott,, 2012; Datta et al.,, 2013). Next
corollary demonstrates this issue.

######  Corollary 5.2.1.

Suppose $\sum_{i=1}^{p}|\beta_{0,j}|^{\frac{1}{2^{\gamma}}}=O(s_{0})$ and
$||\boldsymbol{\beta}_{0}||_{\infty}={E}$, using the empirical Bayesian
estimator (14)
$\hat{b}=\frac{2^{\gamma-1}p}{\sum_{i=1}^{p}|\hat{\beta_{i}}|^{\frac{1}{2^{\gamma}}}}$
for $b$, then the model will fail to recover the signal if
$\lim_{n\rightarrow{\infty}}\frac{p_{n}}{ns_{0}}=\infty$.

Proof. Suppose there exists estimator $\hat{\boldsymbol{\beta}}$ which can
recover the signal, then
$\sum_{i=1}^{p}|\hat{\beta}_{j}|^{\frac{1}{2^{\gamma}}}=\sum_{i=1}^{p}|\beta_{0,j}|^{\frac{1}{2^{\gamma}}}=O(s_{0})$.
By Lemma 5.2,
$h(\tilde{\beta}_{j})=O\left(\left[\frac{p_{n}}{ns_{0}(1+1/2^{\gamma-1}p_{n})}\right]^{\frac{1}{2-\frac{1}{2^{\gamma}}}}\right)$.
Since $\lim_{n\rightarrow{\infty}}\frac{p_{n}}{ns_{0}}=\infty$,
$h(\tilde{\beta}_{j})\rightarrow\infty$ as $n\rightarrow{\infty}$. This
implies that the model will shrink all the elements to zero, which contradicts
the assumption that $\hat{\boldsymbol{\beta}}$ can recover the signal.

This phenomenon is empirically observed in Figure 1.This figure shows the
solution path of $\boldsymbol{\beta}$ in the same setting as Figure 2 in
secetion 5.1. Here, seven of the ten nonzero elements have been estimated to
be zero. One way to remedy this issue is to set $b\propto\frac{\log
p_{n}}{p_{n}}$ (see Figure 2). In this case, the oracle properties of the
$\hat{\beta}$ and $\hat{\sigma}^{2}$ are guaranteed by Theorem 5.5.The exact
value of hyper-parameter $b$(determine the unknown constant) can be determined
by the variable screening strategy as discuss in section 5.4. Here we provide
the details of our CD optimisation algorithm in Algorithm 3:

Input $\epsilon>0$; $\beta^{1}$; $b>0$ ; $\gamma>0$

while  _$||\boldsymbol{\beta}^{(i)}-\boldsymbol{\beta}^{(i-1)}|| >\epsilon$ _
do

for  _$j\leftarrow 1$ to $P$_ do

$z_{j}=\boldsymbol{X}_{j}^{T}(\boldsymbol{Y}-\boldsymbol{X}_{-j}{\boldsymbol{\beta}}^{i}_{-j})$

if  _$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|\leq
u(\beta_{-j})$_ then

$\beta_{j}^{(i+1)}=0$

else

$\beta_{j}^{(i+1,0)}=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|$

$k=1$

while  _$|\beta_{j}^{(i,k)}-\beta_{j}^{(i,k-1)}|
>\epsilon\quad\text{and}\quad\beta_{j}^{(i,k)}>0\quad\text{and}\quad k\leq T$_
do

$C_{2}=||\boldsymbol{\beta}_{-j}^{i}||^{\frac{1}{2\gamma}}_{\frac{1}{2\gamma}}+\frac{1}{b}$

$\beta_{j}^{(i,k+1)}=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|-(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\beta^{(i,k)}_{j}+C_{2}(\beta^{(i,k)}_{j})^{1-\frac{1}{2^{\gamma}}}}$

$k=k+1$

end while

if  _$\beta_{j}^{(i,k)} <0\quad\text{or}\quad k>T\quad\text{or}\quad
L_{-j}(0)>L_{-j}(\operatorname{sign}(z_{j})\beta_{j}^{(i,k)})$_ then

${\beta}^{(i+1)}_{j}=0$;

else

${\beta}^{(i+1)}_{j}=\operatorname{sign}(z_{j})\beta_{j}^{(i,k)}$

end if

end if

end for

$i=i+1$

end while

Result: ${\boldsymbol{\beta}}=\\{\beta^{(i)}_{1},\ldots,\beta^{(i)}_{p}\\}$

Algorithm 3 Coordinate descent (CD) optimization

where $\epsilon$ is some suitable small error tolerance, $T$ is the terminal
number for the fixed point iteration, and $b$ is the value of the hyper-
parameter chosen by the user.
$\boldsymbol{\beta}_{-j}^{i}=(\beta_{1}^{(i+1)},...,\beta_{j-1}^{(i+1)},\beta_{j+1}^{(i)},...,\beta_{p}^{(i)}).$

Figure 1: Solution paths of the CD algorithm. Data is simulated using
$n=100,p=1000,s_{0}=10$ with iteratively update $b$ by
$\hat{b}=\frac{2^{\gamma-1}p}{\sum_{i=1}^{p}|\hat{\beta_{i}}|^{\frac{1}{2^{\gamma}}}}$.
Left panel shows the solution paths for the ten non-zero elements of
$\boldsymbol{\beta}$, middle panel shows the solutions paths for the 990 zero
elements of $\boldsymbol{\beta}$ and the right panel shows the path for the
loss function $L(\boldsymbol{\beta})$.

###  5.2 Convergence analysis

In general, convergence of the loss function alone cannot guarantee the
convergence of $\left\\{\boldsymbol{\beta}^{i}\right\\}_{i\geq 1}$. Mazumder
et al., (2011) showed the convergence of CD algorithms for a subclass of non-
convex non-smooth penalties for which the first and second order derivatives
are bounded. The authors also observed that for the type of log-penalty
proposed by Friedman, (2012), the CD algorithm can produce multiple limit
points without converging.

The derivative of the bridge penalty and its non-separable version are
unbounded at zero, suggesting a local solution at zero to the optimization
problem, this causes discontinuity in the induced threshed operators.
Recently, Marjanovic and Solo, (2014) proved the convergence of the CD
algorithm for the bridge penalty. Here, we show that the CD algorithm with the
proposed non-separable bridge penalty also converges. We provide convergence
analysis in the next two theorems. Details of the proofs can be found in
section 7 of supplementary materials.

######  Theorem 5.3.

Suppose the data $(\boldsymbol{Y},\boldsymbol{X})$ lies on compact set and the
sequence $\left\\{\boldsymbol{\beta}^{i}\right\\}_{i\geq 1}$ is generated by
$\beta_{j}^{(i+1)}=T(\boldsymbol{\beta}_{-j}^{(i)},\beta_{j}^{(i)})$ where the
$T(\cdot)$ map is given by equation (11) and
$\boldsymbol{\beta}_{-j}^{(i)}=(\beta_{1}^{(i+1)},...,\beta_{j-1}^{(i+1)},\beta_{j+1}^{(i)},...,\beta_{p}^{(i)})$,
then $\|\boldsymbol{\beta}^{i}-\boldsymbol{\beta}^{i+1}\|\rightarrow{0}\mbox{
as }i\rightarrow\infty.$

######  Theorem 5.4.

Let $\boldsymbol{\beta}^{\infty}$ be the estimator obtained by the CD
algorithm 3. Suppose $\boldsymbol{X}$ satisfies A4 with
$\|\boldsymbol{\beta}^{\infty}\|_{0}\leq\widetilde{p}$, then $\beta^{\infty}$
is a strict local minimizer of $L(\cdot)$. In other words, for any ${\bf
e}=(e_{1},e_{2},...,e_{p})\in\mathrm{R^{P}}$ and $\|{\bf e}\|=1$

$\lim_{\alpha\downarrow
0+}\inf_{\alpha}\left\\{\frac{L\left(\boldsymbol{\beta}^{\infty}+\alpha\boldsymbol{e}\right)-L\left(\boldsymbol{\beta}^{\infty}\right)}{\alpha}\right\\}>0.$

Proofs for the above theorems are given in Section 7 of supplementary
materials.

Empirically, we observe that convergence of the solution paths are pleasingly
well behaved. Figure 2 demonstrates solution paths of the CD algorithm with
different initialization strategies on a simulated sample of high dimensional
regression problem with 990 zero entries and 10 nonzero entries in $\beta$.
The design matrix $X$ is generated with pairwise correlation between $x_{i}$
and $x_{j}$ equal to $0.5^{|i-j|}$. Both the response and predictors have been
standardized, so $\sigma_{0}^{2}=1$. Data is simulated with $n=100$ and
$p=1000$. From figure 2, we see that two solution paths finally converge to
the same local minima with only two mistakes (they failed to identify one
nonzero element and to exclude one zero element).

Figure 2: Solution paths according to two initialisation strategies,
$\boldsymbol{\beta}_{initial}=0$ (top row) and random starting values
$\beta_{initial}\sim N_{p}(0,I_{p})$ (bottom row). Left panel shows the
solution paths for the (10) non-zero elements of $\boldsymbol{\beta}$, middle
panel shows the solutions paths for the (990) zero elements of
$\boldsymbol{\beta}$ and the right panel shows the path for the loss function
$L(\boldsymbol{\beta})$.

###  5.3 Oracle properties

In this section, we study statistical properties of the non-separable bridge
estimator. We will show that using $b\propto\frac{\log p_{n}}{p_{n}}$, model
consistency and asymptotic normality of $\hat{\sigma}^{2}$ can be achieved.We
state some additional assumptions here:
A6. $\sum_{j\in S_{0}}|\beta_{0j}|^{\frac{1}{2^{\gamma}}}=O(s_{0})$
A7. There exist a constant $0<c<\infty$ such that
$P\left(\left\|\frac{\boldsymbol{X}_{S_{0}^{c}}^{T}\boldsymbol{X}_{S_{0}}}{n}\right\|_{2,\infty}\leq
c\right)\rightarrow 1$ as $n\rightarrow{\infty}$
A8. $\min_{j\in
S_{0}}|\beta_{0,j}|\geq\max\left\\{3\sigma_{0}\sqrt{\frac{2\log
s_{0}}{v_{1}\log p_{n}}},3\sigma_{0}\left(\frac{\log
p_{n}}{n}\right)^{\frac{1}{2-\alpha}}\right\\}$ where
$\alpha\in(\frac{1}{2^{\gamma}},2)$
A9. $\min_{j\in
S_{0}}|\beta_{0,j}|\geq\max\left\\{3\sigma_{0}\sqrt{\frac{2\log
s_{0}}{v_{1}\log
p_{n}}},3\sigma_{0}\left(\frac{p_{n}}{s_{0}n}\right)^{\frac{1}{2-\alpha}}\right\\}$
where $\alpha\in(\frac{1}{2^{\gamma}},2)$

######  Theorem 5.5.

Suppose assumption A1-A4, A6 to A7 are satisfied. Then under the following two
cases,
(1): $\lim_{n\rightarrow{\infty}}\frac{p_{n}}{ns_{0}}=0$, $b\geq 1$ or
$b\propto\frac{\log p_{n}}{p_{n}}$ and assumption A8 holds.
(2): $\lim_{n\rightarrow{\infty}}\frac{p_{n}}{ns_{0}}=\infty$,
$b=O\left(\frac{\log p_{n}}{p_{n}}\right)$ and assumption A9 holds.
we have:
(a) Model consistency There is a strictly local minimizer $\hat{\beta}$ such
that

$\\#\left\\{j:\hat{\beta_{j}}\neq 0\right\\}=S_{0}$

with probability tending to 1 and,

(b) Asymptotic normality With the local estimator $\hat{\boldsymbol{\beta}}$,
the correponding variance estimator $\hat{\sigma}^{2}$ from Equation (6) has
the property that

$(\hat{\sigma}^{2}-\sigma_{0}^{2})\sqrt{n}\rightarrow
N(0,\sigma_{0}^{4}E[\epsilon^{4}]-\sigma_{0}^{4}).$

Proof for the theorem can be found in section 8 of supplementary Materials.

###  5.4 Forward and backward variable screening

There are many variable selection techniques proposed in the literature, under
a variety of assumptions, all of these $\hat{\beta}$ estimator are shown to
have variable selection consistency in high dimensional settings. Based on
these result, it is not difficult to construct the variance estimator, which
also has nice asymptotic properties. However, all the asymptotic arguments
require some crucial assumptions about the sparsity of the underlying true
model. For example, in this paper, it is assumed that the true number of
nonzero $\beta$ is $o(n/\log p_{n})$. In practice, we only have finite sample
size and we are rarely certain about the sparsity level of the true model. It
was found by Reid et al., (2016) that, the variance estimator from Fan et al.,
(2012); Sun and Zhang, (2012); Dicker, (2014) all suffer large biases in
finite samples when the true underlying signals become less sparse (in Section
6.2, you will see that the posterior distribution of the variance is biased
downwards). To overcome this bias, they suggest tuning the hyper-parameter in
LASSO or SCAD by cross validation, which performs well over a broad range of
sparsity and signal strength settings. Here, we propose a cross-validation
free backward screening procedure for moderate to large sample problem and a
forward screening procedure for small samples problem.

This is deployed as a path following warm up strategy. We begin with
small(backward screening starting with big model) or big (forward screening
starting with null model) value of hyper-parameter $b$. The solution is used
as warm-starts for next bigger or smaller value of $b$ and keep going with
this fashion. When sample size is moderate to large, the solution path with
respect to hyper-parameter $b$ from backward screening algorithm will
stabilize automatically. Thus, in this case, it is cross-validation free. When
the sample size is small, the solution path stabilizes slowly, it may
ultimately threshold everything to zero. In this case, we run the forward
screening with cross-validation to identify the hyper-parameter. Details of
the procedures are included in section 4 of supplementary materials.

##  6 Simulation Study

In this section, we present some results from a set of simulations studies.

###  6.1 Comparison of PCG sampler effective sample size

In this section, we benchmark the PCG sampler scheme for $L_{\frac{1}{2}}$
prior with $\gamma=0$ and $\gamma=1$ against the Gibbs sampler scheme of (Park
and Casella,, 2008) for the Bayesian LASSO and the Horseshoe prior (Makalic
and Schmidt,, 2015) in the high dimension sparse regression setting. We
included only these two schemes for comparison since in the high dimensional
setting, very few MCMC schemes have closed form condition posterior, which are
free of Metropolis -hastings proposal.

We simulate two sets of data with moderate correlation between the covariates,
($\rho_{i,j}=0.5^{|i-j|}$) and high ($\rho_{i,j}=0.8^{|i-j|}$) correlation for
the design matrix $\boldsymbol{X}$ with $p=1000$ variables and sample size
$n=100$. The true vector $\boldsymbol{\beta}_{0}$ is constructed by assigning
10 nonzero elements $(3,1.5,2,1,1,0.5,-0.5,2,-1.2,-1)$ to location
$(1,2,5,10,13,19,26,31,46,51)$ and setting all the remaining coefficients to
zero. We set the variance of the error to 1 ($\sigma_{0}^{2}=1$).

For each data set, we ran 10 MCMC chains. Each chain generated 10,000 samples
after an initial burn-in of 10,000. We calculated the effective sample sizes
(ESS) based on the formula from Gelman et al., (2013). In Table 1, we report
the maximum, minimum, median, average and standard deviation of ESS for
$\beta_{j}$ across $p=1000$ dimensions. Since all MCMC schemes considered here
are based on normal mixture representations, they all need to generate from
high dimensional gaussian distributions. All three schemes have the same
$O(np^{2})$ computational complexity (Bhattacharya et al.,, 2016). Thus, if
all their conditional posterior distribution have closed form expressions,
there is not too much difference between their computational times. Thus, we
only consider the indication about convergence rates, via their effective
sample sizes here.

The effective sample size from the table shows that for the Bayesian LASSO,
the PCG sampler ($\gamma=0$) performs similarly to the traditional Gibbs
sampler under all criterion, where the PCG sampler achieves a slightly higher
ESS overall. The performances of both samplers do not appear to be sensitive
to the increase in the correlation of the design matrix.

The PCG sampler for $L_{\frac{1}{2}}$ prior ($\gamma=1$) significantly over
performs the Horseshoe prior (using a Gibbs sampler) with respect to all the
criterions. Here increasing the correlation of the design matrix has had a
large effect on the ESS of both the Horseshoe Gibbs sampler and the PCG
sampler. For the Horseshoe prior, the ESS has decreased by as much as 10 fold
in some cases, and is particularly dramatic in the non-zero elements of
$\boldsymbol{\beta}$. The PCG sampler is also impacted, with the worst case
scenario being a seven fold decrease in ESS. However, comparing PCG to the
Horseshoe prior, it can be seen that the PCG sampler is up to ten times more
efficient in terms of ESS.

Finally, it should be pointed out that the proposed PCG sampler for the
$L_{\frac{1}{2}}$ prior has one limitation in practice. Theoretically, it
works for any $\gamma\in\mathbb{N}^{+}$, but in our python implementation, we
found that when $\gamma=2$, numerical underflow is quite frequently
encountered in the high dimension and very sparse setting. The problem worsens
when $\gamma>2$, occurring even in very simple low dimension problems. When
$\gamma\geq 2$, the $L_{\frac{1}{2}}$ prior places increasingly more
probability mass around zero. In this case, the full conditional posterior of
$\beta$, has a covariance matrix of the form
$\Sigma=\sigma^{2}(\boldsymbol{X}^{T}\boldsymbol{X}+\sigma^{2}\lambda^{2^{\gamma+1}}\boldsymbol{D}_{\tau^{2}}^{-1})^{-1}$,
and is nearly singular. It seems that the 64-bit floating-point numbers in the
numpy package from Python is not enough to keep the precision, therefore we
restrict our numerical experiments to the case $\gamma=1$ only.

| $\rho_{i,j}=0.5^{|i-j|}$ | $\rho_{i,j}=0.8^{|i-j|}$
---|---|---
ESS | Max | Min | Median | Average | SD | Max | Min | Median | Average | SD
| Bayesian LASSO Gibbs Sampler | Bayesian LASSO Gibbs Sampler
$\beta_{j}=0$ | 100053 | 24030 | 90496 | 86022 | 11814 | 100669 | 15986 | 92847 | 89595 | 10356
$\beta_{j}\neq 0$ | 86707 | 9806 | 38738 | 40636 | 21737 | 90138 | 10208 | 24451 | 36125 | 26440
All $\beta_{j}$ | 100053 | 9806 | 90317 | 85568 | 12779 | 100669 | 10208 | 92715 | 89060 | 11894
| PCG sampler $\gamma=0$ | PCG sampler $\gamma=0$
$\beta_{j}=0$ | 100831 | 26760 | 91208 | 87218 | 10949 | 101525 | 21476 | 92851 | 89802 | 10604
$\beta_{j}\neq 0$ | 90860 | 10500 | 54789 | 47604 | 32152 | 91965 | 10399 | 31069 | 37976 | 26254
All $\beta_{j}$ | 100831 | 10500 | 91018 | 86821 | 12023 | 101525 | 10399 | 92805 | 89284 | 12022
| Horseshoe Gibbs Sampler | Horseshoe Gibbs Sampler
$\beta_{j}=0$ | 40536 | 1418 | 18854 | 18178 | 9194 | 53821 | 866 | 23564 | 23551 | 13492
$\beta_{j}\neq 0$ | 6772 | 943 | 3482 | 3554 | 1427 | 5729 | 90 | 2840 | 3151 | 1716
All $\beta_{j}$ | 40536 | 943 | 18530 | 18326 | 9266 | 53821 | 90 | 23331 | 23347 | 13578
| PCG sampler $\gamma=1$ | PCG sampler $\gamma=1$
$\beta_{j}=0$ | 90564 | 12298 | 72514 | 67809 | 15355 | 92276 | 1614 | 69610 | 66323 | 16082
$\beta_{j}\neq 0$ | 37852 | 7166 | 20229 | 20355 | 10002 | 47088 | 1427 | 6768 | 13483 | 14627
All $\beta_{j}$ | 90564 | 7166 | 72093 | 67334 | 16022 | 92276 | 1427 | 69380 | 65794 | 16906

Table 1: Effective sample size of $\beta_{j}$ across $p=1000$ dimensions, for
two data sets generated from different degrees of correlation between
variables ($\rho_{i,j}=0.5^{|i-j|}$ and $\rho_{i,j}=0.8^{|i-j|}$. The final
MCMC sample contains 100,000 samples. The maximum, minimum, median average and
standard deviation of ESS over 10 MCMC runs and $\beta_{j}$ are reported for
PCG, LASSO and Horseshoe priors.

###  6.2 Performance of signal recovery in PCG and CD algorithms

In this section, we examine how $L_{\frac{1}{2}}$ prior recover signal under
high-dimensional settings. We compare both the PCG algorithm and CD
optimization algorithm with two well known Bayesian regression procedures, the
Bayesian LASSO (Park and Casella,, 2008) and the Horseshoe (Carvalho et al.,,
2010). We also consider the penalized likelihood procedures: the MC+ penalty
(Zhang et al.,, 2010) and non-separable spike-and-slab-lasso (NSSL) (Ročková
and George,, 2018).

Similarly with the previous section, we construct the AR(1) correlation
structure for design matrix $\boldsymbol{X}$ with correlation
$\rho_{ij}=0.5^{|i-j|}$. In Section 6.2.1, we test the performance of the
model in both low and high dimensional settings under the following five
setting when the number of nonzero regression coefficients is $s_{0}=10$.

$\displaystyle(1)\,n=500,p=25,\sigma_{0}^{2}=3\quad\,\,\,\,\,(2)\,n=500,p=1000,\sigma_{0}^{2}=1\quad(3)\,n=500,p=1000,\sigma_{0}^{2}=1$
$\displaystyle(4)\,n=100,p=1000,\sigma_{0}^{2}=1\quad(5)\,n=100,p=1000,\sigma_{0}^{2}=3$

For each of the above scenarios, we repeated the experiment 100 times. Each
time, the true vector $\boldsymbol{\beta}_{0}$ is constructed by assigning 10
nonzero elements $(3,1.5,2,1,1,0.5,-0.5,2,-1.2,-1)$ to random locations and
setting all the remaining $\beta$s to zero.Similarly, in Section 6.2.2, we
consider the small sample size and all zero setting i.e., $s_{0}=0$. In this
case, $\boldsymbol{\beta}_{0}=0$ for $p=1000$. We test two scenarios:
$(1)\,n=100,p=1000,\sigma_{0}^{2}=1\quad(2)\,n=100,p=1000,\sigma_{0}^{2}=3$
Finally, we also consider the very challenging case of small sample size and a
less sparse problem (Reid et al.,, 2016) in Section 6.2.3 with $s_{0}=20$,
under the following two scenarios:
$(1)\,n=100,p=1000,\sigma_{0}^{2}=1\quad(2)\,n=100,p=1000,\sigma_{0}^{2}=3$

The true vector $\boldsymbol{\beta}_{0}$ with 20 nonzero elements is
constructed by duplicating the 10 nonzero elements
$(3,1.5,2,1,1,0.5,-0.5,2,-1.2,-1)$ and assigning them to random locations and
setting all the remaining $\beta$s to zero.

For all the Bayesian models, we assign a non-conjugate Jeffrey’s prior to
$\sigma^{2}$. The Bayesian LASSO and $L_{\frac{1}{2}}$ prior use the hyper-
parameter setting as discussed in Section 2.2. For the Horseshoe prior, we
assign the hyper-prior $\tau\sim C^{+}(0,1)$ to the hyper-parameter as
suggested by Carvalho et al., (2010). We use posterior means of
$\boldsymbol{\beta}$ and $\sigma^{2}$ as the estimator. Since the MCMC cannot
provide a sparse solution directly, variables are selected by a hard
thresholding approach where we perform the t-test based on the posterior
samples, under the null hypothesis $\beta_{j}=0$ using $95\%$ as significant
level.

For the optimization procedure, the hyper-parameter $b$ for the NSB penalty is
chosen by using the forward/backward variable screening algorithm described in
Section 5.4. The NSSL(Ročková and George,, 2018) is claimed to be cross-
validation free. The popular MC+ penalty (Zhang et al.,, 2010) has two tuning
parameter $(\lambda,\gamma)$.

We applied the Sparsenet algorithm (Mazumder et al.,, 2011) to fit the MC+
penalty, which performs cross-validation over a two-dimensional grid of values
$(\lambda,\gamma)$. For NSSL, the variance is estimated by the iterative
algorithm from (Moran et al.,, 2018). For NSB and MC+ penalties, we use the
variance estimator from Equation (6).

In Tables 2 \- 5, we report the simulation results by calculating the average
of the following criteria from the 100 experiments:$L_{2}$ error(Root mean
square error),$L_{1}$ error(Mean absolute error),FDR(False discover rate),
FNDR(False non-discover rate),HD(Hamming distance),$\hat{\sigma}^{2}$(variance
estimator),$\hat{s}$(Estimated model size).We also report the standard
derivations of the $L_{2}$ error, $L_{1}$ error, Hamming distance and the
variance estimator $\hat{\sigma}^{2}$ among 100 experiments in the brackets.
We use bold to highlight the best performance for each criteria among the full
Bayesian procedure and optimization procedures respectively.

####  6.2.1 Sparse case with $s_{0}=10$

| Full Bayesian procedure | Optimization procedure
---|---|---
| $L_{\frac{1}{2}}$ | $L_{\frac{1}{2}}$ | Horseshoe | Bayesian | NSB | NSB | NSSL | MC+
| ($\gamma=1$) | ($\gamma=2$) |  | LASSO | ($\gamma=1$) | ($\gamma=3$) |  |
$n=500$, $p=25$, $s_{0}=10$, $\sigma_{0}^{2}=3$ |
L2 | 0.42(0.03) |  0.38(0.02) | 0.49(0.03) | 0.46(0.03) | 0.45(0.04) | 0.40(0.02) | 1.05(0.08) |  0.33(0.02)
L1 | 1.64(0.11) |  1.44(0.09) | 1.92(0.13) | 1.86(0.15) | 1.68(0.20) | 1.37(0.15) | 2.67(0.41) |  0.92(0.16)
FDR | 1.30 | 0.00 | 4.37 | 3.79 | 50.75 | 38.00 | 1.12 | 6.68
FNDR | 0.06 | 6.65 | 0.00 | 0.00 | 0.00 | 0.00 | 11.13 | 0.00
HD |  0.16(0.05) | 1.10(0.11) | 0.50(0.12) | 0.43(0.27) | 10.42(1.41) | 6.41(0.98) | 2.02(0.33) | 1.00(0.10)
$\hat{s}$ | 10.14 | 8.9 | 10.50 | 10.43 | 20.42 | 16.41 | 8.20 | 11.00
$\hat{\sigma}^{2}$ |  3.01(0.35) | 3.13(0.14) | 2.96(0.24) |  3.01(0.40) | 2.96(0.31) | 2.97(0.20) | 1.69(0.40) |  2.99(0.21)
$n=500$, $p=1000$, $s_{0}=10$,$\sigma_{0}^{2}=1$ |
L2 | 0.41(0.04) | – |  0.25(0.04) | 1.15(0.19) | 0.15(0.02) |  0.13(0.01) | 0.14(0.02) | 0.21(0.02)
L1 | 9.15(1.77) | – |  2.90(0.88) | 26.82(2.98) | 0.37(0.06) |  0.35(0.04) | 0.37(0.08) | 0.61(0.16)
FDR | 4.68 | – | 2.16 | 47.67 | 0.73 | 0.00 | 0.00 | 10.00
FNDR | 0.00 | – | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
HD | 0.51(0.10) | – |  0.20(0.06) | 10.4(1.21) | 0.03(0.01) |  0.00(0.00) |  0.00(0.00) | 3.06(0.54)
$\hat{s}$ | 10.51 | – | 10.20 | 20.40 | 10.03 | 10.00 | 10.00 | 13.06
$\hat{\sigma}^{2}$ | 0.70(0.05) | – |  0.81(0.04) | 0.09(0.02) | 1.02(0.10) |  1.01(0.05) | 0.97(0.05) | 1.02(0.19)
$n=500$, $p=1000$, $s_{0}=10$, $\sigma_{0}^{2}=3$ |
L2 | 0.65(0.13) | – |  0.49(0.08) | 1.45(0.31) | 0.81(0.09) | 0.41(0.07) |  0.27(0.08) | 0.34(0.09)
L1 | 13.93(2.01) | – |  4.55(1.25) | 32.17(4.96) | 3.95(0.78) | 1.22(0.35) |  0.65(0.20) | 0.97(0.28)
FDR | 10.73 | – | 9.10 | 59.16 | 63.87 | 14.54 | 0.00 | 10.77
FNDR | 0.00 | – | 0.01 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
HD | 1.29(0.20) | – |  1.01(0.17) | 15.50(2.33) | 18.97(1.52) | 2.66(0.49) |  0.14(0.02) | 2.91(0.53)
$\hat{s}$ | 11.29 | – | 11.01 | 25.50 | 28.97 | 12.66 | 9.86 | 12.91
$\hat{\sigma}^{2}$ | 2.06(0.18) | – |  2.43(0.16) | 1.20(0.15) | 2.39(0.13) | 2.88(0.10) | 2.95(0.01) |  2.97(0.14)

Table 2: Comparison of big sample size ($n=500$) and the sparse case $s_{0}=10$ among the following priors: $L_{\frac{1}{2}}$, Bayesian LASSO, Horseshoe, non-separable bridge penalty (NSB), non-separable spike-and-slab LASSO (NSSL) and MC+.  | Full Bayesian procedure | Optimization procedure
---|---|---
| $L_{\frac{1}{2}}$ | Horseshoe | Bayesian | NSB | NSB | NSSL | MC+
| ($\gamma=1)$ |  | LASSO | ($\gamma=1$) | ($\gamma=3$) |  |
$n=100$, $p=1000$, $s_{0}=10$, $\sigma_{0}^{2}=1$
L2 | 1.19(0.21) |  0.80(0.15) | 3.86(0.66) | 0.81(0.17) | 0.70(0.15) | 0.71(0.19) |  0.64(0.19)
L1 | 14.05(2.57) |  9.40(2.05) | 40.27(5.35) | 2.06(0.73) | 1.65(0.37) |  1.62(0.40) | 1.80(0.41)
FDR | 1.11 | 1.02 | 59.00 | 0.21 | 0.21 | 0.20 | 16.62
FNDR | 0.20 | 0.20 | 0.00 | 0.10 | 0.10 | 0.10 | 0.03
HD | 2.10(0.22) |  1.80 (0.19) | 12.80(2.01) | 1.77(0.15) | 1.56(0.14) |  1.50(0.14) | 4.66(0.57)
$\hat{s}$ | 8.10 | 9.20 | 14.60 | 8.03 | 7.89 | 8.54 | 13.90
$\hat{\sigma}^{2}$ | 0.02(0.01) | 0.10(0.02) |  0.25(0.05) | 1.30(0.10) | 1.44(0.11) |  1.10(0.11) |  1.10(0.14)
$n=100$, $p=1000$, $s_{0}=10$, $\sigma_{0}^{2}=3$
L2 | 2.01(0.25) |  1.60(0.20) | 4.09(0.90) | 1.46(0.22) | 1.35(0.20) | 1.41(0.19) |  1.33(0.22)
L1 | 18.96(2.88) |  13.6(2.44) | 43.33(5.55) | 4.93(1.01) | 3.91(0.41) |  3.52(0.56) | 4.46(0.59)
FDR | 6.03 | 14.37 | 73.11 | 35.99 | 20.00 | 0.80 | 36.88
FNDR | 0.36 | 0.20 | 0.44 | 0.19 | 0.20 | 0.30 | 0.16
HD | 4.04(0.44) |  3.54(0.41) | 21.30(3.54) | 6.92(1.10) | 4.46(0.77) |  3.26(0.67) | 11.63(1.69)
$\hat{s}$ | 6.78 | 9.46 | 22.6 | 13.09 | 10.00 | 6.86 | 18.39
$\hat{\sigma}^{2}$ | 0.44(0.1) |  0.77(0.21) | 0.20(0.05) | 2.49(0.30) |  3.02(0.20) | 3.68(0.41) | 2.91(0.22)
Table 3: Comparison of small sample size ($n=100$) and the sparse case
$s_{0}=10$ among the following priors: $L_{\frac{1}{2}}$, Bayesian LASSO,
Horseshoe, non-separable bridge penalty (NSB), non-separable spike-and-slab
LASSO (NSSL) and MC+.

Tables 2 and 3 summarize the results of the simulation studies for $s_{0}=10$.
In low dimensional problems ($n>p$), all approaches performed well. In high
dimensional problems ($n<p$), all except the Bayesian LASSO performed well.
This is within our expectation as Castillo et al., (2015) showed that the
Bayesian LASSO has poor posterior contraction rate in high dimensional
problems.

One remarkable thing is the under-estimation of the variance by both
$L_{\frac{1}{2}}$ and Horseshoe priors when $p>n$. It was shown by Moran et
al., (2018) that using the conjugate prior to $\sigma^{2}$ can result in the
under-estimation of variance in high dimension settings. However, as shown in
Figure 3, when $p>>n$, using the non-conjugate prior for $\sigma^{2}$ will
also lead to under-estimation. This is because there exists spurious sample
correlation with the realized noises in some predictors or between the
predictor and response, even when they are actually independent (Fan and Lv,,
2008). As a result, the realized noises are explained by the model with extra
irrelevant variables, leading to an underestimate of the variance $\sigma^{2}$
(Fan et al.,, 2012).It can be seen that the marginal posterior distribution
for $\hat{\sigma}^{2}$ has a heavy long tail and is biased downwards. Since
the full conditional posterior of $\sigma^{2}$ has a
$\mathrm{InvGamma}\left(\frac{n}{2},\frac{1}{2}(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta})^{T}(\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta})\right)$
distribution, we conclude that most of the time, the MCMC sampler of $\beta$
overfits the model. Figure 3 also shows that when the sample size becomes
large, this downwards bias becomes smaller.

Figure 3: Boxplots of the marginal posterior samples of the variance parameter
($\hat{\sigma}$), under the $L_{\frac{1}{2}}$ and Horseshoe priors. Horizontal
line indicates the true value of $\sigma^{2}$, $p>n$ cases.

Overall, among all the full Bayesian procedures, the $L_{\frac{1}{2}}$ prior
with $\gamma=2$ has the best performance in terms of $L_{2}$ error and $L_{1}$
error in the low dimensional case (Table 2). For the high dimensional case,
the Horseshoe prior slightly out performs the $L_{\frac{1}{2}}$ prior with
$\gamma=1$. The shrinkage of $L_{\frac{1}{2}}$ prior at $\gamma=1$ is slightly
weaker than the Horseshoe prior, consequently the variance estimates
$\hat{\sigma}^{2}$ from Table 2 and 3 shows that while both procedures
underestimate the variance parameter, the $L_{\frac{1}{2}}$ prior produces
more severe under-estimation. Unfortunately, numerical underflow problem
prevents us from using the $L_{\frac{1}{2}}$ prior with $\gamma=2$ in high
dimensional settings, we expect it could perform better than the Horseshoe and
$L_{\frac{1}{2}}$ with $\gamma=1$ priors, since the $L_{\frac{1}{2}}$ prior
with $\gamma=2$ has a strong shrinkage for values near zero. From Table 2, we
see that this prior tends to select less and has a slightly higher value for
the posterior mean of the variance parameter.

Amongst the optimisation procedures, in the low dimensional case, except for
the NSSL penalty (which slightly under-performs), no procedure clearly
dominated across all situations for all criteria. However, when we look at the
high dimensional problem especially for the small sample size scenario
($n=100$, Table 3), the performances of Bayesian and optimization procedures
are comparable, except for the $L_{1}$ error where optimisation procedure out-
performs the full Bayesian procedure.

All optimization procedure successfully targets the true variance and have
lower $L_{2}$ and $L_{1}$ errors. It should be emphasised again that in the
ultra high dimensional setting, the hyper-parameter $b$ in NSB penalty is
determined by a cross-validation based strategy and we believe this is the key
to success. The hyper-parameter in MC+ penalty is always determined by cross-
validation in all scenarios.

####  6.2.2 No signal case $s_{0}=0$

Table 4 shows the performance of the models under the setting with small
sample size $n=100$, $p=1000$ and no signal ($s_{0}=0$). It is apparently
that, the full Bayesian procedure out-performs the optimization procedure
under criteria which measure the performance of variable selection (i.e.
FDR,FNDR,HD and $\hat{s}$). Both the NSB penalty and MC+ penalties provide
reliable variance estimation, but this is not the case for NSSL penalty. The
NSSL selects too many variables and thus underestimates the variance.

Table 4 provides further evidence that, compared with the Horseshoe, the
$L_{\frac{1}{2}}$ prior with $\gamma=1$ provides relatively weaker shrinkage
at values near zero. We can also see that the $L_{\frac{1}{2}}$ with
$\gamma=1$ again produces smaller variance estimates, although the
underestimations are not severe in this setting. In addition, the NSB penalty
with $\gamma=1$ tends to select more variables than NSB penalty with
$\gamma=3$.

| Full Bayesian procedure | Optimization procedure
---|---|---
| $L_{\frac{1}{2}}$ | Horseshoe | Bayesian | NSB | NSB | NSSL | MC+
| ($\gamma=1)$ |  | LASSO | ($\gamma=1$) | ($\gamma=3$) |  |
$n=100$, $p=1000$, $s_{0}=0$, $\sigma_{0}^{2}=1$
L2 | 0.09(0.01) |  0.08(0.01) | 0.10(0.01) | 0.09(0.01) |  0.07(0.01) | 0.72(0.11) | 0.14(0.02)
L1 | 1.96(0.15) |  0.92(0.10) | 1.41(0.24) |  0.21(0.03) | 0.12(0.02) | 2.78(0.32) | 0.37(0.04)
FDR | 0.00 | 0.00 | 0.00 | 42.00 | 40.00 | 100.00 | 99.00
FNDR | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
HD |  0.00(0.00) |  0.00(0.00) |  0.00(0.00) | 3.18(0.18) |  1.42(0.10) | 13.41(1.17) | 5.82(0.50)
$\hat{s}$ | 0.00 | 0.00 | 0.00 | 3.18 | 1.42 | 13.41 | 5.82
$\hat{\sigma}^{2}$ | 0.75(0.07) |  0.85(0.07) | 0.75(0.07) | 0.92(0.03) |  0.95(0.02) | 0.60(0.03) | 0.90(0.03)
$n=100$, $p=1000$, $s_{0}=0$,$\sigma_{0}^{2}=3$
L2 | 0.10(0.01) | 0.13(0.02) |  0.07(0.01) |  0.13(0.01) | 0.17(0.02) | 0.95(0.15) | 0.20(0.03)
L1 |  1.11(0.12) | 1.21(0.11) | 1.94(0.28) |  0.38(0.04) |  0.38(0.04) | 3.73(0.30) | 0.66(0.05)
FDR | 0.00 | 0.00 | 0.00 | 33.33 | 41.00 | 100.00 | 98.00
FNDR | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
HD |  0.00(0.00) |  0.00(0.00) |  0.00(0.00) | 3.22(0.20) |  2.58(0.23) | 7.89(0.40) | 7.60(0.44)
$\hat{s}$ | 0.00 | 0.00 | 0.00 | 3.22 | 2.58 | 7.89 | 7.60
$\hat{\sigma}^{2}$ | 2.50(0.08) |  2.74(0.10) | 2.45(0.23) |  2.80(0.07) | 2.75(0.06) | 1.20(0.15) | 2.64(0.10)
Table 4: Comparison of high dimension and small sample size problem ($p=1000$,
$n=100$, $s_{0}=0,\sigma_{0}^{2}=1,3$) when there exists no signal, among the
following priors: $L_{\frac{1}{2}}$, Bayesian LASSO, Horseshoe, non-separable
bridge penalty (NSB), non-separable spike-and-slab LASSO (NSSL) and MC+.

####  6.2.3 Less sparse case $s_{0}=20$

Here we consider the less sparse case with $s_{0}=20$ and we set $n=100$,
$p=1000$ (small sample and high dimensions). We see that whereas all the
algorithms in the high dimensional, small sample, but very sparse setting
($s_{0}=10$, Table 3) behave reasonably well and similarly to each other,
their performances deteriorate when the sample size is small and true
underlying model become less sparse($n=100$, $p=1000$, $s_{0}=20$), see Table
5.

Comparing the result from Table 3 and 5, the full Bayesian procedure shows
deterioration of the estimation error (i.e. $L_{1}$,$L_{2}$ and
$\hat{\sigma}^{2}$). The higher Hamming distance(HD) from all the models
indicates that neither the t-statistic test from MCMC nor the thresholding
operator from optimization provides reliable variable selection. In this
difficult situation, the NSB penalty with $\gamma=3$ shows superior
performance in terms of estimation. Our result is consistent with Reid et al.,
(2016), who also investigated the similar less sparse scenario. They
recommended using the cross-validation strategy. We followed their
recommendation and proposed forward variable screening with cross-validation
strategy to determining hyper-parameter $b$.

| Full Bayesian procedure | Optimization procedure
---|---|---
| $L_{\frac{1}{2}}$ | Horseshoe | Bayesian | NSB | NSB | NSSL | MC+
| ($\gamma=1)$ |  | LASSO | ($\gamma=1$) | ($\gamma=3$) |  |
$n=100$, $p=1000$, $s_{0}=20$, $\sigma_{0}^{2}=1$
L2 | 4.03(0.50) |  2.06(0.33) | 6.17(1.03) | 1.74(0.30) | 1.33(0.27) | 2.96(0.33) |  1.32(0.25)
L1 | 39.01(4.05) |  20.14(3.24) | 68.54(6.71) | 8.26(1.60) |  5.70(0.81) | 10.83(1.70) | 6.31(0.99)
FDR | 0.00 | 0.00 | 0.00 | 47.65 | 27.00 | 0.64 | 53.00
FNDR | 1.70 | 0.70 | 1.96 | 0.26 | 0.19 | 1.00 | 0.10
HD | 16.11(2.01) |  6.96(1.05) | 19.60(2.55) | 18.56(1.88) |  8.98 (0.88) | 9.86(0.90) | 28.79(2.03)
$\hat{s}$ | 3.89 | 13.04 | 0.40 | 33.40 | 25.00 | 10.24 | 46.25
$\hat{\sigma}^{2}$ | 0.20(0.01) | 0.31(0.01) |  0.55(0.03) | 1.35(0.10) |  1.29(0.10) | 2.41(0.22) | 1.40(0.13)
$n=100$, $p=1000$, $s_{0}=20$, $\sigma_{0}^{2}=3$
L2 | 4.64(0.60) |  2.93(0.40) | 6.17(1.35) | 3.20(0.60) |  2.79(0.31) | 4.16(0.44) | 3.23(0.30)
L1 | 49.41(5.05) |  24.58(4.00) | 70.1(6.80) | 16.05(2.31) |  12.67(2.01) | 16.06(2.03) | 15.98(2.44)
FDR | 0.00 | 0.83 | 0.00 | 59.06 | 43.00 | 0.00 | 56.00
FNDR | 1.75 | 1.11 | 1.97 | 0.66 | 0.60 | 1.32 | 0.66
HD | 17.46(2.31) |  11.15(1.55) | 19.7(2.60) | 24.46(2.32) | 16.92(1.88) |  13.33(1.67) | 34.54(3.48)
$\hat{s}$ | 2.54 | 8.84 | 0.30 | 33.46 | 24.62 | 6.99 | 41.94
$\hat{\sigma}^{2}$ | 0.25(0.03) | 0.43(0.04) |  0.89(0.20) |  3.32(0.33) | 3.66(0.30) | 3.69(0.35) | 1.05(0.08)

Table 5: Comparison of high dimension and small sample size ($p=1000$,
$n=100$, $s_{0}=20,\sigma_{0}^{2}=1,3$) among the following priors:
$L_{\frac{1}{2}}$, Horseshoe, Bayesian LASSO, non-separable bridge penalty
(NSB), non-separable spike-and-slab LASSO (NSSL) and MC+.

##  7 Conclusion and Discussion

In this article, we proposed an exact scale mixture of normal representation
of the exponential power prior with $\alpha=\frac{1}{2^{\gamma}}$. Based on
this representation, we developed an efficient MCMC sampling scheme with
faster convergence rate than traditional Gibbs sampling scheme and it is
scalable to high dimensional problem. The normal mixture representation allows
the sparse prior to be easily extended to the generalized linear model and
other regression setting with only slight modifications to the MCMC scheme.

In addition, inspired by the full Bayesian approach of the exponential power
prior, we formulated a non-separable bridge penalty. This is done by
integrating out the hyper-parameter $\lambda$ in the bridge penalty with
respect to some suitable choice of hyper-prior. The NSB penalty yields a
combination of global-local adaptive shrinkage and threshold. With the
gaussian likelihood, its KKT condition lends itself to a fast optimization
algorithm via coordinate-descent. Since the NSB penalty is a non-convex, non-
separable and non-lipschitz function, we also provide a theoretical analysis
to guarantee the convergence of the CD algorithm. Both the PCG sampler and CD
optimization algorithms are not limited to the gaussian likelihood, examples
of extensions of our algorithm to more general settings can be found in
Section 5 of Supplementary Materials.

One limitation of our MCMC scheme is the numerical underflow in high
dimensional problem when $\gamma\geq 2$. For $\gamma$ sufficient large, we are
close to the $L_{0}$ penalty. We anticipate that such prior may attain sharper
than current nearly optimal contraction rate. It may also lead to better
empirical performance in high dimensional settings than the existing global-
local shrinkage priors. In fact, we have already observed the superior
performance when $\gamma=2$ in low dimensional problems. We hope that, in the
future, we can overcome the numerical underflow issue at least for the
$\gamma=2$ case. Another future research area is the uncertainty
quantification of the optimization results. One potential direction is to look
at ensemble learning with Bayesian bootstrap, which is processed by sampling
multiple weight vectors using a Dirichlet distribution and fitting multiple
penalized weight least square models. One could study whether the induced
randomness provides good uncertainty quantification.

## Supplementary Materials

Section 1. Normal-mixture representation of exponential power prior
Section 2. Derivations of Partially Collapsed Gibbs Sampler
Section 3. The KKT condition of non-separable bridge penalty
Section 4. Forward and backward variable screening
Section 5. Extension to logistic model
Section 6. Posterior consistency and contraction rate
Section 7. Convergence analysis of the coordinate descent algorithm
Section 8. Oracle property

## 1\. Normal-mixture representation of exponential power prior

Proof of Lemma 2.1

###### Proof.

If $\frac{\lambda^{2}}{4}\exp(-\lambda|\beta|^{\frac{1}{2}})$ is a Laplace
transformation of some non-negative function $f(\cdot)$ evaluated at
$x=|\beta|$, then $f(\cdot)$ is the inverse Laplace transformation

$\frac{\lambda^{2}}{4}\exp(-\lambda|\beta|^{\frac{1}{2}})=\int_{0}^{\infty}f(t)e^{-t|\beta|}dt$

where
$f(t)=\mathcal{L}^{-1}(\frac{\lambda^{2}}{4}\exp(-\lambda\sqrt{x}))(t)=\frac{\lambda^{3}}{8\sqrt{\pi}}t^{-\frac{3}{2}}\exp(-\frac{\lambda^{2}}{4t})$.

Let $\pi(s)=\frac{2}{s}f(\frac{1}{s})$, then
$s\sim{\mathrm{Gamma}(\frac{3}{2},\frac{\lambda^{2}}{4})}$. Following this
rearrangement, we can put a laplace prior on $\beta$ conditionally on $s$,

$\pi(\beta|s)=\frac{1}{2s}e^{-\frac{|\beta|}{s}}.$

We can construct a laplace-gamma mixture representation for exponential power
prior with $\alpha=\frac{1}{2}$

$\frac{\lambda^{2}}{4}\exp(-\lambda|\beta|^{\frac{1}{2}})=\int_{0}^{\infty}\pi(\beta|s)\pi(s)ds.$

To get another parametric form, move the hyper-parameter $\lambda$ outside the
gamma distribution, and apply the change of variable in the right hand side of
the above equation. Let $v=\lambda^{2}s$, then we have

$\frac{\lambda^{2}}{4}\exp(-\lambda|\beta|^{\frac{1}{2}})=\int_{0}^{\infty}\pi(\beta|v)\pi(v)dv$

where $\pi(\beta|v)=\frac{\lambda^{2}}{2v}e^{-\frac{\lambda^{2}|\beta|}{v}}$
and $v\sim\mathrm{Gamma}(\frac{3}{2},\frac{1}{4})$. ∎

Proof of Lemma 2.2

###### Proof.

Write

$\frac{\lambda^{2^{\gamma}}}{2(2^{\gamma}!)}\exp\left(-\lambda\left|\beta\right|^{\frac{1}{2^{\gamma}}}\right)=\int_{0}^{\infty}f(t)e^{-t|\beta|^{\frac{1}{2^{\gamma-1}}}}dt$

then
$\frac{\lambda^{2^{\gamma}}}{2(2^{\gamma}!)}\exp\left(-\lambda\left|\beta\right|^{\frac{1}{2^{\gamma}}}\right)$
is the Laplace transform of of some non-negative function $f(\cdot)$ evaluated
at $x=|\beta|^{\frac{1}{2^{\gamma-1}}}$, where
$f(t)=\frac{1}{\Gamma(\frac{2^{\gamma}+1}{2})}\frac{1}{2^{2^{\gamma}+2}}\frac{\lambda^{2^{\gamma}+1}}{(2^{\gamma-1}!)}t^{-\frac{3}{2}}\exp(-\frac{\lambda^{2}}{4t})$.

Let $\pi(s)=\frac{2(2^{\gamma-1}!)s^{2^{\gamma}-1}}{s}f(\frac{1}{s})$, then
$s\sim{\mathrm{Gamma}(\frac{2^{\gamma}+1}{2},\frac{\lambda^{2}}{4})}$.
Following this rearrangement, we can put a exponential power prior with
$\alpha=\frac{1}{2^{\gamma-1}}$ on $\beta$ conditional on $s_{\gamma}$ then,

$\pi(\beta|s_{\gamma})=\frac{1}{2(2^{\gamma-1}!)s_{\gamma}^{2^{\gamma-1}}}\exp(-\frac{|\beta|^{\frac{1}{2^{\gamma-1}}}}{s_{\gamma}}).$

We can see that the exponential power prior with $\alpha=\frac{1}{2^{\gamma}}$
is the mixture of exponential power prior with $\alpha=\frac{1}{2^{\gamma-1}}$
and $\mathrm{Gamma}(\frac{2^{\gamma}+1}{2},\frac{\lambda^{2}}{4})$

$\frac{\lambda^{2^{\gamma}}}{2(2^{\gamma}!)}\exp\left(-\lambda\left|\beta\right|^{\frac{1}{2^{\gamma}}}\right)=\int_{0}^{\infty}\pi(\beta|s_{\gamma})\pi(s_{\gamma})ds_{\gamma}.$

Using a similar argument as in the proof of Lemma 2.1 moving the hyper-
parameter $\lambda$ outside the gamma distribution by using the change of
variable trick. We let $v_{\gamma}=\lambda^{2}s_{\gamma}$, then

$\pi(\beta|v_{\gamma},\lambda)=\frac{\lambda^{2^{\gamma}}}{2(2^{\gamma-1}!)v_{\gamma}^{2^{\gamma-1}}}\exp(-\frac{\lambda^{2}|\beta|^{\frac{1}{2^{\gamma-1}}}}{v_{\gamma}}),\quad
v_{\gamma}\sim\mathrm{Gamma}(\frac{2^{\gamma}+1}{2},\frac{1}{4}).$

∎

Proof of Theorem 2.3

###### Proof.

Starting with the decomposition given in Lemma 2.2, we make some further
decompositions. Our argument consists of three steps:

Step 1: Set $\lambda^{\prime}=\frac{\lambda^{2}}{v_{\gamma}}$, then

$\pi(\beta|v_{\gamma},\lambda)=\frac{\lambda^{2^{\gamma}}}{2(2^{\gamma-1}!)v_{\gamma}^{2^{\gamma-1}}}\exp(-\frac{\lambda^{2}|\beta|^{\frac{1}{2^{\gamma-1}}}}{v_{\gamma}})=\frac{\lambda^{\prime
2^{\gamma-1}}}{2(2^{\gamma-1}!)}\exp(-\lambda^{\prime}|\beta|^{\frac{1}{2^{\gamma-1}}})$

Step 2: Apply Lemma 2.2 again with $\frac{\lambda^{\prime
2^{\gamma-1}}}{2(2^{\gamma-1}!)}\exp\left(-\lambda\left|\beta\right|^{\frac{1}{2^{\gamma-1}}}\right)=\int_{0}^{\infty}\pi(\beta|s_{\gamma})\pi(s_{\gamma})ds_{\gamma}$,
then

$\pi(\beta|s_{\gamma-1})=\frac{1}{2(2^{\gamma-2}!)s_{\gamma-1}^{2^{\gamma-2}}}\exp(-\frac{|\beta|^{\frac{1}{2^{\gamma-2}}}}{s_{\gamma-1}}),\quad
s_{\gamma-1}\sim\mathrm{Gamma}(\frac{2^{\gamma-1}+1}{2},\frac{\lambda^{\prime
2}}{4})$

Step 3: By change of variables, set $v_{\gamma-1}=\lambda^{4}s_{\gamma-1}$,
then

$\pi(\beta|v_{\gamma-1},\lambda)=\frac{\lambda^{2^{\gamma}}}{2(2^{\gamma-2}!)v_{\gamma-1}^{2^{\gamma-2}}}\exp(-\frac{\lambda^{4}|\beta|^{\frac{1}{2^{\gamma-2}}}}{v_{\gamma-1}}),\quad
v_{\gamma-1}\sim\mathrm{Gamma}(\frac{2^{\gamma-1}+1}{2},\frac{1}{4v_{\gamma}^{2}}).$

Using the above three step argument recursively, then for any
$i=1,2,...,\gamma-1$, we have

$\pi(\beta|v_{\gamma-i},\lambda)=\frac{\lambda^{2^{\gamma}}}{2(2^{\gamma-i-1}!)v_{\gamma-i}^{2^{\gamma-i-1}}}\exp\left(-\frac{\lambda^{2^{i+1}}|\beta|^{\frac{1}{2^{\gamma-i-1}}}}{v_{\gamma-i}}\right),\quad
v_{\gamma-i}\sim\mathrm{Gamma}\left(\frac{2^{\gamma-i}+1}{2},\frac{1}{4v_{\gamma-i+1}^{2}}\right).$

Consequently, the exponential power prior with $\alpha=\frac{1}{2^{\gamma}}$
can be written as

$\beta|v_{1},\lambda\sim\mathrm{DE}(\frac{v_{1}}{\lambda^{2^{\gamma}}}),\quad
v_{i}|v_{i+1}\sim\mathrm{Gamma}(\frac{2^{i}+1}{2},\frac{1}{4v_{i+1}^{2}}),\quad
v_{\gamma}\sim\mathrm{Gamma}(\frac{2^{\gamma}+1}{2},\frac{1}{4})$

for $i=1,\ldots,\gamma-1$. Since it is well know that the laplace distribution
can be represented as a normal exponential mixture. ∎

### Some intuitions for the $L_{\frac{1}{2}}$ prior

We now provide some intuition as to why the $L_{\frac{1}{2}}$ prior is
suitable as a sparsity prior. We carry out the illustrations in terms of the
simpler normal mean model. Our treatment is the same as Carvalho et al.,
(2010). Suppose that we observe data from the probability model
$y_{i}|\beta_{i}\sim N(\beta_{i},1)$ for $i=1,2,...,n$. Here the primary goal
is to estimate the vector of normal means
$\boldsymbol{\beta}=(\beta_{1},\beta_{2},...,\beta_{n})$ and a secondary goal
is to simultaneously test if the $\beta_{i}$’s are zero. To do that, we
consider the Bayesian method by applying the sparsity inducing prior to
$\boldsymbol{\beta}$. We focus on the posterior mean, which is known to be
optimal under quadratic loss. The global-local $L_{\frac{1}{2}}$ shrinkage
prior can be represented as $\beta_{i}\mid\tau_{i}^{2},\lambda\sim
N\left(0,\tau_{i}^{2}\lambda\right)$ and
$\tau_{i}^{2}\sim\pi\left(\tau_{i}^{2}\right).$

For simplicity, we further assume $\lambda=1$. Then the Bayes estimator for
$\beta_{i}$ under quadratic loss is

$\hat{\beta}_{i}\left(\tau_{i}\right)=\mathrm{E}\left(\beta_{i}\mid\mathbf{y},\tau_{i}\right)=\left(1-\frac{1}{1+\tau_{i}^{2}}\right)y_{i},$

hence by Fubini’s theorem, the marginal posterior mean is

$\hat{\beta}_{i}=\mathrm{E}\left(\beta_{i}\mid\mathbf{y}\right)=\int_{0}^{\infty}\left(1-\frac{1}{1+\tau_{i}^{2}}\right)y_{i}\pi\left(\tau_{i}\mid\mathbf{y}\right)\mathrm{d}\tau_{i}=\left[1-\mathrm{E}\left(\frac{1}{1+\tau_{i}^{2}}\mid\mathbf{y}\right)\right]\cdot
y_{i}.$

Reparametrising with $\kappa_{i}=\frac{1}{1+\tau_{i}^{2}}$, then
$\kappa\in[0,1]$ can be viewed as a shrinkage factor, where $\kappa=0$ means
no shrinkage and $\kappa=1$ means total shrinkage to zero.

For the exponential power prior with $\alpha=\frac{1}{2}$, the inducing
density of $\kappa$ is given by

$\pi_{\alpha=\frac{1}{2}}(\kappa)=\frac{1}{8\sqrt{\pi}\kappa^{2}}\int_{0}^{\infty}v^{-1.5}\exp\left(-\frac{1}{4}v+\frac{1}{2v^{2}}-\frac{1}{2v^{2}\kappa}\right)dv.$

Figure 4 demonstrates the behaviour of the $L_{\frac{1}{2}}$ prior with
$\gamma=1$ (left panel), as well as the Horseshoe prior (middle panel) and the
Laplace (Bayesian LASSO) prior (right panel).

The $L_{\frac{1}{2}}$ prior shares much similarity to the Horseshoe prior,
both mimic the behavior of the spike-and-slab prior, putting large probability
mass around $\kappa=0$ and $\kappa=1$. The $L_{\frac{1}{2}}$ prior starts from
zero, then dramatically increase to a first peak at $\kappa_{0}$ near zero and
finally follows a U-shape that is very similar to the Horseshoe prior. In
addition, its density function is almost zero for $\kappa\in[0.4,1)$.

For the more general $\gamma>1$ case, there is no analytic expression for the
density of $\kappa$, and it is difficult to numerically illustrate the
density. Intuitively as $\gamma$ becomes larger, $\alpha$ tends to zero, the
density of $\kappa$ will degenerate to two points $\kappa=0$ and $\kappa=1$,
because the $L_{\frac{1}{2}}$ prior converges to an improper spike-and-slab
prior
$\pi(\beta)\propto\exp(-\|\beta\|_{0})\propto\frac{1}{1+e^{-1}}\delta_{0}(\beta)+\frac{e^{-1}}{1+e^{-1}}\delta_{|\beta|>0}(\beta)$,
induced by the $L_{0}$ norm, and $\delta$ denotes the dirac mass function. In
other words, as $\alpha$ tends to zero the density of $\pi_{\alpha}(\kappa)$
will degenerate to two point mass with $P(\kappa=0)=\frac{e^{-1}}{1+e^{-1}}$
and $P(\kappa=1)=\frac{1}{1+e^{-1}}$. Under the sparse normal mean model, we
have

$\pi(\beta_{i}|y_{i})\propto\frac{1}{1+e^{-1}}\delta_{0}(\beta_{i})+\frac{e^{-1}}{1+e^{-1}}N(y_{i},1).$

and the corresponding posterior mean for $\beta_{i}$ is given by

$\hat{\beta}_{i}=\mathrm{E}\left(\beta_{i}\mid\mathbf{y_{i}}\right)=[1-\mathrm{E}\left(\kappa_{i}\mid\mathbf{y_{i}}\right)]\cdot
y_{i}=\frac{e^{-1}}{1+e^{-1}}y_{i}$

with
$\mathrm{E}\left(\kappa_{i}\mid\mathbf{y_{i}}\right)=P(\kappa=1)=\frac{1}{1+e^{-1}}$.

Figure 4: The density of $\kappa$ for $L_{\frac{1}{2}}$ prior with $\gamma=1$
(left), the Horseshoe prior (middle) and Laplace prior (right). All global
shrinkage parameters are set to 1.

## 2\. Partially collapsed Gibbs sampling

### Steps for PCG sampler

Here we give the details of the three steps that comprise the PCG sampler, for
sampling the parameters
$\boldsymbol{\beta},\sigma^{2},\tau_{1}^{2},\ldots,\tau^{2}_{p},v_{1},\ldots,v_{p},\lambda$
under the $L_{\frac{1}{2}}$ prior with $\gamma=1$.
Step M: Marginalise

$\displaystyle\text{Step 1:
Draw}\;\boldsymbol{\beta}\;\text{from}\;\pi(\boldsymbol{\beta}|\sigma^{2},\boldsymbol{\tau}^{2},\boldsymbol{v},\lambda,b,\boldsymbol{Y})$
$\displaystyle\text{Step 2:
Draw}\;{\boldsymbol{\tau}^{2}}\;\text{from}\;\pi(\boldsymbol{\tau}^{2}|\boldsymbol{\beta},\boldsymbol{v},\lambda,\sigma^{2},b,\boldsymbol{Y})$
$\displaystyle\text{Step 3:
Draw}\;\sigma^{2}\;\text{from}\;\pi(\sigma^{2}|\boldsymbol{\beta},\boldsymbol{\tau}^{2},v,\lambda,b,\boldsymbol{Y})$
$\displaystyle\text{Step 4:
Draw}\;b\;\text{from}\;\pi(b|\boldsymbol{\beta},\boldsymbol{\tau}^{2},v,\lambda,\sigma^{2},\boldsymbol{Y})$
$\displaystyle\text{Step 5:
Draw}\;{\boldsymbol{\tau}^{2}}^{\star},\;\boldsymbol{v}^{\star}\;\text{and}\;\lambda\;\text{from}\;\pi(\boldsymbol{\tau}^{2},v,\lambda|\boldsymbol{\beta},\sigma^{2},b,\boldsymbol{Y})$
$\displaystyle\text{Step 6:
Draw}\;\boldsymbol{\tau}^{2\star}\;\text{and}\;\boldsymbol{v}\;\text{from}\;\pi(\boldsymbol{\tau}^{2},v|\boldsymbol{\beta},\sigma^{2},\lambda,b,\boldsymbol{Y})$

Step P: Permute

$\displaystyle\text{Step 1:
Draw}\;\boldsymbol{\beta}\;\text{from}\;\pi(\boldsymbol{\beta}|\sigma^{2},\boldsymbol{\tau}^{2},v,\lambda,b,\boldsymbol{Y})$
$\displaystyle\text{Step 2:
Draw}\;{\boldsymbol{\tau}^{2}}^{\star},\;\boldsymbol{v}^{\star}\;\text{and}\;\lambda\;\text{from}\;\pi(\boldsymbol{\tau}^{2},v,\lambda|\boldsymbol{\beta},\sigma^{2},b,\boldsymbol{Y})$
$\displaystyle\text{Step 3:
Draw}\;{\boldsymbol{\tau}^{2}}^{\star}\;\text{and}\;\boldsymbol{v}\;\text{from}\;\pi(\boldsymbol{\tau}^{2},\boldsymbol{v}|\boldsymbol{\beta},\sigma^{2},\lambda,b,\boldsymbol{Y})\;$
$\displaystyle\text{Step 4:
Draw}\;\boldsymbol{\tau}^{2}\;\text{from}\;\pi(\boldsymbol{\tau}^{2}|\boldsymbol{\beta},\boldsymbol{v},\lambda,\sigma^{2},b,\boldsymbol{Y})$
$\displaystyle\text{Step 5:
Draw}\;\sigma^{2}\;\text{from}\;\pi(\sigma^{2}|\boldsymbol{\beta},\boldsymbol{\tau}^{2},\boldsymbol{v},\lambda,b,\boldsymbol{Y})$
$\displaystyle\text{Step 6:
Draw}\;b\;\text{from}\;\pi(b|\boldsymbol{\beta},\boldsymbol{\tau}^{2},v,\lambda,\sigma^{2},\boldsymbol{Y})$

Step T: Trim

$\displaystyle\text{Step 1:
Draw}\;\boldsymbol{\beta}\;\text{from}\;\pi(\boldsymbol{\beta}|\sigma^{2},\boldsymbol{\tau}^{2},v,\lambda,\boldsymbol{Y})$
$\displaystyle\text{Step 2:
Draw}\;\lambda\;\text{from}\;\pi(\lambda|\boldsymbol{\beta},\sigma^{2},\boldsymbol{Y})$
$\displaystyle\text{Step 3:
Draw}\;\boldsymbol{v}\;\text{from}\;\pi(\boldsymbol{v}|\boldsymbol{\beta},\sigma^{2},\lambda,\boldsymbol{Y})\;$
$\displaystyle\text{Step 4:
Draw}\;\boldsymbol{\tau}^{2}\;\text{from}\;\pi(\boldsymbol{\tau}^{2}|\boldsymbol{\beta},\boldsymbol{v},\lambda,\sigma^{2},\boldsymbol{Y})$
$\displaystyle\text{Step 5:
Draw}\;\sigma^{2}\;\text{from}\;\pi(\sigma^{2}|\boldsymbol{\beta},\boldsymbol{\tau}^{2},\boldsymbol{v},\lambda,\boldsymbol{Y})$
$\displaystyle\text{Step 6:
Draw}\;b\;\text{from}\;\pi(b|\boldsymbol{\beta},\boldsymbol{\tau}^{2},v,\lambda,\sigma^{2},\boldsymbol{Y})$

We use superscript $\star$ to present intermediate quantities that are sampled
but not retained as part of the output. Step M is generalisation of the
traditional Gibbs sampler with some components being updated multiple times
within each iteration.

In the second step we rearranged the ordering of the update in Step M by
swapping the ordering of steps 5 and 6 with steps 2,3 and 4. This step does
not alter the stationary distribution of the Gibbs sampler. Finally, the
intermediate draws of $\tau^{2\star}$ and $v^{\star}$ are not necessary if we
can sample $\lambda$ and $v$ in steps 2 and 3 directly from the respective
marginal distributions, so trimming the intermediate Steps 2 and 3, we arrive
at the last set of updating procedure.

### Derivations of conditional posteriors

Since

$\pi(\lambda|\boldsymbol{\beta},\boldsymbol{\tau}^{2},\sigma^{2},b,\boldsymbol{Y})\propto\lambda^{2^{\gamma}p}\exp\left(-\lambda\sum_{j=1}^{p}\left|\beta_{j}\right|^{\frac{1}{2^{\gamma}}}\right)\pi(\lambda|b)$

with $\lambda|b\sim\mathrm{Gamma}(\frac{1}{2},\frac{1}{b})$. The conditional
posterior is a gamma distribution with shape parameters $0.5+2^{\gamma}p$ and
$\frac{1}{b}+\sum_{j=1}^{p}\left|\beta_{j}\right|^{\frac{1}{2^{\gamma}}}$.

Next, we derive the conditional posterior for $\boldsymbol{v}_{i}$. Observe
that given $\boldsymbol{v}_{i+1}$, $\boldsymbol{v}_{i}$ is conditionally
independent of the remaining $\boldsymbol{v}_{j}$s. In addition, given
$\boldsymbol{\beta}$, $\boldsymbol{v}_{i}$ is also conditionally independent
of $\sigma^{2}$ and b. Thus,

$\pi(v_{i}|\boldsymbol{\beta},\lambda,v_{i+1},\boldsymbol{Y})=\prod_{j=1}^{p}\pi(v_{ij}|\beta_{j},\lambda,v_{i+1,j},\boldsymbol{Y})\propto\prod_{j=1}^{p}\pi(\beta_{j}|v_{ij},\lambda)\pi(v_{ij}|v_{i+1,j}).$

From Theorem 2.3, we see that

$\pi(\beta_{j}|v_{ij},\lambda)\propto\frac{\lambda^{2^{i}}}{2(2^{i-1}!)v_{ij}^{2^{i-1}}}\exp(-\frac{\lambda^{2}|\beta_{j}|^{\frac{1}{2^{i-1}}}}{v_{ij}}),\quad
v_{ij}|v_{i+1,j}\sim\mathrm{Gamma}(\frac{2^{i}+1}{2},\frac{1}{4v_{i+1,j}^{2}}).$

Thus,

$\pi(\beta_{j}|v_{ij},\lambda)\pi(v_{ij}|v_{i+1,j})\propto
v_{ij}^{-\frac{1}{2}}\exp(-\frac{v_{ij}}{4v_{i+1,j}^{2}}-\frac{\lambda^{2}|\beta_{j}|^{\frac{1}{2^{i-1}}}}{v_{ij}})$
$v_{ij}|\beta_{j},\lambda,v_{i+1,j},\boldsymbol{Y}\sim\mathrm{GIG}(\frac{1}{2v_{i+1,j}^{2}},2\lambda^{2}|\beta_{j}|^{\frac{1}{2^{i-1}}},\frac{1}{2}),$

since sampling the generalized inverse gaussian distribution is not easy, we
can sample $\frac{1}{v_{ij}}$ from the inverse gaussian distribution.

$h_{ij}=\frac{1}{v_{ij}}|\beta_{j},\lambda,v_{i+1,j},\boldsymbol{Y}\sim\mathrm{InvGaussian}(\frac{1}{2\lambda
v_{i+1,j}|\beta_{j}|^{\frac{1}{2^{i}}}},\frac{1}{2v_{i+1,j}^{2}}).$

## 3\. The KKT condition

Proof of Theorem 5.1

###### Proof.

When $\beta_{j}\neq 0$, taking $\frac{\partial L}{\partial\beta_{j}}=0$ yields

$|\beta_{j}|=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|-(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{|\beta_{j}|+C_{2}|\beta_{j}|^{1-\frac{1}{2^{\gamma}}}}$
(15)

and $\mathrm{sign}(z_{j})=\mathrm{sign}(\beta_{j})$. Thus, without loss of
generality, we assume $\beta_{j}\geq 0$.

Step 1: We first show that there exists a hard threshold in the KKT condition.
When $\beta_{j}\neq 0$ is the solution, define

$\psi(\beta_{j})=\beta_{j}+(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\beta_{j}+C_{2}\beta_{j}^{1-\frac{1}{2^{\gamma}}}}.$

Then

$\displaystyle\psi^{{}^{\prime}}(\beta_{j})=$ $\displaystyle
1-(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}C_{1}\frac{1+(1-\frac{1}{2^{\gamma}})C_{2}\beta_{j}^{-\frac{1}{2^{\gamma}}}}{\beta_{j}^{2}\left(1+C_{2}\beta_{j}^{-\frac{1}{2^{\gamma}}}\right)^{2}}$
(16) $\displaystyle\psi^{{}^{\prime\prime}}(\beta_{j})=$
$\displaystyle(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}C_{1}\frac{C_{2}(1-\frac{1}{2^{\gamma}})\frac{1}{2^{\gamma}}(\beta_{j}^{-\frac{1}{2^{\gamma}}}+C_{2}\beta_{j}^{-\frac{1}{2^{\gamma-1}}})+2[1+(1-\frac{1}{2^{\gamma}})C_{2}\beta_{j}^{-\frac{1}{2^{\gamma}}}]^{2}}{(\beta_{j}+C_{2}\beta_{j}^{1-\frac{1}{2^{\gamma}}})^{2}}>0,$

therefore, $\psi(\beta_{j})$ is minimised at $\tilde{\beta}_{j}$ where
$\psi^{{}^{\prime}}(\tilde{\beta}_{j})=0$. i.e., $\tilde{\beta}_{j}$ satisfy
the following equation

$1=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}C_{1}\frac{1+(1-\frac{1}{2^{\gamma}})C_{2}\tilde{\beta}_{j}^{-\frac{1}{2^{\gamma}}}}{\tilde{\beta}_{j}^{2}\left(1+C_{2}\tilde{\beta}_{j}^{-\frac{1}{2^{\gamma}}}\right)^{2}}$
(17)

with minimum value

$h(\tilde{\beta}_{j})=\tilde{\beta}_{j}+(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\tilde{\beta}_{j}+C_{2}\tilde{\beta}_{j}^{(1-\frac{1}{2^{\gamma}})}}.$

Re-writing Equation (15) as

$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|=\beta_{j}+(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\beta_{j}+C_{2}\beta_{j}^{1-\frac{1}{2^{\gamma}}}},$

we see that
$\psi(\beta_{j})=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|$ has
no solution when
$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|<h(\tilde{\beta}_{j})$,
in this case, $\beta_{j}=0$. For
$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|\geq
h(\tilde{\beta}_{j})$, there exists at least one solution. As
$\psi(\beta_{j})$ is strictly increasing when
$\beta_{j}\geq\tilde{\beta}_{j}$, there exists
$\beta_{j}^{{}^{\prime\prime}}\in[\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]$
such that
$\psi(\beta_{j}^{{}^{\prime\prime}})=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|$.
Step 2: Next, we show the convergence of fixed point iteration at
$\beta_{j}^{{}^{\prime\prime}}\in[\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]$.
Let

$\rho(\beta_{j})=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|-(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\beta_{j}+C_{2}\beta_{j}^{1-\frac{1}{2^{\gamma}}}},$

we need to check that

$\rho([\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|])\in[\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|].$

$\rho([\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|])\leq(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|$
is follows from the definition of $\rho(\cdot)$. We see that for
$\beta_{j}\in[\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]$

$h(\tilde{\beta}_{j})=\tilde{\beta}_{j}+(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\tilde{\beta}_{j}+C_{2}\tilde{\beta}_{j}^{(1-\frac{1}{2^{\gamma}})}}\geq\tilde{\beta}_{j}+(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\beta_{j}+C_{2}\beta_{j}^{(1-\frac{1}{2^{\gamma}})}}.$

Using the fact that
$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|>h(\tilde{\beta}_{j})$,
we have

$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|\geq\tilde{\beta}_{j}+(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\beta_{j}+C_{2}\beta_{j}^{(1-\frac{1}{2^{\gamma}})}}.$

Thus,
$\rho([\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|])\geq\tilde{\beta}_{j}$.
Now we show that $\rho(\cdot)$ is a a contraction mapping on
$[\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]$.
We have
$\rho^{{}^{\prime}}(\beta_{j})=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}C_{1}\frac{1+(1-\frac{1}{2^{\gamma}})C_{2}\beta_{j}^{-\frac{1}{2^{\gamma}}}}{\beta_{j}^{2}\left(1+C_{2}\beta_{j}^{-\frac{1}{2^{\gamma}}}\right)^{2}}>0$
and
$\rho^{{}^{\prime\prime}}(\beta_{j})=-\psi^{{}^{\prime\prime}}(\beta_{j})<0$.
Since $\rho^{{}^{\prime}}(\tilde{\beta}_{j})=1$,
$\rho^{{}^{\prime}}(\beta_{j})<1$ for $\beta_{j}>\tilde{\beta}_{j}$.
Therefore, $\rho(\cdot)$ is a a contraction mapping on
$[\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]$.
Since
$[\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]$
is complete, by Banach fixed point theorem, there exists one and only one
fixed point
$\beta_{j}^{{}^{\prime\prime}}\in[\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]$
such that

$\beta_{j}^{{}^{\prime\prime}}=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|-(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\beta_{j}^{{}^{\prime\prime}}+C_{2}\beta_{j}^{\prime\prime
1-\frac{1}{2^{\gamma}}}}$ (18)

Step 3: Up to now, we showed that when
$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|>h(\tilde{\beta}_{j})$,
there exists a
$\beta_{j}^{{}^{\prime\prime}}\in[\tilde{\beta}_{j},(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]$
such that the KKT condition is satisfied. We haven’t precluded the potential
solutions $\beta_{j}^{{}^{\prime}}\in(0,\tilde{\beta})$, which also satisfies
the KKT condition. Now we want to claim that $\beta_{j}^{{}^{\prime\prime}}$
is the only local minimum for
$\beta\in(0,(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]$ given the
parameters in other dimension are fixed. We write the descent function

$\displaystyle\delta_{-j}(|\beta_{j}|)=L_{-j}(\beta_{j})-L_{-j}(0)$
$\displaystyle=\frac{1}{2}\|Y_{-j}-X_{j}\beta_{j}\|^{2}-\frac{1}{2}\|Y_{-j}\|^{2}+(2^{\gamma}p+0.5)\log\left(1+\frac{|\beta_{j}|^{\frac{1}{2^{\gamma}}}}{C_{2}}\right)$
(19)
$\displaystyle=\frac{(X_{j}^{T}X_{j})}{2}\beta_{j}^{2}-|z_{j}||\beta_{j}|+(2^{\gamma}p+0.5)\log\left(1+\frac{|\beta_{j}|^{\frac{1}{2^{\gamma}}}}{C_{2}}\right)$

where $Y_{-j}=Y-X_{-j}\beta_{-j}$. Then it is easy to see that

$\delta_{-j}^{{}^{\prime\prime}}(|\beta_{j}|)=(X_{j}^{T}X_{j})\psi^{{}^{\prime}}(|\beta_{j}|),\quad\delta_{-j}^{{}^{\prime\prime\prime}}(|\beta_{j}|)=(X_{j}^{T}X_{j})\psi^{{}^{\prime\prime}}(|\beta_{j}|)>0$

where the derivative is taken with respect to $|\beta_{j}|$. Since
$\delta_{-j}^{{}^{\prime\prime}}(\tilde{\beta}_{j})=0$ and
$\delta_{-j}^{{}^{\prime\prime\prime}}(|\beta_{j}|)>0$, we see that
$\delta_{-j}^{{}^{\prime\prime}}(|\beta_{j}|)<0$ for
$|\beta_{j}|\in(0,\tilde{\beta}_{j})$ and
$\delta_{-j}^{{}^{\prime\prime}}(|\beta_{j}|)>0$ for
$|\beta_{j}|\in(\tilde{\beta},(X_{j}^{T}X_{j})^{-1}|z_{j}|)$. Thus, if there
exists $\beta_{j}^{{}^{\prime}}\in(0,\tilde{\beta})$ such that
$\delta_{-j}^{{}^{\prime}}(\beta_{j}^{{}^{\prime}})=0$, then
$\beta_{j}^{{}^{\prime}}$ can only be the local maximum.
$\beta_{j}^{{}^{\prime\prime}}$ is the only local minimum.
Step 4. To find the global minimum when
$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|>h(\tilde{\beta}_{j})$,
we need to check the value of descent function
$\delta_{-j}(\beta_{j}^{{}^{\prime\prime}})$ because we couldn’t preclude the
case that $0$ is the global minimum. To see the existence of this potential
case, we provide a constructive proof.
Suppose
$(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|>h(\tilde{\beta}_{j})$,
plug in (18) into (19), then we have

$\delta_{-j}(\beta_{j}^{{}^{\prime\prime}})=-\frac{(X_{j}^{T}X_{j})}{2}{\beta_{j}^{{}^{\prime\prime}}}^{2}-\frac{C_{1}}{1+C_{2}\beta_{j}^{\prime\prime-\frac{1}{2^{\gamma}}}}+(2^{\gamma}p+0.5)\log\left(1+\frac{{\beta_{j}^{{}^{\prime\prime}}}^{\frac{1}{2^{\gamma}}}}{C_{2}}\right)$

we have $\delta_{-j}(\beta_{j}^{{}^{\prime\prime}})>0$ if

$(2^{\gamma}p+0.5)\log\left(1+\frac{{\beta_{j}^{{}^{\prime\prime}}}^{\frac{1}{2^{\gamma}}}}{C_{2}}\right)>\frac{(X_{j}^{T}X_{j})}{2}{\beta_{j}^{{}^{\prime\prime}}}^{2}+\frac{C_{1}}{1+C_{2}\beta_{j}^{\prime\prime-\frac{1}{2^{\gamma}}}}.$

Thus, there exists sufficient small $m$ such that for
$\beta^{{}^{\prime\prime}}<m$, the inequality above will hold. We rewrite
Equation (18) as

$f(\beta_{j}^{{}^{\prime\prime}},|z_{j}|)=\beta_{j}^{{}^{\prime\prime}}-(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|+(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\beta_{j}^{{}^{\prime\prime}}+C_{2}\beta_{j}^{\prime\prime
1-\frac{1}{2^{\gamma}}}}=0.$

Since $\frac{\partial f}{\partial\beta_{j}^{{}^{\prime\prime}}}>0$ for
$\beta_{j}^{{}^{\prime\prime}}\in(\tilde{\beta},(X_{j}^{T}X_{j})^{-1}|z_{j}|)$,
by implicit function theorem, there exists $g\in\mathrm{C}^{1}$ such that
$g(|z_{j}|)=\beta_{j}^{{}^{\prime\prime}}$. In addition, the fixed point
iteration guarantees that $g$ is a one to one map for
$\beta\in(\tilde{\beta},(X_{j}^{T}X_{j})^{-1}|z_{j}|)$. By Equation (18), it
is easy to see that $g$ is a monotonic increasing function. Hence, for
$|z_{j}|<g^{-1}(m)$ the inequality will hold. By Lemma 5.2 in the main paper,
we have $2\tilde{\beta}_{j}\leq h(\tilde{\beta}_{j})\leq 3\tilde{\beta}_{j}$
and

$\tilde{\beta}_{j}^{(2-\frac{1}{2^{\gamma}})}=\frac{C_{1}M^{\prime}(\gamma)}{\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j}}\frac{1}{\tilde{\beta}_{j}^{\frac{1}{2^{\gamma}}}+C_{2}}$

where $M^{\prime}(\gamma)\in(\frac{1}{2},1]$. By setting $b$ sufficient small,
it is possible to construct $h(\tilde{\beta}_{j})$ such that

$h(\tilde{\beta}_{j})<(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|<(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}g^{-1}(m).$

∎

Proof of Lemma 5.2

###### Proof.

From Equation (10) in the main paper, we plug
$\frac{(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}C_{1}}{\tilde{\beta}_{j}\left(1+C_{2}\tilde{\beta}_{j}^{-\frac{1}{2^{\gamma}}}\right)}=\frac{\tilde{\beta}_{j}\left(1+C_{2}\tilde{\beta}_{j}^{-\frac{1}{2^{\gamma}}}\right)}{1+(1-\frac{1}{2^{\gamma}})C_{2}\tilde{\beta}_{j}^{-\frac{1}{2^{\gamma}}}}$
into $h(\tilde{\beta}_{j})$, then

$\displaystyle h(\tilde{\beta}_{j})$
$\displaystyle=\tilde{\beta}_{j}+(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}\frac{C_{1}}{\tilde{\beta}_{j}+C_{2}\tilde{\beta}_{j}^{(1-\frac{1}{2^{\gamma}})}}=\tilde{\beta}_{j}+\frac{\tilde{\beta}_{j}\left(1+C_{2}\tilde{\beta}_{j}^{-\frac{1}{2^{\gamma}}}\right)}{1+(1-\frac{1}{2^{\gamma}})C_{2}\tilde{\beta}_{j}^{-\frac{1}{2^{\gamma}}}}$
$\displaystyle=2\tilde{\beta}_{j}+\frac{\frac{1}{2^{\gamma}}C_{2}}{\tilde{\beta}_{j}^{\frac{1}{2^{\gamma}}}+(1-\frac{1}{2^{\gamma}})C_{2}}\tilde{\beta}_{j}\leq
2\tilde{\beta}_{j}+\frac{\frac{1}{2^{\gamma}}}{(1-\frac{1}{2^{\gamma}})}\tilde{\beta}_{j}\leq
3\tilde{\beta}_{j}.$

By arranging
$1=(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}C_{1}\frac{1+(1-\frac{1}{2^{\gamma}})C_{2}\tilde{\beta}_{j}^{-\frac{1}{2^{\gamma}}}}{\tilde{\beta}_{j}^{2}\left(1+C_{2}\tilde{\beta}_{j}^{-\frac{1}{2^{\gamma}}}\right)^{2}}$,
we have

$\tilde{\beta}_{j}^{2}=\frac{C_{1}}{\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j}}\frac{\tilde{\beta}_{j}^{\frac{1}{2^{\gamma}}}}{\tilde{\beta}_{j}^{\frac{1}{2^{\gamma}}}+C_{2}}\left(1-\frac{\frac{1}{2^{\gamma}}C_{2}}{\tilde{\beta}_{j}^{\frac{1}{2^{\gamma}}}+C_{2}}\right)$

$1-\frac{\frac{1}{2^{\gamma}}C_{2}}{\tilde{\beta}_{j}^{\frac{1}{2^{\gamma}}}+C_{2}}=M^{\prime}(\gamma)$
for $M^{\prime}(\gamma)\in(\frac{1}{2},1]$. Therefore,

$\tilde{\beta}_{j}^{(2-\frac{1}{2^{\gamma}})}=\frac{C_{1}M^{\prime}(\gamma)}{\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j}}\frac{1}{\tilde{\beta}_{j}^{\frac{1}{2^{\gamma}}}+C_{2}}=M^{\prime}(\gamma)\frac{p+\frac{1}{2^{\gamma+1}}}{n\tilde{\beta}_{j}^{\frac{1}{2^{\gamma}}}+nC_{2}}$

Since we have
$\tilde{\beta}_{j}<(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}^{(i)}|$,
it follows directly that

$2\left\\{\frac{C_{1}(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}}{2C_{2}+2[(\boldsymbol{X}_{j}^{T}\boldsymbol{X}_{j})^{-1}|z_{j}|]^{\frac{1}{2^{\gamma}}}}\right\\}^{\frac{1}{2-\frac{1}{2^{\gamma}}}}\leq
h(\tilde{\beta}_{j}).$

∎

## 4\. Forward and backward variable screening

In this section, we provide details of the proposed backward and forward
variable screening algorithms. The backward screening algorithm does not use
cross-validation, and is recommended for when the sample size is moderate to
large. For small sample size, a forward variable screening algorithm cross-
validation is recommended.

Consider a sequence of $L$ parameters
$g=\left\\{g_{1}<g_{2}<...<g_{L}\right\\}$ with $b_{l}=g_{l}*\frac{\log
p}{p}$, the CD optimization begins with hyper-parameter $b_{1}=g_{1}\frac{\log
p}{p}$. In default, we set $g=\left\\{1,2,...,100\right\\}$. The first output
$\hat{\beta}^{(1)}$ is then used as a warm start for the 2nd iteration of the
CD optimisation with $b_{2}=g_{2}\frac{\log p}{p}$. The algorithm continues in
this fashion and we can visualise the entire solution path of
$\boldsymbol{\beta}$ with respect to the hyper-parameter $b$. Since the
function $b\rightarrow T_{b}(\cdot)$ is continuous, this assures that a small
perturbation in $b$ will lead to only small changes in the solution for
$\boldsymbol{\beta}$, this suggests that using the solution of the previous
optimisation for the next one is a reasonable strategy. The backward screening
algorithm is provided in Algorithm 4.

Set $K$, $\gamma$ and $b=\left\\{b_{1}<b_{2}<...<b_{L}\right\\}$

$\boldsymbol{\beta}^{(0)}=0$

for  _$i\leftarrow 1$ to $L$_ do

$\boldsymbol{\beta}^{(i)}=$CD
algorithm($\boldsymbol{Y},\boldsymbol{X},b=b_{i},\gamma,\boldsymbol{\beta}_{initial}=\boldsymbol{\beta}^{(i-1)}$)

$\sigma^{2(i)}=\frac{\|\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{\beta}^{(i)}\|_{2}^{2}}{n-s^{(i)}}$

end for

Result:
$\boldsymbol{\beta}=(\boldsymbol{\beta}^{(1)},...,\boldsymbol{\beta}^{(L)})$,
$\sigma^{2}=(\sigma^{2(1)},...,\sigma^{2(L)})$

Algorithm 4 Backward variable screening

Roughly speaking, the higher the value of $b$, the more the estimated
$\hat{\boldsymbol{\beta}}$ are set to zero. Hence we use the term backward
variable screening. The algorithm starts with a big model and gradually
removes insignificant coefficients. When the sample size is relative large,
the solution path stabilizes very quickly with increasing values of $b$ (See
Figure 5). In fact, we found that, when the sample size is large, we can fit
the model directly by setting $\frac{1}{b}=0.0001$, but we still recommend
backward variable screening because we don’t know how large the sample size
needs to be.

When the sample size is small, the solution path stabilizes slowly, it may
ultimately threshold everything to zero, requiring cross-validation to
identify the best solution. Figure 6 shows the solutions paths from a
simulated example with $n=100,p=1000,s_{0}=100$ and $\sigma^{2}=3$. We can see
that the model quickly shrinks all coefficients to zero. In such cases, we
found that a forward screening strategy, together with cross-validation, works
better. Although cross-validation is computationally intensive, it was found
by Reid et al., (2016) that it guarantees a robust performance for small
sample size and less sparse models. Our numerical experiments are consistent
with the numerical results in their paper.

Figure 5: Solution paths of $\boldsymbol{\beta}$, as $b$ is varied. Nonzero
elements (left); zero elements (middle) and $\hat{s}$ (right). The data is
simulated using $n=500$, $p=1000$, $s_{0}=10$ and $\sigma_{0}^{2}=3$.
Figure 6: Solution paths of $\boldsymbol{\beta}$, as $b$ is varied. Nonzero
elements (left); zero elements (middle) and $\hat{s}$ (right). The data is
simulated using $n=100$, $p=1000$, $s_{0}=10$ and $\sigma_{0}^{2}=3$.

Here we develop a forward variable screening algorithm with K-fold cross-
validation. We first reparametrize the hyper-parameter $b$ to $t=\frac{1}{b}$.
Then we consider a sequence of $L$ increasing parameters
$T=\left\\{t_{1}<t_{2}<...<t_{L}\right\\}$ with $t_{1}=0$ and
$t_{L}=\frac{p}{\log p}$. In default, 100 evenly spaced points have been used.
We use the same re-initialization strategy as backward variable screening. The
optimization output from the previous value of $b$ is used as a warm start for
optimization at current $b$. For each $b$, we obtain $K$ solutions of
$\boldsymbol{\beta}$ from the $K$-fold cross-validation. Once we have the
optimal value of $b$, we fit the model to the entire data set with forward
strategies. Empirically, we found the forward strategies did better than
backward strategy when the sample size is small and the model is less sparse.

The sure independence screening by Fan and Lv, (2008) can be viewed as a
backward strategies. They greatly simplified the original ultra-high
dimensional problem into a low-dimensional one and then fitted the model with
LASSO or SCAD penalty. In another direction, Wang, (2009) investigated the
forward regression in ultra-high dimensional problems. Our forward variable
screening strategies can be viewed as a smooth version of forward regression.
Wang, (2009) used the BIC criterion to select the optimal model size. They
show screening consistency for their approach. Here, we determine the optimal
model by cross-validation. For infinite sample size, the performance of our
forward and backward screening algorithms are guaranteed by the oracle
properties of the non-separable bridge penalty. We provide the description of
forward screening algorithm in Algorithm 5.

Set $K$, $\gamma$ and $T=\left\\{t_{1}<t_{2}<...<t_{L}\right\\}$

Split the data into $K$-fold with
$\boldsymbol{Y}_{validation}=(\boldsymbol{Y}_{-1},...,\boldsymbol{Y}_{-K})$
and
$\boldsymbol{X}_{validation}=(\boldsymbol{X}_{-1},...,\boldsymbol{X}_{-K})$

for  _$k\leftarrow 1$ to $K$_ do

$\boldsymbol{\beta}^{(k,0)}=0$

for  _$i\leftarrow 1$ to $L$_ do

$\boldsymbol{\beta}^{(k,i)}=$CD
algorithm($\boldsymbol{Y}_{-k},\boldsymbol{X}_{-k},b=\frac{1}{t_{i}},\gamma,\boldsymbol{\beta}_{initial}=\boldsymbol{\beta}^{(k,i-1)}$)

end for

end for

for  _$i\leftarrow 1$ to $L$_ do

$Err_{valid}^{(i)}=\frac{1}{N}\sum_{k=1}^{K}\|\boldsymbol{Y}_{-k}-\boldsymbol{X}_{-k}\boldsymbol{\beta}^{(k,i)}\|_{2}^{2}$

end for

$m=\mathbf{argmin}_{i}\,\,Err_{valid}^{(i)}$

$\hat{\boldsymbol{\beta}}=0$

for  _$i\leftarrow 1$ to $m$_ do

$\hat{\boldsymbol{\beta}}=$CD
algorithm($\boldsymbol{Y},\boldsymbol{X},b=\frac{1}{t_{i}},\gamma,\boldsymbol{\beta}_{initial}=\hat{\boldsymbol{\beta}}$)

end for

$\hat{\sigma}^{2}=\frac{\|\boldsymbol{Y}-\boldsymbol{X}\hat{\boldsymbol{\beta}}\|_{2}^{2}}{n-\hat{s}}$

Result: $\hat{\boldsymbol{\beta}}$, $\hat{\sigma}^{2}$

Algorithm 5 Forward variable screening with cross-validation

## 5\. Extension to logistic model

### PCG sampler for Bayesian logistic model

The global-local scale mixture of gaussian representation for the
$L_{\frac{1}{2^{\gamma}}}$ prior allows us to extend the PCG sampler to the
model which can be expressed as gaussian mixtures or the likelihood which are
conjugate with gaussian priors. For example Bayesian quantile regression under
asymmetric Laplace distributions (Yu and Moyeed,, 2001) and logistic
regressions.

Here we only show how to extend the result to the logistic regression models
with binomial likelihood and a logistic link function. We do that by using the
data-augmentation trick in logistic model proposed by Polson et al., (2013).
Applying their technique here, let $y_{i}$ be the number of successes, $n_{i}$
the number of trials and $x_{i}=(x_{i1},...,x_{ip})$ the vector of predictors
for observation $i\in\left\\{1,2,3,...,N\right\\}$. Let
$y_{i}\sim\operatorname{Binom}\left(n_{i},1/\left\\{1+e^{-\psi_{i}}\right\\}\right)$,
where $\psi_{i}=x_{i}^{T}\boldsymbol{\beta}$ is the log odds of success. The
likelihood for observation $i$ is

$L_{i}(\boldsymbol{\beta})=\frac{\left\\{\exp\left(x_{i}^{T}\boldsymbol{\beta}\right)\right\\}^{y_{i}}}{\left\\{1+\exp\left(x_{i}^{T}\boldsymbol{\beta}\right)\right\\}^{n_{i}}}\propto\exp\left(\kappa_{i}x_{i}^{T}\boldsymbol{\beta}\right)\int_{0}^{\infty}\exp\left\\{-\omega_{i}\left(x_{i}^{T}\boldsymbol{\beta}\right)^{2}/2\right\\}p\left(\omega_{i}|n_{i},0\right)dw_{i}$

where $\kappa_{i}=y_{i}-n_{i}/2$ and $p\left(\omega_{i}|n_{i},0\right)$ is the
density of a Polya–Gamma random variable with parameters
$\left(n_{i},0\right)$. Then, conditional on $\omega$, we have a gaussian
likelihood of the form

$\prod_{i=1}^{N}L_{i}(\boldsymbol{\beta}|\omega_{i})\propto\exp\left\\{-\frac{1}{2}(\boldsymbol{z}-\boldsymbol{X}\boldsymbol{\beta})^{T}\Omega(\boldsymbol{z}-\boldsymbol{X}\boldsymbol{\beta})\right\\}$

where
$\boldsymbol{z}=\left(\kappa_{1}/\omega_{1},\ldots,\kappa_{N}/\omega_{N}\right)$
and $\Omega=\operatorname{diag}\left(\omega_{1},\ldots,\omega_{N}\right)$. let
the $L_{\frac{1}{2^{\gamma}}}$ prior be assigned to $\boldsymbol{\beta}$ ,
then we can construct the PCG Sampler scheme:

$\displaystyle{\bf S1.}$ $\displaystyle\mbox{Sample
}\boldsymbol{\beta}^{\prime}\sim\mathrm{N_{p}}\left((\boldsymbol{X}^{T}\Omega\boldsymbol{X}+\lambda^{2^{\gamma+1}}D_{\tau^{2}}^{-1})^{-1})\boldsymbol{X}^{T}\kappa,(\boldsymbol{X}^{T}\Omega\boldsymbol{X}+\lambda^{2^{\gamma+1}}D_{\tau^{2}}^{-1})^{-1}\right)$
$\displaystyle{\bf S2.}$ $\displaystyle\mbox{Sample
}\lambda^{\prime}\sim\mathrm{Gamma}\left(2^{\gamma}p+0.5,\sum_{j=1}^{p}\left|\beta^{\prime}_{j}\right|^{\frac{1}{2^{\gamma}}}+\frac{1}{b}\right)$
$\displaystyle{\bf S3.}$ $\displaystyle\mbox{Sample
}b^{\prime}\sim\mathrm{InvGamma}\left(1,1+\lambda\right)$ $\displaystyle{\bf
S4.}$ $\displaystyle\mbox{Sample
}\frac{1}{v^{\prime}_{\gamma,j}}\sim\mathrm{InvGaussian}\left(\frac{1}{2\lambda|\beta^{\prime}_{j}|^{\frac{1}{2^{\gamma}}}},\frac{1}{2}\right),\quad
j=1,...,p$ $\displaystyle\mbox{For }i=1,\ldots,\gamma-1,$
$\displaystyle\mbox{Sample
}\frac{1}{{v^{\prime}}_{i,j}}\sim\mathrm{InvGaussian}\left(\frac{1}{2{v^{\prime}}_{i+1,j}\lambda|{\beta^{\prime}}_{j}|^{\frac{1}{2^{i}}}},\frac{1}{2{v^{\prime}}_{i+1,j}^{2}}\right),\quad
j=1,...,p$ $\displaystyle{\bf S5.}$ $\displaystyle\mbox{Sample
}\frac{1}{{\tau^{\prime}}_{j}^{2}}\sim\mathrm{InvGaussian}\left(\frac{1}{\lambda^{2^{\gamma}}{v^{\prime}}_{1,j}|\beta^{\prime}_{j}|},\frac{1}{{v^{\prime}}_{1,j}^{2}}\right),\quad
j=1,...,p$ $\displaystyle{\bf S6.}$ $\displaystyle\mbox{Sample
}\omega_{i}\sim\mathrm{PG}(n_{i},x_{i}^{T}\boldsymbol{\beta}),\quad i=1,...,n$
$\displaystyle{\bf S7.}$ $\displaystyle\mbox{Set
}\boldsymbol{\beta},\lambda,\boldsymbol{v},\boldsymbol{\tau}^{2},\sigma^{2}\mbox{
to
}\boldsymbol{\beta}^{\prime},\lambda^{\prime},\boldsymbol{v}^{\prime},\boldsymbol{\tau}^{\prime
2},\sigma^{\prime 2}\mbox{ and increment }t=t+1.$

where $\kappa=(y_{1}-n_{1}/2,...,y_{N}-n_{N}/2)$. Here we have just slightly
modified the conditional posterior for $\boldsymbol{\beta}$ and replaced the
conditional posterior of $\sigma^{2}$ with $\omega_{i}$. The construction of
the PCG sampler is similar to Section 2.

### Coordinate descent with Proximal Newton map

In this section, we extend the coordinate descent algorithm from gaussian
likelihood to the binary logistic likelihood. Our strategy is the same as
Friedman et al., (2010). By coding the response as $0/1$, the likelihood can
be written as

$\displaystyle\ell(\boldsymbol{\beta})$
$\displaystyle=\sum_{i=1}^{N}\left\\{y_{i}\log
p\left(x_{i};\boldsymbol{\beta}\right)+\left(1-y_{i}\right)\log\left(1-p\left(x_{i};\boldsymbol{\beta}\right)\right)\right\\}$
$\displaystyle=\sum_{i=1}^{N}\left\\{y_{i}x_{i}\boldsymbol{\beta}-\log\left(1+e^{x_{i}\boldsymbol{\beta}}\right)\right\\}.$

If we apply the Newton algorithm to fit the (unpenalized) log-likelihood
above, that is equivalent to using the iteratively reweighted least squares
algorithm. For each iteration, we solve

$\boldsymbol{\beta}^{\text{new
}}\leftarrow\arg\min_{\beta}(\mathbf{z}-\mathbf{X}\boldsymbol{\beta})^{T}\mathbf{W}(\mathbf{z}-\mathbf{X}\boldsymbol{\beta})$
(20)

where $\mathbf{z}=\mathbf{X}\boldsymbol{\beta}^{\text{old
}}+\mathbf{W}^{-1}\left(Y-P(Y=1|\boldsymbol{X}\boldsymbol{\beta}^{old})\right)$
and
$W=\mathbf{diag}\left(P(Y=1|X\boldsymbol{\beta}^{old})*(1-P(Y=1|\boldsymbol{X}\boldsymbol{\beta}^{old}))\right)$
is $N\times N$ diagonal matrix of weights.

For logistic regression with non-separable bridge penalty, the fitting
algorithm consist of two steps. First, we create an outer loop, which computes
the quadratic approximation of (unpenalized) log-likelihood as in (20) based
on the old $\boldsymbol{\beta}$ from last iteration. Then in the inner loop,
we use coordinate descent to solve the penalized weighted least-squares
problem

$\boldsymbol{\beta}^{\text{new
}}\leftarrow\arg\min_{\boldsymbol{\beta}}(\mathbf{z}-\mathbf{X}\boldsymbol{\beta})^{T}\mathbf{W}(\mathbf{z}-\mathbf{X}\boldsymbol{\beta})+\mathbf{Pen}(\boldsymbol{\beta})$
(21)

The solution to equation (21) in the inner loop is known as Proximal Newton
Map (Lee et al.,, 2014). By slightly modifying the coordinate descent
algorithm for penalized least square problem in Section 5.2, we can get the
coordinate descent algorithm for penalized weighted least square problem.

## 6\. Posterior consistency and contraction rate

Assuming that the observations $\boldsymbol{Y}=(Y_{1},\ldots,Y_{n})$
consisting of $n$ independent observations with product measure
$P_{\theta}^{n}=P_{\theta,1}\times P_{\theta,2}\ldots\times P_{\theta,n}$,
where $\theta\in\Theta$. If the measure $P_{\theta,i}$ with density
$f_{\theta,X_{i}}(Y_{i})$ is considered random and distributed according to
$\Pi_{n}$, as it is the case in Bayesian inference, then the posterior
distribution $\Pi_{n}\left(B\mid Y_{1},\ldots,Y_{n}\right)$ is the conditional
distribution of $P_{\theta}^{n}$ based on the prior $\Pi_{n}$

$\Pi_{n}\left(B\mid
Y_{1},\ldots,Y_{n}\right)=\frac{\int_{B}\prod_{i=1}^{n}f_{\theta,X_{i}}\left(Y_{i}\right)d\Pi_{n}(\theta)}{\int_{\Theta}\prod_{i=1}^{n}f_{\theta,X_{i}}\left(Y_{i}\right)d\Pi_{n}(\theta)},$

where $B\in\mathcal{B}$ and $\mathcal{B}$ is a $\sigma$-field on $\Theta$.

The posterior is said to be weakly consistent if it concentrates on
arbitrarily small neighborhoods of the true $P_{\theta_{0}}^{n}$,
$\theta_{0}\in\Theta$, in $P_{\theta_{0}}^{n}$-probability tending to 1 as
$n\rightarrow{\infty}$ or strongly consistent if this convergence is in the
almost-sure sense. We study the rate at which the neighborhoods of
$P_{\theta_{0}}^{n}$ may decrease to a small ball while still capturing most
of the posterior mass, in the weak sense (See Chapter 6 of Ghosal and Van der
Vaart, (2017)). Note that in high dimensional regression setting our prior
depends on the number of predictors, and hence sample size dependent,
therefore we use $\Pi_{n}$ to denote this dependence.

The following notation will be used in this section and in proving the oracle
property in a later section: we rewrite the dimension $p$ of the model by
$p_{n}$ to indicate that the number of predictors can increase with the sample
size $n$. We use $\beta_{0}$ and $\sigma_{0}$ to indicate true regression
coefficients and the true standard deviation and
$\theta_{0}=(\beta_{0},\sigma_{0})$. Let $S_{0}$ be the set containing the
indices of the true nonzero coefficients, where
$S_{0}\subseteq\left\\{1,..,p_{n}\right\\}$, then $s_{0}=|S_{0}|$ denotes the
size of the true model. The $\varepsilon$-covering number of the space
$\mathcal{F}$ with respect to a metric $d$ is denoted by
$N(\varepsilon,\mathcal{F},d)$, which is the minimal number of $d$-balls of
radius $\varepsilon$ needed to cover $\mathcal{F}$. We further let
$\boldsymbol{X}_{s}$ denote the design matrix corresponding to the model $S$,
and use $Pf$ to abbreviate $\int fdP$, finally we use $\gtrsim$ and $\lesssim$
to denote inequality up to a constant.

We make the following assumptions on contraction theorems:

A1.

The dimensional is high, with $p_{n}>>n$ and $\log(p_{n})=o(n)$.

A2.

The true number of nonzero $\beta_{0}$ is $s_{0}$, where $s_{0}=o(n/\log
p_{n})$.

A3.

All the covariates are uniformly bounded. In other words, there exists a
constant $k>0$ such that $\lambda_{\max}(X^{T}X)\leq{kn^{\alpha}}$ for some
$\alpha\in{[1,+\infty)}$.

A4.

There exist constants $v_{1}>0$, $v_{2}>0$ and an integer $\widetilde{p}$
satisfying $s_{0}=o(\widetilde{p})$ and $\widetilde{p}=o(s_{0}\log n)$, so
that
$nv_{1}\leq{\lambda_{min}(X_{S}^{T}X_{S})}\leq{\lambda_{max}(X_{S}^{T}X_{S})}\leq{nv_{2}}$
for any model of size $|s|\leq{\widetilde{p}}$.

A5.

$||\beta_{0}||_{\infty}={E}$ and E is some positive number independent of $n$.

### Some preliminary results

Before proving the posterior contraction rates, we first state and prove some
auxiliary results.

######  Lemma 7.1.

Suppose that assumptions A1, A2 and A5 hold, the exponential power prior
$\Pi_{n}(\beta)=\frac{\lambda_{n}^{\frac{1}{\alpha}}}{2\Gamma(1+\frac{1}{\alpha})}\exp{\left(-\lambda_{n}|\beta|^{\alpha}\right)}$
with $0<\alpha\leq\frac{1}{2}$ and $\lambda_{n}\propto\log p_{n}$, satisfies

$\int_{|\beta|\geq{a_{n}}}\Pi_{n}(\beta)d\beta\leq p_{n}^{-(1+u)}$ (22)
$-\log\left(\inf_{\beta\in\left[-E,E\right]}\Pi_{n}(\beta)\right)=O\left(\log
p_{n}\right)$ (23)

where $u,E$ are positive constants and
$a_{n}=\frac{\sigma_{0}}{\sqrt{2k}}\sqrt{s\log p_{n}/n}/p_{n}$.

###### Proof.

It is easy to see that

$-\log\left(\inf_{\beta\in\left[-E,E\right]}\Pi_{n}(\beta)\right)=\lambda_{n}E^{\alpha}-\frac{1}{\alpha}\log\lambda_{n}+\log
2+\log\Gamma(1+\frac{1}{\alpha})$