# Horseshoe shrinkage methods for Bayesian fusion estimation

Sayantan Banerjee 111Corresponding author. Address: IIM Indore, Rau-Pithampur
Road, Indore, M.P. - 453 556, India. e-mail<EMAIL_ADDRESS>Operations
Management and Quantitative Techniques Area
Indian Institute of Management Indore
Indore, M.P., India

###### Abstract

We consider the problem of estimation and structure learning of high
dimensional signals via a normal sequence model, where the underlying
parameter vector is piecewise constant, or has a block structure. We develop a
Bayesian fusion estimation method by using the Horseshoe prior to induce a
strong shrinkage effect on successive differences in the mean parameters,
simultaneously imposing sufficient prior concentration for non-zero values of
the same. The proposed method thus facilitates consistent estimation and
structure recovery of the signal pieces. We provide theoretical justifications
of our approach by deriving posterior convergence rates and establishing
selection consistency under suitable assumptions. We also extend our proposed
method to signal de-noising over arbitrary graphs and develop efficient
computational methods along with providing theoretical guarantees. We
demonstrate the superior performance of the Horseshoe based Bayesian fusion
estimation method through extensive simulations and two real-life examples on
signal de-noising in biological and geophysical applications. We also
demonstrate the estimation performance of our method on a real-world large
network for the graph signal de-noising problem.

> Keywords: Bayesian shrinkage; Fusion estimation; Horseshoe prior; piecewise
> constant function; posterior convergence rate.

##  1 Introduction

With modern technological advances and massive data storage capabilities,
large datasets are becoming increasingly common in a plethora of areas
including finance, econometrics, bioinformatics, engineering, signal-
processing, among others. Flexible modeling of such datasets often require
parameters whose dimension exceed the available sample size. Plausible
inference is possible in such scenarios by finding a lower dimensional
embedding of the high dimensional parameter.

Sparsity plays an important role in statistical learning in a high-dimensional
scenario, where the underlying parameters are ‘nearly black’ (Donoho et al., ,
1992), implying that only a small subset of the same is non-zero. Several
regularization based methods have been proposed in this regard in the
literature, that induce sparsity in the models. From a frequentist
perspective, such methods include penalization based methods like ridge
regression (Hoerl and Kennard, , 1970), lasso (Tibshirani, , 1996), elastic
net (Zou and Hastie, , 2005), SCAD (Fan and Li, , 2001), fused lasso
(Tibshirani et al., , 2005; Rinaldo, , 2009), graphical lasso (Friedman et
al., , 2008), among others; see Bühlmann and Van De Geer,  (2011) for
frequentist methods in high dimensions. Bayesian methods in high dimensions
have been developed more recently, where sparsity is induced via suitable
prior distributions on the parameters. These include spike-and-slab priors
(Mitchell and Beauchamp, , 1988; Ishwaran and Rao, , 2005), Bayesian lasso
(Park and Casella, , 2008), Bayesian graphical lasso (Wang, , 2012), and other
shrinkage priors like spike-and-slab lasso (Ročková and George, , 2018),
Normal Exponential Gamma (Griffin and Brown, , 2010), Horseshoe (Carvalho et
al., , 2009, 2010) and other variants like Horseshoe+ (Bhadra et al., , 2017)
and Horseshoe-like (Bhadra et al., , 2019), Dirichlet-Laplace (Bhattacharya et
al., , 2015), R2-D2 (Zhang et al., , 2016), among others. Theoretical
guarantees of such Bayesian procedures have been developed recently as well.
We refer the readers to Banerjee et al.,  (2021) for a comprehensive overview
of Bayesian methods in a high-dimensional framework.

Piecewise constant signal estimation and structure recovery is an important
and widely studied problem in this regard. Such signals occur in numerous
applications, including bioinformatics, remote sensing, digital image
processing, finance, and geophysics. We refer the readers to Little and Jones,
(2011) for a review on generalized methods for noise removal in piecewise
constant signals, along with potential applications. In this paper, we
consider a noisy high-dimensional signal modeled via the normal sequence model
given by

$y_{i}=\theta_{i}+\epsilon_{i},\;i=1,\ldots,n,$ (1)

where $\epsilon_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\sigma^{2}),$
$\sigma^{2}$ being the noise (or, error) variance,
$y=(y_{1},\ldots,y_{n})^{T}$ is the observed noisy signal, and
$\theta=(\theta_{1},\ldots,\theta_{n})^{T}\in\mathbb{R}^{n}$ is the vector of
mean parameters. We further assume that the true parameter
$\theta_{0}=(\theta_{0,1},\ldots,\theta_{0,n})^{T}\in\mathbb{R}^{n}$ is
piecewise constant, or has an underlying block structure, in the sense that
the transformed parameters $\eta_{0,j}=\theta_{0,j}-\theta_{0,j-1},\,2\leq
j\leq n$ belong to a nearly-black class
$l_{0}[s]=\\{\eta\in\mathbb{R}^{n-1}:\\#\\{j:\eta_{j}\neq 0,\,2\leq j\leq
n\\}\leq s,\,0\leq s=s(n)\leq n\\}$. Here $\\#$ denotes the cardinality of a
finite set. The number of true blocks in the parameter vector is
$s_{0}:=\\#\\{j:\eta_{0,j}\neq 0,\,2\leq j\leq n\\}.$ Our fundamental goals
are two-fold – (i) estimating the mean parameter $\theta$, and (ii) recovering
the underlying piecewise constant structure.

Frequentist methods in addressing the piecewise constant signal de-noising
problem include penalization methods like the fused lasso method (Tibshirani
et al., , 2005) and the $L_{1}$-fusion method (Rinaldo, , 2009). The Bayesian
equivalent of the fused lasso penalty is using a Laplace shrinkage prior on
the successive differences $\eta$ (Kyung et al., , 2010). However, the Laplace
prior leads to posterior inconsistency (Castillo et al., , 2015). To overcome
this problem, Song and Liang,  (2017) used a heavy-tailed shrinkage prior for
the coefficients in a linear regression framework. Motivated by this, Song and
Cheng,  (2020) proposed to use a $t$-shrinkage prior for Bayesian fusion
estimation. Shimamura et al.,  (2019) use a Normal Exponential Gamma (NEG)
prior in this context, and make inference based on the posterior mode.

Motivated by strong theoretical guarantees regarding estimation and structure
learning induced by the Horseshoe prior in a regression model (Datta and
Ghosh, , 2013; Van Der Pas et al., , 2014; van der Pas et al., , 2017), we
propose to investigate the performance of the same in our fusion estimation
framework. The Horseshoe prior induces sparsity via an infinite spike at zero,
and also possess a heavy tail to ensure consistent selection of the underlying
pieces or blocks. Furthermore, the global scale parameter that controls the
level of sparsity in the model automatically adapts to the actual level of
sparsity in the true model, as opposed to choosing the scale parameter in the
$t$-shrinkage prior on the basis of the underlying dimension $n$. The
conjugate structure of the posterior distribution arising out of modeling via
a Horseshoe prior allows fast computation along with working with the full
posterior distribution, so that the samples obtained via Markov Chain Monte
Carlo (MCMC) can further be used for uncertainty quantification.

The piecewise constant normal sequence model can be further extended to
piecewise constant signals over arbitrary graphs. Note that the former is a
special case of signals over graphs, when the underlying graph is a linear
chain. We consider the normal sequence model as given by (1), but in this
case, the signal $\theta=(\theta_{1},\ldots,\theta_{n})^{T}$ is defined over
an undirected graph $G=(V,E)$, where $V=\\{1,\ldots,n\\}$ is the vertex set of
the graph and $E$ is the corresponding edge-set. The true signal is then
assumed to be piecewise constant over the graph $G$, so that the parameter
space is now given by $l_{0}[G,s]:=\\{\theta\in\mathbb{R}^{n}:\\#\\{(i,j)\in
E:\theta_{i}\neq\theta_{j}\\}\leq s,\,0\leq s(n)=s\leq n\\}$. The resulting
graph de-noising problem has several potential applications, including
multiple change-point detection in linear chain graphs, image segmentation,
and anomaly detection in large networks. We refer the readers to Fan and Guan,
(2018) and references therein for related work in this field. In this work, we
propose a Bayesian fusion method via Horseshoe prior specification on adjacent
edges of a transformed graph obtained from the original graph and illustrate
its excellent inference and posterior concentration properties.

The paper is organized as follows. In the next section, we specify the
Bayesian model along with the prior specifications, followed by posterior
computations in Section 3. In Section 4, we provide theoretical guarantees of
our proposed method via determining posterior convergence rates and
establishing posterior selection consistency of the signal pieces. We
demonstrate the numerical performance of our method along with other competing
methods via simulations in Section 5, followed by real-life applications in
two different areas – DNA copy number analysis and Solar X-Ray flux data
analysis. The methodological extension of our proposed method to the graph
signal de-noising problem, along with theoretical guarantees and real-life
illustration are given in Section 7. We conclude our paper with a brief
discussion of our proposed fusion estimation method along with future
directions for research. Additional lemmas, proofs of main results, and a
practical solution for block structure recovery along with numerical results
for the same are provided in the Appendix.

The notations used in the paper are as follows. For real-valued sequences
$\\{a_{n}\\}$ and $\\{b_{n}\\}$, $a_{n}=O(b_{n})$ implies that $a_{n}/b_{n}$
is bounded, and $a_{n}=o(b_{n})$ implies that $a_{n}/b_{n}\rightarrow 0$ as
$n\rightarrow\infty.$ By $a_{n}\lesssim b_{n}$, we mean that $a_{n}=O(b_{n})$,
while $a_{n}\asymp b_{n}$ means that both $a_{n}\lesssim b_{n}$ and
$b_{n}\lesssim a_{n}$ hold. $a_{n}\prec b_{n}$ means $a_{n}=o(b_{n}).$ For a
real vector $x=(x_{1},\ldots,x_{n})^{T}$, the $L_{r}$-norm of $x$ for $r>0$ is
defined as $\|x\|_{r}=\left(\sum_{i=1}^{n}|x_{i}|^{r}\right)^{1/r}.$ We denote
the cardinality of a finite set $S$ as $\\#S$. The indicator function is
denoted by $1\\!\mathrm{l}.$

##  2 Bayesian modeling and prior specification

We consider the normal sequence model (1) and assume that the successive
differences of the means belong to a nearly-black class $l_{0}[s]$. As
discussed earlier, frequentist procedures induce sparsity in the model via
suitable regularization based procedures like penalization of the underlying
parameters, whereas Bayesian methods usually induce sparsity via imposing
suitable prior distributions on the same. For example, for learning a high-
dimensional parameter $\theta\in\mathbb{R}^{n}$, a general version of a
penalized optimization procedure can be written as
$\arg\min_{\theta\in\mathbb{R}^{n}}\\{l(\theta;y)+\pi(\theta)\\},$ where
$l(\theta;y)$ is the empirical risk and $\pi(\theta)$ is the penalty function.
If $l(\theta;y)$ is defined as the negative log-likelihood (upto a constant)
of the observations $y$, the above optimization problem becomes equivalent to
finding the mode of the posterior distribution $p(\theta\mid y)$,
corresponding to the prior density
$p(\theta)\propto\exp\left(-\pi(\theta)\right).$

In our context, Tibshirani et al.,  (2005) proposed the _fused lasso_
estimator $\hat{\theta}^{FL}$ that induces sparsity on both $\theta$ and
$\eta$, defined as

$\hat{\theta}^{FL}=\arg\min_{\theta\in\mathbb{R}^{n}}\left\\{\dfrac{1}{2}\|y-\theta\|_{2}^{2}+\lambda_{1}\|\theta\|_{1}+\lambda_{2}\|\eta\|_{1}\right\\},$
(2)

for suitable penalty parameters $\lambda_{1}$ and $\lambda_{2}$. Rinaldo,
(2009) considered the $L_{1}$-fusion estimator $\hat{\theta}^{F}$ with
penalization of the successive differences only, given by,

$\hat{\theta}^{F}=\arg\min_{\theta\in\mathbb{R}^{n}}\left\\{\dfrac{1}{2}\|y-\theta\|_{2}^{2}+\lambda\|\eta\|_{1}\right\\},$
(3)

where $\lambda$ is the corresponding penalty (tuning) parameter. A Bayesian
equivalent of the fused lasso estimator (3) can be obtained by putting a
$Laplace(\lambda/\sigma)$ prior on the successive differences $\eta$ and
finding the corresponding posterior mode. As in a normal regression model with
Bayesian lasso (Park and Casella, , 2008), the Bayesian fused lasso estimator
given by the posterior mode will converge to the true $\eta_{0}$ at a nearly
optimal rate. However, the induced posterior distribution has sub-optimal
contraction rate (Castillo et al., , 2015), owing to insufficient prior
concentration near zero.

The posterior inconsistency of the Bayesian fused lasso method motivates us to
explore other approaches that would mitigate the problems leading to the
undesirable behavior of the posterior distribution. Shrinkage priors qualify
as naturally good choices for our problem, as they can address the dual issue
of shrinking true zero parameters to zero, and retrieving the ‘boundaries’ of
the blocks effectively by placing sufficient mass on the non-zero values of
successive differences of the normal means. In the normal sequence model,
optimal posterior concentration has been achieved via using shrinkage prior
distributions with polynomially decaying tails (Song and Liang, , 2017). This
is in contrast to the Laplace prior, that has exponentially light tails.

The Horseshoe prior is a widely acclaimed choice as a shrinkage prior owing to
its infinite spike at 0 and simultaneously possessing a thick tail. The tails
decay like a second-order polynomial, and hence the corresponding penalty
function behaves like a logarithmic penalty, and is non-convex (see Carvalho
et al.,  (2010); Bhadra et al.,  (2019) for more details). The Horseshoe prior
can be expressed as a scale mixture of normals with half-Cauchy prior, thus
acting as a global-local shrinkage prior. We put a Horseshoe prior on the
pairwise differences in the parameters
$\eta_{i}=\theta_{i}-\theta_{i-1},\;i=2,\ldots,n.$ We also need to put
suitable priors on the mean parameter $\theta_{1}$ and the error variance
$\sigma^{2}$. The full prior specification is given by,

$\displaystyle\theta_{1}\mid\lambda_{1}^{2},\sigma^{2}\sim
N(0,\lambda_{1}^{2}\sigma^{2}),\;\eta_{i}\mid\lambda_{i}^{2},\tau^{2},\sigma^{2}\stackrel{{\scriptstyle
ind}}{{\sim}}N(0,\lambda_{i}^{2}\tau^{2}\sigma^{2}),\;2\leq i\leq n,$
$\displaystyle\lambda_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}C^{+}(0,1),\;2\leq i\leq n,\;\tau\sim C^{+}(0,1),\;\sigma^{2}\sim
IG(a_{\sigma},b_{\sigma}).$ (4)

Here $C^{+}(\cdot,\cdot)$ and $IG(\cdot,\cdot)$ respectively denote the half
Cauchy density and Inverse Gamma density. The level of sparsity induced in the
model is controlled by the global scale parameter $\tau$, and choosing the
same is a non-trivial problem. Using a plug-in estimate for $\tau$ based on
empirical Bayes method suffers from a potential danger of resulting in a
degenerate solution resulting in a heavily sparse model. There are several
works (Carvalho et al., , 2010; Piironen and Vehtari, 2017b, ; Piironen and
Vehtari, 2017a, ) that suggest effective methods regarding the choice of
$\tau$. In this paper, we have proposed to take a fully Bayesian approach as
suggested in Carvalho et al.,  (2010); Piironen and Vehtari, 2017b and use a
half-Cauchy prior.

The half-Cauchy distribution can further be written as a scale-mixture of
Inverse-Gamma distributions. For a random variable $X\sim C^{+}(0,\psi)$, we
can write,

$X^{2}\mid\phi\sim IG(1/2,1/\phi),\;\phi\sim IG(1/2,1/\psi^{2}).$

Thus, the full hierarchical prior specification for our model is given by,

$\displaystyle\theta_{1}\mid\lambda_{1}^{2},\sigma^{2}\sim
N(0,\lambda_{1}^{2}\sigma^{2}),\;\eta_{i}\mid\lambda_{i}^{2},\tau^{2},\sigma^{2}\stackrel{{\scriptstyle
ind}}{{\sim}}N(0,\lambda_{i}^{2}\tau^{2}\sigma^{2}),\;2\leq i\leq n,$
$\displaystyle\lambda_{i}^{2}\mid\nu_{i}\stackrel{{\scriptstyle
ind}}{{\sim}}IG(1/2,1/\nu_{i}),\;2\leq i\leq n,\;\tau^{2}\mid\xi\sim
IG(1/2,1/\xi),$
$\displaystyle\nu_{2},\ldots,\nu_{n},\xi\stackrel{{\scriptstyle
iid}}{{\sim}}IG(1/2,1),\;\sigma^{2}\sim IG(a_{\sigma},b_{\sigma})$ (5)

The hyperparameters $a_{\sigma}$ and $b_{\sigma}$ may be chosen in such a way
that the corresponding prior becomes non-informative. The local scale
parameter $\lambda_{1}$ is considered to be fixed as well.

##  3 Posterior computation

The conditional posterior distributions of the underlying parameters can be
explicitly derived via exploring the conjugate structure. Hence, the posterior
computations can be accomplished easily via Gibbs sampling. We present the
conditional posterior distributions of the parameters below.

The normal means have the conditional posterior distribution

$\theta_{i}\mid\cdots\sim N(\mu_{i},\zeta_{i}),\;1\leq i\leq n,$ (6)

where $\mu_{i}$ and $\zeta_{i}$ are given by,

$\zeta_{1}^{-1}=\dfrac{1}{\sigma^{2}}\left(1+\dfrac{1}{\lambda_{2}^{2}\tau^{2}}+\dfrac{1}{\lambda_{1}^{2}}\right),\;\mu_{1}=\dfrac{\zeta_{1}}{\sigma^{2}}\left(y_{1}+\dfrac{\theta_{2}}{\lambda_{2}^{2}\tau^{2}}\right),$
$\zeta_{i}^{-1}=\dfrac{1}{\sigma^{2}}\left(1+\dfrac{1}{\lambda_{i+1}^{2}\tau^{2}}+\dfrac{1}{\lambda_{i}^{2}\tau^{2}}\right),\;\mu_{i}=\dfrac{\zeta_{i}}{\sigma^{2}}\left(y_{i}+\dfrac{\theta_{i+1}}{\lambda_{i+1}^{2}\tau^{2}}+\dfrac{\theta_{i-1}}{\lambda_{i}^{2}\tau^{2}}\right),\;2\leq
i\leq n.$

Here $\lambda_{n+1}$ is considered to be infinity. The conditional posteriors
for the rest of the parameters are given by,

$\displaystyle\lambda_{i}^{2}\mid\cdots$ $\displaystyle\sim$ $\displaystyle
IG\left(1,\dfrac{1}{\nu_{i}}+\dfrac{(\theta_{i}-\theta_{i-1})^{2}}{2\tau^{2}\sigma^{2}}\right),\;2\leq
i\leq n,$ $\displaystyle\sigma^{2}\mid\cdots$ $\displaystyle\sim$
$\displaystyle
IG\left(n+a_{\sigma},b_{\sigma}+\dfrac{1}{2}\left[\sum_{i=1}^{n}(y_{i}-\theta_{i})^{2}+\dfrac{1}{\tau^{2}}\sum_{i=2}^{n}\dfrac{(\theta_{i}-\theta_{i-1})^{2}}{\lambda_{i}^{2}}+\dfrac{\theta_{1}^{2}}{\lambda_{1}^{2}}\right]\right),$
$\displaystyle\tau^{2}\mid\cdots$ $\displaystyle\sim$ $\displaystyle
IG\left(\dfrac{n}{2},\dfrac{1}{\xi}+\dfrac{1}{2\sigma^{2}}\sum_{i=2}^{n}\dfrac{(\theta_{i}-\theta_{i-1})^{2}}{\lambda_{i}^{2}}\right),$
$\displaystyle\nu_{i}\mid\cdots$ $\displaystyle\sim$ $\displaystyle
IG\left(1,1+\dfrac{1}{\lambda_{i}^{2}}\right),\;2\leq i\leq n,$
$\displaystyle\xi\mid\cdots$ $\displaystyle\sim$ $\displaystyle
IG\left(1,1+\dfrac{1}{\tau^{2}}\right).$ (7)

##  4 Theoretical results

In this section, we present the theoretical validations of using a Horseshoe
prior for fusion estimation. We first discuss the result involving posterior
convergence rate of the mean parameter $\theta$ under certain assumptions and
then proceed to discuss the result on posterior selection of the underlying
true block structure of the mean parameter.

######  Assumption 1.

The number of true blocks in the model satisfies $s_{0}\prec n/\log n$.

######  Assumption 2.

The true mean parameter vector
$\theta_{0}=(\theta_{0,1},\ldots,\theta_{0,n})^{T}$ and the true error
variance $\sigma_{0}^{2}$ satisfy the following conditions:

  1. (i)

Define $\eta_{0,j}=\theta_{0,j}-\theta_{0,j-1},\,2\leq j\leq n.$ Then,
$\max_{j}|\eta_{0,j}/\sigma_{0}|<L,$ where $\log L=O(\log n).$

  2. (ii)

$\theta_{0,1}/(\lambda_{1}^{2}\sigma_{0}^{2})+2\log\lambda_{1}=O(\log n),$
where $\lambda_{1}$ is the prior hyperparameter appearing in the prior for
$\theta_{1}$ in (4).

######  Assumption 3.

The global scale parameter $\tau$ in the prior specification (4) satisfies
$\tau<n^{-(2+b)}$ for some constant $b>0$, and $-\log\tau=O(\log n)$.

The first assumption above involves the true block size, that is ubiquitous in
structure recovery problems in high-dimensions. In Assumption 2, we have
considered $L$ as a bound on the maximum value of $|\eta_{0,j}/\sigma_{0}|$.
Such a restriction is necessary for block recovery at a desired contraction
rate. We shall see (in Lemma 16) that the aforesaid condition on $L$, along
with one of the conditions on the global scale parameter $\tau$ as in
Assumption 3 would guarantee that the tail of the prior is not too sharp. This
would ensure a minimum prior concentration for large non-zero values of the
successive difference in the means. An equivalent prior condition cannot hold
uniformly over an unbounded parameter space, thus prompting a restriction on
the successive differences, and also on $\theta_{0,1}$. Similar restrictions
in the context of linear and generalized regression models have been
considered in the literature; for example, see Wei and Ghosal,  (2020); Song,
(2020). Assumption 3 involves an upper bound on the global scale parameter
$\tau$, that would be necessary for ensuring that the prior puts sufficient
mass around zero for the successive differences in means, thus leading to
effective fusion estimation via imposing sparsity. However, $\tau$ should not
be too small, otherwise it would lead to degeneracy in inference by picking up
too sparse a model. The lower bound on $\tau$, along with the condition on $L$
in Assumption 2 guarantees that the Horseshoe prior is ‘thick’ at non-zero
parameter values, so that it is not too sharp. We now present our main result
on posterior contraction rate.

######  Theorem 4.

Consider the Gaussian means model (1) with prior specification as in (4), and
suppose that assumptions 1, 2 and 3 hold. Then the posterior distribution of
$\theta$, given by $\Pi^{n}(\cdot\mid y)$, satisfies

$\Pi^{n}(\|\theta-\theta_{0}\|_{2}/\sqrt{n}\geq M\sigma_{0}\epsilon_{n}\mid
y)\rightarrow 0,\;\mathrm{as}\,n\rightarrow\infty,$ (8)

in probability or in $L_{1}$ wrt the probability measure of $y$, for
$\epsilon_{n}\asymp\sqrt{s_{0}\log n/n}$ and a constant $M>0$.

The above result implies that the posterior convergence rate for
$\|\theta-\theta_{0}\|_{2}/\sqrt{n}$ is of the order
$\sigma_{0}\sqrt{s_{0}\log n/n}$. When the exact piecewise constant structure
is known, the proposed Bayesian fusion estimation method achieves the optimal
convergence rate $O(\sigma_{0}\sqrt{s/n})$ upto a logarithmic term in $n$. The
posterior convergence rate also adapts to the (unknown) size of the pieces.
The rate directly compares with the Bayesian fusion estimation method as
proposed in Song and Cheng,  (2020).

The Horseshoe prior (and other global-local shrinkage priors) are continuous
shrinkage priors, and hence exact block structure recovery is not possible.
However, we can consider discretization of the posterior samples via the
posterior projection of the samples $(\theta,\sigma)$ to a discrete set
$S(\theta,\sigma)=\\{2\leq j\leq
n:|\theta_{j}-\theta_{j-1}|/\sigma<\epsilon_{n}/n\\}.$ The number of false-
positives resulting from such a discretization can be expressed as the
cardinality of the set $A(\theta,\sigma)=S^{c}(\theta,\sigma)-\\{2\leq j\leq
n:\theta_{0,j}-\theta_{0,j-1}\neq 0\\}.$ The induced posterior distribution of
$S(\theta,\sigma)$ (and hence that of $A(\theta,\sigma)$) can be shown to be
‘selection consistent’, in the sense that the number of false-positives as
defined above is bounded in probability. We formally present this result
below.

######  Theorem 5.

Under the assumptions of Theorem 4, the posterior distribution of
$A(\theta,\sigma)$ satisfies

$\Pi^{n}(\\#A(\theta,\sigma)>Ks_{0}\mid y)\rightarrow 0,$ (9)

in probability or in $L_{1}$ wrt the measure of $y$ for some fixed constant
$K>0$.

Note that the thresholding rule depends on the posterior contraction rate,
that involves the number $(s_{0})$ of true blocks. However, in practical
scenarios, knowledge of $s_{0}$ may not be readily available. To tackle such
situations, we propose to use a threshold analogous to the concept of
shrinkage weights as in the sparse normal sequence model. We outline the
details in the appendix.

##  5 Simulation studies

To evaluate the performance of our method and compare with other competing
approaches, we carry out simulation studies for varying signal and noise
levels. Three different types of piecewise constant functions of length
$n=100$ are considered – (i) 10 evenly spaced signal pieces, with each of the
pieces having length 10, (ii) 10 unevenly spaced signal pieces, with the
shorter pieces each having length 5, and (iii) 10 very unevenly spaced signal
pieces, with the shorter pieces each having length 2. The response variable
$y$ is generated from an $n$-dimensional Gaussian distribution with mean
$\theta_{0}$ and variance $\sigma^{2}$, where $\theta_{0}$ is the true signal
vector, and $\sigma\in\\{0.1,0.3,0.5\\}.$

We estimate the true signal using our Horseshoe prior based fusion estimation
approach, and compare with three other approaches – (a) Bayesian $t$-fusion
method, as proposed in Song and Cheng,  (2020), (b) Bayesian fusion approach
based on a Laplace prior, and (c) the $L_{1}$-fusion method. For the Bayesian
methods, we consider 5000 MCMC samples, with initial 500 samples as burn-in.
The MCMC-details including Gibbs sampler updates for the $t$-fusion and
Laplace fusion approaches, and the choice of the corresponding scale
parameters for the above priors are taken as suggested in Song and Cheng,
(2020). The hyperparameter values for the prior on error variance are taken as
$a_{\sigma}=b_{\sigma}=0.5$ across all the Bayesian methods, with the local
scale parameter $\lambda_{1}=5$. For the frequentist fused lasso method based
on $L_{1}$-fusion penalty, we use the genlasso package in R, and choose the
fusion penalty parameter using a 5-fold cross-validation approach. To evaluate
the performance of the estimators, we use the Mean Squared Error (MSE) and the
adjusted MSE, respectively defined as,

$\displaystyle\mathrm{MSE}$ $\displaystyle=$
$\displaystyle\|\hat{\theta}-\theta_{0}\|_{2}^{2}/n,$
$\displaystyle\mathrm{adj.MSE}$ $\displaystyle=$
$\displaystyle\|\hat{\theta}-\theta_{0}\|_{2}^{2}/\|\theta_{0}\|_{2}^{2},$

where $\hat{\theta}$ is the estimated signal, given by the posterior mean in
case of the Bayesian methods and the fused lasso estimate for the frequentist
method. All the computations were performed in R on a laptop having an
Intel(R) Core(TM) i7-10750H CPU @ 2.60GHz with 16GB RAM and a 64-bit OS,
x64-based processor. The R codes to implement our proposed method are
available at https://github.com/sayantanbiimi/HSfusion.

The summary measures for the performance of the four methods using MSE,
adjusted MSE and their associated standard errors (in brackets) based on 100
Monte-Carlo replications are presented in Table 1. We find that our proposed
Horseshoe based fusion estimation method has excellent signal estimation
performance across all the different types of signal pieces and noise levels.
Within a Bayesian framework, the performance of the Laplace based fusion
estimation method is not at all promising. Though the estimation performances
for the Horseshoe-fusion, $t$-fusion and fused lasso are comparable in case of
low noise levels, our proposed method performs much better in the higher noise
scenario. Additionally, the Bayesian methods provide credible bands as well
that can be utilized further for uncertainty quantification. In that regard,
we observe from Figures 1, 2 and 3 that the Horseshoe based estimates have
comparatively narrower credible bands as compared to other Bayesian competing
methods considered here. Overall, we could successfully demonstrate the
superiority of using a heavy-tailed prior for successive differences in the
signals in contrast to an exponentially lighter one. In addition to that,
using a prior with a comparatively sharper spike at zero results in better
signal and structure recovery, especially in situations where the noise level
is higher.

|  | HS-fusion |  $t$-fusion | Laplace fusion |  $L_{1}$ fusion
---|---|---|---|---|---
Signal | $\sigma$ | MSE | adj MSE | MSE | adj MSE | MSE | adj MSE | MSE | adj MSE
| 0.1 | 0.002 (0.000) | 0.000 (0.000) | 0.002 (0.000) | 0.001 (0.000) | 0.504 (0.002) | 0.114 (0.004) | 0.005 (0.003) | 0.001 (0.001)
Even | 0.3 | 0.020 (0.001) | 0.004 (0.000) | 0.030 (0.002) | 0.007 (0.001) | 0.516 (0.002) | 0.116 (0.004) | 0.035 (0.019) | 0.008 (0.004)
| 0.5 | 0.062 (0.002) | 0.014 (0.000) | 0.133 (0.009) | 0.030 (0.002) | 0.539 (0.003) | 0.121 (0.004) | 0.087 (0.039) | 0.019 (0.009)
| 0.1 | 0.002 (0.000) | 0.001 (0.000) | 0.002 (0.000) | 0.001 (0.000) | 0.566 (0.002) | 0.249 (0.003) | 0.005 (0.003) | 0.002 (0.001)
Uneven | 0.3 | 0.019 (0.001) | 0.008 (0.000) | 0.036 (0.002) | 0.016 (0.001) | 0.575 (0.002) | 0.253 (0.004) | 0.037 (0.019) | 0.016 (0.008)
| 0.5 | 0.060 (0.002) | 0.026 (0.001) | 0.212 (0.012) | 0.094 (0.006) | 0.596 (0.003) | 0.263 (0.004) | 0.090 (0.044) | 0.039 (0.019)
| 0.1 | 0.002 (0.000) | 0.002 (0.000) | 0.005 (0.000) | 0.005 (0.000) | 0.479 (0.001) | 0.520 (0.001) | 0.009 (0.002) | 0.010 (0.002)
V. Uneven | 0.3 | 0.020 (0.001) | 0.022 (0.001) | 0.114 (0.007) | 0.124 (0.007) | 0.489 (0.001) | 0.531 (0.001) | 0.066 (0.024) | 0.072 (0.026)
| 0.5 | 0.064 (0.002) | 0.070 (0.002) | 0.581 (0.020) | 0.631 (0.021) | 0.507 (0.002) | 0.551 (0.002) | 0.143 (0.057) | 0.155 (0.062)
Table 1: MSE and adjusted MSE values (with associated standard errors in
parentheses) for Horseshoe-fusion, $t$-fusion, Laplace fusion, and fused lasso
method, when the true signal is evenly spaced (“Even”), unevenly spaced
(“Uneven”), and very unevenly spaced (“V. Uneven”).

| Evenly spaced pieces | Unevenly spaced pieces | Very unevenly spaced pieces
---|---|---|---
True |  |  |
|  |  |
HS |  |  |
|  |  |
$t$ |  |  |
|  |  |
Laplace |  |  |
|  |  |
Fused |  |  |
|  |  |
Figure 1: Fusion estimation performance with differently spaced signals and error sd $\sigma=0.1.$ Observations are represented in red dots, true signals in blue dots, point estimates in black dots, and 95% credible bands of the Bayesian procedures in green. | Evenly spaced pieces | Unevenly spaced pieces | Very unevenly spaced pieces
---|---|---|---
True |  |  |
|  |  |
HS |  |  |
|  |  |
$t$ |  |  |
|  |  |
Laplace |  |  |
|  |  |
Fused |  |  |
|  |  |
Figure 2: Fusion estimation performance with differently spaced signals and error sd $\sigma=0.3.$ Observations are represented in red dots, true signals in blue dots, point estimates in black dots, and 95% credible bands of the Bayesian procedures in green. | Evenly spaced pieces | Unevenly spaced pieces | Very unevenly spaced pieces
---|---|---|---
True |  |  |
|  |  |
HS |  |  |
|  |  |
$t$ |  |  |
|  |  |
Laplace |  |  |
|  |  |
Fused |  |  |
|  |  |
Figure 3: Fusion estimation performance with differently spaced signals and
error sd $\sigma=0.5.$ Observations are represented in red dots, true signals
in blue dots, point estimates in black dots, and 95% credible bands of the
Bayesian procedures in green.

##  6 Real-data examples

We check the performance of our proposed method on two real-world data
examples where the observed data can be modeled using piecewise constant
functions. The first example involves DNA copy-number analysis of Glioblastoma
multiforme (GBM) data, and the second one involves analysis of Solar X-Ray
flux data for the period of October-November 2003 that encompasses the
Halloween 2003 Solar storm events.

###  6.1 DNA copy-number analysis of Glioblastoma multiforme (GBM)

Array Comparative Genomic Hybridization (aCGH) is a high-throughput technique
used to identify chromosomal abnormalities in the genomic DNA. Identification
and characterization of chromosomal aberrations are extremely important for
pathogenesis of various diseases, including cancer. In cancerous cells, genes
in a particular chromosome may get amplified or deleted owing to mutations,
thus resulting in gains or losses in DNA copies of the genes. Array CGH data
measure the $\log_{2}$-ratio between the DNA copy number of genes in tumor
cells and those in control cells. The aCGH data thus may be considered as a
piecewise constant signal with certain non-zero blocks corresponding to
aberrations owing to additions/deletions at the DNA level.

Glioblastoma multiforme (GBM) is the most common malignant brain tumors in
adults (Holland, , 2000). Lai et al.,  (2005) analyzed array CGH data in GBM
using glioma data from Bredel et al.,  (2005) along with control samples.
Tibshirani and Wang,  (2008) used the data from Lai et al.,  (2005), and
created aCGH data corresponding to a pseudo-chromosome so that the resulting
genome sequence have shorter regions of high amplitude gains (copy number
additions) as well as a broad region of low amplitude loss (copy number
deletions). To be precise, the combined array consists of aCGH data from
chromosome 7 in GBM29 and chromosome 13 in GBM31. The aCGH data is presented
in Figure 4, and is available from the archived version of the cghFlasso
package in R.

We apply the Horseshoe prior based Bayesian fusion estimation method and also
compare the performance with the fused lasso method. We post-process the
Bayesian estimates by discretization via the thresholding method as outlined
in the appendix. The results are displayed in Figure 5. We find that our
proposed Bayesian fusion estimation method has successfully detected the high
amplitude regions of copy number additions, as well as the broad low amplitude
region of copy number deletions. The results are very much comparable with the
fused lasso estimates, and the findings are similar to that in Lai et al.,
(2005). The Bayesian credible intervals obtained using our proposed method may
further be utilized in addressing future research problems, like controlling
false discovery rates.

Figure 4: Figure showing array CGH data for Chromosome 7 in GBM29 and
chromosome 13 in GBM31.
---
(a) Bayesian Horseshoe fusion method
(b) $L_{1}$ fused lasso method
Figure 5: Figure showing estimation performances of (a) the proposed Bayesian
Horseshoe prior based fusion method, and (b) the $L_{1}$ fused lasso method.
The blue lines represent the estimated signals, and the red dots represent the
$\log_{2}$ ratios between the DNA copy numbers of genes in tumor cells and
those in control cells.

###  6.2 Solar X-Ray flux analysis: Halloween Solar Storms 2003

The Earth’s ionosphere is responsible for blocking high frequency radio waves
on the side facing the Sun. Large solar flares affect the ionosphere, thus
interrupting satellite communications, and posing threats of radiation hazards
to astronauts and spacecrafts. These solar flares are also associated with
Coronal Mass Ejections (CMEs) that have the potential to trigger geomagnetic
storms, which in turn, have the ability to disrupt satellite and radio
communications, cause massive black-outs via electric grid failures, and
affect global navigation systems. The Geostationary Operational Environmental
Satellite (GOES), operated by the United States’ National Oceanic and
Atmospheric Administration (NOAA), records data on solar flares. Measured in
Watts-per-square-meter ($W/m^{2}$) unit, solar X-Ray flux can be classified in
five categories, namely, A, B, C, M, and X, depending on the approximate peak
flux with wavelengths in the range of $1-8$ Angstroms (long wave) and $0.5-4$
Angstroms (short wave). Solar flares in category A have the least peak flux
($<10^{-7}\,W/m^{2}$), and those in caregory X are the most powerful ones
($>10^{-4}\,W/m^{2}$).

Around Halloween in 2003, a series of strong solar storms affected our planet,
leading to disruptions in satellite communication, damages to spacecrafts, and
massive power outages in Sweden. The increased levels of solar activity made
it possible to observe the Northern lights (Aurora Borealis) even from far
south latitudes in Texas, USA. We consider the short wave solar X-ray flux
data for the two months of October and November 2003. The data have been
accessed from NOAA’s website
(https://www.ngdc.noaa.gov/stp/satellite/goes/dataaccess.html). We use the
5-minute average flux data for all the days in the two months, and extract the
information for the flux level for shortwave X-rays.

The data contain some missing values, which we impute via a linear
interpolation method. Also, we further use 1-hour averages of the X-ray flux
data, so as to arrive at 24 data points for each day, thus leading to a total
of 1464 data points. Figure 6 shows the plot of the logarithm of the flux
levels for the two months. As mentioned in Little and Jones,  (2011), we can
approximate this data as a piecewise constant signal, and use fusion-based
estimation methods to identify days of extreme solar activities. We apply our
Horseshoe prior based Bayesian fusion estimation method, and post-process the
estimates by using the thresholding approach as mentioned in the appendix. We
also apply the $L_{1}$-fusion estimation method. The fitted estimates are
shown in Figure 7. We find that our Horseshoe prior based Bayesian approach
has been able to detect the periods of increased solar activity around October
28 and also on November 04, 2003, besides giving piecewise constant estimates.
Our findings match with related analyses; see Pulkkinen et al.,  (2005). On
the other hand, the fused lasso method suffers from serious over-fitting
problems.

Figure 6: Figure showing scatter plot of logarithm of solar X-ray flare flux
values for the months of October and November 2003, recorded by GOES.
---
(a) Bayesian Horseshoe fusion method
(b) $L_{1}$ fused lasso method
Figure 7: Figure showing estimation performances of (a) the proposed Bayesian
Horseshoe prior based fusion method, and (b) the $L_{1}$ fused lasso method.
The blue lines represent the estimated signals, and the red dots represent the
logarithm of the shortwave X-ray solar flare flux values as recorded by GOES.

##  7 Signal de-noising over arbitrary graphs

In this section, we present an extension of our method for signal de-noising
over arbitrary graphs. As mentioned earlier, we consider the normal sequence
model as given by (1), where the signal
$\theta=(\theta_{1},\ldots,\theta_{n})^{T}$ is defined over an undirected
graph $G=(V,E)$. The true signal is then assumed to be piecewise constant over
$G$, and the parameter space is given by
$l_{0}[G,s]:=\\{\theta\in\mathbb{R}^{n}:\\#\\{(i,j)\in
E:\theta_{i}\neq\theta_{j}\\}\leq s,\,0\leq s(n)=s\leq n\\}$.

Motivated by the fused lasso approach of Padilla et al.,  (2017) to deal with
the graph de-noising problem, Banerjee and Shen,  (2020) recently proposed a
Bayesian fusion estimation method using $t$-shrinkage priors. In these two
works, the actual graph $G$ is first reduced to a linear chain graph via the
Depth-First-Search algorithm (Tarjan, , 1972). Padilla et al.,  (2017)
provided theoretical justifications of reducing an arbitrary graph $G$ to a
DFS-chain graph $G_{C}=(V,E_{C})$, by showing that the total variation of the
signal defined over $G_{C}$ achieves an upper bound of at most twice its total
variation wrt the $L_{1}$-norm over $G$, that is, $\sum_{(i,j)\in
E_{C}}|\theta_{i}-\theta_{j}|\leq 2\sum_{(i,j)\in E}|\theta_{i}-\theta_{j}|.$
Hence, if the true signal is smooth over the graph $G$, so that
$\sum_{(i,j)\in E}|\theta_{i}-\theta_{j}|\leq\delta$, then, it is reasonable
to perform de-noising over $G_{C}$ as well, in the sense that the total
variation over $G_{C}$ would satisfy $\sum_{(i,j)\in
E_{C}}|\theta_{i}-\theta_{j}|\leq 2\delta.$ In addition to this, a similar
result holds for the total variation in the $L_{0}$-norm, so that,
$\sum_{(i,j)\in E_{C}}1\\!\mathrm{l}\\{\theta_{i}\neq\theta_{j}\\}\leq
2\sum_{(i,j)\in E}1\\!\mathrm{l}\\{\theta_{i}\neq\theta_{j}\\}.$ A signal
defined over a linear chain graph is then identical to the sequence model
considered earlier in this paper, upto a permutation of the indices of the
signal.

###  7.1 Graph de-noising via Horseshoe fusion

We use our proposed Horseshoe fusion estimation method for graph signal de-
noising. Note that reducing a graph to a DFS-chain graph is not unique, and
hence would result in different chain graphs. We need to first put a prior on
the root vertex, say $\theta_{r},\,r\in V$, and then put suitable priors on
the successive differences $\theta_{i}-\theta_{j}$, for $(i,j)\in E_{C}$. For
inference, we propose to use multiple root nodes to arrive at different DFS-
chain graphs, and then combine the posterior samples to arrive at a final
estimate. Thus, our proposed Bayesian approach can account for the uncertainty
in fusion estimation as well as that in the DFS method. The prior
specification is given by,

$\displaystyle\theta_{r}\mid\lambda_{0}^{2},\sigma^{2}\sim
N(0,\lambda_{0}^{2}\sigma^{2}),\;\theta_{i}-\theta_{j}\mid\lambda_{k(i,j)}^{2},\tau^{2},\sigma^{2}\stackrel{{\scriptstyle
ind}}{{\sim}}N(0,\lambda_{k(i,j)}^{2}\tau^{2}\sigma^{2}),\;(i,j)\in E_{C},$
$\displaystyle\lambda_{i}^{2}\mid\nu_{i}\stackrel{{\scriptstyle
ind}}{{\sim}}IG(1/2,1/\nu_{i}),\;1\leq i\leq n-1,\;\tau^{2}\mid\xi\sim
IG(1/2,1/\xi),$
$\displaystyle\nu_{1},\ldots,\nu_{n-1},\xi\stackrel{{\scriptstyle
iid}}{{\sim}}IG(1/2,1),\;\sigma^{2}\sim IG(a_{\sigma},b_{\sigma}).$ (10)

Here $k(i,j)$ refers to the index corresponding to the edge in $E_{C}$ that
connects vertices $i$ and $j$. The hyperparameters in the above prior
specification are identical to that in (5). The conditional posterior
distributions of the parameters and the hyperparameters are given as follows.
The normal means have the conditional posterior distribution

$\theta_{i}\mid\cdots\sim N(\mu_{i},\zeta_{i}),\;1\leq i\leq n,$ (11)

where $\mu_{i}$ and $\zeta_{i}$ are given by,

$\zeta_{r}^{-1}=\dfrac{1}{\sigma^{2}}\left(1+\sum_{j=1}^{n}\dfrac{1\\!\mathrm{l}\\{(r,j)\in
E_{C}\\}}{\lambda_{k(r,j)}^{2}\tau^{2}}+\dfrac{1}{\lambda_{0}^{2}}\right),\;\mu_{r}=\dfrac{\zeta_{r}}{\sigma^{2}}\left(y_{r}+\sum_{j=1}^{n}\dfrac{\theta_{j}1\\!\mathrm{l}\\{(r,j)\in
E_{C}\\}}{\lambda_{k(r,j)}^{2}\tau^{2}}\right),$
$\zeta_{i}^{-1}=\dfrac{1}{\sigma^{2}}\left(1+\sum_{j=1}^{n}\dfrac{1\\!\mathrm{l}\\{(i,j)\in
E_{C}\\}}{\lambda_{k(i,j)}^{2}\tau^{2}}+\sum_{j=1}^{n}\dfrac{1\\!\mathrm{l}\\{(j,i)\in
E_{C}\\}}{\lambda_{k(j,i)}^{2}\tau^{2}}\right),\;i\in V\backslash\\{r\\},$
$\mu_{i}=\dfrac{\zeta_{i}}{\sigma^{2}}\left(y_{i}+\sum_{j=1}^{n}\dfrac{\theta_{j}1\\!\mathrm{l}\\{(i,j)\in
E_{C}\\}}{\lambda_{k(i,j)}^{2}\tau^{2}}+\sum_{j=1}^{n}\dfrac{\theta_{j}1\\!\mathrm{l}\\{(j,i)\in
E_{C}\\}}{\lambda_{k(j,i)}^{2}\tau^{2}}\right),\;i\in V\backslash\\{r\\}.$

The conditional posteriors for the rest of the parameters are given by,

$\displaystyle\lambda_{k(i,j)}^{2}\mid\cdots$ $\displaystyle\sim$
$\displaystyle
IG\left(1,\dfrac{1}{\nu_{k(i,j)}}+\dfrac{(\theta_{i}-\theta_{j})^{2}}{2\tau^{2}\sigma^{2}}\right),\;(i,j)\in
E_{C},$ $\displaystyle\sigma^{2}\mid\cdots$ $\displaystyle\sim$ $\displaystyle
IG\left(n+a_{\sigma},b_{\sigma}+\dfrac{1}{2}\left[\sum_{i=1}^{n}(y_{i}-\theta_{i})^{2}+\dfrac{1}{\tau^{2}}\sum_{(i,j)\in
E_{C}}\dfrac{(\theta_{i}-\theta_{j})^{2}}{\lambda_{k(i,j)}^{2}}+\dfrac{\theta_{r}^{2}}{\lambda_{0}^{2}}\right]\right),$
$\displaystyle\tau^{2}\mid\cdots$ $\displaystyle\sim$ $\displaystyle
IG\left(\dfrac{n}{2},\dfrac{1}{\xi}+\dfrac{1}{2\sigma^{2}}\sum_{(i,j)\in
E_{C}}^{n}\dfrac{(\theta_{i}-\theta_{j})^{2}}{\lambda_{k(i,j)}^{2}}\right),$
$\displaystyle\nu_{i}\mid\cdots$ $\displaystyle\sim$ $\displaystyle
IG\left(1,1+\dfrac{1}{\lambda_{i}^{2}}\right),\;1\leq i\leq n-1,$
$\displaystyle\xi\mid\cdots$ $\displaystyle\sim$ $\displaystyle
IG\left(1,1+\dfrac{1}{\tau^{2}}\right).$ (12)

###  7.2 Theoretical guarantees

We now present the results on posterior concentration of the estimates
obtained using the proposed Horseshoe based method. The assumptions on the
true model and prior parameters (and, hyperparameters) are almost identical to
that in Section 3, except for some modifications owing to the presence of the
graph $G$. We present the assumptions in detail for the sake of completeness.

######  Assumption 6.

The number $s_{0}$ of true subgraphs of $G$ for which the components of the
signal $\theta_{0}$ are piecewise constant in the model, satisfies $s_{0}\prec
n/\log n$.

######  Assumption 7.

The true mean parameter vector
$\theta_{0}=(\theta_{0,1},\ldots,\theta_{0,n})^{T}$ and the true error
variance $\sigma_{0}^{2}$ satisfy the following conditions:

  1. (i)

$\max_{(i,j)\in E}|\theta_{0,i}-\theta_{0,j}|/\sigma_{0}<L,$ where $\log
L=O(\log n).$

  2. (ii)

For root node $r\in V$,
$\theta_{0,r}/(\lambda_{0}^{2}\sigma_{0}^{2})+2\log\lambda_{0}=O(\log n),$
where $\lambda_{0}$ is the prior hyperparameter appearing in the prior for
$\theta_{r}$ in (10).

######  Assumption 8.

The global scale parameter $\tau$ in the prior specification (10) satisfies
$\tau<n^{-(2+b)}$ for some constant $b>0$, and $-\log\tau=O(\log n)$.

Under the above assumptions, we present the result on posterior concentration
of the proposed estimator.

######  Theorem 9.

Consider the Gaussian means model (1) over a graph $G=(V,E)$ with prior
specification as in (10), and suppose that assumptions 6, 7 and 8 hold. Then
the posterior distribution of $\theta$, given by $\Pi^{n}(\cdot\mid y)$,
satisfies

$\Pi^{n}(\|\theta-\theta_{0}\|_{2}/\sqrt{n}\geq M\sigma_{0}\epsilon_{n}\mid
y)\rightarrow 0,\;\mathrm{as}\,n\rightarrow\infty,$ (13)

in probability or in $L_{1}$ wrt the probability measure of $y$, for
$\epsilon_{n}\asymp\sqrt{s_{0}\log n/n}$ and a constant $M>0$.

The above theorem implies that the mean squared error of $\theta$ is of the
order $s_{0}\log n/n$, which differs from the oracle contraction rate
$s_{0}/n$ under known graph sparsity by only a logarithmic factor in $n$. Note
that our prior construction does not require any knowledge regarding the
sparsity pattern of the graph, and hence the proposed Bayesian estimation
procedure can adapt to the total variation metric(wrt the $L_{0}$-norm) and
acheive the desired optimality in estimating the signal over the graph $G$. We
discuss the proof of the above result in the appendix.

We now apply our method on a real-world network and check the performance of
signal de-noising over the same.

###  7.3 Application: Signal de-noising over Euroroad network

We consider signals defined on a real-world graph, given by the Euroroad
network (Rossi and Ahmed, , 2015), obtained from the repository
http://networkrepository.com/inf-euroroad.php. The Euroraod network is a real
road network connecting various cities in Europe, having a total of around
1200 vertices and about 1400 edges. We first extract the largest connected
component of the network, which results in 1039 vertices and 1305 edges. Next,
to generate true piecewise constant signals on this network, we extract non-
overlapping communities from the network using the Louvain community detection
algorithm (Blondel et al., , 2008). The Louvain algorithm is one of the most
popular community detection algorithms available in the literature which
detects underlying communities in a network via optimising the modularity of
the same. We use the cluster_louvain function available in the igraph package
in R to extract the communities. The method results in 25 connected
communities, with the smallest community having 11 vertices and the largest
community having 82 vertices as members. The true signal on vertex $j$ with
community label $k$ is given by
$\theta_{0,j}=k1\\!\mathrm{l}\\{\mathrm{vertex}\,j\in\mathrm{community}\,k\\},\,1\leq
j\leq n.$ Finally we generate a noisy version of the signal on the network
given by $y_{j}\sim N(\theta_{0,j},0.3^{2}),\,1\leq j\leq n.$

We apply our proposed Horseshoe fusion method on the above network by first
reducing the same to a DFS-chain graph. We choose three root vertices at
random from $V$, and construct three different DFS-chain graphs, followed by
drawing MCMC samples from the posterior independently for the three chains.
Finally, we merge the three sets of posterior samples after discarding the
burn-in samples in each of the cases and find the posterior mean of the
combined sample to arrive at the estimates. The MSE in estimation came out to
be $0.022$ and the adjusted MSE is $0.0001$. The network visualization for the
graph with true signal and the de-noised graph are provided in Figure 8, where
the vertices of the former are color-coded according to their respective
communities. Note that since the Horseshoe prior is a continuous shrinkage
prior, no two estimated signals are of equal magnitude. Hence while coloring
the vertices of the de-noised graph, we classify the estimated values in 25
bins of equal length, and color the vertices so that vertices having estimated
signals falling in the same bin have identical colors. Additionally, in Figure
9, we provide visualizations of the true signal and the corresponding
estimates as a step-function via permutation of the graph vertices to
demonstrate the estimation performance of our proposed Bayesian fusion method.
The figure re-iterates the excellent performance of our method to detect the
piecewise constant signals over different communities in the actual graph. It
is worthwhile to note in this regard that we work only with a transformed
version of the original graph structure, yet we are able to learn the signals
over the communities in an efficient manner, adapting to the sparsity level of
the graphical structure.

---
(a) Euroroad network with true signal values
(b) Euroroad network with Bayesian Horseshoe fusion estimation-based values
Figure 8: Figure showing estimation performances of the proposed Bayesian
Horseshoe prior based fusion method. Figure 9: Step plot showing the Bayesian
Horseshoe fusion performance for graph signal de-noising. The true signals are
in blue and the estimated signals are in black.

##  8 Discussion

In this paper, we have addressed the problem of fusion estimation in a normal
sequence model where the underlying mean parameter vector is piecewise
constant. Noise removal and block structure learning is an important and
widely studied problem with applications in a plethora of areas. We developed
a Horseshoe prior based fusion estimation method and provided theoretical
guarantees via establishing posterior convergence and selection consistency.
Numerical studies demonstrated excellent performance of our method as compared
to other competing Bayesian and frequentist approaches.

In many applications, we observe noisy signals defined over an undirected
graph. The excellent performance of the Horseshoe prior based learning
motivated us to extend the same to the graph signal denoising problem as well.
We constructed a Gibbs sampler for efficient learning of graph signals via
transformation of an arbitrary graph to a DFS-chain graph, and provided
theoretical guarantees of the procedure.

There are several directions for extending our work. One important area is to
move beyond normal sequence models and consider multiple graphical models. The
graph fused lasso method for estimation and structure learning in multiple
Gaussian graphical models (Danaher et al., , 2014) has got much attention in
recent years. The Horseshoe fusion prior may be utilized in this instance for
precision matrix estimation over several groups and structure learning of
respective graphical models. Additionally, our proposed approach may be
extended to isotonic regression models, where there is a natural ordering in
the means of the observations. We propose to explore these ideas as future
work.

## Data availability and software

Data for the DNA copynumber analysis are available for free from the archived
version of the cghFlasso package in R. Data on solar X-ray flux are available
for free from the National Oceanic and Atmospheric Administration (NOAA)’s
website at https://www.ngdc.noaa.gov/stp/satellite/goes/dataaccess.html. The
Euroroad network data are available for free from the network repository
http://networkrepository.com/inf-euroroad.php.

R codes to implement the fusion estimation procedure for the normal sequence
model and for the graph signal de-noising problem, along with the codes to
transform an arbitrary graph to a DFS-chain graph are available at
https://github.com/sayantanbiimi/HSfusion.

## Acknowledgements

The author is supported by DST INSPIRE Faculty Award Grant No. 04/2015/002165,
and also by IIM Indore Young Faculty Research Chair Award grant.

## Appendix

### Additional lemmas and proofs of results

######  Lemma 10.

Under the assumptions on the true parameter $\theta_{0}$ and the prior
distribution on $\theta$ as outlined in Theorem 4, for some $b^{\prime}>0$, we
have,

$\int_{-s_{0}\log n/n^{2}}^{s_{0}\log n/n^{2}}p_{HS}(\eta;\tau)\,d\eta\geq
1-n^{-b^{\prime}},$ (14)

where $p_{HS}(\eta;\tau)$ denotes the Horseshoe prior density on $\eta$ with
hyperparameter $\tau>0$.

###### Proof.

Define $a_{n}=s_{0}\log n/n^{2}$. We will show that
$\int_{-a_{n}}^{a_{n}}p_{HS}(\eta;\tau)\,d\eta\geq 1-n^{-b}$ for some $b>0$
under the assumptions outlined in Theorem 4. The Horseshoe prior density does
not assume a closed form; however, Carvalho et al.,  (2010) showed that the
density admits tight lower and upper bounds given by,

$\dfrac{1}{\tau(2\pi)^{3/2}}\log\left(1+\dfrac{4\tau^{2}}{\theta^{2}}\right)<p_{HS}(\theta\mid\tau)<\dfrac{2}{\tau(2\pi)^{3/2}}\log\left(1+\dfrac{2\tau^{2}}{\theta^{2}}\right).$
(15)

Hence, we have,

$\displaystyle
1-\int_{-a_{n}}^{a_{n}}p_{HS}(\eta;\tau)\,d\eta=2\int_{a_{n}}^{\infty}p_{HS}(\eta;\tau)\,d\eta$
$\displaystyle<\left(\dfrac{2}{\pi^{3}}\right)^{1/2}\dfrac{1}{\tau}\int_{a_{n}}^{\infty}\log\left(1+\dfrac{2\tau^{2}}{\eta^{2}}\right)\,d\eta$
$\displaystyle<\left(\dfrac{2}{\pi^{3}}\right)^{1/2}\dfrac{1}{\tau}\int_{a_{n}}^{\infty}\dfrac{2\tau^{2}}{\eta^{2}}\,d\eta$
$\displaystyle\lesssim\dfrac{\tau}{a_{n}}\leq n^{-b^{\prime}},$

where $0<b^{\prime}<b$, for $\tau\leq a_{n}n^{-b}\asymp n^{-(2+b)}.$ This
completes the proof.

∎

As mentioned earlier, the above result guarantees that the Horseshoe prior
puts sufficient mass around zero for the successive differences $\eta_{i}$ so
as to facilitate Bayesian shrinkage for the block-structured mean parameters.
However, the prior should be able to retrieve the blocks effectively as well,
so that successive differences that are non-zero should not be shrunk too
much. The following result below guarantees that the Horseshoe prior is
‘thick’ at non-zero parameter values, so that it is not too sharp.

######  Lemma 11.

Consider the prior structure and the assumptions as mentioned in Theorem 4.
Then, we have,

$-\log\left(\inf_{\eta/\sigma\in[-L,L]}p_{HS}(\eta;\tau)\right)=O(\log n).$
(16)

###### Proof.

We utilise the tight lower bound for the Horseshoe prior density as in 15. We
have,

$\displaystyle\inf_{\eta/\sigma\in[-L,L]}p_{HS}(\eta;\tau)$ $\displaystyle>$
$\displaystyle\inf_{\eta/\sigma\in[-L,L]}\dfrac{1}{\tau(2\pi)^{3/2}}\log\left(1+\dfrac{4\tau^{2}}{\eta^{2}}\right)$
$\displaystyle\geq$
$\displaystyle\dfrac{1}{\tau(2\pi)^{3/2}}\log\left(1+\dfrac{4\tau^{2}}{L^{2}}\right)\asymp\dfrac{\tau}{L^{2}}$
$\displaystyle=$ $\displaystyle
O(n^{-b^{\prime\prime}}),\;b^{\prime\prime}>0,$

for $-\log(\tau)=O(\log n)$ and $\log L=O(\log n).$ Thus, we get,
$-\log\left(\inf_{\eta/\sigma\in[-L,L]}p_{HS}(\eta;\tau)\right)=O(\log n),$
hence completing the proof.

∎

We now present the proofs of the main results in our paper.

###### Proofs of Theorem 4 and Theorem 5.

The proof readily follows from Theorems 2.1 and 2.2 in Song and Cheng,
(2020), if we can verify the conditions (see display (2.5) in the
aforementioned paper) for posterior convergence specified therein. Lemma 10
and Lemma 16 imply that the first two conditions are satisfied. The third
condition is satisfied for a Normal prior distribution on $\theta_{1}$ and for
a fixed hyperparameter $\lambda_{1}>0.$ The fourth condition is trivially
satisfied for fixed choices of (non-zero) hyperparamters $a_{\sigma}$ and
$b_{\sigma}$, and for some fixed (unknown) true error variance
$\sigma_{0}^{2}.$ ∎

###### Proof of Theorem 9.

Note that, by Lemma 1 in Padilla et al.,  (2017), the total variation of
$\theta$ wrt to the $L_{0}$-norm in the DFS-chain graph $G_{C}$ is bounded by
at most twice that in $G$, that is, $\theta\in l_{0}[G,s]$ always implies that
$\theta\in l_{0}[2s].$ Also, assumptions 6 \- 8 for an arbitrary graph $G$
implies assumptions 1 \- 3 for a linear chain graph. Hence the result readily
follows from that of Theorem 4. ∎

### A practical solution to block structure recovery

The Horseshoe prior is a continuous shrinkage prior, and hence block structure
recovery is not straight-forward. In Bayesian fusion estimation with Laplace
shrinkage prior or with $t$-shrinkage prior, one may discretize the samples
obtained via MCMC using a suitable threshold. Song and Cheng,  (2020)
recommended using the $1/2n$-th quantile of the corresponding prior for
discretization of the scaled samples. In Section 4, we proposed a threshold
based on the posterior contraction rate. However, that would also require an
idea of the true block size, which may not be available always, or may be
difficult to ascertain beforehand. We provide a practical approach for
discretization of the samples. Note that in a sparse normal sequence problem,
the posterior mean of the mean parameter $\theta_{i}$ is given by
$\kappa_{i}y_{i}$, where $\kappa_{i}$ is the shrinkage weight. These shrinkage
weights mimic the posterior selection probability of the means, and hence
thresholding the weights to 0.5 provide excellent selection performance, which
is justified numerically (Carvalho et al., , 2010) and later theoretically
(Datta and Ghosh, , 2013). Motivated by the same, we propose to use the
following thresholding rule for selection of the blocks for our proposed
method: $\theta_{j_{1}}$ and $\theta_{j_{2}}$ are equal if
$|\hat{\theta}_{j_{1}}-\hat{\theta}_{j_{2}}|<0.5|y_{j_{1}}-y_{j_{2}}|$, where
$\theta_{j_{1}},\theta_{j_{2}}$ are the posterior means of $\theta_{j_{1}}$
and $\theta_{j_{2}}$respectively, for $1\leq j_{1}\neq j_{2}\leq n.$ We thus
estimate the block structure indicator
$s_{j_{1}j_{2}}=1\\!\mathrm{l}\\{\theta_{j_{1}}=\theta_{j_{2}}\\}$ by
$\hat{s}_{j_{1}j_{2}}=1\\!\mathrm{l}\\{|\hat{\theta}_{j_{1}}-\hat{\theta}_{j_{2}}|<0.5|y_{j_{1}}-y_{j_{2}}|\\}$.
To evaluate the performance of the structure recovery method using the
proposed thresholding approach, we define the metrics $W$ and $B$, where $W$
is the average within-block variation defined as
$W:=\mathrm{mean}_{\\{s_{j_{1}j_{2}}\neq
0\\}}|\hat{\theta}_{j_{1}}-\hat{\theta}_{j_{2}}|$, and $B$ is the between-
block separation defined as
$B:=\min_{\\{s_{j_{1}j_{2}}=0\\}}|\hat{\theta}_{j_{1}}-\hat{\theta}_{j_{2}}|.$

This practical thresholding approach resulted in excellent block structure
recovery performance with respect to the metrics $W$ and $B$. We compare our
results with the other Bayesian fusion estimation methods and also with the
frequentist fused lasso method. The fused lasso method results in exact
structure recovery, and hence no thresholding is required. As mentioned
earlier, the thresholding rules for the competing Bayesian methods are taken
as suggested in Song and Cheng,  (2020). We notice that for low error
variances, the within-block average variation is lower in case of the
Horseshoe fusion and $t$-fusion priors, with the Horseshoe fusion method
having larger between-block speration especially in case of higher noise, thus
indicating superior structure learning capabilities. The results are
summarized in Table 2.

Other possible directions for coming up with a threshold include using a
multiple hypothesis testing approach for successive differences in the means,
and also using simultaneous credible intervals. We leave the theoretical
treatment of using (or improving, if possible) our practical thresholding
approach and exploring other proposed methods as a future work.

|  | HS-fusion |  $t$-fusion | Laplace fusion |  $L_{1}$ fusion
---|---|---|---|---|---
Signal | $\sigma$ | W | B | W | B | W | B | W | B
| 0.1 | 0.039 (0.001) | 0.846 (0.005) | 0.035 (0.001) | 0.892 (0.010) | 0.678 (0.002) | 0.003 (0.000) | 0.069 (0.003) | 0.739 (0.011)
Even | 0.3 | 0.140 (0.003) | 0.425 (0.020) | 0.121 (0.005) | 0.381 (0.035) | 0.687 (0.003) | 0.003 (0.000) | 0.181 (0.007) | 0.333 (0.023)
| 0.5 | 0.250 (0.005) | 0.139 (0.016) | 0.221 (0.007) | 0.079 (0.017) | 0.704 (0.004) | 0.003 (0.000) | 0.285 (0.009) | 0.119 (0.016)
| 0.1 | 0.037 (0.001) | 0.846 (0.005) | 0.034 (0.001) | 0.864 (0.010) | 0.537 (0.002) | 0.006 (0.001) | 0.070 (0.003) | 0.751 (0.010)
Uneven | 0.3 | 0.134 (0.003) | 0.457 (0.020) | 0.106 (0.004) | 0.266 (0.031) | 0.548 (0.003) | 0.005 (0.001) | 0.182 (0.006) | 0.348 (0.021)
| 0.5 | 0.238 (0.005) | 0.146 (0.015) | 0.173 (0.006) | 0.022 (0.007) | 0.570 (0.004) | 0.004 (0.000) | 0.284 (0.009) | 0.131 (0.015)
| 0.1 | 0.036 (0.001) | 0.800 (0.007) | 0.034 (0.001) | 0.670 (0.022) | 0.325 (0.002) | 0.008 (0.001) | 0.106 (0.001) | 0.710 (0.010)
V. Uneven | 0.3 | 0.128 (0.002) | 0.401 (0.025) | 0.092 (0.003) | 0.013 (0.004) | 0.344 (0.003) | 0.007 (0.001) | 0.276 (0.006) | 0.269 (0.023)
| 0.5 | 0.222 (0.004) | 0.149 (0.019) | 0.124 (0.004) | 0.001 (0.000) | 0.376 (0.004) | 0.006 (0.001) | 0.391 (0.010) | 0.105 (0.017)
Table 2: Within (W) and between (B) block average variation (with associated
standard errors in parentheses) for different fusion estimation methods. For
good structure recovery, W values should be low and B values should be high.

## References

  * Banerjee et al., (2021) Banerjee, S., Castillo, I., and Ghosal, S. (2021).  Bayesian inference in high-dimensional models.  arXiv preprint arXiv:2101.04491.
  * Banerjee and Shen, (2020) Banerjee, S. and Shen, W. (2020).  Graph signal denoising using $t$-shrinkage priors.  arXiv preprint arXiv:2012.13696.
  * Bhadra et al., (2017) Bhadra, A., Datta, J., Polson, N. G., and Willard, B. (2017).  The horseshoe+ estimator of ultra-sparse signals.  Bayesian Analysis, 12(4):1105–1131.
  * Bhadra et al., (2019) Bhadra, A., Datta, J., Polson, N. G., and Willard, B. T. (2019).  The horseshoe-like regularization for feature subset selection.  Sankhya B, pages 1–30.
  * Bhattacharya et al., (2015) Bhattacharya, A., Pati, D., Pillai, N. S., and Dunson, D. B. (2015).  Dirichlet–laplace priors for optimal shrinkage.  Journal of the American Statistical Association, 110(512):1479–1490.
  * Blondel et al., (2008) Blondel, V. D., Guillaume, J.-L., Lambiotte, R., and Lefebvre, E. (2008).  Fast unfolding of communities in large networks.  Journal of statistical mechanics: theory and experiment, 2008(10):P10008.
  * Bredel et al., (2005) Bredel, M., Bredel, C., Juric, D., Harsh, G. R., Vogel, H., Recht, L. D., and Sikic, B. I. (2005).  High-resolution genome-wide mapping of genetic alterations in human glial brain tumors.  Cancer Research, 65(10):4088–4096.
  * Bühlmann and Van De Geer, (2011) Bühlmann, P. and Van De Geer, S. (2011).  Statistics for high-dimensional data: methods, theory and applications.  Springer Science & Business Media.
  * Carvalho et al., (2009) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2009).  Handling sparsity via the horseshoe.  In Artificial Intelligence and Statistics, pages 73–80. PMLR.
  * Carvalho et al., (2010) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2010).  The horseshoe estimator for sparse signals.  Biometrika, 97(2):465–480.
  * Castillo et al., (2015) Castillo, I., Schmidt-Hieber, J., and Van der Vaart, A. (2015).  Bayesian linear regression with sparse priors.  Annals of Statistics, 43(5):1986–2018.
  * Danaher et al., (2014) Danaher, P., Wang, P., and Witten, D. M. (2014).  The joint graphical lasso for inverse covariance estimation across multiple classes.  Journal of the Royal Statistical Society. Series B, Statistical methodology, 76(2):373.
  * Datta and Ghosh, (2013) Datta, J. and Ghosh, J. K. (2013).  Asymptotic properties of bayes risk for the horseshoe prior.  Bayesian Analysis, 8(1):111–132.
  * Donoho et al., (1992) Donoho, D. L., Johnstone, I. M., Hoch, J. C., and Stern, A. S. (1992).  Maximum entropy and the nearly black object.  Journal of the Royal Statistical Society: Series B (Methodological), 54(1):41–67.
  * Fan and Li, (2001) Fan, J. and Li, R. (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American Statistical Association, 96(456):1348–1360.
  * Fan and Guan, (2018) Fan, Z. and Guan, L. (2018).  Approximate $\ell_{0}$-penalized estimation of piecewise-constant signals on graphs.  The Annals of Statistics, 46(6B):3217–3245.
  * Friedman et al., (2008) Friedman, J., Hastie, T., and Tibshirani, R. (2008).  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics, 9(3):432–441.
  * Griffin and Brown, (2010) Griffin, J. E. and Brown, P. J. (2010).  Inference with normal-gamma prior distributions in regression problems.  Bayesian Analysis, 5(1):171–188.
  * Hoerl and Kennard, (1970) Hoerl, A. E. and Kennard, R. W. (1970).  Ridge regression: applications to nonorthogonal problems.  Technometrics, 12(1):69–82.
  * Holland, (2000) Holland, E. C. (2000).  Glioblastoma multiforme: the terminator.  Proceedings of the National Academy of Sciences, 97(12):6242–6244.
  * Ishwaran and Rao, (2005) Ishwaran, H. and Rao, J. S. (2005).  Spike and slab variable selection: frequentist and bayesian strategies.  Annals of Statistics, 33(2):730–773.
  * Kyung et al., (2010) Kyung, M., Gill, J., Ghosh, M., and Casella, G. (2010).  Penalized regression, standard errors, and bayesian lassos.  Bayesian Analysis, 5(2):369–411.
  * Lai et al., (2005) Lai, W. R., Johnson, M. D., Kucherlapati, R., and Park, P. J. (2005).  Comparative analysis of algorithms for identifying amplifications and deletions in array cgh data.  Bioinformatics, 21(19):3763–3770.
  * Little and Jones, (2011) Little, M. A. and Jones, N. S. (2011).  Generalized methods and solvers for noise removal from piecewise constant signals. i. background theory.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 467(2135):3088–3114.
  * Mitchell and Beauchamp, (1988) Mitchell, T. J. and Beauchamp, J. J. (1988).  Bayesian variable selection in linear regression.  Journal of the American Statistical Association, 83(404):1023–1032.
  * Padilla et al., (2017) Padilla, O., Sharpnack, J., and Scott, J. (2017).  The dfs fused lasso: Linear-time denoising over general graphs.  The Journal of Machine Learning Research, 18(1):6410–6445.
  * Park and Casella, (2008) Park, T. and Casella, G. (2008).  The bayesian lasso.  Journal of the American Statistical Association, 103(482):681–686.
  * (28) Piironen, J. and Vehtari, A. (2017a).  On the hyperprior choice for the global shrinkage parameter in the horseshoe prior.  In Artificial Intelligence and Statistics, pages 905–913. PMLR.
  * (29) Piironen, J. and Vehtari, A. (2017b).  Sparsity information and regularization in the horseshoe and other shrinkage priors.  Electronic Journal of Statistics, 11(2):5018–5051.
  * Pulkkinen et al., (2005) Pulkkinen, A., Lindahl, S., Viljanen, A., and Pirjola, R. (2005).  Geomagnetic storm of 29–31 october 2003: Geomagnetically induced currents and their relation to problems in the swedish high-voltage power transmission system.  Space Weather, 3(8).
  * Rinaldo, (2009) Rinaldo, A. (2009).  Properties and refinements of the fused lasso.  Annals of Statistics, 37(5B):2922–2952.
  * Ročková and George, (2018) Ročková, V. and George, E. I. (2018).  The spike-and-slab lasso.  Journal of the American Statistical Association, 113(521):431–444.
  * Rossi and Ahmed, (2015) Rossi, R. A. and Ahmed, N. K. (2015).  The network data repository with interactive graph analytics and visualization.  In AAAI.
  * Shimamura et al., (2019) Shimamura, K., Ueki, M., Kawano, S., and Konishi, S. (2019).  Bayesian generalized fused lasso modeling via neg distribution.  Communications in Statistics-Theory and Methods, 48(16):4132–4153.
  * Song, (2020) Song, Q. (2020).  Bayesian shrinkage towards sharp minimaxity.  Electronic Journal of Statistics, 14(2):2714–2741.
  * Song and Cheng, (2020) Song, Q. and Cheng, G. (2020).  Bayesian fusion estimation via t-shrinkage.  Sankhya A, 82(2):353–385.
  * Song and Liang, (2017) Song, Q. and Liang, F. (2017).  Nearly optimal bayesian shrinkage for high dimensional regression.  arXiv preprint arXiv:1712.08964.
  * Tarjan, (1972) Tarjan, R. (1972).  Depth-first search and linear graph algorithms.  SIAM journal on computing, 1(2):146–160.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288.
  * Tibshirani et al., (2005) Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., and Knight, K. (2005).  Sparsity and smoothness via the fused lasso.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(1):91–108.
  * Tibshirani and Wang, (2008) Tibshirani, R. and Wang, P. (2008).  Spatial smoothing and hot spot detection for cgh data using the fused lasso.  Biostatistics, 9(1):18–29.
  * van der Pas et al., (2017) van der Pas, S., Szabó, B., and van der Vaart, A. (2017).  Adaptive posterior contraction rates for the horseshoe.  Electronic Journal of Statistics, 11(2):3196–3225.
  * Van Der Pas et al., (2014) Van Der Pas, S. L., Kleijn, B. J., and Van Der Vaart, A. W. (2014).  The horseshoe estimator: Posterior concentration around nearly black vectors.  Electronic Journal of Statistics, 8(2):2585–2618.
  * Wang, (2012) Wang, H. (2012).  Bayesian graphical lasso models and efficient posterior computation.  Bayesian Analysis, 7(4):867–886.
  * Wei and Ghosal, (2020) Wei, R. and Ghosal, S. (2020).  Contraction properties of shrinkage priors in logistic regression.  Journal of Statistical Planning and Inference, 207:215–229.
  * Zhang et al., (2016) Zhang, Y., Reich, B. J., and Bondell, H. D. (2016).  High dimensional linear regression via the r2-d2 shrinkage prior.  arXiv preprint arXiv:1609.00046.
  * Zou and Hastie, (2005) Zou, H. and Hastie, T. (2005).  Regularization and variable selection via the elastic net.  Journal of the Royal Statistical Society: Series B (Statistical methodology), 67(2):301–320.