{\bf b}egin{document}

{\bf t}itle{A Laplace Mixture Representation of the Horseshoe and Some Implications}

{{\bf b}f a}uthor{Ksheera Sagar and Anindya Bhadra
{\bf t}hanks{Submitted on 11/28/2022. Bhadra is supported by US National Science Foundation Grant DMS-2014371.}
{\bf t}hanks{The authors are with the Department of Statistics, Purdue University, IN 47906, United States of America  (e-mails: kkeralap@purdue.edu; bhadra@purdue.edu).}}
{{\bf b}f m}aketitle

{\bf b}egin{abstract}
The horseshoe prior, defined as a half Cauchy scale mixture of normal, provides a state of the art approach to Bayesian sparse signal recovery. We provide a new representation of the horseshoe density as a scale mixture of the Laplace density, explicitly identifying the mixing measure. Using the celebrated Bernstein--Widder theorem and a result due to Bochner, our representation immediately establishes the complete monotonicity of the horseshoe density and strong concavity of the corresponding penalty. Consequently, the equivalence between local linear approximation and expectation--maximization algorithms for finding the posterior mode under the horseshoe penalized regression is established. Further, the resultant estimate is shown to be sparse.
{{\bf r}m e}nd{abstract}

{\bf b}egin{IEEEkeywords}
Bayesian sparse signal recovery, Bernstein--Widder, completely monotone function, local linear approximation
{{\bf r}m e}nd{IEEEkeywords}

{{\bf b}f I}EEEpeerreviewmaketitle

{{\bf b}f s}ection{Introduction}
{{\bf b}f I}EEEPARstart{F}{or} response $y$, covariates ${\bf P}hi$ and regression coefficients $x$, consider the high-dimensional linear regression model:
{\bf b}egin{eqnarray}
y={\bf P}hi x+ {{\bf r}m e}psilon,\label{regression_model}
{{\bf r}m e}nd{eqnarray}
where $y=\{y_j\}\in {{\bf b}f m}athbb{R}^n, {\bf P}hi=\{{{\bf b}f p}hi_{ji}\}\in {{\bf b}f m}athbb{R}^{n{\bf t}imes p}, x=\{x_i\}\in{{\bf b}f m}athbb{R}^{p}$ with $p{{\bf b}f g}g n$ and ${{\bf r}m e}psilon {{\bf b}f s}im {{\bf b}f m}athcal{N}(0,{{\bf b}f s}igma^2{{\bf b}f m}athbf{I}_n)$.
The horseshoe prior hierarchy {\bf c}ite{carvalho2010horseshoe} is defined for $i=i,\ldots,p$ as a half Cauchy scale mixture of normal as:
{\bf b}egin{equation}
x_i {{\bf b}f m}id \lambda_i, {\bf t}au {{\bf b}f s}tackrel{ind}{{\bf b}f s}im {{\bf b}f m}athcal{N}(0,\lambda_i^2{\bf t}au^2),\; \lambda_i{{\bf b}f s}tackrel{ind}{{\bf b}f s}im{{\bf b}f m}athcal{C}^{+}(0,1), \; {\bf t}au>0, \label{eq:prior}
{{\bf r}m e}nd{equation}
where ${{\bf b}f m}athcal{C}^{+}$ denotes a standard half Cauchy random variable with density $p(\lambda_i) = (2/{{\bf b}f p}i)(1+ \lambda_i^2)^{-1},\; \lambda_i>0$ and ${\bf t}au$ is a hyperparameter. This prior has found numerous applications in sparse signal recovery{\bf c}ite{carvalho2009handling, carvalho2010horseshoe, li2019graphical, yu2019variational, bhadra2019lasso, bhadra2020horseshoe, lee2020variation, bhadra2021horseshoe, li2021joint}. Although there is no closed form to the marginal horseshoe density, obtained by integrating out the latent $\lambda_i$s, the main benefits of the horseshoe stem from two crucial properties of the marginal prior. The first is infinite prior mass at the origin, allowing strong shrinkage towards zero in a sparse regime; and the second is heavy, polynomially decaying tails, aiding detection of signals {\bf c}ite{carvalho2010horseshoe}. 

Since the tails of the horseshoe prior follow a quadratic rate of decay, it is conceivable to expect the prior density to admit a mixture representation with respect to the double exponential (Laplace) density as well. Yet, to our knowledge, no such explicit representation is available in the literature. We close this gap by showing the horseshoe density can also be represented as a mixture of double exponential with respect to a random variable whose density is proportional to the Dawson function. This representation, in turn, allows us to appeal to the celebrated Bernstein--Widder theorem {\bf c}ite{widder2015laplace} and a result due to Bochner {\bf c}ite{bochner2005harmonic} to establish complete monotonicity of the horseshoe density and  strong concavity of the horseshoe penalty, both new results. Further, appealing to the same mixture representation, we are able to establish an equivalence between local linear approximation or LLA {\bf c}ite{zou2008one,polson2016mixtures} and expectation--maximization or EM {\bf c}ite{dempster1977maximum} schemes for finding the maximum a posteriori (MAP) estimate under the horseshoe prior for sparse signal recovery.

{{\bf b}f s}ection{A Laplace Mixture Representation of the Horseshoe}
{\bf b}egin{proposition}
\label{prop_equiv_HS_Laplace}
The marginal horseshoe density for a scalar random variable ($p=1$) admits the following representation as a Laplace mixture:
{\bf b}egin{eqnarray*}
p_{HS}(x) &=& {{\bf b}f f}rac{{2}}{{{\bf b}f p}i{{\bf b}f s}qrt{{{\bf b}f p}i}{\bf t}au} \int_0^\infty {{\bf r}m e}xp\left(-u{{\bf b}f f}rac{|x|}{{\bf t}au}{\bf r}ight) D_{+}\left({{\bf b}f f}rac{u}{{{\bf b}f s}qrt{2}}{\bf r}ight) du,
{{\bf r}m e}nd{eqnarray*}
where $D_{+}(z)={{\bf r}m e}xp(-z^2)\int_0^z{{\bf r}m e}xp(t^2)dt,\; z>0,$ is the Dawson function.
{{\bf r}m e}nd{proposition}

{\bf b}egin{proof}
We have using the hierarchy of Equation~{{\bf r}m e}qref{eq:prior},
{\bf b}egin{equation*}
{\bf b}egin{split}
p_{HS}(x)=& \int_0^\infty {{\bf b}f f}rac{1}{{{\bf b}f s}qrt{2{{\bf b}f p}i\lambda^2{\bf t}au^2}}  {{\bf r}m e}xp\left(-{{\bf b}f f}rac{x^2}{2\lambda^2 {\bf t}au^2}{\bf r}ight) {{\bf b}f f}rac{2}{{{\bf b}f p}i} {{\bf b}f f}rac{1}{1+\lambda^2} d\lambda
{{\bf r}m e}nd{split}
{{\bf r}m e}nd{equation*}

{\bf b}egin{equation*}
{\bf b}egin{split}
=& {{\bf b}f f}rac{{{\bf b}f s}qrt{2}}{{{\bf b}f p}i{{\bf b}f s}qrt{{{\bf b}f p}i}{\bf t}au} \int_0^\infty  {{\bf b}f f}rac{1}{\lambda} {{\bf r}m e}xp\left(-{{\bf b}f f}rac{x^2}{2\lambda^2{\bf t}au^2}{\bf r}ight)  \int_0^\infty {{\bf r}m e}xp(-(\lambda^2 +1) t) dt d\lambda\\
=& {{\bf b}f f}rac{{{\bf b}f s}qrt{2}}{{{\bf b}f p}i{{\bf b}f s}qrt{{{\bf b}f p}i}{\bf t}au} \int_0^\infty \int_0^\infty  {{\bf b}f f}rac{1}{\lambda} {{\bf r}m e}xp\left(-{{\bf b}f f}rac{x^2}{2\lambda^2{\bf t}au^2} - \lambda^2 t {\bf r}ight)  {{\bf r}m e}xp(-t) dt d\lambda\\
=& {{\bf b}f f}rac{{{\bf b}f s}qrt{2}}{{{\bf b}f p}i{{\bf b}f s}qrt{{{\bf b}f p}i}{\bf t}au} \int_0^\infty \left\{\int_0^\infty  {{\bf r}m e}xp\left(-{{\bf b}f f}rac{x^2}{2\omega{\bf t}au^2} - \omega t {\bf r}ight){{\bf b}f f}rac{1}{2} \omega^{-1} d\omega{\bf r}ight\} {{\bf r}m e}xp(-t) dt\\
& ({\bf t}ext{substituting }\lambda^2 = \omega, {\bf t}ext{and using Fubini's theorem})\\
=& {{\bf b}f f}rac{1}{{{\bf b}f p}i{{\bf b}f s}qrt{2{{\bf b}f p}i}{\bf t}au} \int_0^{\infty} 2 K_0\left({{\bf b}f s}qrt{2t} {{\bf b}f f}rac{|x|}{{\bf t}au}{\bf r}ight) {{\bf r}m e}xp(-t) dt\\
&{\bf t}ext{(generalized inverse Gaussian integral {\bf c}ite{barndorff1977exponentially};}\\
&{\bf t}ext{$K_0=$ modified Bessel function of the second kind of order 0)}\\
=&  {{\bf b}f f}rac{{{\bf b}f s}qrt{2}}{{{\bf b}f p}i{{\bf b}f s}qrt{{{\bf b}f p}i}{\bf t}au} \int_0^{\infty} \left\{\int_0^{\infty} {{\bf r}m e}xp\left(-{{\bf b}f f}rac{|x|}{{\bf t}au} {{\bf b}f s}qrt{2t} {\bf c}osh {{\bf b}f z}eta{\bf r}ight) d{{\bf b}f z}eta{\bf r}ight\} {{\bf r}m e}xp(-t) dt\\
&{\bf t}ext{(using identity 9.6.24 in {\bf c}ite{abramowitz1964handbook})}\\
=&  {{\bf b}f f}rac{{{\bf b}f s}qrt{2}}{{{\bf b}f p}i{{\bf b}f s}qrt{{{\bf b}f p}i}{\bf t}au} \int_0^{\infty} \int_0^{\infty} {{\bf r}m e}xp\left(-{{\bf b}f f}rac{|x|}{{\bf t}au} {{\bf b}f s}qrt{2t} {\bf c}osh {{\bf b}f z}eta -t {\bf r}ight) d{{\bf b}f z}eta dt.
{{\bf r}m e}nd{split} 
{{\bf r}m e}nd{equation*}

{{\bf b}f n}oindent Let ${{\bf b}f s}qrt{2t} {\bf c}osh {{\bf b}f z}eta = u$ and $t=v$. The inverse transformations are, ${\bf c}osh^{-1} (u/{{\bf b}f s}qrt{2v}) = {{\bf b}f z}eta$ and $v=t$; with Jacobian:
{\bf b}egin{eqnarray*}
J&=& {\bf b}egin{vmatrix}
    {{\bf b}f f}rac{{{\bf b}f p}artial {{\bf b}f z}eta} {{{\bf b}f p}artial u}  & {{\bf b}f f}rac{{{\bf b}f p}artial {{\bf b}f z}eta} {{{\bf b}f p}artial v} \\
    & \\
    {{\bf b}f f}rac{{{\bf b}f p}artial t} {{{\bf b}f p}artial u}& {{\bf b}f f}rac{{{\bf b}f p}artial t} {{{\bf b}f p}artial v} 
  {{\bf r}m e}nd{vmatrix}= {\bf b}egin{vmatrix}
    {{\bf b}f f}rac{{{\bf b}f p}artial {{\bf b}f z}eta} {{{\bf b}f p}artial u}  & {{\bf b}f f}rac{{{\bf b}f p}artial {{\bf b}f z}eta} {{{\bf b}f p}artial v} \\
    & \\
    0 & 1
  {{\bf r}m e}nd{vmatrix}= \left|{{\bf b}f f}rac{{{\bf b}f p}artial{{\bf b}f z}eta} {{{\bf b}f p}artial u}{\bf r}ight| =  {{\bf b}f f}rac{1}{{{\bf b}f s}qrt{u^2 - 2v}}.
{{\bf r}m e}nd{eqnarray*}
Thus, $p_{HS}(x)$ equals,
{\bf b}egin{eqnarray*}
&&{{\bf b}f f}rac{{{\bf b}f s}qrt{2}}{{{\bf b}f p}i{{\bf b}f s}qrt{{{\bf b}f p}i}{\bf t}au} \int_0^{\infty} \int_0^{u^2/2} {{\bf r}m e}xp\left(-{{\bf b}f f}rac{|x|}{{\bf t}au}u -v {\bf r}ight) {{\bf b}f f}rac{1}{{{\bf b}f s}qrt{u^2 - 2v}} dv du\\
&=& {{\bf b}f f}rac{{{\bf b}f s}qrt{2}}{{{\bf b}f p}i{{\bf b}f s}qrt{{{\bf b}f p}i}{\bf t}au} \int_0^{\infty} {{\bf r}m e}xp\left(-{{\bf b}f f}rac{|x|}{{\bf t}au}u{\bf r}ight)   \left\{\int_0^{u^2/2}  {{\bf b}f f}rac{{{\bf r}m e}xp(-v )}{{{\bf b}f s}qrt{u^2 - 2v}} dv {\bf r}ight\} du\\
&=& {{\bf b}f f}rac{{2}}{{{\bf b}f p}i{{\bf b}f s}qrt{{{\bf b}f p}i}{\bf t}au} \int_0^{\infty} {{\bf r}m e}xp\left(-{{\bf b}f f}rac{|x|}{{\bf t}au}u{\bf r}ight) \left\{  e^{-u^2/2}\int_0^{{{\bf b}f f}rac{u}{{{\bf b}f s}qrt{2}}}e^{t^2}dt {\bf r}ight\} du\\
&&({\bf t}ext{substituting }u^2-2v=2t^2)\\
&=& {{\bf b}f f}rac{{2}}{{{\bf b}f p}i{{\bf b}f s}qrt{{{\bf b}f p}i}{\bf t}au} \int_0^{\infty} {{\bf r}m e}xp\left(-{{\bf b}f f}rac{|x|}{{\bf t}au}u{\bf r}ight) D_{+}\left({{\bf b}f f}rac{u}{{{\bf b}f s}qrt{2}}{\bf r}ight) du. \qedhere\\
{{\bf r}m e}nd{eqnarray*}
{{\bf r}m e}nd{proof}
Proposition~{\bf r}ef{prop_equiv_HS_Laplace} is our main result and the remainder of the paper discusses its implications in sparse signal recovery. Numerical validation of Proposition~{\bf r}ef{prop_equiv_HS_Laplace} is presented in Fig.~{\bf r}ef{fig:validation} by evaluating the horseshoe density using both the usual normal scale mixture (cf. Equation~{{\bf r}m e}qref{eq:prior}) and the newly derived Laplace scale mixture representations over $[-2,2]$. The point-wise maximum difference in density evaluation (in logarithmic scale) between the two representations is less than $10^{-10}$.   
{\bf b}egin{figure*}[tb]
{\bf c}entering
\includegraphics[width={\bf t}extwidth, height =0.18{\bf t}extheight]{Plot_prop_verification.pdf}
{\bf c}aption{
\label{fig:validation}
Numerical validation of Proposition~{\bf r}ef{prop_equiv_HS_Laplace}. Left panel: logarithm of horseshoe density, $\log p_{HS}(x)$, evaluated numerically, using the Laplace scale mixture in Proposition~{{\bf r}m e}qref{prop_equiv_HS_Laplace}. Middle panel: logarithm of horseshoe density, $\log p_{HS}(x)$, evaluated numerically, using the normal scale mixture in Equation~{{\bf r}m e}qref{eq:prior}. Right panel: Difference between the logarithms of densities evaluated using normal and Laplace mixtures.} 
{{\bf r}m e}nd{figure*}

{{\bf b}f s}ection{Main Implications of Proposition~{\bf r}ef{prop_equiv_HS_Laplace}: Complete Monotonicity of 
\label{main_implications}
the Horseshoe Density and Strong Concavity of the Penalty}
We begin by stating the definition of a completely monotone function and the celebrated Bernstein--Widder theorem {\bf c}ite[Chapter 4]{widder2015laplace}.
{\bf b}egin{definition}
A function $g({\bf c}dot)$ is said to be completely monotone if $(-1)^k g^{(k)}(x){{\bf b}f g}e 0$ for $x{{\bf b}f g}e0,\; k=0,1,\ldots,$ with superscripts $(k)$ denoting the order of derivatives and $g^{(0)}{{\bf r}m e}quiv g$.
{{\bf r}m e}nd{definition}
{\bf b}egin{theorem}\label{th:bernstein} (Bernstein--Widder). The function $g({\bf c}dot)$ is completely monotone if and only if it is the Laplace transform of a positive measure ${{\bf b}f g}amma$, i.e., admits the representation: $g(x) = \int_{0}^{\infty} {{\bf r}m e}xp(-xu) d{{\bf b}f g}amma(u),\; x{{\bf b}f g}e 0$.
{{\bf r}m e}nd{theorem}
Since the Dawson function is positive over $(0,\infty)$, and the horseshoe density is symmetric around zero, Proposition~{\bf r}ef{prop_equiv_HS_Laplace} in conjunction with Theorem~{\bf r}ef{th:bernstein} immediately establishes the complete monotonicity of $p_{HS}(|x|)$.

In a penalized likelihood setting, it is customary to view the negative log likelihood as the loss function and the negative log prior density as the penalty in an optimization problem. The following result, due to Bochner {\bf c}ite[Theorem 4.1.5]{bochner2005harmonic} establishes the strong concavity of the horseshoe penalty.
{\bf b}egin{theorem}\label{lem:bochner} (Bochner). The function $g(x) = {{\bf r}m e}xp(-a {{\bf b}f p}hi(x)), x{{\bf b}f g}e 0,$ is completely monotone for every $a>0$ if and only if ${{\bf b}f p}hi'(x)={{\bf b}f p}hi^{(1)}(x)$ is completely monotone.
{{\bf r}m e}nd{theorem}
Since $p_{HS}(|x|)$ is completely monotone, Theorem~{\bf r}ef{lem:bochner} implies the corresponding penalty function ${{\bf b}f m}athrm{pen}_{HS}(|x|)=-\log p_{HS}(|x|)$  is a Bernstein function (integral of a completely monotone function) and thus strongly concave. In other words, ${{\bf b}f m}athrm{pen'}_{HS}(|x|) = -(p'_{HS}(|x|))/(p_{HS}(|x|))$ is completely monotone. This property of the penalty has important implications on the convergence of optimization algorithms for penalized likelihood approaches, which we discuss next.

{{\bf b}f s}ection{Sparse Signal Recovery via LLA}
{{\bf b}f s}ubsection{Local Linear Approximation of Non-convex Penalties}
Consider the penalized estimation problem:
{\bf b}egin{equation}
\operatornamewithlimits{argmin}_{x}\left({{\bf r}m e}ll(x; y) + {{\bf b}f m}athrm{pen}(|x|){\bf r}ight), \label{eq:pen}
{{\bf r}m e}nd{equation}
where ${{\bf r}m e}ll(x;y)= -\log {{\bf b}f m}athcal{L}(x;y)$ is a loss function providing some measure of model fit to the data $y$, usually a negative log likelihood function and ${{\bf b}f m}athrm{pen}(|x|)$ is a penalty symmetric around zero, with $|{\bf c}dot|$ denoting the vector ${{\bf r}m e}ll_1$ norm. Although convex penalties are simpler to study, both theoretically and computationally, concave penalties enjoy several desirable statistical properties such as near-unbiasedness, minimax optimal properties and sparsity~{\bf c}ite{fan2001variable, fan2014strong}. Nevertheless, the optimization problem becomes much more challenging, and usually, it is only possible to find a local optimum. Required conditions for global optimum under concave penalties are discussed by~{\bf c}ite{mazumder2011sparsenet}.

The local linear approximation (LLA) algorithm of {\bf c}ite{zou2008one} provides one approach for solving the problem in Equation~{{\bf r}m e}qref{eq:pen} for concave penalties. LLA iterations proceed via a first order Taylor approximation of the penalty function, as follows
{\bf b}egin{equation*}
   {{\bf b}f m}athrm{pen}(|x|) {{\bf b}f a}pprox   {{\bf b}f m}athrm{pen}(|x_{(k)}|) + {{\bf b}f m}athrm{pen'}(|x_{(k)}|)(|x|-|x_{(k)}|),
{{\bf r}m e}nd{equation*}
{{\bf b}f n}oindent where $x_{(k)}$ is the current parameter value. With this approximation, and assuming the loss (but not necessarily the penalty) is convex, the non-convex optimization problem in Equation~{{\bf r}m e}qref{eq:pen} reduces to the convex problem
{\bf b}egin{equation}
\label{skeleton_beta_update}
    x_{(k+1)} = {\bf u}nderset{x}{{\bf t}ext{argmin}}\left({{\bf r}m e}ll(x;y)+\lambda_{(k)}|x|{\bf r}ight),
{{\bf r}m e}nd{equation}
where $\lambda_{(k)}={{\bf b}f m}athrm{pen'}(|x_{(k)}|)$. With the common squared error loss function, ${{\bf r}m e}ll(x;y) = |x-y|^2$, Equation~{{\bf r}m e}qref{skeleton_beta_update} defines a lasso problem {\bf c}ite{tibshirani1996regression}, which can be solved very efficiently. In addition, the resultant solution is naturally sparse.

{{\bf b}f s}ubsection{Implications of Complete Monotonicity on LLA}
Numerical convergence and statistical properties of LLA estimates are non-trivial objects of study for general non-convex penalties and properties of LLA for general penalties are largely unknown. Nevertheless, in the special case when the penalty under consideration is a Bernstein function, certain simplifications ensue due to the following result of {\bf c}ite{zou2008one}.
{\bf b}egin{theorem} (Zou and Li). 
Consider a penalty ${{\bf b}f m}athrm{pen}(|x|)$ such that ${{\bf r}m e}xp(-{{\bf b}f m}athrm{pen}(|x|))$ is completely monotone. Then, under a Gaussian likelihood, the LLA algorithm is equivalent to an EM algorithm for MAP estimation under a prior, not necessarily proper, with density $p_{X}(x){{\bf b}f p}ropto {{\bf r}m e}xp(-{{\bf b}f m}athrm{pen}(x))$.
{{\bf r}m e}nd{theorem}
Thus, an application of Proposition~{\bf r}ef{prop_equiv_HS_Laplace} immediately establishes the LLA algorithm as an EM algorithm for horseshoe penalties. This result is of practical significance, since the convergence properties of LLA can then be studied using the same tools for EM, which are better understood~{\bf c}ite{wu1983convergence}. This equivalence is also established by {\bf c}ite{polson2016mixtures}, who provide a variational interpretation. Thus, it is established that the MAP estimate under the horseshoe prior is sparse; due to the soft thresholding operator in LLA resulting in exact zeros. Further, to our knowledge, this is the first viable algorithm for finding the MAP estimate under the horseshoe.

{{\bf b}f s}ubsection{Implementing LLA under the Horseshoe Penalty}
Consider now the practicalities of LLA under the horseshoe penalty. For implementation, one needs to compute $\lambda_{(k)}={{\bf b}f m}athrm{pen'}_{HS}(|x_{(k)}|)$. Differentiating the expression of Proposition~{\bf r}ef{prop_equiv_HS_Laplace} under the integral sign, we have:
{\bf b}egin{equation}
\label{derivative_peanlty}
   \lambda_{(k)} = {{\bf b}f m}athrm{pen'}_{HS}(|x_{(k)}|) = {{\bf b}f f}rac{\int_{0}^\infty {{\bf b}f f}rac{u}{{\bf t}au} {{\bf r}m e}xp\left(-u{{\bf b}f f}rac{|x_{(k)}|}{{\bf t}au}{\bf r}ight)D_+\left({{\bf b}f f}rac{u}{{{\bf b}f s}qrt{2}}{\bf r}ight)du}{\int_{0}^\infty{{\bf r}m e}xp\left(-u{{\bf b}f f}rac{|x_{(k)}|}{{\bf t}au}{\bf r}ight)D_+\left({{\bf b}f f}rac{u}{{{\bf b}f s}qrt{2}}{\bf r}ight)du}.
{{\bf r}m e}nd{equation}
We provide numerical details of evaluating $\lambda_{(k)}$ in the Appendix. Subsequent to this, we use efficient solvers based on pathwise coordinate optimization~{\bf c}ite{friedman2007pathwise}, for solving the lasso problem in Equation~{{\bf r}m e}qref{skeleton_beta_update}. We note here that using the expression from Equation~{{\bf r}m e}qref{derivative_peanlty} is not necessary for evaluating ${{\bf b}f m}athrm{pen'}_{HS}$, strictly speaking, and other alternative numerical methods can be used. However, in our experience, Equation~{{\bf r}m e}qref{derivative_peanlty} is numerically very stable.

{{\bf b}f s}ection{Numerical Experiments}\label{sec:num}

{{\bf b}f s}ubsection{Sparse Normal Means Model}
\label{normal_means_model}

Revisiting the regression model in Equation~{{\bf r}m e}qref{regression_model}, and setting $p=n,\,{\bf P}hi ={{\bf b}f m}athbf{I}_n$, we recover the normal means model $y_i{{\bf b}f s}tackrel{ind}{{\bf b}f s}im{{\bf b}f m}athcal{N}(x_i,{{\bf b}f s}igma^2)$. In our experiments, we set $n=50,\,{{\bf b}f s}igma^2=1$ and the components in $x$ as, $x_1=x_2=3,\,x_3,\ldots,x_{50}=0$; $y_i{\bf t}ext{ is sampled from }{{\bf b}f m}athcal{N}(x_i,{{\bf b}f s}igma^2)$. With the generated data, optimizing Equation~{{\bf r}m e}qref{skeleton_beta_update} with ${{\bf r}m e}ll(x;y) = {{\bf b}f f}rac{1}{2}{{\bf b}f s}um_{i=1}^n(y_i - x_i)^2$, with pathwise coordinate optimization~{\bf c}ite{friedman2007pathwise}, we have :
{\bf b}egin{equation*}
    {\bf b}egin{split}
         x_{(k+1)} &= {\bf u}nderset{x}{{\bf t}ext{argmin}}\left({{\bf b}f f}rac{1}{2}{{\bf b}f s}um_{i=1}^n(y_i-x_i)^2+{{\bf b}f s}um_{i=1}^n\lambda_{i(k)}|x_i|{\bf r}ight),\\
    {\bf t}ext{yielding, }  x_{i(k+1)} &= S(y_i,\lambda_{i(k)}),{\bf t}ext{ where},
    {{\bf r}m e}nd{split}
{{\bf r}m e}nd{equation*}
{\bf b}egin{equation}
    \label{soft_threshold}
    S({\bf b}eta,{{\bf b}f g}amma) = {\bf b}egin{cases}{\bf b}eta - {{\bf b}f g}amma,&{\bf t}ext{ if }{\bf b}eta>0{\bf t}ext{ and } {{\bf b}f g}amma<|{\bf b}eta|,\\
    {\bf b}eta + {{\bf b}f g}amma,&{\bf t}ext{ if }{\bf b}eta<0{\bf t}ext{ and } {{\bf b}f g}amma<|{\bf b}eta|,\\
    0, &{\bf t}ext{ if }{{\bf b}f g}amma{{\bf b}f g}eq|{\bf b}eta|,
    {{\bf r}m e}nd{cases}
{{\bf r}m e}nd{equation}
is the soft thresholding operator. In the above, $x_{i(k+1)}$ denotes the value of $i^{th}$ component of $x$, at $(k+1)^{th}$ iteration. Similarly, $\lambda_{i(k)}$ denotes the derivative of the penalty in Equation~{{\bf r}m e}qref{derivative_peanlty}, evaluated at $x_{i(k)}$. We estimate the normal means, ${{\bf b}f h}at{x}=\{{{\bf b}f h}at{x}_i\}$, under the horseshoe penalty using LLA; and compare with the posterior mean of  $x$ under the horseshoe prior obtained by Markov chain Monte Carlo (HS-MCMC), the minimax concave penalty or MCP {\bf c}ite{zhang2010nearly}, the smoothly clipped absolute deviation penalty or SCAD {\bf c}ite{fan2001variable} and the lasso {\bf c}ite{tibshirani1996regression} penalties, using the {\bf t}exttt{R} packages {\bf t}exttt{horseshoe}~{\bf c}ite{van2016horseshoe}, {\bf t}exttt{sparsenet}~{\bf c}ite{mazumder2011sparsenet}, {\bf t}exttt{ncvreg}~{\bf c}ite{ncvreg} and {\bf t}exttt{glmnet}~{\bf c}ite{glmnet} respectively. MCP are SCAD are non-convex penalties, while lasso is convex. Under all the point estimation methods, the respective tuning parameters are chosen using 10-fold cross validation; whereas for HS-MCMC, we consider ${\bf t}au{{\bf b}f s}im{{\bf b}f m}athcal{C}^+(0,1)$, a total of $1.5{\bf t}imes 10^4$ posterior samples with 5000 samples as burnin, and 50\% credible intervals for variable selection~{\bf c}ite{li2019graphical}. For the LLA procedure, we start with $x_{i(0)}=1,\,{{\bf b}f f}orall  i$ and set the stopping criterion as ${{{\bf b}f s}um_{i=1}^n\left({{\bf b}f h}at{x}_{i(k+1)}-{{\bf b}f h}at{x}_{i(k)}{\bf r}ight)^2}<10^{-6}$. To compare the estimates, we compute the in sample sum of squared errors or ${\bf t}ext{SSE}={{\bf b}f s}um_{i=1}^n({{\bf b}f h}at{x}_i-x_i)^2$, out of sample prediction SSE (pSSE) $={{\bf b}f n}orm{y_{{\bf b}f m}athrm{out}-{\bf P}hi{{\bf b}f h}at{x}}^2$; for $y_{{\bf b}f m}athrm{out}$ generated similarly as $y$, true negative (positive) rates, for zeros (non-zeros) in $x$ identified correctly as TNR, TPR  respectively and the computational time in seconds. We estimate ${{\bf b}f h}at{x}$ for 50 replications, and present the mean and standard deviation (in parentheses) of the performance measures in Table~{\bf r}ef{normal_means_result}. We observe that HS-LLA has the best pSSE and its SSE and TNR are competitive with the resultant MCMC estimate. The good statistical performance of the posterior mean estimate obtained by MCMC is not surprising, since it is the optimal Bayes estimate under ${{\bf r}m e}ll_2$ loss. However, the MCMC estimate is about 3 times slower than than HS-LLA, for negligible gain in statistical performance.  
{\bf b}egin{table}[!htb]
{\bf c}entering
{\bf c}aption{Comparison of results for competing procedures, in a sparse normal means model.}
\label{normal_means_result}
{\bf b}egin{tabular}{|c|ccccc|}
{{\bf b}f h}line
         & HS-LLA           & HS-MCMC                                     & SCAD                                                      & MCP                                                       & lasso                                                     \\ {{\bf b}f h}line
SSE      & {\bf b}egin{tabular}[c]{@{}c@{}}10.051\\ (6.721){{\bf r}m e}nd{tabular}   & {\bf b}egin{tabular}[c]{@{}c@{}}9.757\\ (3.803){{\bf r}m e}nd{tabular} & {\bf b}egin{tabular}[c]{@{}c@{}}26.22\\ (18.629){{\bf r}m e}nd{tabular}  & {\bf b}egin{tabular}[c]{@{}c@{}}32.724\\ (17.794){{\bf r}m e}nd{tabular} & {\bf b}egin{tabular}[c]{@{}c@{}}28.121\\ (18.177){{\bf r}m e}nd{tabular}  \\
pSSE     & {\bf b}egin{tabular}[c]{@{}c@{}}59.238\\ (16.771){{\bf r}m e}nd{tabular}  & {\bf b}egin{tabular}[c]{@{}c@{}}61.381\\ (12.506){{\bf r}m e}nd{tabular}& {\bf b}egin{tabular}[c]{@{}c@{}}75.641\\ (25.609){{\bf r}m e}nd{tabular} & {\bf b}egin{tabular}[c]{@{}c@{}}80.358\\ (22.146){{\bf r}m e}nd{tabular} & {\bf b}egin{tabular}[c]{@{}c@{}}76.082\\ (23.228){{\bf r}m e}nd{tabular} \\
TNR      & {\bf b}egin{tabular}[c]{@{}c@{}}0.982\\ (0.023){{\bf r}m e}nd{tabular}    & {\bf b}egin{tabular}[c]{@{}c@{}}0.992\\ (0.013){{\bf r}m e}nd{tabular}& {\bf b}egin{tabular}[c]{@{}c@{}}0.633\\ (0.336){{\bf r}m e}nd{tabular}   & {\bf b}egin{tabular}[c]{@{}c@{}}0.65\\ (0.325){{\bf r}m e}nd{tabular}    & {\bf b}egin{tabular}[c]{@{}c@{}}0.464\\ (0.366){{\bf r}m e}nd{tabular}   \\
TPR      & {\bf b}egin{tabular}[c]{@{}c@{}}0.64\\ (0.351){{\bf r}m e}nd{tabular}     & {\bf b}egin{tabular}[c]{@{}c@{}}0.53\\ (0.409){{\bf r}m e}nd{tabular}& {\bf b}egin{tabular}[c]{@{}c@{}}0.86\\ (0.268){{\bf r}m e}nd{tabular}    & {\bf b}egin{tabular}[c]{@{}c@{}}0.84\\ (0.31){{\bf r}m e}nd{tabular}     & {\bf b}egin{tabular}[c]{@{}c@{}}0.93\\ (0.202){{\bf r}m e}nd{tabular}     \\
Time (s) & 0.351  & 1.07 & 0.268 & 0.31 & 0.202\\ {{\bf b}f h}line
{{\bf r}m e}nd{tabular}
{{\bf r}m e}nd{table}

{{\bf b}f s}ubsection{Sparse Linear Regression}
\label{linear_regression}
 We set $n=50,\,p=100$ and the components $x$ as $x_1 =3,\,x_2 = 1.5,\,x_3=2,\,x_4,\ldots,x_{100}=0$. The entries in columns of ${\bf P}hi$, where ${\bf P}hi=[{\bf P}hi_1,\ldots,{\bf P}hi_p]$ are generated from a multivariate Gaussian, with zero mean, such that the covariance between ${\bf P}hi_k$ and ${\bf P}hi_l$, for $1\leq k,\,l\leq p$ is $0.5^{|k-l|}$.  The observations $y_j$ are generated as $y_j{{\bf b}f s}im{{\bf b}f m}athcal{N}({\bf P}hi x,1)$, for $j\in\{1,\ldots,n\}$. With the generated data, optimizing~{{\bf r}m e}qref{skeleton_beta_update} with ${{\bf r}m e}ll(x;y) = {{\bf b}f f}rac{1}{2}{{\bf b}f s}um_{i=1}^n(y_i -{\bf P}hi x)^2$, with pathwise coordinate optimization~{\bf c}ite{friedman2007pathwise}, we have :
{\bf b}egin{equation*}
    {\bf b}egin{split}
         x_{(k+1)} =& {\bf u}nderset{x}{{\bf t}ext{argmin}}\left({{\bf b}f f}rac{1}{2}{{\bf b}f s}um_{j=1}^n(y_i-{\bf P}hi x)^2+{{\bf b}f s}um_{i=1}^p\lambda_{i(k)}|x_i|{\bf r}ight),\\
    {\bf t}ext{yielding, }  x_{i(k+1)} &= {{\bf b}f f}rac{1}{{{\bf b}f s}um_{j=1}^n {{\bf b}f p}hi_{ji}^2}S\left({{\bf b}f s}um_{j=1}^n {{\bf b}f p}hi_{ji}\left(y_j - {{\bf b}f w}idetilde{y}_{ji}{\bf r}ight),\lambda_{i(k)}{\bf r}ight),\\
   {\bf t}ext{where},\, {{\bf b}f w}idetilde{y}_{ji}& = {\bf u}nderset{l{{\bf b}f n}eq i}{{{\bf b}f s}um}{{\bf b}f p}hi_{jl}{{\bf b}f w}idetilde{x}_{l}\,,\,{{\bf b}f w}idetilde{x}_{l}={\bf b}egin{cases}x_{l(k+1)},&{\bf t}ext{ if }i>l,\\ x_{l(k)},& {\bf t}ext{ else.}{{\bf r}m e}nd{cases}
    {{\bf r}m e}nd{split}
{{\bf r}m e}nd{equation*}
and the definition of $S({\bf c}dot,\,{\bf c}dot)$ follows from Equation~{{\bf r}m e}qref{soft_threshold}. For the LLA procedure, we start with $x_{i(0)}=0.1,\,{{\bf b}f f}orall  i$ and set the stopping criterion as mentioned in Section~{\bf r}ef{normal_means_model}. We estimate ${{\bf b}f h}at{x}$ for 50 replications, and present the mean and standard deviation (in parentheses) of the performance measures in Table~{\bf r}ef{linear_regression_results}. From these results, we observe that HS-MCMC has the best TNR, followed by HS-LLA and other measures under HS-LLA are comparable to the best performer, which is SCAD for this setting. 
{\bf b}egin{table}[!b]
{\bf c}entering
{\bf c}aption{Comparison of results for competing procedures, in sparse linear regression.}
\label{linear_regression_results}
{\bf b}egin{tabular}{|c|ccccc|}
{{\bf b}f h}line
         & HS-LLA    & HS-MCMC                                                & SCAD                                                      & MCP                                                       & lasso                                                     \\ {{\bf b}f h}line
SSE      & {\bf b}egin{tabular}[c]{@{}c@{}}0.207\\ (0.2){{\bf r}m e}nd{tabular}      & {\bf b}egin{tabular}[c]{@{}c@{}}0.183\\ (0.175){{\bf r}m e}nd{tabular}& {\bf b}egin{tabular}[c]{@{}c@{}}0.118\\ (0.108){{\bf r}m e}nd{tabular}   & {\bf b}egin{tabular}[c]{@{}c@{}}0.182\\ (0.166){{\bf r}m e}nd{tabular}   & {\bf b}egin{tabular}[c]{@{}c@{}}0.281\\ (0.165){{\bf r}m e}nd{tabular}   \\
pSSE     & {\bf b}egin{tabular}[c]{@{}c@{}}59.612\\ (12.918){{\bf r}m e}nd{tabular}  & {\bf b}egin{tabular}[c]{@{}c@{}}58.038\\ (11.45){{\bf r}m e}nd{tabular}& {\bf b}egin{tabular}[c]{@{}c@{}}56.607\\ (11.753){{\bf r}m e}nd{tabular} & {\bf b}egin{tabular}[c]{@{}c@{}}59.233\\ (14.665){{\bf r}m e}nd{tabular} & {\bf b}egin{tabular}[c]{@{}c@{}}65.817\\ (16.688){{\bf r}m e}nd{tabular} \\
TNR      & {\bf b}egin{tabular}[c]{@{}c@{}}0.987\\ (0.019){{\bf r}m e}nd{tabular}    & {\bf b}egin{tabular}[c]{@{}c@{}}1\\ (0){{\bf r}m e}nd{tabular}& {\bf b}egin{tabular}[c]{@{}c@{}}0.975\\ (0.026){{\bf r}m e}nd{tabular}   & {\bf b}egin{tabular}[c]{@{}c@{}}0.975\\ (0.064){{\bf r}m e}nd{tabular}   & {\bf b}egin{tabular}[c]{@{}c@{}}0.92\\ (0.077){{\bf r}m e}nd{tabular}    \\
TPR      & {\bf b}egin{tabular}[c]{@{}c@{}}1\\ (0){{\bf r}m e}nd{tabular}            & {\bf b}egin{tabular}[c]{@{}c@{}}1\\ (0){{\bf r}m e}nd{tabular}& {\bf b}egin{tabular}[c]{@{}c@{}}1\\ (0){{\bf r}m e}nd{tabular}           & {\bf b}egin{tabular}[c]{@{}c@{}}1\\ (0){{\bf r}m e}nd{tabular}           & {\bf b}egin{tabular}[c]{@{}c@{}}1\\ (0){{\bf r}m e}nd{tabular}           \\
Time (s) & 0.26             &   4.23                                      & 0.05                                                      & 0.29                                                      & 0.1                                                       \\ {{\bf b}f h}line
{{\bf r}m e}nd{tabular}
{{\bf r}m e}nd{table}

We again observe from the results presented in Table~{\bf r}ef{linear_regression_results} that the performance measures are comparable for HS-MCMC and HS-LLA. However, HS-MCMC is drastically slower and also requires a choice for the credible interval width for variable selection; which is not the case under HS-LLA. This underlines the scalability of our approach for Bayesian sparse signal recovery.

{{\bf b}f s}ection{Concluding Remarks}
We provide a new Laplace mixture representation of the popular horseshoe prior. While this can be used as an alternative generative hierarchy for the horseshoe, similar to the usual representation of half Cauchy scale mixture of normals, our results show horseshoe prior enjoys some particularly nice properties due to the resultant complete monotonicity of the density function. This is relevant because not all normal scale mixtures are also automatically Laplace mixtures. For example, the class of power exponential densities, $p_{X}(x) {{\bf b}f p}ropto {{\bf r}m e}xp(-|x|^{{\bf b}f a}lpha)$, can be written as a normal scale mixture with respect to a positive ${{\bf b}f a}lpha/2$ stable variable for ${{\bf b}f a}lpha \in (1,2)$ {\bf c}ite{west1987scale}, but these densities are not Laplace mixtures. In addition, our results imply a numerically stable LLA approach for penalized likelihood estimation problem under the horseshoe, and establishes LLA as  special case of EM, thus proving LLA in this case results in the MAP estimate, which is sparse, due to the soft thresholding step in LLA. Deeper studies of these connections should be considered future work.

{{\bf b}f a}ppendix
{{\bf r}m e}mph{The hyperparameter ${\bf t}au$.}  Proposition~{\bf r}ef{prop_equiv_HS_Laplace} is presented for a fixed ${\bf t}au>0$. In both fully Bayesian and penalized likelihood settings, ${\bf t}au$ needs to be tuned for practical implementation. We used cross validation for our numerical results in Section~{\bf r}ef{sec:num}. Fully Bayesian approaches require a hyperprior on ${\bf t}au$, popular among them the half-Cauchy~{\bf c}ite{polson2012half}. Other techniques for choosing ${\bf t}au$ include Akaike information criterion or AIC~{\bf c}ite{lingjaerde2022scalable}, marginal likelihood~{\bf c}ite{bhadra2022graphical} and effective model size~{\bf c}ite{piironen2017sparsity}. For a detailed survey, see {\bf c}ite[Section 5]{bhadra2019lasso}.

{{\bf r}m e}mph{Fast computation of $\lambda_{(k)}$.} Computing the ratio of integrals in Equation~{{\bf r}m e}qref{derivative_peanlty}, at every iteration of the optimization procedure might be perceived as a potential drawback, for two possible reasons: (a) closed form of the integrals do not exist and (b) numerical integration might be a time intensive exercise. Luckily, near-minimax rational approximation of the Dawson function exists~{\bf c}ite{lether1997constrained}, and is given by:
{\bf b}egin{equation*}
    {{\bf b}f f}rac{D_+(u)}{u} {{\bf b}f a}pprox G(u) = {{\bf b}f f}rac{1+{{\bf b}f f}rac{33}{232}u^2+{{\bf b}f f}rac{19}{632}u^4+{{\bf b}f f}rac{23}{1471}u^6}{1+{{\bf b}f f}rac{517}{646}u^2+{{\bf b}f f}rac{58}{173}u^4+{{\bf b}f f}rac{11}{262}u^6+{{\bf b}f f}rac{46}{1471}u^8}{\bf c}dot
{{\bf r}m e}nd{equation*}
{{\bf b}f n}oindent This approximation has a maximum relative error of $6.1{\bf t}imes 10^{-4}$ for $|u|<\infty$, is exact at $u=0$ and has the asymptotic property, $uG(u){{\bf b}f s}im D_+(u)$, as $u{\bf r}ightarrow \infty$~{\bf c}ite{lether1997constrained}. Therefore, having the values of $D_+(u)$ pre-computed and stored on a large enough and fine grid of $u$, a simple Riemann sum can be used to evaluate the expressions in Equation~{{\bf r}m e}qref{derivative_peanlty}. Lipschitz continuity of the rational approximation can also be established with some further algebra, aiding a selection of grid size.

{{\bf b}f I}EEEtriggeratref{17}
{\bf b}ibliographystyle{IEEEtran}
{\bf b}ibliography{ref}

{{\bf r}m e}nd{document}