The Horseshoe at Scale — Selection when p ≫ n

Overview

Part 8 of the variable-selection arc: selection when $p\gg n$. Every dataset so far had a handful of predictors; real variable selection often looks the opposite — thousands of candidate predictors, tens of samples. Genomics is the canonical case, and exactly where the earlier machinery breaks down: enumeration would need $2^{4088}$ models, and discrete spike-and-slab mixes painfully with 4,088 indicators. The tool built for this regime is the continuous global-local shrinkage prior — above all the horseshoe (Carvalho–Polson–Scott 2010). The benchmark is the riboflavin dataset (Bühlmann–Kalisch–Meier 2014): $n=71$ Bacillus subtilis cultures, $p=4088$ gene-expression measurements, and a continuous response — the log vitamin-B2 production rate. The question: which handful of genes drives production? Fit from scratch, cross-checked against PyMC (same prior, NUTS) and R (glmnet, the lasso — the frequentist workhorse for $p\gg n$).

Global-local shrinkage, and the sampler that makes it scale

The horseshoe gives each coefficient its own local scale on top of a global one: $\beta_j\sim\mathcal N(0,\lambda_j^2\tau^2\sigma^2)$, with $\lambda_j\sim C^+(0,1)$ and $\tau\sim C^+(0,1)$. The global $\tau$ pulls everything toward zero (essential when $p\gg n$); the heavy-tailed local $\lambda_j$ lets a few coefficients escape. It is essentially parameter-free — the half-Cauchys need no tuning — and its action is summarised by the shrinkage weight $\kappa_j=1/(1+\lambda_j^2\tau^2)\in[0,1]$ ($\kappa_j\approx1$: shrunk to zero; $\kappa_j\approx0$: kept). Two ingredients make it tractable. (1) Makalic–Schmidt (2016) auxiliary variables turn each half-Cauchy into a pair of conjugate inverse-gamma draws, so the whole sampler is Gibbs. (2) The Bhattacharya–Chakraborty–Mallick (2016) sampler for the coefficient block: the naive draw needs the $p\times p$ precision — $O(p^3)=O(4088^3)$, hopeless — but their exact algorithm never forms it, solving an $n\times n$ system instead at $O(n^2p)$ cost. With $n=71$ that turns an impossible step into a millisecond one — the entire run is ~9 seconds. There are no free hyperparameters (all half-Cauchy scales are 1, the standard horseshoe). Engine: a self-contained hshd.py.

Results

A needle in 4,088 haystacks. The horseshoe drives essentially every coefficient to zero — the median shrinkage weight is 0.998 — and lets just a couple escape: YOAB and YXLE stand out sharply ($\kappa\approx 0.25, 0.29$), with a weaker second tier (ARGB, YISU). This is what selection must look like when $p\gg n$: aggressive global shrinkage, a few local exceptions — and the signals are biologically sensible (the YXL-region genes and YOAB recur across analyses of this benchmark), recovered without any pre-screening. Both Bayesian engines agree on what is real, not on what is noise: with priors matched, the from-scratch Gibbs and PyMC's NUTS pick the same dominant gene (YOAB; YXLE the runner-up in the from-scratch Gibbs) and R's lasso puts YOAB first too — but across all 4,088 coefficients the two samplers correlate only $r\approx 0.80$. That is the posterior being honest: the ~4,086 near-zero coefficients are $\approx 0$ with large uncertainty, so their exact values are not meant to be reproduced — the reliable outputs of a $p\gg n$ fit are the shrinkage weights and the top-gene ranking, not individual near-zero estimates. The sampler matters enormously at scale: the specialised Gibbs finishes in ~10 seconds; general-purpose NUTS takes ~2.7 minutes (~16× longer) and still mixes poorly against the horseshoe funnel in 4,088 dimensions. When $p\gg n$, the algorithm, not just the prior, is the enabling idea. The horseshoe answers the same question as SSVS — but by shrinking rather than switching: no $2^p$ model space, no binary indicators, just a heavy-tailed prior that lets a few coefficients through, which is why global-local priors dominate the $p\gg n$ regime.

Notebooks

Downloads

References

• Carvalho, C. M., Polson, N. G. & Scott, J. G. (2010). The horseshoe estimator for sparse signals. Biometrika 97(2), 465–480. — the horseshoe prior
• Bhattacharya, A., Chakraborty, A. & Mallick, B. K. (2016). Fast sampling with Gaussian scale-mixture priors in high-dimensional regression. Biometrika 103(4), 985–991. — the O(n²p) coefficient sampler that makes p ≫ n tractable
• Makalic, E. & Schmidt, D. F. (2016). A simple sampler for the horseshoe estimator. IEEE Signal Processing Letters 23(1), 179–182. — the inverse-gamma augmentation giving a fully-conjugate Gibbs sampler
• Polson, N. G. & Scott, J. G. (2011). Shrink globally, act locally: sparse Bayesian regularization and prediction. In Bayesian Statistics 9, 501–538. Oxford University Press. — the global-local shrinkage framework
• Bühlmann, P., Kalisch, M. & Meier, L. (2014). High-dimensional statistics with a view toward applications in biology. Annual Review of Statistics and Its Application 1, 255–278. — the riboflavin benchmark dataset (hdi::riboflavin)
• Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society B 58(1), 267–288. — the lasso, the R glmnet cross-check
• Piironen, J. & Vehtari, A. (2017). Sparsity information and regularization in the horseshoe and other shrinkage priors. Electronic Journal of Statistics 11(2), 5018–5051. — shrinkage weights and effective number of nonzero coefficients