T-LoHo: A Bayesian Regularization Model for Structured Sparsity and Smoothness on Graphs
Year: 2022

Abstract:
Graphs have been commonly used to represent complex data structures. In models dealing with graph-structured data, multivariate parameters may not only exhibit sparse patterns but have structured sparsity and smoothness in the sense that both zero and non-zero parameters tend to cluster together. We propose a new prior for high-dimensional parameters with graphical relations, referred to as the Treebased Low-rank Horseshoe (T-LoHo) model, that generalizes the popular univariate Bayesian horseshoe shrinkage prior to the multivariate setting to detect structured sparsity and smoothness simultaneously. The T-LoHo prior can be embedded in many high-dimensional hierarchical models. To illustrate its utility, we apply it to regularize a Bayesian high-dimensional regression problem where the regression coefficients are linked by a graph, so that the resulting clusters have flexible shapes and satisfy the cluster contiguity constraint with respect to the graph. We design an efficient Markov chain Monte Carlo algorithm that delivers full Bayesian inference with uncertainty measures for model parameters such as the number of clusters. We offer theoretical investigations of the clustering effects and posterior concentration results. Finally, we illustrate the performance of the model with simulation studies and a real data application for anomaly detection on a road network. The results indicate substantial improvements over other competing methods such as the sparse fused lasso. Code is available at https://github.com/changwoo-lee/TLOHO.

heading: 1 Introduction
text: In high-dimensional models such as linear regressions where the number of parameters p may exceed the sample size n, it is often assumed that the parameter vector β ∈ Rp has many zero components, namely the sparsity assumption. This sparsity assumption allows the unknown parameter β to lie on a low-dimensional subspace of Rp, which addresses overfitting and leads to improved predictions and easier interpretations (Hastie et al., 2015). In many real-life applications, however, the parameter of interest β should be understood in a certain context where specific data structures exist, such as in time series, spatial, image, and network data analysis. In many such cases, it is desirable to consider another type of low-dimensional structure, where β is assumed to have clustered patterns and many possibly clustered zeros, which we shall call sparse homogeneity assumption.
One of the most popular models which assumes sparse homogeneity is the sparse fused lasso (Tibshirani et al., 2005) in linear regression settings. 1-D sparse fused lasso imposes an `1 penalty on the differences of time-neighboring coefficients as well as individual coefficients. The notion of ‘time-neighboring’ can be further generalized by considering a graph G= (V,E) with |V | = p so that G represents the general context or structure in which β should be interpreted. It leads to the generalized lasso (Tibshirani et al., 2011) where the penalty term now involves m := |E| number of pairwise differences of neighboring vertices with respect to G which facilitates βi = βj , where
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
(i, j) ∈ E. In other words, the generalized lasso now allows neighboring coefficients with respect to the graph G to be clustered together. Other similar graph-based regularization approaches include graph OSCAR (Yang et al., 2012), grouping pursuit (Zhu et al., 2013) and graph trend filtering (Wang et al., 2016). However, when the number of edges, m, is large relative to p, regularization over the whole graph structure often faces severe computational burden. Ke et al. (2015) and Tang and Song (2016) addressed this problem by constructing a preliminary ranking of coefficients and performing a segmentation using 1-D fused lasso. Padilla et al. (2017) used depth-first search (DFS) algorithm to establish a chain graph order for constructing generalized lasso. In a spatial setting, Li and Sang (2019) used a Euclidean minimum spanning tree to form the fused lasso penalty which also enjoys computational benefits. Nevertheless, these fixed chain or tree orders may not be compatible with the true cluster with respect to its context G and hence may lead to over-clustering. Besides, these optimization-based penalized estimators do not usually come with uncertainty measures.
Bayesian regularization methods have gained great popularity in high-dimensional models due to its flexibility and convenience in quantifying estimation and prediction uncertainties. There is rich literature in Bayesian high-dimensional regression models under the sparsity assumption, which can be roughly summarized into two categories: (i) spike-and-slab priors (George and McCulloch, 1993, 1997) and (ii) global-local shrinkage priors (see Bhadra et al., 2019, and references therein). A particular shrinkage prior, called the horseshoe prior (Carvalho et al., 2010), has gained a lot of attraction due to its tail-robustness and super-efficiency as well as substantial computational benefits of the global-local shrinkage prior family (Polson and Scott, 2010). However, relatively few works have been done in the area of Bayesian high-dimensional regression models under the sparse homogeneity assumption. Most approaches put sparsity-inducing priors directly on the pairwise differences of neighboring vertices of a graph to achieve the graph-structured homogeneity of β (Kyung et al., 2010; Shimamura et al., 2019; Song and Cheng, 2020; Kim and Gao, 2020; Banerjee, 2021). But similar to the generalized lasso, unless resorting to approximate methods such as EM algorithm, these approaches become computationally expensive to get posterior samples for general graphs with a large number of edges. Another more subtle limitation of these methods is that they fail to incorporate structural assumptions among local shrinkage parameters.
In this paper, we propose a Bayesian Tree-based Low-rank Horseshoe (T-LoHo) model to identify structured sparsity and smoothness of parameters whose prior structure is represented by a graph. Unlike other existing Bayesian methods which put independent priors on pairwise differences, TLoHo extends the univariate horseshoe shrinkage prior to a multivariate setting where β and its local shrinkage parameters are assumed to be piecewise constants on a graph. This low-rank structured model for local shrinkage parameters allows both clustering and sparsity effects to have strong local adaptivities. A random spanning forest (RSF) based graph clustering prior is introduced to adaptively learn a compatible neighboring order to model graph partitions, which extends the recently developed random spanning tree partition models (Teixeira et al., 2019; Luo et al., 2021) to allow for possibly unconnected graphs. We show that this RSF-based prior retains model richness as its support is flexible enough to accommodate all possible contiguous partitions. We introduce T-LoHo prior to model high-dimensional linear regression coefficients, although it can be flexibly embedded in other high-dimensional models. The resulting cluster estimates not only provide uncertainty measures but also have strong flexibility to accommodate sharp discontinuities. With state of the art computational strategies, we provide a highly efficient MCMC algorithm for posterior inference utilizing tree structures. We also study theoretical aspects of our model including posterior consistency results and its effect on clustering. We demonstrate the efficacy of T-LoHo model by using synthetic data and a real data analysis for anomaly detection on the Manhattan road network.

heading: 2 T-LoHo: Tree-based Low-rank Horseshoe Model
text: Consider a graph G = (V,E) with |V | = p and |E| = m which represents the pre-known structure of our parameter of interest β ∈ Rp. For example, when β has a certain order (e.g. time), then we set G as a linear chain graph. When (vectorized) β lies on a 2-D image, we let G be a 2-D lattice. A contiguous graph partition Π={C1,. . ., CK} is a disjoint partition of V such that each Ck induces a connected subgraph of G (the term graph partition will always refer to a contiguous graph partition hereafter). Under the sparse homogeneity assumption on β, the goal is to find a graph partition Π with unknown size K and corresponding estimate of β which contains many zeros. We introduce the model for β conditional on a partition Π in Section 2.1 and describe the model for Π in Section 2.2.

heading: 2.1 Low-Rank Multivariate Horseshoe Prior
text: First we review the horseshoe prior (Carvalho et al., 2010), which provides a sparse estimate of β by shrinking small effects towards zero while maintaining robustness to large signals:
β|σ2, τ2, {λj}pj=1 ∼ Np(0, σ 2τ2diag(λ21, . . . , λ 2 p)),
λj iid∼ C+(0, 1), τ ∼ C+(0, τ0), p(σ2) ∝ 1/σ2,
where Np(m,V ) denotes a p-dimensional multivariate normal distribution with meanm and covariance V , and C+(0, γ) denotes a half-Cauchy distribution with density p(λ) = 2/[πγ(1 + λ2/γ2)] for λ > 0. Here τ is a global shrinkage parameter with hyperparameter τ0 to enforce global shrinkage towards zero, {λj} are local shrinkage parameters with heavy tails which allow some of the βj’s to escape the shrinkage, and σ2 is a scaling factor which is often assumed to be the noise variance.
Given a graph partition Π of G, we seek a new joint prior distribution of β ∈ Rp to impose sparse homogeneity assumption while incorporating the clustering structure of Π . For example when p = 5 and Π = {{1}, {2, 4}, {3, 5}} as shown in Figure 1(a), we seek a distribution which lies on a 3-dimensional subspace {β ∈ R5 : β2 = β4, β3 = β5}. In general, given a graph partition Π with size K, we can construct a K × p matrix Φ from Π:
Φkj = 1/ √ |Ck| if j ∈ Ck and 0 otherwise, k = 1, . . . ,K, j = 1, . . . , p. (1)
This K × p matrix Φ represents a projection from a p-dimensional space to a reduced K-dimensional space. Note that rows of Φ are mutually orthonormal so that ΦΦ> = IK . Using Φ obtained from Π, we introduce the low-rank horseshoe (LoHo) prior to impose sparse homogeneity assumption on β:
β |σ2, τ2,Λ,Π ∼ Np(0, σ2τ2Φ>ΛΦ), Λ := diag(λ21, . . . , λ2K), (2)
λk iid∼ C+(0, 1), τ ∼ C+(0, τ0), p(σ2) ∝ 1/σ2. (3)
See fig. 1(a) for an example of Φ and a covariance matrix Φ>ΛΦ. Note LoHo introduces a covariance matrix capturing clustering dependence among β and its marginal local shrinkage for simultaneous sparsity and fusion, a key distinction from most existing methods (Carvalho et al., 2010; Kyung et al., 2010; Shimamura et al., 2019) that assume independence or only dependence among β but not local shrinkage parameters. Since K is assumed to be small relative to p, the covariance matrix Φ>ΛΦ has low rank (i.e., rank(Φ>ΛΦ) = K << p), and thus LoHo does not have a density with respect to Lebesgue measure on Rp. Instead, its distribution lies on the row space of Φ with dim(rowsp(Φ)) = K, and we can consider the transformation β̃ = Φβ so that β̃ has a NK(0, σ2τ2Λ) density with respect to Lebesgue measure on rowsp(Φ). Observe that Φ> is the Moore–Penrose pseudoinverse of Φ, which implies Φ>Φ is a projection onto rowsp(Φ) so that we can recover β = Φ>β̃.
By assigning a half-Cauchy prior, global shrinkage parameter τ creates a strong pull towards zero while clusterwise local shrinkage parameters {λk}Kk=1 allow some of the clusterwise β̃k’s to escape the shrinkage. Assuming Φ having orthonormal rows is important in the sense that the effect of local shrinkage parameters {λk}Kk=1 remains same across the clusters with varying size. Thus using the projection onto the low-dimensional subspace, LoHo gives more parsimonious estimate of β by forming clusters of zero and non-zero parameters.
LoHo can be naturally incorporated into a linear model. With response vector y ∈ Rn and columnstandardized design matrix X ∈ Rn×p so that each column has unit `2 norm, we can write
y = Xβ + , ∼ Nn(0, σ2In)
Under this formulation, LoHo has a close connection with Bayesian compressed regression (Guhaniyogi and Dunson, 2015). It randomly projects predictors Xi 7→ ΦXi with a certain family of matrix Φ and performs model averaging, at the cost of losing interpretability of β. But LoHo directly introduces a prior on β using the projection matrix Φ defined as (1) so that it induces clustered coefficient while maintaining interpretability of β. Also, functional horseshoe (Shin et al., 2020) shrinks Xβ towards the subspace of colsp(X) while LoHo shrinks β towards 0 along the rowsp(Φ).
We remark that, although it appears (τ,Λ) and Φ handle sparsity and homogeneity separately, shrinkage component in LoHo also affects clustering by facilitating cluster fusion when signal is low and improving cluster identification when signal is high. More details are discussed in Section 3.2.

heading: 2.2 Tree-based Graph Partition Prior
text: Now we describe how we model the unknown partition Π. A graph partition can be defined as a collection of disjoint connected subgraphs such that the union of vertices is V . To bypass the need to handle a complex combinatorial graph partition problem, we consider an equivalent formulation of graph partition through edge cuts of spanning forests of G. Prop. 1 guarantees that for any choice of partition Π, there exist a spanning forest F and a corresponding set of cut-edges EC that induce Π, i.e., some edges in F can be removed so that vertices connected to each other form a cluster. Proposition 1. Let G=(V,E) be a graph with nc connected components and Π={C1, . . . , CK} be a graph partition of G. Then there exists a spanning forest F=(V,EF ) with |EF |= |V | − nc, and a set of cut-edges EC ⊂ EF with |EC | = K − nc such that F and EC induce Π.
Proof is deferred to Appendix A2. We will say a spanning forest F is compatible with Π if we can construct Π by cutting some of its edges. See Figure 1(b) for two examples of compatible F and cut-edge(s). Prop.1 suggests that to induce a graph partition prior with full support, it amounts to first assigning a prior model on all possible spanning forests of G, and then assigning a prior model on all possible cut-edges sets conditional on a spanning forest. When G is a connected spatial or spatial temporal graph, Teixeira et al. (2019) considered a discrete uniform prior over all possible spanning trees. But in this case, approximate method is required to sample a spanning tree from its full conditional distribution due to serious inefficiency of the rejection sampler. In contrast, Luo et al. (2021) considered a random minimum spanning tree approach by assigning iid uniform random weights to the edges, so it can generate any arbitrarily given spanning tree of G and enable an exact and efficient posterior conditional sampling algorithm. Thus we follow a similar approach as in Luo et al. (2021) which leads to the following random minimum spanning forest prior on F with full support,
F = MSF(G,W ), Wij iid∼ Unif(0, 1), (i, j) ∈ E, (4)
where MSF(G,W ) denotes the minimum spanning forest of the graph G with edge weightsW .
AfterF is given, selectingK−nc cut-edges forms a partition Π with sizeK. Following approaches of Knorr-Held and Raßer (2000), Feng et al. (2016), and Luo et al. (2021), we introduce a geometrically decaying prior distribution on K, and then select cut-edges uniformly at random given (F ,K). The following prior specification completes the T-LoHo model:
Pr(K = k) ∝ (1− c)k, k = nc, nc + 1, · · · , p, c ∈ [0, 1) (5) p(Π | F ,K) ∝ 1(F is compatible with Π and |Π| = K). (6)
T-LoHo involves two hyperparameters related to model complexity penalization. One is c in (5) which, if selected to be closer to 1, strongly penalizes models with larger numbers of clusters. Another is τ0 in (3) controlling the strength of global shrinkage. As τ0 reduces to 0, the posterior distribution of β̃ tends to concentrate more at zero. More detailed hyperparameter sensitivity analysis and selection criteria are deferred to Appendix A5.

heading: 3 Posterior Inference and Theoretical Properties

heading: 3.1 Posterior Sampler and Computational Strategies
text: Here we briefly describe a reversible-jump Markov chain Monte Carlo algorithm (RJMCMC) (Green, 1995) and discuss computational strategies therein. Denote Θ := (β̃, σ2,Λ, τ,Π,K,F) be the set of parameters, and also denote X̃ := XΦ> so that X̃β̃ = XΦ>Φβ = Xβ since β ∈ rowsp(Φ).
The posterior p(Θ|y) is
p(Θ|y) ∝ Nn(y|X̃β̃, σ2In)×NK(β̃|0, σ2τ2Λ)× 1/σ2 (7) × (1 + τ2)−1 ∏K k=1(1 + λ 2 k) −1 × ( p−nc K−nc )−1 × (1− c)K × 1 (8) where line (8) is the product of priors p(τ) ∏K k=1 p(λk)p(Π|K,F)p(K)p(W ). We draw posterior samples of Θ|y using a collapsed RJMCMC posterior sampler as described in Algorithm 1.
Algorithm 1: One full iteration of the RJMCMC posterior sampler
Step 1. Update Π,K,F using collapsed conditional [Π,K,F|Λ, τ,y] where β̃, σ2 are integrated out.† With probabilities (pa, pb, pc, pd) summing up to 1, perform one of the following substeps:
1-a. (split) Propose (Π?,K?=K+1) compatible with F , and accept with probability min{1,Aa · Pa · La}, where Aa is prior ratio, Pa is proposal ratio, La is likelihood ratio. 1-b. (merge) Propose (Π?,K?=K−1) compatible with F , and accept w.p. min{1,Ab ·Pb ·Lb}. 1-c. (change) Propose (Π?,K?=K) compatible with F , and accept w.p. min{1,Ac · Pc · Lc}. 1-d. (hyper) Update F compatible with current Π.
Step 2. Jointly update (τ, σ2, β̃) from [τ, σ2, β̃ |Λ,Π,K,F ,y], by performing: 2-1. Update τ from [τ |Λ,Π,K,F ,y] using a Metropolis-Hastings sampler, 2-2. Update σ2 from [σ2 | τ,Λ,Π,K,F ,y] with an inverse gamma distribution, 2-3. Update β̃ from [β̃ |σ2, τ,Λ,Π,K,F ,y] with a multivariate normal distribution.
Step 3. Update Λ from [Λ | τ, σ2, β̃,Π,K,F ,y] using a slice sampler. † When X = In (i.e., normal means model), it is possible to integrate out Λ instead of σ2; see Appendix A1.
Step 1 updates cluster assignment by proposing one of the four moves (fig. 2); see Appendix A1 for details on acceptance probabilities. Here instead of using the full conditional distribution, we use the collapsed conditional distribution [Π,K,F|Λ, τ,y] where (β, σ2) are integrated out to calculate the likelihood ratio L which significantly improves mixing. Specifically, the collapsed conditional [Π,K,F|Λ, τ,y] is proportional to |Σ|−1/2(y>Σ−1y/2)−n/2 where Σ = In + τ2X̃ΛX̃> so that
Likelihood ratio L = |Σ ?|−1/2(y>Σ?−1y/2)−n/2
|Σ|−1/2(y>Σ−1y/2)−n/2 . (9)
Here superscript ? indicates the proposed parameters. Step 2 jointly updates [τ, σ2, β̃|−] = [τ |−]× [σ2|τ,−] × [β̃|τ, σ2,−] following the approach of Johndrow et al. (2020). Finally, step 3 updates local shrinkage parameters Λ using the slice sampler (Neal, 2003; Polson et al., 2014).
One demanding computational bottleneck is the likelihood ratio calculation that involves Σ−1 and |Σ|. A notable advantage of T-LoHo is that it projects X to a low-dimensional space under the sparse homogeneity assumption, leading to an n by K transformed design matrix X̃ = XΦ> where typically K << n. Applying Sherman-Woodbury-Morrison formula, we can easily calculate
Σ−1 = (In + τ 2X̃ΛX̃>)−1 = In − X̃(τ−2Λ−1 + X̃>X̃)−1X̃>
by reducing the rank of the inverting matrix from n to K. Furthermore, we utilize Cholesky decomposition R>R = τ−2Λ−1 + X̃>X̃ to efficiently calculate Σ−1 and |Σ|. Once we have right triangular R, calculation of Σ−1 only involves backward substitution, and also we can calculate |Σ| at no extra cost using the matrix determinant lemma. A naive implementation of the above likelihood
calculation would require to repeatedly update τ−2Λ−1 + X̃>X̃ and its Cholesky factor R for each MCMC step, except in steps (2-2) and (2-3) where we use Σ−1 and R calculated from step (2-1) respectively via Rue (2001)’s algorithm. For instance, step 1 has a new X̃?=XΦ?> and steps (2-1) and 3 have a new τ? or Λ?. Alternatively, we propose to directly update R from its previous value. For step 1 where X̃>X̃ changes, we use Cholesky rank-1 update/downdate (Golub and Van Loan, 2013, Sec. 6.5.4) which reduces computational cost from O(max{nK2,K3}) to O(K2). For steps (2-1) and 3 where the diagonal part changes, we present a direct algorithm (see Alg. 2 in Appendix A1) for updating the corresponding Cholesky factor. These Cholesky updating schemes prove to be simple yet powerful computational strategies as MCMC typically requires many iterations.
Excluding step (1-d), the total computational cost of Alg. 1 is O ( max{nK,K3} ) , compared to the direct computation of Σ−1 and |Σ| which leads to O(n3). Step (1-d) takes O(m log p) to construct MSF where m is the number of edges in G. Note that step 1-d (hyper) is selected with probability pd at each iteration. We suggest a small value of pd such as 0.05, so that it would save computation time and give RJMCMC enough iterations to explore the partition compatible with the current MSF F .

heading: 3.2 Clustering Effect of T-LoHo
text: In this section, we investigate how T-LoHo prior differs from the Gaussian prior in terms of the effect on clustering. We focus on the sparsity of edge differences because of its important role in clustering. Examples include `1 (Tibshirani et al., 2011) and `0 penalties (Fan and Guan, 2018) on {βi − βj : (i, j) ∈ E}. Existing Bayesian methods seek sparsity of βi − βj by putting suitable prior distributions on edge differences, such as Laplace (Kyung et al., 2010), normal-exponential-gamma (Shimamura et al., 2019), student’s t (Song and Cheng, 2020; Banerjee and Shen, 2020), spikeand-slab (Kim and Gao, 2020), and horseshoe (Banerjee, 2021). But the aforementioned Bayesian methods have several limitations: (i) an additional post-processing step is often required due to the lack of explicit cluster estimates; (ii) sparsity assumption cannot be easily incorporated into their methods; and (iii) posterior inference method is either inefficient or inflexible (e.g., relying on a single F) when the underlying graph G has many edges. In contrast, although T-LoHo does not put sparsity-inducing priors directly on edge differences, it overcomes all these limitations and effectively finds clusters by its flexible low-rank structure. Below, we show that the use of horseshoe prior not only introduces shrinkage of β but also has a less obvious but profound clustering effect to facilitate homogeneity pursuit.
To analyze the effect of T-LoHo prior on clustering, we focus on the simple case when X = In so that observations yi are independent and normally distributed, i.e., yi|βi, σ2 ∼ N (βi, σ2), i = 1, . . . , n. Without loss of generality, consider the merge step (1-b) where the comparison between the proposed merged modelM1 := (Π?,K − 1) by combining two existing clusters C1, C2 and the current model M2 := (Π,K) is made: M1 : mean of {yi}i∈C1 = µ1 = µ2 = mean of {yi}i∈C2 v.s. M2 : µ1 6= µ2 We analyze the acceptance probability min{1,A · P · L} in step (1-b) of Alg. 1 because it is crucial in the clustering mechanism. The term A · P is nothing but 1/(1− c) which reflects the model size penalty imposed by p(K). The key part is the likelihood ratio L, where the different choice of prior on β (equivalently cluster mean µ) leads to the different L under the same data. This likelihood ratio corresponds to the Bayes factor (Kass and Raftery, 1995) of the Bayesian two-sample t test (Gönen et al., 2005). Thus, here we compare the Bayes factor BF12 = p(data|M1)/p(data|M2) under the normal and T-LoHo prior respectively to analyze their effects on clustering.
Following the formulation of Bayesian two-sample t test, priors are reparametrized as (δ, µ̄, σ2) := ((µ1 − µ2)/σ, (µ1 + µ2)/2, σ2) where the standardized difference δ=(µ1 − µ2)/σ is the parameter of interest. Here we assume a noninformative prior on nuisance parameters p(µ̄, σ2) ∝ 1/σ2, since otherwise Bayes factor is no longer a function of two-sample t statistics t and (n1, n2):
t = (ȳ1 − ȳ2)/(sp/ √ nδ), s 2 p = ( (n1 − 1)s21 + (n2 − 1)s22 ) /ν, nδ = (n −1 1 + n −1 2 ) −1
where nk = |Ck|, ν = n1 +n2−2, and ȳk, s2k are the sample mean and variance of group k, k = 1, 2.
It is obvious that independent normal priors on µk (with variance scaled with σ2) leads to a normal prior on δ as well. When δ ∼ N (0, 1), the Bayes factor BFn12 (Gönen et al., 2005) is
BFn12 = ∫ p(data|δ = 0, µ̄, σ2)p(µ̄, σ2)d(µ̄, σ2)∫ p(data|δ, µ̄, σ2)p(δ, µ̄, σ2)d(δ, µ̄, σ2) =
( 1 + t2/ν )−(ν+1)/2 (1 + nδ)−1/2 {1 + t2/ [ν(1 + nδ)]}−(ν+1)/2 .
Now if we change the priors on µ1 and µ2 from independent normal to independent horseshoe distributions, it induces the heavy-tailed prior π∆ on δ which is a convolution of two horseshoe priors: Proposition 2. Let πHS(µ|σ, τ) = ∫∞ 0 N (µ|0, σ2τ2λ2)C+(λ|0, 1)dλ be a horseshoe prior. If µ1 ∼ πHS(σ, τ1) and µ2 ∼ πHS(σ, τ2) independently, then it induces a distribution of the standardized difference δ := (µ1 − µ2)/σ given τ1, τ2, denoted as π∆(δ|τ1, τ2), which can be written as a scale mixture of normal with mixing distribution fW (w) where w > 0:
π∆(δ|τ1, τ2) = ∫ ∞
0
N (δ|0, w)fW (w)dw, fW (w) = 1
π
τ1 √ w + τ21 + τ2 √ w + τ22√
w + τ21 √ w + τ22 (w + τ 2 1 + τ 2 2 )
Proof is deferred to Appendix A3. See left panel of fig. 3 for the graphical illustration of π∆(δ|1, 1). Its tail behaves similarly with the Strawderman-Berger prior (Strawderman, 1971; Berger et al., 1980). By Proposition 2, the Bayes factor BFhs12 under the prior induced by horseshoe δ ∼ π∆(δ|τ1, τ2) is
BFhs12 =
( 1 + t2/ν )−(ν+1)/2∫∞ 0 (1 + nδw)−1/2 {1 + t2/[ν(1 + nδw)]}−(ν+1)/2 fW (w)dw .
Under the normal and π∆ priors on δ, we compare BF12 as a function of |t| when the group sizes are (i) balanced, n1 = n2; (ii) unbalanced, n1 : n2 = 9 : 1 with ν ∈ {5, 10, 20}, and display results at the right panel of fig. 3. Since an arbitrary choice of scale leads to the different BF, we set δ ∼ N (0, 1) which corresponds to the unit information prior (Kass and Wasserman, 1995) and choose the scale of δ ∼ π∆ such that median of scale mixture distribution fW (w) becomes 1. See Appendix A3 for details of specific choices of (τ1, τ2) under different settings of (n1, n2). Although this is not the only possible scale matching criterion, it is a reasonable choice for a fair BF comparison so that we can focus on different origin/tail behavior of π∆ compared to the normal distribution.
When |t| is small to moderate, we can see that BFhs12 is generally greater than BFn12, except for the case when nδ is large (balanced, large ν) and |t| is small (in which case large BF12 in both models leads to accepting merge proposals anyways). This implies the prior δ ∼ π∆ more strongly favors the one-groupM1 over the two-groupM2 compared to the normal prior δ ∼ N (0, 1). This can be explained by the fact that, underM2, the heavy-tailed prior π∆ anticipates a large effect size a priori compared to the normal prior. Thus T-LoHo facilitates cluster fusion when |t| is small to moderate. Now when |t| is large, we can see that π∆ more strongly favorsM2 compared to the normal prior, and the difference becomes more noticeable when group sizes are small and unbalanced. In fact, when (n1, n2) are fixed, BFhs12 converges to 0 as |t| → ∞ whenever ν ≥ 1 but BFn12 never converges to 0 and is lower bounded by (1 + nδ)−ν/2 > 0, which is also known as the information paradox (Liang et al., 2008). This finite sample consistency is important since it allows the model to identify the (small, unbalanced) cluster with a high signal difference, which might not be possible under the usual normal prior because the penalty term p(K) ∝ (1− c)K may overwhelm it. In summary, compared to the normal prior, T-LoHo tends to reduce redundant cluster representations (when |t| is small) while better capturing the highly significant cluster differences (when |t| is large).

heading: 3.3 Posterior Consistency Results
text: Notations. Let (β∗, β̃∗, σ∗) denote the true β, β̃ and σ, respectively. Let ξ∗ = {j ∈ V : β∗j 6= 0} denote the true active set of indices. Let Π̌ denote an arbitrary partition of V = {1, · · · , p} whose corresponding partition of ξ∗ is determined by removing those edges with β∗i − β∗j > 0 from the subgraph of any F compatible with Π̌ at vertex set ξ∗. We define g∗n = maxΠ̌ |Π̌(ξ∗)| among all possible Π̌. Let Pn denote all unique partitions that have at most g∗n(1 + cδ) clusters and their corresponding partitions of ξ∗ are nested in the true partition of ξ∗ for some constant cδ > 0.
Below, we consider the case when p can be much larger than n and establish posterior concentration results for the T-LoHo model as n goes to infinity. Our results rely on the following assumptions:
(A-1) The graph satisfies g∗n ≺ n/ log p, nc = o(g∗n), and log |Pn| = O(g∗n log p). (A-2) All the covariates are uniformly bounded. There exists some fixed constant λ0 > 0, such
that λmin(X̃T X̃) ≥ nλ0 for any partition in Pn.
(A-3) maxj |β̃∗j |/σ∗ < L, where log(L) = O(log p). (A-4) − log τ = O(log p), τ < p−(2+cτ ) √ g∗n log p/n , 1 − c ≥ p−cα , and
minσ2∈[σ∗2, σ∗2(1+cσε2n)] π(σ 2) > 0 for some positive constants cτ , cα and cσ .
Assumption (A-1) is a regularity condition on G such that the resulting space of spanning forests is not too large. Assumption (A-2) is a commonly adopted condition on design matrix in highdimensional linear regressions. Assumption (A-3) bounds the growth rate of the standardized true coefficients. Assumption (A-4) is a condition on the prior distributions of τ and σ2 as well as choice of hyperparameter c. The proof of Theorem 1 is given in the Appendix A4.
Theorem 1. (Posterior contraction) Under Assumptions (A-1) to (A-4), there exists a large enough constant M1 > 0 and εn √ g∗n log p/n such that the posterior distribution satisfies πn (‖β − β∗‖2 ≥M1σ∗εn | y) ≤ exp(−c1nε2n) with probability 1− exp(−c2nε2n) for some constants c1 > 0 and c2 > 0.

heading: 4 Numerical Examples

heading: 4.1 Simulation Studies
text: We conduct a simulation study to demonstrate the utility of our model. Motivated by the scalaron-image regression problem, we consider a similar setting as Kang et al. (2018). We construct a 30 × 30 lattice graph which represents a structure of 2-D image. Column-standardized image predictors Xi ∈ R900, i = 1, . . . , 100 are generated from mean zero Gaussian process with kernel K(xj , xl) = exp(−djl/ϑ) where djl is the distance between pixels j and l, ϑ is the range parameter with ϑ = 0 indicating no dependence. True coefficient β ∈ R900 is sparse (84% are zero) and has irregular cluster shapes with sharp discontinuities as shown in fig. 4(a). We let SNR ∈ {2, 4} to set error variance σ2 = Var(Xβ)/SNR and generate scalar responses y ∈ R100 by y ∼ N (Xβ, σ2I). We compared our model with the soft-thresholded Gaussian process (STGP, Kang et al., 2018), sparse fused lasso on graph (FL, Tibshirani et al., 2011; Wang et al., 2016), graph OSCAR (GOSCAR, Yang et al., 2012), and Bayesian graph Laplacian (BGL, Chakraborty and Lozano, 2019). We used mean squared prediction error (MSPE) with test set size 1000 to measure the predictive power and Rand index (RI) (Rand, 1971) to measure the clustering accuracy. For T-LoHo and STGP, we collected 4,000 posterior samples after 10,000 burn-in with 10 thin-in rate. For T-LoHo, we calculated the posterior median estimate of β and the cluster point estimate of Π using Dahl (2006)’s method. For STGP, posterior mean estimate of β is used, and Rand index is calculated based on the binary classification (zero/non-zero) using posterior thresholding probabilities. For FL, the two tuning parameters are selected among the fixed candidate set of tuning parameter ratio γFL ∈ {0.2, 1, 5} using the Bayesian information criterion. The T-LoHo and FL results presented here are when (τ0, c) = (1, 0.5) and γFL = 0.2 respectively; the hyperparameter sensitivity analysis of T-LoHo and the detailed settings of other models are available in Appendix A5. All computations were performed on Intel E5-2690 v3 CPU with 128GB of memory.
From fig. 4, we can see that T-LoHo successfully captures the irregular shape of cluster boundaries and sharp discontinuities. STGP gives a much more smoothed estimate (which is expected because it does not assume homogeneity), and FL estimate contains many falsely identified non-zero clusters. The superior performance of T-LoHo over FL is partly attributed to (i) the use of horseshoe that reduces bias in FL, and (ii) the use of RSF-based prior which more efficiently searches non-zero edge differences from a spanning forest instead of a full graph G used in FL. Table 1 shows that in all (ϑ,SNR) settings, T-LoHo indeed outperforms other models in terms of both predictive and clustering accuracy. GOSCAR and BGL give very similar results as their penalty functions have similar octagonal shapes, and both perform poorly in prediction and result in partitions with (nearly) all singletons. It is partly because GOSCAR and BGL allow the coefficients within the same cluster having similar magnitudes but with different signs. Additional simulation studies with more non-zero clusters are available in Appendix A5.

heading: 4.2 Anomaly Detection in Road Networks
text: We applied T-LoHo to the problem of detecting anomalies in road graphs. NYC Pride March is an annual event held in every June in Manhattan, New York. As the march causes traffic congestion along the route, Wang et al. (2016) considered the problem of detecting clusters on the Manhattan road network which have different taxi pickup/dropoff patterns from usual. We constructed the road network graph using GRASS GIS (GRASS Development Team, 2020) and got 3748 nodes (junctions) and 8474 edges (sections of roads). Following Wang et al. (2016), we considered the event held in
2011 (12:00–2:00 pm, June 26th) and processed the number of taxi pickup/dropoff counts data1to the nearest nodes with log transformation. A log baseline seasonal average was calculated from the same time block 12:00–2:00 pm on the same day of the week across the nearest 8 weeks. Then it is plausible to assume that the difference between log counts on event day and log baseline seasonal average has many zeros and clustered patterns over the graph.
We fit T-LoHo and FL and compare the results. For T-LoHo, (τ0, c) = (1, 0.8) were used to collect 5000 posterior samples after 1.5× 105 burn-in with 10 thin-in rate. For the FL estimate, we followed Wang et al. (2016)’s specification where the regularizaton parameter in the fused term was chosen such that the model has 200 degrees of freedom and the one in the sparsity term is set as 0.2.
Figure 5 shows the fitted results of two models. Both successfully capture the decreased taxi activity along the parade route and slightly increased taxi activity around the starting/ending point of the parade. From the subplots, we can see that FL is biased because of soft-thresholding (Rinaldo et al., 2009) while T-LoHo appears to be less biased and capable of producing reasonable uncertainty measures. A notable difference is T-LoHo can capture the decreased taxi activity around the lower Manhattan area while FL cannot due to the bias. In summary, T-LoHo gives better insight on how taxi activity changes when such event occurs.

heading: 5 Concluding Remarks
text: We propose a Tree-based Low-rank Horseshoe (T-LoHo) model to carry out Bayesian inference for graph-structured parameter which is assumed to be sparse and smooth. Accompanied with theoretical grounds and computational strategies, our simulation studies and real data example demonstrate that T-LoHo outperforms other competing methods such as fused lasso in a high-dimensional linear regression context. Extensions to other types of high-dimensional models are possible. In addition, following a similar model construction spirit as T-LoHo, we can build a general class of tree-based low-rank sparse homogeneity model extending other global-local shrinkage priors (Polson and Scott, 2010). Another scenario not addressed in this paper is when we have a weighted graph G = (V,E,W0) as a parameter structure. In this case, incorporating a priori edge weight W0 to T-LoHo is a nontrivial but interesting future research question which might be useful for many possible real data applications. This work does not present any foreseeable societal consequence, but users must be fully aware of the context represented as a graph G when giving interpretation on clustered parameters to avoid any misleading conclusions.
1Publicly available from NYC Taxi & Limousine Commission website, CC0.
Acknowledgments and Disclosure of Funding
The research was partially supported by NSF DMS-1854655, NSF CCF-1934904 and NIH R01AG064010. The authors thank the referees and the area chair for valuable comments.

Review summaries:
The authors generalize the univariate Bayesian horseshoe shrinkage prior to the multivariate setting and propose a new prior that detects sparsity and encourages smoothness within clusters of variables. The paper illustrates the application of the prior in Bayesian regression and present an MCMC algorithm. Numerical experiments with synthetic data and results with one real dataset illustrate the potential advantages over alternative methods.

This paper propose a new prior referred to as Tree-based Low-rank Horseshoe (T-Loho), which can incorporate relationship among variables as a graph. Theoretical characteristics of the clustering effects and posterior concentrations are investigated. In computational experiments based on both synthetic and real-world data, T-Loho turned out to be superior to fused lasso and soft-thresholded Gaussian process in terms of both regression performance and support recovery.

The paper proposes a new prior for high dimensional parameters with graphical relations. It can be regarded as an extension of univariate horseshoe shrinkage prior to multivariate setting by borrowing ideas from studies on homogenous sparse priors and tree-based graph partition prior. The paper describes posterior sampling strategy with efficient implementation details and theoretical guarantees. In the end, the paper demonstrates its methods via both simulated and real data experiments.