Dirichlet-Laplace priors for optimal shrinkage
Anirban Bhattacharya, Debdeep Pati, Natesh S. Pillai, David B. Dunson
Introduction
High-dimensional data have become commonplace in broad application areas, and there is an exponentially increasing literature on statistical and computational methods for big data. In such settings, it is well known that classical methods such as maximum likelihood estimation break down, motivating a rich variety of alternatives based on penalization and thresholding. There is a rich theoretical literature justifying the optimality properties of such penalization approaches , with fast algorithms and compelling applied results leading to routine use of regularization in particular.
The overwhelming emphasis in this literature has been on rapidly producing a point estimate with good empirical and theoretical properties. However, in many applications, it is crucial to obtain a realistic characterization of uncertainty in the estimates of parameters and functions of the parameters, and in predictions of future outcomes. Usual frequentist approaches to characterize uncertainty, such as constructing asymptotic confidence regions or using the bootstrap, can break down in high-dimensional settings. For example, in regression when the number of subjects is much less than the number of predictors , one cannot naively appeal to asymptotic normality and resampling from the data may not provide an adequate characterization of uncertainty.
Most penalization approaches have a Bayesian interpretation as corresponding to the mode of a posterior distribution obtained under a shrinkage prior. For example, the wildly popular Lasso/ regularization approach to regression is equivalent to maximum a posteriori (MAP) estimation under a Gaussian linear regression model having a double exponential (Laplace) prior on the coefficients. Given this connection, it is natural to ask whether we can use the entire posterior distribution to provide a probabilistic measure of uncertainty. In addition to providing a characterization of uncertainty, a Bayesian perspective has distinct advantages in terms of tuning parameter choice, allowing key penalty parameters to be marginalized over the posterior distribution instead of relying on cross-validation. In addition, by inducing penalties through shrinkage priors, important new classes of penalties can be discovered that may outperform usual -type choices.
From a frequentist perspective, we would like to be able to choose a default shrinkage prior that leads to similar optimality properties to those shown for penalization and other approaches. However, instead of showing that a particular penalty leads to a point estimator having a minimax optimal rate of convergence under sparsity assumptions, we would rather like to show that the entire posterior distribution concentrates at the optimal rate, i.e., the posterior probability assigned to a shrinking neighborhood of the true parameter value converges to one, with the neighborhood size proportional to the frequentist minimax rate.
An amazing variety of shrinkage priors have been proposed in the Bayesian literature; however with essentially no theoretical justification for the performance of these priors in the high-dimensional settings for which they were designed. and provided conditions on the prior for asymptotic normality of linear regression coefficients allowing the number of predictors to increase with sample size , with requiring a very slow rate of growth and assuming . These results required the prior to be sufficiently flat in a neighborhood of the true parameter value, essentially ruling out shrinkage priors. considered shrinkage priors in providing simple sufficient conditions for posterior consistency in linear regression where the number of variables grows slower than the sample size, though no rate of contraction was provided.
In studying posterior contraction in high-dimensional settings, it becomes clear that it is critical to understand several aspects of the prior distribution on the high-dimensional space, including (but not limited to) the prior concentration around sparse vectors and the implied dimensionality of the prior. Specifically, studying the reduction in dimension induced by shrinkage priors is challenging due to the lack of exact zeros, with the prior draws being sparse in only an approximate sense. This substantial technical hurdle has prevented any previous results (to our knowledge) on posterior concentration in high-dimensional settings for shrinkage priors. In fact, investigating these properties is critically important not just in studying frequentist optimality properties of Bayesian procedures but for Bayesians in obtaining a better understanding of the behavior of their priors and choosing associated hyperparameters. Without such technical handle, it becomes an art to use intuition and practical experience to indirectly induce a shrinkage prior, while focusing on Gaussian scale families for computational tractability. Some beautiful classes of priors have been proposed by among others, with showing that essentially all existing shrinkage priors fall within the Gaussian global-local scale mixture family. One of our primary goals is to obtain theory that can allow evaluation of existing priors and design of novel priors, which are appealing from a Bayesian perspective in allowing incorporation of prior knowledge and from a frequentist perspective in leading to minimax optimality under weak sparsity assumptions.
A new class of shrinkage priors
For concreteness, we focus on the widely studied normal means problem (see, for example, and references therein); although most of the ideas developed in this paper generalize directly to high-dimensional linear and generalized linear models. In the normal means setting, one aims to estimate a -dimensional meanFollowing standard practice in this literature, we use to denote the dimensionality and it should not be confused with the sample size. based on a single observation corrupted with i.i.d. standard normal noise:
In a fully Bayesian framework, it is common to place a beta prior on , leading to a beta-Bernoulli prior on the model size, which conveys an automatic multiplicity adjustment . In a beautiful recent paper, established that prior (3) with an appropriate beta prior on and suitable tail conditions on leads to a minimax optimal rate of posterior contraction, i.e., the posterior concentrates most of its mass on a ball around of squared radius of the order of :
where is a constant and . obtained consistency in model selection using point-mass mixture priors with appropriate data-driven hyperparameters.
2 Global-local shrinkage rules
Although point mass mixture priors are intuitively appealing and possess attractive theoretical properties, posterior sampling requires a stochastic search over an enormous space, leading to slow mixing and convergence . Computational issues and consideration that many of the s may be small but not exactly zero has motivated a rich literature on continuous shrinkage priors; for some flavor of the vast literature refer to . noted that essentially all such shrinkage priors can be represented as global-local (GL) mixtures of Gaussians,
where controls global shrinkage towards the origin while the local scales allow deviations in the degree of shrinkage. If puts sufficient mass near zero and is appropriately chosen, GL priors in (5) can intuitively approximate (3) but through a continuous density concentrated near zero with heavy tails.
GL priors potentially have substantial computational advantages over point mass priors, since the normal scale mixture representation allows for conjugate updating of and in a block. Moreover, a number of frequentist regularization procedures such as ridge, lasso, bridge and elastic net correspond to posterior modes under GL priors with appropriate choices of and . For example, one obtains a double-exponential prior corresponding to the popular or lasso penalty if is an exponential distribution. However, unlike point mass priors (3), many aspects of shrinkage priors are poorly understood, with the lack of exact zeroes compounding the difficulty in studying basic properties, such as prior expectation, tail bounds for the number of large signals, and prior concentration around sparse vectors. Hence, subjective Bayesians face difficulties in incorporating prior information regarding sparsity, and frequentists tend to be skeptical due to the lack of theoretical justification.
This skepticism is somewhat warranted, as it is clearly the case that reasonable seeming priors can have poor performance in high-dimensional settings. For example, choosing in prior (3) leads to an exponentially small prior probability of assigned to the null model, so that it becomes literally impossible to override that prior informativeness with the information in the data to pick the null model. However, with a beta prior on , this problem can be avoided . In the same vein, if one places i.i.d. priors on the entries of , then the induced prior on is highly concentrated around leading to misleading inferences on almost everywhere. Although these are simple examples, similar multiplicity problems can transpire more subtly in cases where complicated models/priors are involved and hence it is fundamentally important to understand properties of the prior and the posterior in the setting of (1).
3 Dirichlet-kernel priors
Let us revisit the global-local specification (5). Integrating out the local scales ’s, (5) can be equivalently represented as a global scale mixture of a kernel ,
These traditional choices lead to a kernel which is bounded in a neighborhood of zero. However, if one instead uses a half Cauchy prior , then the resulting horseshoe kernel is unbounded with a singularity at zero. This phenomenon coupled with tail robustness properties leads to excellent empirical performance of the horseshoe. However, the joint distribution of under a horseshoe prior is understudied and further theoretical investigation is required to understand its operating characteristics. One can imagine that it concentrates more along sparse regions of the parameter space compared to common shrinkage priors since the singularity at zero potentially allows most of the entries to be concentrated around zero with the heavy tails ensuring concentration around the relatively small number of signals.
In (7), is any symmetric (about zero) unimodal density with exponential or heavier tails; for computational purposes, we shall restrict attention to the class of kernels that can be represented as scale mixture of normals . While previous shrinkage priors in the literature obtain marginal behavior similar to the point mass mixture priors (3), our construction aims at resembling the joint distribution of under a two-component mixture prior. Constraining on restrains the degrees of freedom of the ’s, offering better control on the number of dominant entries in . In particular, letting for a suitably chosen allows (7) to behave like (3) jointly, forcing a large subset of to be simultaneously close to zero with high probability.
We focus on the Laplace kernel from now on for concreteness, noting that all the results stated below can be generalized to other choices. The corresponding hierarchical prior given ,
is referred to as a Dirichlet–Laplace prior, denoted .
To understand the role of , we undertake a study of the marginal properties of conditional on , integrating out . The results are summarized in Proposition 2.1 below.
Thus, marginalizing over , we obtain an unbounded kernel , so that the marginal density of has a singularity at 0 while retaining exponential tails. A proof of Proposition 2.1 can be found in the appendix.
The parameter plays a critical role in determining the tails of the marginal distribution of ’s. We consider a fully Bayesian framework where is assigned a prior on the positive real line and learnt from the data through the posterior. Specifically, we assume a prior on with . We continue to refer to the induced prior on implied by the hierarchical structure,
as a Dirichlet–Laplace prior, denoted .
There is a recent frequentist literature on including a local penalty specific to each coefficient. The adaptive Lasso relies on empirically estimated weights that are plugged in. instead propose to sample the penalty parameters from a posterior, with a sparse point estimate obtained for each draw. These approaches do not produce a full posterior distribution but focus on sparse point estimates.
4 Posterior computation
The proposed class of DL priors leads to straightforward posterior computation via an efficient data augmented Gibbs sampler. Note that the prior (9) can be equivalently represented as
We detail the steps in the normal means setting noting that the algorithm is trivially modified to accommodate normal linear regression, robust regression with heavy tailed residuals, probit models, logistic regression, factor models and other hierarchical Gaussian cases. To reduce auto-correlation, we rely on marginalization and blocking as much as possible. Our sampler cycles through (i) , (ii) , (iii) and (iv) . We use the fact that the joint posterior of is conditionally independent of given . Steps (ii) - (iv) together give us a draw from the conditional distribution of , since
Steps (i) – (iii) are standard and hence not derived. Step (iv) is non-trivial and we develop an efficient sampling algorithm for jointly sampling . Usual one at a time updates of a Dirichlet vector leads to tremendously slow mixing and convergence, and hence the joint update in Theorem 2.2 is an important feature of our proposed prior; a proof can be found in the Appendix. Consider the following parametrization for the three-parameter generalized inverse Gaussian (giG) distribution: if for .
The joint posterior of has the same distribution as , where are independently distributed according to a distribution, and .
The summary of each step are finally provided below.
To sample , draw independently from a distribution with
The conditional posterior of can be sampled efficiently in a block by independently sampling from an inverse-Gaussian distribution with .
Sample the conditional posterior of from a distribution.
To sample , draw independently with and set with .
Concentration properties of Dirchlet–Laplace priors
The formulation (10) is analytically convenient since the joint distribution factors as a product of marginals and the marginal density can be obtained analytically in Proposition 3.1 below. The proof follows from standard properties of the modified Bessel function ; a proof is sketched in the Appendix.
The marginal density of for any is given by
is the modified Bessel function of the second kind.
Figure 1 plots the marginal density (11) to compare with other common shrinkage priors.
If instead, then (12) holds when .
A proof of Theorem 3.1 can be found in Section 6. To best of our knowledge, Theorem 3.1 is the first result obtaining posterior contraction rates for a continuous shrinkage prior in the normal means setting or the closely related high-dimensional regression problem. Theorem 3.1 posits that when the parameter in the Dirichlet–Laplace prior is chosen, depending on the sample size, to be for any small, the resulting posterior contracts at the minimax rate (2), provided . Using the Cauchy–Schwartz inequality, for and the bound on implies that . Hence, the condition in Theorem 3.1 permits each non-zero signal to grow at a rate, which is a fairly mild assumption. Moreover, in a recent technical report, the authors showed that a large subclass of global-local priors (5) including the Bayesian lasso lead to a sub-optimal rate of posterior convergence; i.e., the expression in (12) converges to whenever . Therefore, Theorem 3.1 indeed provides a substantial improvement over a large class of GL priors.
The choice will be evident from the various auxiliary results in Section 3.1, specifically Lemma 3.3 and Theorem 3.4. The conclusion of Theorem 3.1 continues to hold when under an additional mild assumption on the sparsity . In Table 1 of Section 4, detailed empirical results are provided with as a default choice.
The lower bound result alluded to in the previous paragraph precludes GL priors with polynomial tails, such as the horseshoe. We hope to address the polynomial tails case elsewhere, though based on strong empirical performance, we conjecture that the horseshoe leads to the optimal posterior contraction in a much broader domain compared to the Bayesian lasso and other common shrinkage priors.
In this section, we state a number of properties of the DL prior which provide a better understanding of the joint prior structure and also crucially help us in proving Theorem 3.1.
We first provide useful bounds on the joint density of the DL prior in Lemma 3.2 below; a proof can be found in the Appendix.
If for small, then
If for large, then
It is evident from Figure 1 that the univariate marginal density has an infinite spike near zero. We quantify the probability assigned to a small -neighborhood of the origin in Lemma 3.3 below.
A proof of Lemma 3.3 can be found in the Appendix.
for some constant . If instead, then (15) holds when .
The posterior compressibility property in Theorem 3.4 ensures that the dimensionality of the posterior distribution of (in an approximate sense) doesn’t substantially overshoot the true dimensionality of , which together with the bounds on the joint prior density near zero and infinity in Lemma 3.2 delivers the minimax rate in Theorem 3.1.
Simulation Study
The squared error loss corresponding to the posterior median averaged across simulation replicates is provided in Table 1. To offer further grounds for comparison, we have also tabulated the results for Lasso (LS), Empirical Bayes median (EBMed) as in The EBMed procedure was implemented using the package . , posterior median with a point mass prior (PM) as in and the posterior median corresponding to the horseshoe prior . For the fully Bayesian analysis using point mass mixture priors, we use a complexity prior on the subset-size, with and independent standard Laplace priors for the non-zero entries as in .Given a draw for , a subset of size is drawn uniformly. Set for all and draw i.i.d. from standard Laplace. The beta-bernoulli priors in (3) induce a similar prior on the subset size.
Even in this succinct summary of the results, a wide difference between the Bayesian Lasso and the proposed is observed in Table 1, vindicating our theoretical results. The horseshoe performs similarly as the . The superior performance of the prior can be attributed to its strong concentration around the origin. However, in cases where there are several relatively small signals, the prior can shrink all of them towards zero. In such settings, depending on the practitioner’s utility function, the singularity at zero can be softened using a prior for a larger value of . In the next set of simulations, we report results for , whence computational gains arise as the distribution of in (iv) turns out to be inverse-Gaussian (iG), for which exact samplers are available. In practice, one could as well use a discrete uniform prior on as in Section 5.
For visual illustration and comparison, we finally present the results from a single replicate in the first simulation setting with , and in Figure 2 & 3. The blue circles indicate the entries of , while the red circles correspond to the posterior median of . The shaded region corresponds to a point wise credible interval for .
Prostate data application
We consider a popular dataset from a microarray experiment consisting of expression levels for genes for normal control subjects and patients diagnosed with prostate cancer. The data takes the form of a matrix with the th entry corresponding to the expression level for gene on patient ; the first columns correspond to the normal control subjects with the remaining for the cancer patients. The goal of the study is to discover genes whose expression levels differ between the prostate cancer patients (treatment) and normal subjects (control). A two sample -test with degrees of freedom was implemented for each gene and the resulting t-statistic was converted to a -statistic . Under the null hypothesis of no difference in expression levels between the treatment and control group for the th gene, the null distribution of is . Figure 4 shows a histogram of the -values, comparing it to a density with a multiplier chosen to make the curve integrate to the same area as the histogram. The shape of the histogram suggests the presence of certain interesting genes .
The classical Bonferroni correction for multiple testing flags only genes as significant, while the two-group empirical Bayes method of found significant genes, being much less conservative. The local Bayes false discovery rate (fdr) control method identified genes as non-null. For detailed analysis of this dataset using existing methods, refer to .
To apply our method, we set up a normal means model and assign a prior. Instead of fixing , we use a discrete uniform prior on supported on the interval , with the support points of the form . Such a fully Bayesian approach allows the data to dictate the choice of the tuning parameter which is only specified up to a constant by the theory and also avoids potential numerical issues arising from fixing when is large. Updating is straightforward since the full conditional distribution of is again a discrete distribution on the chosen support points.
We implemented the Gibbs sampler in Section 2.4 for 10,000 draws discarding a burn-in of 5,000. Mixing and convergence of the Gibbs sampler was satisfactory based on examination of trace plots, with the 5,000 retained samples having an effective sample size of 2369.2 averaged across the ’s. The computational time per iteration scaled approximately linearly with the dimension. The posterior mode of was at .
In this application, we expect there to be two clusters of s, with one concentrated closely near zero corresponding to genes that are effectively not differentially expressed and another away from zero corresponding to interesting genes for further study. As a simple automated approach, we cluster s at each MCMC iteration using kmeans with 2 clusters. For each iteration, the number of non-zero signals is then estimated by the smaller cluster size out of the two clusters. A final estimate () of the number of non-zero signals is obtained by taking the mode over all the MCMC iterations. The largest (in absolute magnitude) entries of the posterior median are identified as the non-zero signals.
Using the above selection scheme, our method declared 128 genes as non-null. Interestingly, out of the 128 genes, 100 are common with the ones selected by EBMed. Also all the 54 genes obtained using FDR control form a subset of the selected genes. Horseshoe is overly conservative; it selected only gene (index: 610) using the same clustering procedure; the selected gene was the one with the largest effect size (refer to Table 11.2 in ).
Proof of Theorem 3.1
For a sequence of positive real numbers to be chosen later, let . Define . Let
Let . The proof of Theorem 3.1 is completed by deriving an upper bound to in the following Lemma 6.1 akin to Lemma 5.4 in .
.
Let denote the -dimensional Euclidean ball of radius centered at zero and denote its volume. For the sake of brevity, denote , so that . The numerator of (18) can be clearly bounded above by . Since the set is contained in the ball and , the denominator of (18) can be bounded below by , where . Putting together these inequalities and invoking Lemma 3.2, we have
Using (see Lemma 5.3 in ) and , we can bound from above by . Therefore, we have
Now, since and , we have and hence . Substituting in (20), we have
Substituting the upper bound for obtained in Lemma 6.1, and noting that and , the expression in the left hand side of (16) can be bounded above by
When , the conclusion of Lemma 6.1 remains unchanged and the proof of Theorem 3.4 does not require . The rest of the proof remains exactly the same.
Appendix
When , marginally. Hence, the marginal distribution of given is proportional to
Substituting so that , the above integral reduces to
In the general case, marginally. Substituting as before, the marginal density of is proportional to
The above integral can clearly be bounded below by a constant multiple of
The above expression clearly diverges to infinity as by the monotone convergence theorem.
Proof of Theorem 2.2
Integrating out , the joint posterior of has the form
We now state a result from the theory of normalized random measures (see, for example, (36) in ). Suppose are independent random variables with having a density on . Let with . Then, the joint density of supported on the simplex has the form
where . Setting in (22), we get
We aim to equate the expression in (23) with the expression in (21). Comparing the exponent of gives us . The other requirement is also satisfied, since . The proof is completed by observing that corresponds to a when .
Proof of Proposition 3.1
Proof of Lemma 3.2
Letting , we have .
We first prove (13). Since , and hence , is monotonically decreasing in , and for all , we have . Using for small and for small, we have from (11) that and hence .
We next prove (14). Noting that for large (section 9.7 of ), we have from (11) that for large. Using Cauchy–Schwartz inequality twice, we have , which implies .
Proof of Lemma 3.3
where is a constant independent of . Using a bound for the incomplete gamma function from Theorem 2 of ,
for small. The proof is completed by noting that is bounded above by a constant and for small enough.
Proof of Theorem 3.4
where and respectively denote the numerator and denominator of the expression in (26). Observe that
where is a subset of as in Lemma 5.2 of (replacing by and by ) defined as
with for some sequence of positive real numbers . We set here. With this choice, from (27),