Fast sampling with Gaussian scale-mixture priors in high-dimensional regression

Introduction

Continuous shrinkage priors have recently received significant attention as a mechanism to induce approximate sparsity in high-dimensional parameters. Such prior distributions can mostly be expressed as global-local scale mixtures of Gaussians (Polson & Scott, 2010; Bhattacharya et al., 2015). These global-local priors (Polson & Scott, 2010) aim to shrink noise coefficients while retaining any signal, thereby providing an approximation to the operating characteristics of discrete mixture priors (George & McCulloch, 1997; Scott & Berger, 2010), which allow a subset of the parameters to be exactly zero.

A major attraction of global-local priors has been computational efficiency and simplicity. Posterior inference poses a stiff challenge for discrete mixture priors in moderate to high-dimensional settings, but the scale-mixture representation of global-local priors allows parameters to be updated in blocks via a fairly automatic Gibbs sampler in a wide variety of problems. These include regression (Caron & Doucet, 2008; Armagan et al., 2013), variable selection (Hahn & Carvalho, 2015), wavelet denoising (Polson & Scott, 2010), factor models and covariance estimation (Bhattacharya & Dunson, 2011; Pati et al., 2014), and time series (Durante et al., 2014). Rapid mixing and convergence of the resulting Gibbs sampler for specific classes of priors has been recently established in the high-dimensional regression context by Khare & Hobert (2013) and Pal & Khare (2014). Moreover, recent results suggest that a subclass of global-local priors can achieve the same minimax rates of posterior concentration as the discrete mixture priors in high-dimensional estimation problems (Bhattacharya et al., 2015; van der Pas et al., 2014; Pati et al., 2014).

In this article, we focus on computational aspects of global-local priors in the high-dimensional linear regression setting

where $X\in\Re^{n\times p}$ is a $n\times p$ matrix of covariates with the number of variables $p$ potentially much larger than the sample size $n$. A global-local prior on $\beta$ assumes that

where $f,g$ and $h$ are densities supported on $(0,\infty)$. The $\lambda_{j}$s are usually referred to as local scale parameters while $\tau$ is a global scale parameter. Different choices of $f$ and $g$ lead to different classes of priors. For instance, a half-Cauchy distribution for $f$ and $g$ leads to the horseshoe prior of Carvalho et al. (2010). In the $p\gg n$ setting where most entries of $\beta$ are assumed to be zero or close to zero, the choices of $f$ and $g$ play a key role in controlling the effective sparsity and concentration of the prior and posterior (Polson & Scott, 2010; Pati et al., 2014).

The algorithm

Frequentist operating characteristics in high dimensions

The proposed algorithm provides an opportunity to compare the frequentist operating characteristics of shrinkage priors in high-dimensional regression problems. We compare various aspects of the horseshoe prior (Carvalho et al., 2010) to frequentist procedures and obtain highly promising results. We expect similar results for the Dirichlet–Laplace (Bhattacharya et al., 2015), normal-gamma (Griffin & Brown, 2010) and generalized double-Pareto (Armagan et al., 2013) priors, which we hope to report elsewhere.

While there is now a huge literature on penalized point estimation, uncertainty characterization in $p>n$ settings has received attention only recently (Zhang & Zhang, 2014; van de Geer et al., 2014; Javanmard & Montanari, 2014). Although Bayesian procedures provide an automatic characterization of uncertainty, the resulting credible intervals may not possess the correct frequentist coverage in nonparametric/high-dimensional problems (Szabó et al., 2015). This led us to investigate the frequentist coverage of shrinkage priors in $p>n$ settings; it is trivial to obtain element-wise credible intervals for the $\beta_{j}$s from the posterior samples. We compared the horseshoe prior with van de Geer et al. (2014) and Javanmard & Montanari (2014), which can be used to obtain asymptotically optimal element wise confidence intervals for the $\beta_{j}$s. We considered a similar simulation scenario as before. We let $p\in\{500,1000\}$, and considered a Toeplitz structure, $\Sigma_{jj^{\prime}}=0.9^{|j-j^{\prime}|}$, for the covariate design (van de Geer et al., 2014) in addition to the independent and compound symmetry cases stated already. The first two rows of Table 1 report the average coverage percentages and $100\times$lengths of confidence intervals over $100$ simulation replicates, averaged over the $5$ signal variables. The last two rows report the same averaged over the $(p-5)$ noise variables.

Table 1 shows that the horseshoe has a superior performance. An attractive adaptive property of shrinkage priors emerges, where the lengths of the intervals automatically adapt between the signal and noise variables, maintaining the nominal coverage. The frequentist procedures seem to yield approximately equal sized intervals for the signals and noise variables. The default choice of the tuning parameter $\lambda\asymp(\log p/n)^{1/2}$ suggested in van de Geer et al. (2014) seemed to provide substantially poorer coverage for the signal variables at the cost of improved coverage for the noise, and substantial tuning was required to arrive at the coverage probabilities reported. The default approach of Javanmard & Montanari (2014) produced better coverages for the signals compared to van de Geer et al. (2014). The horseshoe and other shrinkage priors on the other hand are free of tuning parameters. The same procedure used for estimation automatically provides valid frequentist uncertainty characterization.

Discussion

Our numerical results warrant additional numerical and theoretical investigations into properties of shrinkage priors in high dimensions. The proposed algorithm can be used for essentially all the shrinkage priors in the literature and should prove useful in an exhaustive comparison of existing priors. Its scope extends well beyond linear regression. For example, extensions to logistic and probit regression are immediate using standard data augmentation tricks (Albert & Chib, 1993; Holmes & Held, 2006). Multivariate regression problems where one has a matrix of regression coefficients can be handled by block updating the vectorized coefficient matrix ; even if $p<n$, the number of regression coefficients may be large if the dimension of the response if moderate. Shrinkage priors have been used as a prior of factor loadings in Bhattacharya & Dunson (2011). While Bhattacharya & Dunson (2011) update the $p>n$ rows of the factor loadings independently, exploiting the assumption of independence in the idiyosyncratic components, their algorithm does not extend to approximate factor models, where the idiyosyncratic errors are dependent. The proposed algorithm can be adapted to such situations by block updating the vectorized loadings. Finally, we envision applications in high-dimensional additive models where each of a large number of functions is expanded in a basis, and the basis coefficients are updated in a block.

Appendix

Here we give a constructive argument which leads to the poroposed algorithm. By the Sherman–Morrison–Woodbury formula (Hager, 1989) and some algebra we have,

References