Fast Markov chain Monte Carlo for high dimensional Bayesian regression models with shrinkage priors

TL;DR

This paper introduces two-block Gibbs samplers for high-dimensional Bayesian regression with shrinkage priors, demonstrating faster convergence and lower computational costs compared to traditional methods, enabling practical analysis of large datasets.

Contribution

The paper proposes and theoretically justifies two-block Gibbs samplers that outperform existing methods in efficiency and scalability for Bayesian shrinkage models in high-dimensional settings.

Findings

2BG converges faster than 3BG in simulations and real data.

2BG has comparable or lower per-iteration costs than 3BG.

2BG outperforms HMC in efficiency for large p datasets.

Abstract

In the past decade, many Bayesian shrinkage models have been developed for linear regression problems where the number of covariates, $p$, is large. Computing the intractable posterior are often done with three-block Gibbs samplers (3BG), based on representing the shrinkage priors as scale mixtures of Normal distributions. An alternative computing tool is a state of the art Hamiltonian Monte Carlo (HMC) method, which can be easily implemented in the Stan software. However, we found both existing methods to be inefficient and often impractical for large $p$ problems. Following the general idea of Rajaratnam et al. (2018), we propose two-block Gibbs samplers (2BG) for three commonly used shrinkage models, namely, the Bayesian group lasso, the Bayesian sparse group lasso and the Bayesian fused lasso models. We demonstrate with simulated and real data examples that the Markov chains underlying 2BG's converge much faster than that of 3BG's, and no worse than that of HMC. At the same time, the computing costs of 2BG's per iteration are as low as that of 3BG's, and can be several orders of magnitude lower than that of HMC. As a result, the newly proposed 2BG is the only practical computing solution to do Bayesian shrinkage analysis for datasets with large $p$. Further, we provide theoretical justifications for the superior performance of 2BG's. We establish geometric ergodicity (GE) of Markov chains associated with the 2BG for each of the three Bayesian shrinkage models. We also prove, for most cases of the Bayesian group lasso and the Bayesian sparse group lasso model, the Markov operators for the 2BG chains are trace-class. Whereas for all cases of all three Bayesian shrinkage models, the Markov operator for the 3BG chains are not even Hilbert-Schmidt.

Figures

Tables

Table 1. Table 1: Lag-one Autocorrelation and effective sample size per second for the σ 2 superscript 𝜎 2 \sigma^{2} -component of 2BG, 3BG and HMC as applied to the three datasets in section  5 .

			Autocorrelation				$\frac{N_{eff}}{T}$
Data Set	n	p	2BG	3BG	Stan		2BG	3BG	Stan
Gene Expression Data	120	100	0.057	0.40	0.19		578	245	0.08
Economic Data	72	51	0.014	0.42	0.029		1213	494	8.63
CGH Data	200	200	0.288	0.637	0.232		49	19	1.38

Equations

$$\begin{array}{rl}\mathbf{Y}|\mathbf{\beta },{\sigma }^2& \sim {\mathcal{N}}_n\left(\mathbf{X}\mathbf{\beta },{\sigma }^2{\mathbf{I}}_n\right),\\ \mathbf{\beta }|\mathbf{\eta },{\sigma }^2& \sim {\mathcal{N}}_p\left(0,{\sigma }^2{\Sigma }_{\mathbf{\eta }}\right),\\ \mathbf{\eta }& \sim p(\mathbf{\eta }),\\ {\sigma }^2& \sim \text{Inverse-Gamma}(\alpha ,\xi ).\end{array}$$

$$\begin{array}{rl}\mathbf{\eta }|\mathbf{\beta },{\sigma }^2,\mathbf{Y}& \sim \pi \left(\mathbf{\eta }|\mathbf{\beta },{\sigma }^2,\mathbf{Y}\right),\\ {\sigma }^2|\mathbf{\beta },\mathbf{\eta },\mathbf{Y}& \sim \text{Inverse-Gamma}\left(\frac{n+p+2\alpha }{2},\frac{\Vert \mathbf{Y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2+{\mathbf{\beta }}^T{\Sigma }_{\mathbf{\eta }}^{-1}\mathbf{\beta }+2\xi }{2}\right),\\ \mathbf{\beta }|{\sigma }^2,\mathbf{\eta },\mathbf{Y}& \sim {\mathcal{N}}_p\left({\mathbf{A}}_{\mathbf{\eta }}^{-1}{\mathbf{X}}^T\mathbf{Y},{\sigma }^2{\mathbf{A}}_{\mathbf{\eta }}^{-1}\right),\end{array}$$

$$\displaystyle{\color[rgb]{0,0,0}{\tilde{k}\left[(\tilde{\boldsymbol{\beta}}_{0},\tilde{\sigma}_{0}^{2}),(\tilde{\boldsymbol{\beta}}_{1},\tilde{\sigma}_{1}^{2})\right]=\int_{\mathbb{R}^{s}_{+}}\pi(\tilde{\boldsymbol{\beta}}_{1}|\tilde{\boldsymbol{\eta}},\tilde{\sigma}_{1}^{2},\boldsymbol{Y})\pi(\tilde{\sigma}_{1}^{2}|\tilde{\boldsymbol{\beta}}_{0},\tilde{\boldsymbol{\eta}},\boldsymbol{Y})\pi(\tilde{\boldsymbol{\eta}}|\tilde{\boldsymbol{\beta}}_{0},\tilde{\sigma}_{0}^{2},\boldsymbol{Y})d\tilde{\boldsymbol{\eta}}.}}$$

$${\hat{\mathbf{\beta }}}_{\text{group}}=\arg \underset{\mathbf{\beta }}{\min }{\Vert \mathbf{Y}-\sum _{k=1}^K{\mathbf{X}}_{G_k}{\mathbf{\beta }}_{G_k}\Vert }_2^2+\lambda \sum _{k=1}^K{\Vert {\mathbf{\beta }}_{G_k}\Vert }_2,$$

$$\pi (\mathbf{\beta }|{\sigma }^2)\propto \exp \left(-\frac{\lambda }{\sigma }\sum _{k=1}^K\Vert {\mathbf{\beta }}_{G_k}{\Vert }_2\right).$$

$$\begin{array}{rl}\mathbf{Y}|\mathbf{\beta },{\sigma }^2& \sim {\mathcal{N}}_n(\mathbf{X}\mathbf{\beta },{\sigma }^2{\mathbf{I}}_n),\\ {\mathbf{\beta }}_{G_k}|{\tau }_k^2,{\sigma }^2& \stackrel{ind}{\sim }\mathcal{N}(0_{m_k},{\sigma }^2{\tau }_k^2{\mathbf{I}}_{m_k}),k=1,\dots ,K,\\ {\tau }_k^2& \stackrel{ind}{\sim }\text{Gamma}\left(\frac{m_k+1}{2},\frac{{\lambda }^2}{2}\right),k=1,\dots ,K,\\ {\sigma }^2& \sim \text{Inverse-Gamma}(\alpha ,\xi ).\end{array}$$

$${\mathbf{D}}_{\mathbf{\tau }}=\text{diag}(\underset{m_1}{\underbrace{{\tau }_1^2,\dots ,{\tau }_1^2}},\underset{m_2}{\underbrace{{\tau }_2^2,\dots ,{\tau }_2^2}},\dots ,\underset{m_K}{\underbrace{{\tau }_K^2,\dots ,{\tau }_K^2}}).$$

$${\tilde{k}}_{gl}\left[({\tilde{\mathbf{\beta }}}_0,{\tilde{\sigma }}_0^2),({\tilde{\mathbf{\beta }}}_1,{\tilde{\sigma }}_1^2)\right]={\int }_{{\mathbb{R}}_{+}^K}\pi ({\tilde{\mathbf{\beta }}}_1|{\tilde{\sigma }}_1^2,{\tilde{\mathbf{\tau }}}^2,\mathbf{Y})\pi ({\tilde{\sigma }}_1^2|{\tilde{\mathbf{\beta }}}_0,{\tilde{\mathbf{\tau }}}^2,\mathbf{Y})\pi ({\tilde{\mathbf{\tau }}}^2|{\tilde{\mathbf{\beta }}}_0,{\tilde{\sigma }}_0^2,\mathbf{Y})d{\tilde{\mathbf{\tau }}}^2.$$

$${\hat{\mathbf{\beta }}}_{\text{sgroup}}=\arg \underset{\mathbf{\beta }}{\min }{\Vert \mathbf{Y}-\sum _{k=1}^K{\mathbf{X}}_{G_k}{\mathbf{\beta }}_{G_k}\Vert }_2^2+{\lambda }_1\sum _{k=1}^K{\Vert {\mathbf{\beta }}_{G_k}\Vert }_2+{\lambda }_2{\Vert \mathbf{\beta }\Vert }_1,$$

$$\pi (\mathbf{\beta }|{\sigma }^2)\propto \exp \left(-\frac{{\lambda }_1}{\sigma }\sum _{k=1}^K\Vert {\mathbf{\beta }}_{G_k}{\Vert }_2-\frac{{\lambda }_2}{\sigma }\Vert \mathbf{\beta }{\Vert }_1\right).$$

$${\mathbf{V}}_k=\text{diag}\left\{{\left(\frac{1}{{\tau }_k^2}+\frac{1}{{\gamma }_{k,j}^2}\right)}^{-1}\text{vspace}.2cm:j=1,\dots ,m_k\right\},$$

$${\pi }_k\propto \prod _{j=1}^{m_k}\left[({\gamma }_{k,j}^2{)}^{-\frac{1}{2}}{\left(\frac{1}{{\gamma }_{k,j}^2}+\frac{1}{{\tau }_k^2}\right)}^{-\frac{1}{2}}\right]({\tau }_k^2{)}^{-\frac{1}{2}}\exp \left(-\frac{{\lambda }_2^2}{2}\sum _{j=1}^{m_k}{\gamma }_{k,j}^2-\frac{{\lambda }_1^2}{2}{\tau }_k^2\right).$$

$$\begin{array}{rl}\mathbf{Y}|\mathbf{\beta },{\sigma }^2& \sim {\mathcal{N}}_n(\mathbf{X}\mathbf{\beta },{\sigma }^2{\mathbf{I}}_n),\\ {\mathbf{\beta }}_{G_k}|{\mathbf{\tau }}^2,{\mathbf{\gamma }}^2,{\sigma }^2& \stackrel{ind}{\sim }\mathcal{N}(0_{m_k},{\sigma }^2{\mathbf{V}}_k),k=1,\dots ,K,\\ ({\gamma }_{k,1}^2,\dots ,{\gamma }_{k,m_k}^2,{\tau }_k^2)& \stackrel{ind}{\sim }{\pi }_k,k=1,\dots ,K,\\ {\sigma }^2& \sim \text{Inverse-Gamma}(\alpha ,\xi ),\end{array}$$

$${\tilde{k}}_{sgl}\left[({\tilde{\mathbf{\beta }}}_0,{\tilde{\sigma }}_0^2),({\tilde{\mathbf{\beta }}}_1,{\tilde{\sigma }}_1^2)\right]={\int }_{{\mathbb{R}}_{+}^K}{\int }_{{\mathbb{R}}_{+}^p}\pi ({\tilde{\mathbf{\beta }}}_1|{\tilde{\sigma }}_1^2,{\tilde{\mathbf{\tau }}}^2,{\tilde{\mathbf{\gamma }}}^2,\mathbf{Y})\pi ({\tilde{\sigma }}_1^2|{\tilde{\mathbf{\beta }}}_0,{\tilde{\mathbf{\tau }}}^2,{\tilde{\mathbf{\gamma }}}^2,\mathbf{Y})\pi ({\tilde{\mathbf{\tau }}}^2,{\tilde{\mathbf{\gamma }}}^2|{\tilde{\mathbf{\beta }}}_0,{\tilde{\sigma }}_0^2,\mathbf{Y})d{\tilde{\mathbf{\gamma }}}^{2}d{\tilde{\mathbf{\tau }}}^2.$$

$${\hat{\mathbf{\beta }}}_{\text{fused}}=\arg \underset{\mathbf{\beta }}{\min }\Vert \mathbf{Y}-\mathbf{X}\mathbf{\beta }{\Vert }_2^2+{\lambda }_1\Vert \mathbf{\beta }{\Vert }_1+{\lambda }_2\sum _{j=1}^{p-1}|{\beta }_{j+1}-{\beta }_j|.$$

$$\pi (\mathbf{\beta }|{\sigma }^2)\propto \exp \left(-\frac{{\lambda }_1}{\sigma }\Vert \mathbf{\beta }{\Vert }_1-\frac{{\lambda }_2}{\sigma }\sum _{j=1}^{p-1}|{\beta }_{j+1}-{\beta }_j|\right).$$

$${\mathbf{\Sigma }}_{\tau ,\omega }^{-1}=\left[\begin{array}{cccc}\frac{1}{{\tau }_1^2}+\frac{1}{{\omega }_1^2}& -\frac{1}{{\omega }_1^2}& \cdots & 0\\ -\frac{1}{{\omega }_1^2}& \frac{1}{{\tau }_2^2}+\frac{1}{{\omega }_2^2}+\frac{1}{{\omega }_1^2}& \cdots & 0\\ 0& -\frac{1}{{\omega }_2^2}& \cdots & 0\\ \cdots & \cdots & \ddots & \cdots \\ 0& 0& \frac{1}{{\tau }_{p-1}^2}+\frac{1}{{\omega }_{p-2}^2}+\frac{1}{{\omega }_{p-1}^2}& -\frac{1}{{\omega }_{p-1}^2}\\ 0& 0& -\frac{1}{{\omega }_{p-1}^2}& \frac{1}{{\tau }_p^2}+\frac{1}{{\omega }_{p-1}^2}\end{array}\right].$$

$$\begin{split}\boldsymbol{Y}|\boldsymbol{\beta},\sigma^{2}&\sim{\cal N}_{n}(\boldsymbol{X}\boldsymbol{\beta},\sigma^{2}\boldsymbol{I}_{n}),\\ \boldsymbol{\beta}|\boldsymbol{\tau}^{2},\boldsymbol{\omega}^{2},\sigma^{2}&\sim{\cal N}_{p}(\mathbf{0}_{p},\sigma^{2}\boldsymbol{\Sigma}_{\tau,\omega}),\\ {\color[rgb]{0,0,0}{\pi(\boldsymbol{\tau}^{2},\boldsymbol{\omega}^{2})}}&{\color[rgb]{0,0,0}{\propto}}{\color[rgb]{0,0,0}{\text{ det}(\boldsymbol{\Sigma}_{\tau,\omega})^{1/2}}}{\color[rgb]{0,0,0}{\left(\prod_{i=1}^{p}(\tau_{i}^{2})^{-1/2}e^{-\lambda_{1}\tau_{i}^{2}/2}\right)\left(\prod_{i=1}^{p-1}(\omega_{i}^{2})^{-1/2}e^{-\lambda_{2}\omega_{i}^{2}/2}\right)}}\\ \sigma^{2}&\sim\text{Inverse- Gamma}(\alpha,\xi),\end{split}$$

$${\tilde{k}}_{fl}\left[({\tilde{\mathbf{\beta }}}_0,{\tilde{\sigma }}_0^2),({\tilde{\mathbf{\beta }}}_1,{\tilde{\sigma }}_1^2)\right]={\int }_{{\mathbb{R}}_{+}^p}{\int }_{{\mathbb{R}}_{+}^{p-1}}\pi ({\tilde{\mathbf{\beta }}}_1|{\tilde{\sigma }}_1^2,{\tilde{\mathbf{\tau }}}^2,{\tilde{\mathbf{\omega }}}^2,\mathbf{Y})\pi ({\tilde{\sigma }}_1^2|{\tilde{\mathbf{\beta }}}_0,{\tilde{\mathbf{\tau }}}^2,{\tilde{\mathbf{\omega }}}^2,\mathbf{Y})\pi ({\tilde{\mathbf{\tau }}}^2,{\tilde{\mathbf{\omega }}}^2|{\tilde{\mathbf{\beta }}}_0,{\tilde{\sigma }}_0^2,\mathbf{Y})d{\tilde{\mathbf{\omega }}}^{2}d{\tilde{\mathbf{\tau }}}^2.$$

$$\begin{split}\boldsymbol{\eta}|\boldsymbol{\beta},\sigma^{2},\boldsymbol{Y}&\sim\pi(\boldsymbol{\eta}|\boldsymbol{\beta},\sigma^{2},\boldsymbol{Y}),\\ (\sigma^{2},\boldsymbol{\beta})|{\boldsymbol{\eta}},\boldsymbol{Y}&\begin{cases}\sigma^{2}|\boldsymbol{\eta},\boldsymbol{Y}\sim\text{Inverse-Gamma}\left(\frac{n}{2}+\alpha,\frac{\boldsymbol{Y}^{T}(\mathbf{I}_{n}-\boldsymbol{X}(\boldsymbol{X}^{T}\boldsymbol{X}+\boldsymbol{{\color[rgb]{0,0,0}{\Sigma}}}_{\eta}^{-1})^{-1}\boldsymbol{X}^{T})\boldsymbol{Y}}{2}+\xi\right),\\ &\\ \boldsymbol{\beta}|\sigma^{2},\boldsymbol{\eta},\boldsymbol{Y}\sim{\cal N}_{p}\left(\left(\boldsymbol{X}^{T}\boldsymbol{X}+\boldsymbol{{\color[rgb]{0,0,0}{\Sigma}}}_{\eta}^{-1}\right)^{-1}\boldsymbol{X}^{T}\boldsymbol{Y},\sigma^{2}\left(\boldsymbol{X}^{T}\boldsymbol{X}+\boldsymbol{{\color[rgb]{0,0,0}{\Sigma}}}_{\eta}^{-1}\right)^{-1}\right).\end{cases}\end{split}$$

$$k\left[({\mathbf{\beta }}_0,{\sigma }_0^2),({\mathbf{\beta }}_1,{\sigma }_1^2)\right]={\int }_{{\mathbb{R}}_{+}^s}\pi ({\mathbf{\beta }}_1|{\sigma }_1^2,\mathbf{\eta },\mathbf{Y})\pi ({\sigma }_1^2|\mathbf{\eta },\mathbf{Y})\pi (\mathbf{\eta }|{\mathbf{\beta }}_0,{\sigma }_0^2,\mathbf{Y})d\mathbf{\eta }.$$

$$\begin{array}{rl}\frac{1}{{\tau }_k^2}|\mathbf{\beta },{\sigma }^2,\mathbf{Y}& \stackrel{ind}{\sim }\text{Inverse-Gaussian}\left(\sqrt{\frac{{\lambda }^2{\sigma }^2}{\Vert {\mathbf{\beta }}_{G_K}{\Vert }_2^2}},{\lambda }^2\right),\text{ for }k=1,\dots ,K,\\ ({\sigma }^2,\mathbf{\beta })|{\mathbf{\tau }}^2,\mathbf{Y}& \{\begin{array}{l}{\sigma }^2|{\mathbf{\tau }}^2,\mathbf{Y}\sim \text{Inverse-Gamma}\left(\frac{n}{2}+\alpha ,\frac{{\mathbf{Y}}^T({\mathbf{I}}_n-\mathbf{X}({\mathbf{X}}^T\mathbf{X}+{\mathbf{D}}_{\tau }^{-1}{)}^{-1}{\mathbf{X}}^T)\mathbf{Y}}{2}+\xi \right),\\ \\ \mathbf{\beta }|{\sigma }^2,{\mathbf{\tau }}^2,\mathbf{Y}\sim {\mathcal{N}}_p\left({\left({\mathbf{X}}^T\mathbf{X}+{\mathbf{D}}_{\tau }^{-1}\right)}^{-1}{\mathbf{X}}^T\mathbf{Y},{\sigma }^2{\left({\mathbf{X}}^T\mathbf{X}+{\mathbf{D}}_{\tau }^{-1}\right)}^{-1}\right).\end{array}\end{array}$$

$$k_{gl}\left[({\mathbf{\beta }}_0,{\sigma }_0^2),({\mathbf{\beta }}_1,{\sigma }_1^2)\right]={\int }_{{\mathbb{R}}_{+}^K}\pi ({\mathbf{\beta }}_1|{\sigma }_1^2,{\mathbf{\tau }}^2,\mathbf{Y})\pi ({\sigma }_1^2|{\mathbf{\tau }}^2,\mathbf{Y})\pi ({\mathbf{\tau }}^2|{\mathbf{\beta }}_0,{\sigma }_0^2,\mathbf{Y})d{\mathbf{\tau }}^2.$$

$$\begin{array}{rl}\frac{1}{{\tau }_k^2}|\mathbf{\beta },{\sigma }^2,\mathbf{Y}& \stackrel{ind}{\sim }\text{Inverse-Gaussian}\left(\sqrt{\frac{{\lambda }_1^2{\sigma }^2}{\Vert {\mathbf{\beta }}_{G_K}{\Vert }_2^2}},{\lambda }_1^2\right),\text{ for k=1,…,K, and independently, }\\ \frac{1}{{\gamma }_{k,j}^2}|\mathbf{\beta },{\sigma }^2,\mathbf{Y}& \stackrel{ind}{\sim }\text{Inverse-Gaussian}\left(\sqrt{\frac{{\lambda }_2^2{\sigma }^2}{{\beta }_{k,j}^2}},{\lambda }_2^2\right),\text{ for }k=1,\dots ,K\text{ and }j=1,\dots ,m_k,\\ ({\sigma }^2,\mathbf{\beta })|{\mathbf{\tau }}^2,{\mathbf{\gamma }}^2,\mathbf{Y}& \sim \{\begin{array}{ll}{\sigma }^2|{\mathbf{\tau }}^2,{\mathbf{\gamma }}^2,\mathbf{Y}\sim \text{Inverse-Gamma}\left(\frac{n}{2}+\alpha ,\frac{{\mathbf{Y}}^T({\mathbf{I}}_n-\mathbf{X}({\mathbf{X}}^T\mathbf{X}+{\mathbf{V}}_{\tau ,\gamma }^{-1}{)}^{-1}{\mathbf{X}}^T)\mathbf{Y}}{2}+\xi \right),& \\ & \\ \mathbf{\beta }|{\sigma }^2,{\mathbf{\tau }}^2,{\mathbf{\gamma }}^2,\mathbf{Y}\sim {\mathcal{N}}_p\left({\left({\mathbf{X}}^T\mathbf{X}+{\mathbf{V}}_{\tau ,\gamma }^{-1}\right)}^{-1}{\mathbf{X}}^T\mathbf{Y},{\sigma }^2{\left({\mathbf{X}}^T\mathbf{X}+{\mathbf{V}}_{\tau ,\gamma }^{-1}\right)}^{-1}\right).& \end{array}\end{array}$$

$$k_{sgl}\left[({\mathbf{\beta }}_0,{\sigma }_0^2),({\mathbf{\beta }}_1,{\sigma }_1^2)\right]={\int }_{{\mathbb{R}}_{+}^K}{\int }_{{\mathbb{R}}_{+}^p}\pi ({\mathbf{\beta }}_1|{\sigma }_1^2,{\mathbf{\tau }}^2,{\mathbf{\gamma }}^2,\mathbf{Y})\pi ({\sigma }_1^2|{\mathbf{\tau }}^2,{\mathbf{\gamma }}^2,\mathbf{Y})\pi ({\mathbf{\tau }}^2,{\mathbf{\gamma }}^2|{\mathbf{\beta }}_0,{\sigma }_0^2,\mathbf{Y})d{\mathbf{\tau }}^{2}d{\mathbf{\gamma }}^2.$$

$$\begin{array}{rl}\frac{1}{{\tau }_j^2}|\mathbf{\beta },{\sigma }^2,\mathbf{Y}& \stackrel{ind}{\sim }\text{Inverse-Gaussian}\left(\sqrt{\frac{{\lambda }_1^2{\sigma }^2}{{\beta }_j^2}},{\lambda }_1^2\right),\text{ for j=1,…,p, and independently,}\\ \frac{1}{{\omega }_j^2}|\mathbf{\beta },{\sigma }^2,\mathbf{Y}& \stackrel{ind}{\sim }\text{Inverse-Gaussian}\left(\sqrt{\frac{{\lambda }_2^2{\sigma }^2}{({\beta }_{j+1}-{\beta }_j{)}^2}},{\lambda }_2^2\right),\text{ for }j=1,\dots ,p-1\\ ({\sigma }^2,\mathbf{\beta })|{\mathbf{\tau }}^2,{\mathbf{\omega }}^2,\mathbf{Y}& \{\begin{array}{ll}{\sigma }^2|{\mathbf{\tau }}^2,{\mathbf{\omega }}^2,\mathbf{Y}\sim \text{Inverse-Gamma}\left(\frac{n}{2}+\alpha ,\frac{{\mathbf{Y}}^T({\mathbf{I}}_n-\mathbf{X}({\mathbf{X}}^T\mathbf{X}+{\mathbf{\Sigma }}_{\tau ,\omega }^{-1}{)}^{-1}{\mathbf{X}}^T)\mathbf{Y}}{2}+\xi \right),& \\ & \\ \mathbf{\beta }|{\sigma }^2,{\mathbf{\tau }}^2,{\mathbf{\omega }}^2,\mathbf{Y}\sim {\mathcal{N}}_p\left({\left({\mathbf{X}}^T\mathbf{X}+{\mathbf{\Sigma }}_{\tau ,\omega }^{-1}\right)}^{-1}{\mathbf{X}}^T\mathbf{Y},{\sigma }^2{\left({\mathbf{X}}^T\mathbf{X}+{\mathbf{\Sigma }}_{\tau ,\omega }^{-1}\right)}^{-1}\right),& \end{array}\end{array}$$

$$k_{fl}\left[({\mathbf{\beta }}_0,{\sigma }_0^2),({\mathbf{\beta }}_1,{\sigma }_1^2)\right]={\int }_{{\mathbb{R}}_{+}^p}{\int }_{{\mathbb{R}}_{+}^{p-1}}\pi ({\mathbf{\beta }}_1|{\sigma }_1^2,{\mathbf{\tau }}^2,{\mathbf{\omega }}^2,\mathbf{Y})\pi ({\sigma }_1^2|{\mathbf{\tau }}^2,{\mathbf{\omega }}^2,\mathbf{Y})\pi ({\mathbf{\tau }}^2,{\mathbf{\omega }}^2|{\mathbf{\beta }}_0,{\sigma }_0^2,\mathbf{Y})d{\mathbf{\omega }}^{2}d{\mathbf{\tau }}^2.$$

$$N_{\text{eff}}=\frac{N}{1+2{\sum }_{k=1}^{\infty }{\rho }_k}.$$

$$\mathbf{Y}=\mathbf{X}{\mathbf{\beta }}_{\ast }+\mathbf{ϵ},$$

Full text

Key words and phrases:  Bayesian fused lasso, Bayesian group lasso, Geometric ergodicity, Gibbs sampler, Hamilton Monte Carlo, trace class.

1 Introduction

A linear regression model describes how a response variable $y$ is related to potential explanatory variables $x_{1},\cdots,x_{p}$ through a vector of regression coefficients, $\boldsymbol{\beta}=(\beta_{1},\cdots,\beta_{p})$. Sparse estimates of $\boldsymbol{\beta}$ are often desired to achieve simple models and stable predictions, especially when the number of explanatory variables $p$ is larger than the sample size $n$. One popular set of tools is the penalized regression method, such that, with proper choice of the penalty function, some $\beta_{j}$’s will be estimated to be zero. See, e.g., Fan and Li (2001). However, it’s generally challenging to quantify the uncertainty of the coefficient estimates and that of the resulting sparse models, though there are important recent works in this direction. See, for e.g., Chatterjee and Lahiri (2011); Lockhart et al. (2014); Tibshirani et al. (2016).

Alternatively, Bayesian models quantify uncertainties of parameters and latent variables through posterior distributions. Two major types of Bayesian regression models are used to induce sparsity. First, the spike-and-slab models assign mixtures of the degenerate distribution at zero and continuous distributions as priors for the $\beta_{j}$’s, to explore exactly sparse solutions (Mitchell and Beauchamp, 1988; George and McCulloch, 1993). Such models are natural for sparse problems and enjoy optimality in various senses (Johnstone and Silverman, 2004; Castillo et al., 2015). But both exact and approximate computing are notoriously difficult for even a moderate number of predictors, as an exploration of $2^{p}$ sub-models is needed essentially. Models of a second kind, which we will call the Bayesian shrinkage models, assign the $\beta_{j}$’s continuous prior distributions. Bayesian shrinkage models have gained popularity as their approximate computation are more easily achieved, and many such models lead to good inference results. Some notable examples include the horseshoe priors (Carvalho et al., 2010), the generalized double Pareto priors (Armagan et al., 2013), and the Dirichlet-Laplace priors (Bhattacharya et al., 2015).

Bayesian shrinkage methods have been tailored to different regression setups. More often than not, predictors have known structures, such as they form groups, or are ordered in some meaningful way. The Bayesian group lasso and the Bayesian sparse group lasso models are designed for the first situation, while the Bayesian fused lasso model is for the latter. See Xu and Ghosh (2015); Lee and Chen (2015); Kyung et al. (2010) for their introduction into the Bayesian framework, and see Yuan and Lin (2006); Huang et al. (2012); Simon et al. (2013); Tibshirani et al. (2005) for their origins in the penalized regression literature. These Bayesian shrinkage models have been applied to diverse areas such as agriculture, genetics, image processing and ecology. Examples can be found in Hefley et al. (2017); Li et al. (2015); Raman et al. (2009); Fan and Peng (2016); Hooten and Hobbs (2015); Greenlaw et al. (2017) and Zhang et al. (2014). Successful applications of the Bayesian group lasso, the Bayesian sparse group lasso and the Baysian fused lasso models depend on fast and reliable computing, which is the pursuit of this paper.

Since there are usually no closed form expressions for posterior quantities of Bayesian shrinkage models, one would resort to Markov chain Monte Carlo (MCMC) as a most reliable method for approximations. Many shrinkage priors, including the three that we will focus on, enjoy a hierarchical representation. This allows an MCMC scheme called the three-block Gibbs sampler (3BG) to explore the posterior. However, 3BG’s often suffer from slow convergence due to high correlation between components of the different blocks, especially that between the regression coefficients in one block and the variance parameters in another (Rajaratnam and Sparks, 2015). An alternative MCMC solution is to use Hamiltonian Monte Carlo (HMC), a general sampling method for exploring continuous distributions. We are not aware of previous usage of HMC to compute the aforementioned Bayesian shrinkage models, but this is easy to set up using the Stan (Carpenter et al., 2017) software, with no need of the hierarchical representation of the prior. Stan implements a certain leapfrog variant of HMC, and is widely-known for producing nearly independent samples in numerous empirical examples, for which standard MCMC algorithms like random walk Metropolis do not mix well. In contrast to the empirical success of HMC, there is limited understanding of its theoretical properties, such as the rate of convergence to the target distribution, and whether the resulting Monte Carlo estimators have finite asymptotic variances (Diaconis et al., 2014). Indeed, performance of the HMC heavily depends on the choice of tuning parameters including the number of leapfrog steps and the step sizes. A rare attempt to identify conditions under which certain variants of HMC converge at a geometric rate is Livingstone et al. (2016), but the conditions therein are non-trivial to check and do not concern exactly the type of HMC implemented in Stan. Despite the lack of theoretical guarantees, we run HMC for the Bayesian shrinkage models in sections 4 and 5. It turns out that, in high-dimensional regression problems with correlated covariates, the computing cost of each iteration of HMC can be very high, while the resulting chains still exhibit moderate to high autocorrelations.

So far, we mentioned that for the three Bayesian shrinkage models of interest, both existing computing solutions, 3BG and HMC, are far from satisfactory. We develop and study two-block Gibbs samplers (2BG) that are reliable and fast in high dimensional cases, where these models are most needed. Empirically, we demonstrate that the proposed 2BG greatly outperforms both 3BG and HMC in computing efficiency, notably with multiple orders of magnitude improvements in challenging cases where $p$ exceeds $n$. Here, computing efficiency is measured as the effective sample size per unit time, where the effective sample size is inversely proportional to the square root of the asymptotic variance of the Monte Carlo estimator of interest. Looking further into what leads to the improvement, we see that the proposed 2BG mixes much better than 3BG in terms of having lag-one autocorrelations closer to 0, with similar computing cost per iteration; and 2BG mixes similarly to or better than HMC, but costs only a small fraction of the latter per iteration.

Beyond showing the empirical success of 2BG’s for computing the three Bayesian shrinkage models of interest, we prove in theory that they are reliable, and have superior spectral properties than the corresponding 3BG’s. In more details, there are two commonly used criteria to evaluate Markov chains that explore a target distribution, disregarding computing cost: the rate of convergence (in total variation distance) to the target, and the asymptotic variance of Monte Carlo estimators of interest. Concerning convergence rates, it is desirable that a Markov chain has geometric ergodicity (GE), a most regularly used condition to establish the central limit theorem (CLT), which allows standard error evaluations of Monte Carlo estimates. For the proposed 2BG for each of the three variants of the Bayesian lasso model, we establish GE of the respective Markov chain. In addition, both the convergence rate of a Markov chain and asymptotic variances of Monte Carlo estimators depend on the spectrum of the corresponding Markov operator. See, e.g., Rosenthal et al. (2015) and Mira (2001). For the Bayesian group lasso and the Bayesian sparse group lasso models, we are able to prove for almost all cases that the Markov operators associated with 2BG are trace-class whereas that of 3BG are not even Hilbert-Schmidt. According to the definition of these properties that we provide in section 2 of the supplement, this means that the (self-adjoint) 2BG Markov operator has eigenvalues that are absolutely summable, and hence square-summable (as the eigenvalues are bounded above by $1$ in absolute value). While for the 3BG chain, the absolute value of its operator has non-square-summable eigenvalues. In other words, the 2BG Markov operators have qualitative smaller spectrum than that of 3BG, and the former are expected to enjoy faster convergence rate, which agree with our empirical observations.

Another benefit of showing the trace class results for 2BG is to provide theoretical foundations for further acceleration of these algorithms. Indeed, 2BG is a kind of Data Augmentation (DA) algorithm. There is a well-studied “sandwich” trick (Liu and Wu, 1999; Meng and Van Dyk, 1999; Hobert and Marchev, 2008) that allows many ways to insert an extra step in each iteration of a DA, to achieve potentially much faster convergence rate. The sandwich trick will improve the computational efficiency of DA provided the extra step is inexpensive, and the trick has seen many successful implementations (Liu and Wu, 1999; Meng and Van Dyk, 1999; Van Dyk and Meng, 2001; Marchev and Hobert, 2004; Hobert et al., 2011; Ghosh and Tan, 2015). The design of good sandwich algorithms for the Bayesian shrinkage models is beyond the scope of this paper. Nevertheless, our results that certain 2BG samplers are trace class imply that any sandwich improvements of them are also trace-class and have better spectral properties, in the sense that eigenvalues of the improved algorithms are dominated by that of the basic 2BG (Khare and Hobert, 2011).

The 2BG algorithms proposed and analyzed in this paper stem from the work of Rajaratnam et al. (2018), which introduced the two-block strategy for a class of Bayesian models for regression. In the case of the Bayesian lasso model, they analyzed the convergence rates and the spectral properties for the corresponding 2BG and 3BG, and acknowledged that such a theoretical study is “a big and challenging undertaking”. A main contribution of our work is to develop these theoretical results for the aforementioned three variants of the Bayesian lasso models, which complement and support the empirical advantage we observed of using 2BG over other existing computing solutions.

The rest of the paper is organized as follows. In section 2, we present the three Bayesian shrinkage models and review the existing 3BG algorithms that explore their posterior distributions. In section 3, we propose the 2BG algorithms for these models, establish geometric ergodicity of the underlying Markov chains and show they have superior spectral properties than their 3BG counterparts. In sections 4 and 5, we perform empirical comparisons among HMC, 3BG and 2BG in simulation studies and real data applications.

2 Bayesian Shrinkage Models and Standard Three-block Gibbs Samplers

Let $\boldsymbol{Y}\in\mathbb{R}^{n}$ be the response vector, $\boldsymbol{X}$ be the $n\times p$ design matrix, $\boldsymbol{\beta}\in\mathbb{R}^{p}$ be the vector of regression coefficients and $\sigma^{2}>0$ be the residual variance in a regression problem. We consider a general framework of Bayesian models as follows:

[TABLE]

Here, $\Sigma_{\boldsymbol{\eta}}$ is a $p\times p$ positive definite matrix indexed by $\boldsymbol{\eta}$, with a proper prior $p(\boldsymbol{\eta})$ on $\boldsymbol{\eta}$. The hyperparameters $\alpha$ and $\xi$ are non-negative real numbers specified by users. Many Bayesian shrinkage models can be represented using (2.1), e.g., those under the lasso prior (Park and Casella, 2008), the enet prior (Li and Lin, 2010), the horseshoe prior (Carvalho et al., 2010), and the Dirichlet–Laplace prior (Bhattacharya et al., 2015). Bayesian inference of the model is based on the joint posterior distribution of $(\boldsymbol{\beta},\sigma^{2},\boldsymbol{\eta})$, especially the $(\boldsymbol{\beta},\sigma^{2})-$marginal and various summaries of it. These distributions are analytically intractable and difficult to draw independent samples from. The focus of this paper is to find efficient MCMC algorithms to sample from the joint posterior distribution of $(\boldsymbol{\beta},\sigma^{2})$ for the model in (2.1). An alternative solution through variational approximation is usually of much lower accuracy (see, e.g., Neville et al. (2014)), and not discussed here.

Note that the model in (2.1) yields the following conditional distributions:

[TABLE]

where $\boldsymbol{A}_{\boldsymbol{\eta}}=\boldsymbol{X}^{T}\boldsymbol{X}+\Sigma_{\boldsymbol{\eta}}^{-1}$. Hence, one can use a Gibbs sampler to update the three blocks of components, $\boldsymbol{\eta}$, $\sigma^{2}$ and $\boldsymbol{\beta}$ in turn, as long as $\pi\left(\boldsymbol{\eta}|\boldsymbol{\beta},\sigma^{2},\boldsymbol{Y}\right)$ is easy to sample from. Such Gibbs samplers have been adopted in Kyung et al. (2010); Griffin and Brown (2010); Pal et al. (2017) and Li and Lin (2010) to study special cases of model (2.1). Note that $(\boldsymbol{\beta},\sigma^{2})$ are the parameters of main interest in the Bayesian regression models, and the $(\boldsymbol{\beta},\sigma^{2})$-marginal of the Gibbs sampler, $\tilde{\boldsymbol{\Phi}}:=\left\{(\tilde{\boldsymbol{\beta}}_{m},\tilde{\sigma}_{m}^{2})\right\}_{m=0}^{\infty}$ is itself a Markov chain on the state space $\mathbb{R}^{p}\times\mathbb{R}_{+}$. This Markov chain is what we will refer to as the 3BG in this paper, with Markov transition density:

[TABLE]

We next review three special cases of model (2.1) that are useful for different regression setups.

2.1 The Bayesian Group Lasso

In regression problems, explanatory variables often occur in groups. Grouping structures can be implied by knowledge of the area of application. For instance, in modeling health of individuals, variables may belong to different categories such as social and economic environment, physical environment, and individual characteristics. Also, quantitative variables may enter a model as groups of basis functions, and qualitative variables with multiple levels could be represented by groups of dummy variables. See Huang et al. (2012) for an extensive review of problems that involve groups of variables. Suppose there are $K$ groups of explanatory variables. For $k=1,\ldots,K$, let $m_{k}$ denote the size of the $k$th group, $\boldsymbol{\beta}_{G_{k}}$ the $m_{k}$-dim vector of regression coefficients, and $\boldsymbol{X}_{G_{k}}$ the corresponding sub-design matrices. Yuan and Lin (2006) presented the group lasso method for estimating $\boldsymbol{\beta}$ as follows,

[TABLE]

where the penalty parameter $\lambda$ decides the intensity of the shrinkage and is often tuned by cross-validation. An obvious Bayesian analogue of the above method can be obtained by combining the first line of (2.1) and the prior

[TABLE]

Kyung et al. (2010) showed a hierarchical representation of (2.5), and the model is completed with a prior on $\sigma^{2}$:

[TABLE]

Here, $\alpha,\xi\geq 0$ and $\lambda>0$ are user specified hyperparameters. Note that (2.6) is indeed a special case of (2.1), with $\boldsymbol{\eta}$ specified to be $\boldsymbol{\tau}=(\tau_{1},\dots,\tau_{K})$, and $\boldsymbol{\Sigma_{\eta}}$ specified to be

[TABLE]

For the model in (2.6), derivations of the full conditionals of $\boldsymbol{\tau}^{2}$, $\sigma^{2}$ and $\boldsymbol{\beta}$ are in the supplementary material. These conditionals enable a Gibbs sampler to explore the corresponding posterior, which updates the three blocks $\boldsymbol{\tau}^{2}$, $\sigma^{2}$ and $\boldsymbol{\beta}$ in that order. The $(\sigma^{2},\boldsymbol{\beta})$-marginal of this Gibbs sampler is a Markov chain that we call the 3BG chain for the Bayesian group lasso model. Its transition density is given by

[TABLE]

2.2 The Bayesian Sparse Group Lasso

In the same setup of the previous section where the explanatory variables form $K$ groups, Simon et al. (2013) introduced the sparse group lasso method that promotes sparsity in regression coefficients both among groups and within each group. For penalty parameters $\lambda_{1},\lambda_{2}>0$, the sparse group lasso estimator is defined as

[TABLE]

where $\lambda_{1}$ and $\lambda_{2}$ are often tuned by cross validation. A Bayesian analogue of (2.8) can be obtained by specifying

[TABLE]

We will refer to the model that uses (2.9) and completed with an Inverse-Gamma prior on $\sigma^{2}$ as the Bayesian sparse group lasso model (Xu and Ghosh, 2015). For $k=1,\dots,K$, let

[TABLE]

and

[TABLE]

The Bayesian sparse group lasso model is a special case of the Bayesian shrinkage model (2.1), as it has the following hierarchical representation:

[TABLE]

where $\boldsymbol{\gamma}^{2}=(\gamma_{1,1}^{2},\dots,\gamma_{1,m_{1}}^{2},\dots,\gamma_{K,1}^{2},\dots,\gamma_{K,m_{K}}^{2})$, $\boldsymbol{\tau}^{2}=(\tau_{1}^{2},\dots,\tau_{K}^{2})$, $\alpha,\xi\geq 0$, and $\lambda_{1},\lambda_{2}>0$ are user specified hyperparameters.

For the model in (2.10), the full conditionals of $(\boldsymbol{\tau}^{2},\boldsymbol{\gamma}^{2})$, $\boldsymbol{\beta}$ and $\sigma^{2}$ are available in the supplementary material. Hence a Gibbs sampler can be run that updates the three blocks $(\boldsymbol{\tau}^{2},\boldsymbol{\gamma}^{2})$, $\sigma^{2}$, and $\boldsymbol{\beta}$, in that order. The $(\sigma^{2},\boldsymbol{\beta})$-marginal forms a Markov chain that we call the 3BG chain for the Bayesian sparse group lasso model. Its transition density is given by

[TABLE]

2.3 The Bayesian Fused Lasso

For regression problems in which explanatory variables can be ordered in some meaningful way, Tibshirani et al. (2005) developed the fused lasso method, that penalizes the $\ell_{1}$ norm of both the coefficients and their successive differences. For fixed tuning parameters $\lambda_{1},\lambda_{2}>0$, the fused lasso estimator is defined as

[TABLE]

We now present the Bayesian counterpart of the fused lasso method. Let the conditional prior of $\boldsymbol{\beta}|\sigma^{2}$ be

[TABLE]

Let $\boldsymbol{\Sigma}_{\tau,\omega}$ be a matrix such that

[TABLE]

A full Bayesian fused lasso model based on a hierarchical representation of (2.13) and a prior on $\sigma^{2}$ was presented in Vats (2017), as shown below:

[TABLE]

where $\boldsymbol{\tau}=(\tau_{1},\dots,\tau_{p})$ and $\boldsymbol{\omega}=(\omega_{1},\dots,\omega_{p-1})$ play the role of $\eta$ of (2.1), and $\alpha,\xi\geq 0$, and $\lambda_{1},\lambda_{2}>0$ are hyperparameters specified by the user. The full conditionals of $(\boldsymbol{\tau}^{2},\boldsymbol{\omega}^{2})$, $\boldsymbol{\beta}$ and $\sigma^{2}$ can be derived from (2.14) and are included in the supplementary material. These conditionals enable a Gibbs sampler that updates the three blocks $(\boldsymbol{\tau}^{2},\boldsymbol{\omega})$, $\sigma^{2}$ and $\boldsymbol{\beta}$ in that order. The $(\sigma^{2},\boldsymbol{\beta})$-marginal forms a Markov chain with transition density

[TABLE]

We refer to the above as the 3BG chain for the Bayesian fused lasso model.

3 Fast MCMC for Bayesian Shrinkage Models

The 3BG algorithms in section 2 are the most commonly used computing solution for the aforementioned Bayesian shrinkage models. Despite straightforward implementations, 3BG converges slowly in high-dimensional settings. For a special case of the Bayesian shrinkage model in (2.1) with $\alpha=\xi=0$, which assigns a “non-informative” improper prior on $\sigma^{2}$, Rajaratnam et al. (2018) noted that the conditional distribution $\pi(\sigma^{2}|\boldsymbol{\eta},\boldsymbol{Y})$ is also tractable, and inexpensive to sample from. Hence, $(\boldsymbol{\beta},\sigma^{2})$ can be updated in one block by drawing from $\pi(\sigma^{2}|\boldsymbol{\eta},\boldsymbol{Y})$ and then $\pi(\boldsymbol{\beta}|\sigma^{2},\boldsymbol{\eta},\boldsymbol{Y})$. We now slightly generalize Lemma 1 of Rajaratnam et al. (2018) by developing a tractable form of the conditional posterior distribution $\pi(\sigma^{2}|\boldsymbol{\eta},\boldsymbol{Y})$ for the general Bayesian shrinkage model in (2.1), where the prior can be either improper or proper. The proof can be found in the supplementary materials.

Lemma 1

For the Bayesian shrinkage model in (2.1), suppose that the prior on $\sigma^{2}$ is $\text{Inverse-Gamma}(\alpha,\xi)$ for $\alpha,\xi\geq 0$. Then, $\sigma^{2}|\boldsymbol{\eta},\boldsymbol{Y}$ follows an Inverse-Gamma distribution with shape parameter $\frac{n}{2}+\alpha$ and scale parameter $\boldsymbol{Y}^{T}(\boldsymbol{I}-\boldsymbol{X}(\boldsymbol{X}^{T}\boldsymbol{X}+\boldsymbol{\Sigma}_{\boldsymbol{\eta}}^{-1})^{-1}\boldsymbol{X}^{T})\boldsymbol{Y}/2+\xi$.

Lemma 1 leads to a Gibbs sampler that updates the two blocks $\boldsymbol{\eta}$ and $(\sigma^{2},\boldsymbol{\beta})$ in turn:

[TABLE]

Comparing (3.1) to (2.2), note that the only technical difference between 2BG and its 3BG counterpart lies in the step of drawing $\sigma^{2}$, which has similar cost for the two algorithms. Note that despite a seemingly harder-to-evaluate conditional distribution of $\sigma^{2}$ in 2BG that requires the inversion of a $p\times p$ matrix $\boldsymbol{A}_{\boldsymbol{\eta}}=\boldsymbol{X}^{T}\boldsymbol{X}+\Sigma_{\boldsymbol{\eta}}^{-1}$, this quantity is actually needed in other steps of both algorithms, hence not an extra burden for 2BG.

Convergence properties of the Markov chain that underlies this Gibbs sampler is characterized by its $\left(\boldsymbol{\beta},\sigma^{2}\right)$-marginal, $\boldsymbol{\Phi}:=\left\{({\boldsymbol{\beta}}_{m},{\sigma}_{m}^{2})\right\}_{m=0}^{\infty}$. We call this Markov chain on the state space $\mathbb{R}^{p}\times\mathbb{R}_{+}$ the 2BG chain. Its transition density is given by

[TABLE]

In the subsections below, we provide details of the new 2BG samplers for the three aforementioned Bayesian shrinkage models, each followed by a theorem on their theoretical properties and that of their 3BG counterparts. Proofs of the three Theorems are in the supplementary materials, and we give a brief overview of them here. For Bayesian group lasso models and Bayesian sparse group lasso models with hyperparameters $\alpha\geq 0$ and $\xi>0$, we show that 2BG chains are trace-class, which imply they are Hilbert-Schmidt and have GE. For these two Bayesian models with hyperparameters $\alpha\geq 0$ and $\xi=0$, and for the Bayesian fused lasso models with $\alpha,\xi\geq 0$ we do not yet know if their 2BG chains are trace-class, due to technical difficulties; instead, we establish their GE directly, using the classical drift and minorization condition (Rosenthal, 1995).

Remark 1

It’s straightforward from their definitions that a self-adjoint trace-class operator is also Hilbert-Schmidt. Further, a Hilbert-Schmidt operator is always compact. For a Harris ergodic Markov chain, the compactness implies that the spectral radius of the Markov operator is strictly less than $1$, hence geometrically ergodic. See, e.g., Chan and Geyer (1994) and Roberts and Rosenthal (1997). Note that the Markov transition densities (3.4), (3.6) and (3.8) are indeed strictly positive, it follows that the corresponding 2BG chains are Harris ergodic (Tan and Hobert, 2009, Lemma 1). Therefore, the Theorems proved in this paper show that all 2BG chains for the three types of Bayesian shrinkage models have GE, while many of them are also shown to have the stronger trace-class property. In comparison, 3BG chains for the same Bayesian shrinkage models were shown to have GE by Vats (2017), but our Theorems show that none of them are Hilbert-Schmidt.

3.1 The Two-block Gibbs Sampler for the Bayesian Group Lasso

For the Bayesian group lasso model, the 2BG updates described in (3.1) become

[TABLE]

The one-step transition density of the $\left(\boldsymbol{\beta},\sigma^{2}\right)$-chain is given by

[TABLE]

Theorem 1

For the Bayesian group lasso model with hyperparameters $\alpha\geq 0$ and $\xi>0$, the Markov operator of the 2BG chain with transition density $k_{gl}$ is trace-class, and hence is Hilbert-Schmidt; while that of the 3BG chain with transition density $\tilde{k}_{gl}$ is not Hilbert-Schmidt for any $\alpha,\xi\geq 0$. In the case of $\alpha\geq 0$ and $\xi=0$, the 2BG chain is geometrically ergodic.

3.2 The Two-block Gibbs Sampler for the Bayesian Sparse Group Lasso

For the Bayesian sparse group lasso model, the 2BG updates the two blocks $(\boldsymbol{\tau},\boldsymbol{\gamma})$ and $(\sigma^{2},\boldsymbol{\beta})$ in each iteration as follows:

[TABLE]

The one-step transition density of the $\left(\boldsymbol{\beta},\sigma^{2}\right)$-chain is given by

[TABLE]

Theorem 2

For the Bayesian sparse group lasso model with hyperparameters $\alpha\geq 0$ and $\xi>0$, the Markov operator of the 2BG chain with transition density $k_{sgl}$ is trace-class, and hence is Hilbert-Schmidt; while that of the 3BG chain with transition density $\tilde{k}_{sgl}$ is not Hilbert-Schmidt for any $\alpha,\xi\geq 0$. In the case of $\alpha\geq 0$ and $\xi=0$, the 2BG chain is geometrically ergodic.

3.3 The Two-block Gibbs Sampler for the Bayesian Fused Lasso

For the Bayesian fused lasso model, the 2BG updates the two blocks $(\boldsymbol{\tau},\boldsymbol{\omega})$ and $(\sigma^{2},\boldsymbol{\beta})$ in each iteration as follows:

[TABLE]

The one-step transition density of the $\left(\boldsymbol{\beta},\sigma^{2}\right)$-chain is given by

[TABLE]

Theorem 3

For any $\alpha,\xi\geq 0$, the 2BG chain for the Bayesian fused lasso model is geometrically ergodic. The Markov operator of the 3BG chain with transition density $\tilde{k}_{fl}$ is not Hilbert-Schmidt for any $\alpha,\xi\geq 0$.

4 Simulation Studies

We use simulation studies to demonstrate the advantage of the newly developed 2BG’s over existing MCMC computing solutions for the three Bayesian shrinkage models of interest. Existing methods include 3BG and HMC. Briefly, HMC is a state of the art generic MCMC technique that has shown remarkable empirical success in exploring many difficult Bayesian posteriors. See Neal (2011) for an introduction. All HMC samplers in our study are run using the Stan software, which simply requires specification of the Bayesian models in Stan language, then the software will automatically derive the log posterior density, and tune and run HMC algorithms accordingly. Note that, similar to that of Gibbs samplers, different representations or parameterizations of the same model will lead to different HMC algorithms. We experimented two ways to implement HMC for the Bayesian group lasso model. The first one uses the conditional prior $\pi(\boldsymbol{\beta}|\sigma^{2})$ in (2.5) directly, while the second one uses its hierarchical representation in (2.6). For the sparse group lasso and the fused lasso models, we implemented them in Stan where the conditional priors of $\boldsymbol{\beta}|\sigma^{2}$ are directly specified as (2.9) and (2.13), respectively. Their hierarchical versions are not implemented because they are highly inconvenient to put down in the language of Stan, and are expected to be much more costly per iteration. The 2BG and 3BG samplers are coded in $\text{C}_{++}$, then called from within R using the package Rcpp. The code can be found in an online supplement to this paper.

For the Markov chain underlying each algorithm, let $\rho_{k}$ denote the lag-$k$ autocorrelation of the $\sigma^{2}-$marginal of the chain (under stationarity). Rajaratnam and Sparks (2015, sec 3) provided comprehensive reasons why $\rho_{1}$ is a good proxy for the mixing rate of the joint Markov chain, more so than that of the autocorrelations of the regression coefficients. For the empirical evaluation of an MCMC algorithm, the closer the estimates of $\rho_{1}$ is to [math], the better the chain mixes. Another criteria we used to compare the chains is the per unit time effective sample size of the $\sigma^{2}$ component. Specifically, $N$ iterations of the Markov chain produce an effective sample size (ESS) of

[TABLE]

We estimate $N_{\text{eff}}$ of the HMC, the 3BG and the 2BG samples using the effectiveSize function from the R package coda (Plummer et al., 2006). Then the estimates of $N_{\text{eff}}$ are divided by $T$, their respective running time in seconds, for the eventual comparisons reported in this section.

Remark 2

Besides $\rho_{1}$ and $N_{eff}/T$, we tried to compare the different MCMC algorithms using the multivariate ESS (Wilks, 1932; Vats et al., 2019), carried out using the function multiESS from the R package mcmcse (Flegal et al., 2017). Under our simulation setups, the estimates of the multivariate ESS for the full vector of interest $(\sigma^{2},\boldsymbol{\beta})$ do not appear stable enough, potentially due to the very high dimension of this vector. Hence, we opt to report the estimated multivariate ESS of the 3-dimensional vector $(\sigma^{2},\beta_{1},\beta_{2})$ for our simulation study, which turns out to largely agree with the comparison results based on the (univariate) ESS for the $\sigma^{2}$ component, $N_{\text{eff}}$, defined above. The results for multivariate ESS can be found in the supplementary material.

4.1 Comparing MCMC Algorithms for the Three Bayesian Shrinkage Models

Data in this study are simulated using the following model:

[TABLE]

where $\boldsymbol{\epsilon}$ is a $n\times 1$ vector of independent standard normal random variables, and $\boldsymbol{\beta}_{*}$ is a $p\times 1$ vector consisting of the true regression coefficients.

• •

Scenario 1. Grouped explanatory variables. Denote $\mathcal{K}=\{5,6,7,8,9,10,20,30,40,50\}$. For $n=50$, let $K\in\mathcal{K}$. For $n=100$, let $K\in 2\mathcal{K}$. We generate $K$ variables independently from the multivariate normal distribution with mean vector [math] and the identity variance-covariance matrix. We build order-5 polynomials based on each variable to get a $n\times p$ design matrix $\boldsymbol{X}$, where $p=5K$ and every 5 consecutive columns belong to a group. For each $p$ we experiment with, $\boldsymbol{\beta}_{*}$ is such that only its first $p/5$ elements are nonzero, and are drawn independently from the $t_{2}$ distribution.

• •

Scenario 2. Similar coefficients for adjacent explanatory variables. At $n=100$ and $200$, consider $\frac{p}{n}=\{0.5,0.6,0.7,0.8,0.9,1,2,3,4,5\}$. To generate the $n\times p$ design matrix $\boldsymbol{X}$, the rows are independently drawn, each from the $p$-dim multivariate normal distribution, where the mean vector is [math], and the variance-covariance matrix has $1$ on the diagonal, and $0.2$ everywhere else. Then, the columns of $\boldsymbol{X}$ are standardized to have mean zero and squared Euclidean norm $n$. The first and the third $p/10$ elements of $\boldsymbol{\beta}_{*}$ are drawn independently from ${\cal{N}}(1,0.1^{2})$, and all other elements are set to zero.

In analyzing the simulated data, both the Bayesian group lasso and the Bayesian sparse group lasso are applied to scenario 1, and the Bayesian fused lasso to scenario 2. Take the comparison of MCMC algorithms for the Bayesian group lasso model for example. At each $(n,p)$ combination, $100$ datasets are generated from scenario 1. For every data set, each kind of MCMC algorithm is run for $10,000$ iterations, with the first $1,000$ discarded as burn in. Each point in Figures 1, 3 and 5 represents the average lag-one autocorrelation estimates from $100$ Markov chains of the same kind, over different datasets. Similarly, each point in Figures 2, 4 and 6 represents the decadic logarithm of the average of $N_{eff}/T$ over the $100$ Markov chains. To show that the comparison results in the Figures are meaningful and not due to randomness, we include summaries of the standard error of these averages in the main text, and include boxplots for the 100 estimates of $\rho_{1}$ and that of $N_{eff}/T$ in the supplementary material. The series of experiments are run on the High Performance Argon Cluster at University of Iowa, on a standard machine with 64 GB RAM, 16 cores and 2.6GHz processor.

Remark 3

Note that we choose the relatively short chain length of $10,000$ for each chain when comparing all the methods because the HMC samplers are very costly. We are actually able to easily run each of the 2BG and 3BG samplers for $100,000$ iterations in all scenarios. Estimates of $\rho_{1}$ and $N_{\text{eff}}/T$ based on the longer chains are only included in the supplementary material, and they are numerically similar to that of the shorter chains reported here.

Remark 4

The main goal of the our simulation study is to learn the behavior of the different algorithms, not to see which models or hyperparameter values fit each scenario the best. Hence for each dataset, we have only considered the most reasonable model with their default hyperparameter values. We always set $\alpha=\xi=0$, which lead to improper priors on $\sigma^{2}$. Also, we set $\lambda=1$ for the Bayesian group lasso, and $\lambda_{1}=\lambda_{2}=1$ for the other two Bayesian models. The resulting posteriors are proper (Kyung et al., 2010; Xu and Ghosh, 2015). Our future work includes investigating the impact of hyperparameter values through Bayesian cross-validation and Bayesian sensitivity analysis, or doing fully Bayesian analysis that employs hyperpriors. All such endeavors will heavily depend on the fast algorithms developed in this article.

4.1.1 Results for the Bayesian Group Lasso

Figure 1 displays the average lag-one autocorrelation of the $\sigma^{2}$-chain. We can see that the proposed 2BG has dramatically smaller autocorrelations than the classical 3BG across all setups of $(n,p)$, indicating better mixing of 2BG. Both HMC samplers have comparable mixing rates to that of 2BG when $p$ is no bigger than $n$. But with larger $p$, the lag-one autocorrelation of the HMCs deteriorates, while that of 2BG stays below $0.4$. The standard error of the average $\rho_{1}$ at all setups are calculated, which never exceed $0.0065$ and are typically between $0.00073$ and $0.0064$. Further, Figure 2 takes computing time into consideration by examining the average effective sample size per second, $N_{eff}/{T}$, in the log scale of each sampler. The standard error of the average $N_{eff}/{T}$ (in base-$10$ log scale) at all setups are less than $0.021$ with typical values between $0.0047$ and $0.021$. These observations show that 2BG is computationally more efficient than all others. For example, at $n=p=50$, 2BG produces about twice as many effective samples per second as that of 3BG, $430$ times that of “Stan” and $766$ times that of “Stan_hier’. At $n=50,p=250$, the relative efficiency of 2BG over the other three methods increases to $7$, $2172$ and $71640$, respectively. Note that “Stan_hier” is computationally much less efficient than “Stan”, as a result of similar mixing rate of the underlying Markov chains, but that the hierarchical version involves a higher-dimensional posterior that greatly increases the computational cost per iteration.

4.1.2 Results for the Bayesian Sparse Group Lasso

Figure 3 displays the lag-one autocorrelation of the $\sigma^{2}$-chain of the 2BG, the 3BG and the HMC sampler for the Bayesian sparse group lasso model, which shows a similar pattern as that of the group lasso model. Again, 2BG always mixes much faster than 3BG. HMC mixes faster than 2BG when $p\leq n$, but 2BG shows an increasingly large advantage over HMC as $p$ grows larger than $n$. After taking computing time into account, in Figure 4, we see that 2BG produces the most effective samples per second in all setups. For example, at $n=p=50$, 2BG produces around $3$ and $43$ times as many effective samples per second as that of 3BG and HMC, respectively. And at $n=50,p=250$, the relative efficiency of 2BG to the other two samplers grow to about $9$ and $273$, respectively. As for the standard error associated with the averages in Figures 3 and 4, that of the $\rho_{1}$ estimates never exceed $0.0087$ with typical values between $0.00081$ and $0.0077$, and that of the $N_{eff}/{T}$ estimates (in base-$10$ log scale) never exceed $0.015$, with typical values between $0.0047$ and $0.014$. Boxplots of the estimates from all the replications are also available in the supplementary material. These suggest that the patterns observed in Figures 3 and  4 are meaningful.

4.1.3 Results for the Bayesian Fused Lasso

Figure 5 displays the autocorrelation of the $\sigma^{2}$-chain of the three samplers for the Bayesian fused lasso model. 2BG and HMC always mix much faster than 3BG, while 2BG has an advantage over HMC for large values when $p$ exceeds $n$. Taking computing time into account, Figure 6 shows that 2BG is the most efficient sampler in all situations. For example, at $n=p=100$, the 2BG produces $3$ times as many effective samples as that of 3BG and $60$ times that of “Stan”. When $p$ increases to $500$, 2BG is $8$ times as efficient as 3BG and twice that of HMC, though all three generate less than $1$ effective sample per second. The standard error of the averages are small, that for the $\rho_{1}$ estimates never exceed $0.0046$ and are typically between $0.00058$ and $0.0044$, and that for the $N_{eff}/{T}$ estimates (in the log scale) never exceed $0.018$ and are typically between $0.0032$ and $0.018$. Boxplots of the estimates from all the replications are also available in the supplementary material. These suggest that the patterns observed in Figures 5 and  6 are meaningful.

4.2 Analyzing Extra Wide and Extra Tall Datasets with Bayesian Group Lasso Models

We perform additional studies following directions suggested by a reviewer, including analyzing datasets that are extra wide, extra tall and have non-normal noises. The first two are presented below, and the third in the supplement.

First, we consider extra wide datasets and restrict the study to Bayesian group lasso models due to space limit. We form datasets with $n=50$ and $\frac{p}{n}=(10,20,30,\cdots,100)$, respectively. For each $p$, the design matrix $\boldsymbol{X}$ is generated using the same scheme as that of scenario 1 in section 4.1 and $\boldsymbol{\beta}_{*}$ is such that only its first $5$ elements are nonzero, and are drawn independently from the $t_{2}$ distribution. Note that we fix the number of true covariates to achieve more and more sparse signals as $p$ increases, which forms more extreme setups than the ones in section 4.1. At each $(n,p)$ combination, $100$ datasets are generated, and for every data set, each kind of MCMC algorithm is run for $10,000$ iterations, with the first $1,000$ discarded as burn in. Figure 7 displays the average lag-one autocorrelations and the average estimates of $N_{eff}/T$ for the $\sigma^{2}$-chain, and corresponding boxplots can be found in the supplement. We see that, as $p$ grows much larger than $n$, the mixing rates of 3BG and HMC deteriorate, while that of the 2BG actually improves.

Next, we consider Bayesian group lasso models applied to tall datasets. Specifically, for $p=25$ and $n\in\{50,100,150,\cdots,500\}$, $100$ datasets are generated for each $(n,p)$ combination using the same generation scheme as above. MCMC samples of the three algorithms are generated in the same manner as before. The left panel of Figure 8 shows that all Markov chains mix quite well for tall datasets, where HMC mixes the best and often enjoys negative autocorrelations, and 2BG generates almost independent samples. The right panel of Figure 8 shows that, once taking computing time into consideration, 2BG is still more efficient than HMC. It’s also interesting to see that the ESS per second of 3BG is comparable to that of 2BG. Indeed, despite having slightly worse mixing rates at large $n$, 3BG chains cost a little less per iteration than 2BG in these setups.

4.3 Summary of Simulation Results

For all three Bayesian shrinkage models, 2BG is almost always the most efficient among existing algorithms. In the most demanding, high-dimensional regression problems where the number of predictors $p$ is larger than the data size $n$, 2BG produces about $10$ – $10^{4}$ times as many effective samples per second as that of its competitors. Contributing to the advantage of 2BG are (1) the fast mixing rate of 2BG, which is similar to that of the state of the art method HMC for small $p$, and even faster than that of the HMC for larger $p$, and (2) the low computing cost per iteration of 2BG, which is similar to that of 3BG, and much lower than that of HMC.

Besides the theoretical results on 2BG and 3BG chains presented in section 3, our intuition for why 2BG mixes much better than HMC and 3BG in our studies are as follows. For Bayesian shrinkage models, Rajaratnam and Sparks (2015) established that $\boldsymbol{\beta}$ and $\sigma^{2}$ has high posteriori dependence when $p>>n$. For HMC in general, Livingstone and Zanella (2019) argued that HMC is highly sensitive to tuning when the posterior features strong dependence among components, partly because the spectral gap of its Markov operator decays to zero very fast when the mismatch between the proposal and target scales increases in some directions. Hence, for Bayesian shrinkage models with large $p$, HMC can be extremely hard to tune to overcome the high posterior dependence given the current techniques and softwares, as we see in Figure 7 for the example of extra wide datasets. Similarly, 3BG involves drawing from the conditionals $\boldsymbol{\beta}|\sigma^{2},\tau$ and $\sigma^{2}|\boldsymbol{\beta},\tau$, which produces high auto-correlations between successive $\boldsymbol{\beta}$ and $\sigma^{2}$ samples in the presence of high posteriori dependency between these two components. In contrast, the newly proposed 2BG avoid the high dependence problem by updating $\boldsymbol{\beta}$ and $\sigma^{2}$ within the same block, and hence explore the posterior distribution of high-dimensional Bayeisan shrinkage models more efficiently. See Liu et al. (1994) for some general theory on the benefit of blocking in Markov chains.

5 Real Data Analysis

In this section, we observe superior performance of the proposed 2BG over 3BG and HMC in real data applications. Three data sets are analyzed, each using a different Bayesian shrinkage model that is appropriate for the structure of the covariates. For each model, all three computing methods are used and evaluated in terms of mixing rate and computational efficiency, as summarized in Table 1.

Each algorithm is run for $20,000$ iterations with the first one-tenth discarded as “burn-in”. The experiments are run on a 2.2 GHz machine with 16GB RAM and MacOS Sierra system. More details of the data and the comparisons in computation are given in the subsections below.

5.1 Gene Expression Data

We first evaluate the performance of 2BG, 3BG and HMC for a Bayesian group lasso model applied to a gene expression data of Scheetz et al. (2006), which is publicly available as the data set bardet in the R package gglasso (Yang and Zou, 2015). The response variable is the expression level of gene TRIM32, which causes Bardet-Biedl syndrome. The expression level of $20$ other genes are available and expanded using $5$ basis B-splines. The resulting design matrix $\boldsymbol{X}$ contains $p=100$ columns where each consecutive $5$ correspond to the same gene. Each column is then standardized to have mean zero and squared Euclidean norm $n$. The Bayesian group lasso model was implemented with regularization parameter $\lambda=0.06$, based on a cross-validation step of its frequentist counterpart. Details can be found in the supplementary materials.

The 2BG chain mixes very fast with an average lag-one autocorrelation of $0.057$ for the $\sigma^{2}$-chain, compared to that of $0.19$ of HMC, and $0.40$ of 3BG. Next, consider the number of effective samples of the $\sigma^{2}$-component generated per second, 2BG comes first at around $579$, followed by 3BG at $245$, and HMC at only $0.08$.

5.2 Economic Data

We use the Bayesian sparse group lasso model to analyze an economic data from Rose and Spiegel (2010), which is publicly available as the data set Crisis2008BalancedData in the R package BSGS (Lee and Chen, 2015). The response variable is the growth rate of a crisis measure from 2008 to 2009. The design matrix $\boldsymbol{X}$ contains $51$ explanatory variables for the 2008 financial crisis for $72$ countries. These $51$ variables can be classified into $9$ groups. All the columns of $\boldsymbol{X}$ are standardized except for those corresponding to the dummy variables. The regularization parameters were taken to be $\lambda_{1}=0.104$ and $\lambda_{2}=0.082$.

Again, using the lag-one autocorrelation of the chain $\sigma^{2}$, 2BG comes in first at $0.014$, followed by HMC at $0.029$ and 3BG at a substantially larger correlation at $0.42$. In terms of the number of effective samples per second, $N_{eff}/T$, 2BG is the best at $1,213$, followed by $494$ of 3BG, and $9$ of HMC.

5.3 Comparative Genomic Hybridization (CGH) Data

We apply the Bayesian fused lasso model to a subset of the CGH data of Tibshirani and Wang (2007), which is publicly available as the data set CGH in the R package cghFLasso. Tibshirani and Wang (2007) analyzed this data with fused lasso model. The response contains $200$ samples of CGH array. The design matrix is the identity matrix. The regularization parameters were taken to be $\lambda_{1}=0.129$ and $\lambda_{2}=0.962$.

Concerning the lag-one autocorrelation of the $\sigma^{2}$-chain, 2BG and the HMC are comparable, at $0.29$ and $0.23$, respectively. While 3BG falls behind at $0.64$. After taking time into account, 2BG is clearly the winner by producing $48.86$ effective samples per second, compared to $19.14$ from 3BG, and $1.38$ from HMC.

6 Discussion

As we have seen above, various heuristic arguments and theoretical results suggest that 2BG is expected to be more efficient than 3BG, and the empirical evidence in sections 4 and 5 show that the improvements are substantial. There are currently no theoretical result that directly compares 2BG, which is a simple Gibbs type algorithm, to HMC, which is a gradient based algorithm. It may be surprising to many that iteration-wise, 2BG mixes better than HMC in the most challenging high-dimensional cases, while performing similarly in other cases. And the advantage of 2BG is more striking when computing cost is taken into account. For example, to obtain $10,000$ effective samples for the Bayesian group lasso model with a data of size $p=5n=500$, it roughly takes 2BG 1 minute, 3BG 10 minutes, and HMC (via stan) 13.5 hours of computing time. That is, when using these Bayesian shrinkage models in some moderate to high-dimensional setups, 2BG could be the only viable MCMC solution.

Another interesting observation common to all experiments in section 4 is that, while the autocorrelations of 3BG and HMC chains increase quickly towards $1$ as $p/n$ increases, that of 2BG chains stay bounded and actually decrease. One of our on-going work is performing complexity analysis of the 2BG algorithms, that is, studying their asymptotic geometric convergence rate as the dimension of the data increases. Although such theoretical results are extremely hard to obtain, and the proofs need to be tailored to each model specification, there have been several successful examples along with new analysis methods. See, e.g., Qin and Hobert (2018, 2019); Yang and Rosenthal (2017).

An important practical issue of applying any Bayesian models for data analysis is assigning priors or specifying values for the hyperparameters. This topic is beyond the scope of the current paper but deserves further research. Bayesian cross-validation and Bayesian sensitivity analysis provide ways to calibrate and evaluate different choices of priors and hyperparameters. The fast algorithms developed here serve as critical building blocks to perform these analysis that is computationally more demanding.

7 Supplementary Materials

Additional Results

A document containing background knowledge and proofs of theoretical results, details on the application of the models, and additional simulation results. (pdf)

R code

Code for the examples in this paper. (zip)