The Horseshoe Estimator: Posterior Concentration around Nearly Black Vectors

Introduction

for independent normal random variables $\varepsilon_{i}$ with mean zero and variance $\sigma^{2}$. The vector $\theta=(\theta_{1},\ldots,\theta_{n})$ is assumed to be sparse, in the ‘nearly black’ sense that the number of nonzero means

In general one would use a posterior distribution both for recovery and for uncertainty quantification. For the first, a measure of centre, such as a median or mode, suffices. For the second, one typically employs a credible set, which is defined as a central set of prescribed posterior probability. For realistic uncertainty quantification it is necessary that the posterior contracts to its center at the same rate as the posterior median or mode approaches the true parameter.

In this paper we study the posterior distribution resulting from the horseshoe prior, which is a one-component prior, introduced in (Carvalho, Polson and Scott, 2010, 2009) and expanded upon in (Scott, 2011; Polson and Scott, 2012a, b). It combines a pole at zero with Cauchy-like tails. The corresponding estimator does not face the computational issues of the point mass mixture models. Carvalho, Polson and Scott (2010) already showed good behaviour of the horseshoe estimator in terms of Kullback-Leibler risk when the true mean is zero. Datta and Ghosh (2013) proved some optimality properties of a multiple testing rule induced by the horseshoe estimator. In this paper, we prove that the horseshoe estimator achieves the minimax quadratic risk, possibly up to a multiplicative constant. We furthermore prove that the posterior variance is of the order of the minimax risk, and thus the posterior contracts at the minimax rate around the underlying mean vector. These results are proven under the assumption that the number $p_{n}$ of nonzero parameters is known. However, we also provide conditions under which the horseshoe estimator combined with an empirical Bayes estimator still attains the minimax rate, when $p_{n}$ is unknown.

The horseshoe prior

In this section, we give an overview of some known properties of the horseshoe estimator which will be relevant to the remainder of our discussion. The horseshoe prior for a parameter $\theta$ modelling an observation $Y\sim\mathcal{N}(\theta,\sigma^{2}I_{n})$ is defined hierarchically (Carvalho, Polson and Scott, 2010):

for $i=1,\ldots,n$, where $C^{+}(0,1)$ is a standard half-Cauchy distribution. The parameter $\tau$ is assumed to be fixed in this paper, rendering the $\theta_{i}$ independent a priori. The corresponding density $p_{\tau}$ increases logarithmically around zero, while its tails decay quadratically. The posterior density of $\theta_{i}$ given $\lambda_{i}$ and $\tau$ is normal with mean $(1-\kappa_{i})y_{i}$, where $\kappa_{i}=\frac{1}{1+\tau^{2}\lambda_{i}^{2}}$. Hence, by Fubini’s theorem:

If $\tau=1$, this reduces to a Be$(\tfrac{1}{2},\tfrac{1}{2})$ distribution, which looks like a horseshoe. As illustrated in Figure 1, decreasing $\tau$ skews the prior distribution on $\kappa_{i}$ towards one, corresponding to more mass near zero in the prior on $\theta_{i}$ and a stronger shrinkage effect in $T_{\tau}(y)$.

where $\Phi_{1}(\alpha,\beta,\gamma;x,y)$ denotes the degenerate hypergeometric function of two variables (Gradshteyn and Ryzhik, 1965).

An unanswered question so far has been how $\tau$ should be chosen. Intuitively, $\tau$ should be small if the mean vector is very sparse, as the horseshoe prior will then place more of its mass near zero. By approximating the posterior distribution of $\tau^{2}$ given $\kappa=(\kappa_{1},\ldots,\kappa_{n})$ in case a prior on $\tau$ is used, Carvalho, Polson and Scott (2010) show that if most observations are shrunk near zero, $\tau$ will be very small with high probability. They suggest a half-Cauchy prior on $\tau$. Datta and Ghosh (2013) implemented this prior on $\tau$ and their plots of posterior draws for $\tau$ at various sparsity levels indicate the expected relationship between $\tau$ and the sparsity level: the posterior distribution of $\tau$ tends to concentrate around smaller values when the underlying mean vector is sparser. As will be discussed further in the next section, the value $\tau=\frac{p_{n}}{n}$ (up to a log factor) is optimal in terms of mean square error and posterior contraction rates.

In case $\tau$ is estimated empirically, as will be considered in section 4, the horseshoe estimator can be computed by plugging this estimate into expression (1), thereby avoiding the use of MCMC. Other aspects of the posterior, such as the posterior variance, can be computed using such a plug-in procedure as well. Polson and Scott (2012a) and Polson and Scott (2012b) consider computation of the horseshoe estimator based on the representation in terms of degenerate hypergeometric functions, as these can be efficiently computed using converging series of confluent hypergeometric functions. They report unproblematic computations for $\tau^{2}$ between $\frac{1}{1000}$ and $1000$. A second option is to apply a quadrature routine to the integral representation in (1). As the continuity and symmetry of $T_{\tau}(y)$ in $y$ can be taken advantage of when computing the horseshoe estimator for a large number of observations, the complexity of these computations mostly depends on the value of $\tau$. Both approaches will be slower for smaller values of $\tau$. Hence, if we use the (estimated) sparsity level $\frac{p_{n}}{n}$ (up to a log factor) for $\tau$, the computation of the horseshoe estimator will be slower if there are fewer nonzero parameters. As noted by Scott (2010), problems arise in Gibbs sampling precisely when $\tau$ is small as well. Hence care needs to be taken with any computational approach if $\frac{p_{n}}{n}$ is suspected to be very close to zero.

The performance of the horseshoe prior, with additional priors on $\tau$ and $\sigma^{2}$, in various simulation studies has been very promising. Carvalho, Polson and Scott (2010) simulated sparse data where the nonzero components were drawn from a Student-$t$ density and found that the horseshoe estimator systematically beat the MLE, the double-exponential (DE) and normal-exponential-gamma (NEG) priors, and the empirical Bayes model due to Johnstone and Silverman (2004) in terms of square error loss. Only when the signal was neither sparse nor heavy-tailed did the MLE, DE and NEG priors have an edge over the horseshoe estimator. In similar experiments in (Carvalho, Polson and Scott, 2009; Polson and Scott, 2012a) the horseshoe prior outperformed the DE prior, while behaving similarly to a heavy-tailed discrete mixture. In a wavelet-denoising experiment under several noise levels and loss functions, the horseshoe estimator compared favorably to the discrete wavelet transform and the empirical Bayes model (Polson and Scott, 2010). Bhattacharya et al. (2012) applied several shrinkage priors to data with the underlying mean vector consisting of zeroes and fixed nonzero values and found the posterior median of the horseshoe prior performing better in terms of squared error than the Bayesian Lasso (BL), the Lasso, the posterior median of a point mass mixture prior as in (Castillo and Van der Vaart, 2012) and the empirical Bayes model proposed by Johnstone and Silverman (2004), and comparable to their proposed Dirichlet-Laplace (DL) prior with parameter $\frac{1}{n}$. Results in (Armagan, Dunson and Lee, 2013) are similar. In a second simulation setting, Bhattacharya et al. (2012) generated data of length $n=1000$, with the first ten means equal to 10, the next 90 equal to a number $A\in\{2,\ldots,7\}$ and the remainder equal to zero. In this simulation, the horseshoe prior beat the BL (except when $A=2$) and the DL prior with parameter $\frac{1}{n}$ (except when $A=7$), while performing similarly to the DL prior with parameter $\frac{1}{2}$. It is worthy of note that Koenker (2014) generated data according to the same scheme and applied the empirical Bayes procedures due to Martin and Walker (2014) (EBMW) and Koenker and Mizera (2014) (EBKM) to it. The MSE of EBMW was lower than that of the horseshoe prior for $A\in\{5,6,7\}$, while that of EBKM was much lower in all cases.

Mean square error and bounds on the posterior variance

Suppose $Y\sim\mathcal{N}(\theta_{0},\sigma^{2}I_{n})$. Then the estimator $T_{\tau}(y)$ satisfies

for $\tau\to 0$, as $n,p_{n}\to\infty$ and $p_{n}=o(n)$

By the minimax risk result in (Donoho et al., 1992), we also have a lower bound:

as $n,p_{n}\to\infty$ and $p_{n}=o(n)$. The choice $\tau=\left(\frac{p_{n}}{n}\right)^{\alpha}$, for $\alpha\geq 1$, leads to an upper bound (2) of order $p_{n}\log(n/p_{n})$, with (as can be seen from the proof) a multiplicative constant of at most $4\alpha\sigma^{2}$. Thus, for this choice of $\tau$, we have:

The horseshoe estimator therefore performs well as a point estimator, as it attains the minimax risk (possibly up to a multiplicative constant of at most 2 for $\alpha=1$). This may seem surprising, as the prior does not include a point mass at zero to account for the assumed sparsity in the underlying mean vector. Theorem 3.1 shows that the pole at zero of the horseshoe prior mimics the point mass well enough, while the heavy tails ensure that large observations are not shrunk too much.

An upper bound on the rate of contraction of the posterior can be obtained through an upper bound on the posterior variance. The posterior variance can be expressed as:

Suppose $Y\sim\mathcal{N}(\theta_{0},\sigma^{2}I_{n})$. Then the variance of the posterior distribution corresponding to the horseshoe prior satisfies

for $\tau\to 0$, as $n,p_{n}\to\infty$ and $p_{n}=o(n)$.

Again, the choice $\tau=\left(\frac{p_{n}}{n}\right)^{\alpha}$, for $\alpha\geq 1$ leads to an upper bound (3) of the order of the minimax risk. This result indicates that the posterior contracts fast enough to be able to provide a measure of uncertainty of adequate size around the point estimate. Theorems 3.1 and 3.2 combined allow us to find an upper bound on the rate of contraction of the full posterior distribution, both around the underlying mean vector and around the horseshoe estimator.

Under the assumptions of Theorem 3.1, with $\tau=\left(\frac{p_{n}}{n}\right)^{\alpha}$, $\alpha\geq 1$:

for every $M_{n}\to\infty$ as $n\to\infty$.

Combine Markov’s inequality with the results of Theorems 3.1 and 3.2 for (4), and only with the result of Theorem 3.2 for (5). ∎

for $\tau\to 0$ and $p_{n}=o(n)$, as $n\to\infty$. This lower bound is sharp for vectors $\theta_{0,n}$ with $p_{n}$ entries equal to $a_{n}$ and the remaining entries equal to zero, if $a_{n}$ is such that $|a_{n}|\lesssim 1/\sqrt{\log(1/\tau)}$.

for $\tau\to 0$ and $p_{n}=o(n)$, as $n\to\infty$.

Consider $\tau=\left(\frac{p_{n}}{n}\right)^{\alpha}$. Three cases can be discerned:

$\alpha=1$. Bounds (3) and (6) differ by a factor $\sqrt{\log(n/p_{n})}$, as do (7) and (8). The gap can be closed by choosing $\tau=\frac{p_{n}}{n}\sqrt{\log\frac{n}{p_{n}}}$.

These observations suggest that $\tau=\frac{p_{n}}{n}\sqrt{\log(n/p_{n})}$ is a good choice, because then (2), (3), (6), (7), (8) are all of the order $p_{n}\log(n/p_{n})$, suggesting that the posterior contracts at the minimax rate around both the truth and the horseshoe estimator.

Empirical Bayes estimation of τ𝜏\tau

A natural follow-up question is how to choose $\tau$ in practice, when $p_{n}$ is unknown. As discussed in section 2, the full Bayesian approach suggested by Carvalho, Polson and Scott (2010) performs well in simulations. The analysis of such a hierarchical prior would however require different tools than the ones we have used so far. An empirical Bayes estimate of $\tau$ would be a natural solution, and allows us in practice to use one of the representations in (1) for computations, instead of an MCMC-type algorithm.

Suppose we observe an $n$-dimensional vector $Y\sim\mathcal{N}(\theta_{0},\sigma^{2}I_{n})$ and we use $T_{\widehat{\tau}}(y)$ as our estimator of $\theta_{0}$. If $\widehat{\tau}\in(0,1)$ satisfies the following two conditions for $\tau=\frac{p_{n}}{n}$ or $\tau=\frac{p_{n}}{n}\sqrt{\log(n/p_{n})}$:

The first condition requires that $\widehat{\tau}$ does not overestimate the fraction $\frac{p_{n}}{n}$ of nonzero means (up to a log factor) too much or by a too large probability. If $p_{n}\geq 1$, as we have assumed, then it is satisfied already by $\widehat{\tau}=\frac{1}{n}$ (and $c=1$). According to the last assertion of the theorem, this ‘universal threshold’ yields the rate $p_{n}\log n$ (possibly up to a multiplicative constant). This is equal to the rate of the Lasso estimator with the usual choice of $\lambda=2\sqrt{2\sigma^{2}\log n}$ (Bickel, Ritov and Tsybakov, 2009). However, in the framework where $p_{n}\to\infty$, the estimator $\widehat{\tau}=\frac{1}{n}$ will certainly underestimate the sparsity level. A more natural estimator of $\frac{p_{n}}{n}$ is:

where $c_{1}$ and $c_{2}$ are positive constants. By Lemma A.7, this estimator satisfies the first condition for $\tau=\frac{p_{n}}{n}$ and $\tau=\frac{p_{n}}{n}\sqrt{\log(n/p_{n})}$ if $c_{1}>2,c_{2}>1$ and $p_{n}\to\infty$ or $c_{1}=2,c_{2}>1$ and $p_{n}\gtrsim\log n$. Thus $\max\left\{\widehat{\tau},\frac{1}{n}\right\}$ will also lead to a rate of at most order $p_{n}\log n$ under these conditions. Its behaviour will be explored further in section 5.

The rate can be improved to $p_{n}\log(n/p_{n})$ if the second condition is met as well, which ensures that the sparsity level is not underestimated too much or by a too large probability. As we are not aware of any estimators meeting this condition for all $\theta_{0}$, this condition is currently mostly of theoretical interest. If the true mean vector is very sparse, in the sense that there are relatively few nonzero means or the nonzero means are close to zero, there is not much to be gained in terms of rates by meeting this condition. The extra occurrence of $p_{n}$ relative to the rate $p_{n}\log n$ is of interest only if $p_{n}$ is relatively large. For instance, if $p_{n}\asymp n^{\alpha}$ for $\alpha\in(0,1)$, then $p_{n}\log(n/p_{n})=(1-\alpha)p_{n}\log n$, which suggests a decrease of the proportionality constant in (9), particularly if $\alpha$ is close to one. Furthermore, when $p_{n}$ is large, the constant in (9) may be sensitive to the fine properties of $\widehat{\tau}$, as it depends on $g(n,p_{n})$ (as can be seen in the proof). If $\widehat{\tau}$ seriously underestimates the sparsity level, the corresponding value of $g(n,p_{n})$ from the second condition may be so small that the upper bound on the multiplicative constant before (9) becomes very large. Hence in this case, $\widehat{\tau}$ is required to be close to the proportion $\frac{p_{n}}{n}$ (up to a log factor) with large probability in order to get an optimal rate.

Datta and Ghosh (2013) warned against the use of an empirical Bayes estimate of $\tau$ for the horseshoe prior, because the estimate might collapse to zero. Their references for this statement, Scott and Berger (2010) and Bogdan, Ghosh and Tokdar (2008), indicate that they are thinking of a marginal maximum likelihood estimate of $\tau$. However, an empirical Bayes estimate of $\tau$ does not need to be based on this principle. Furthermore, an estimator that satisfies the second condition from Theorem 4.1 or that is truncated from below by $\frac{1}{n}$, would not be susceptible to this potential problem.

Simulation study

A simulation study provides more insight into the behaviour of the horseshoe estimator, both when using an empirical Bayes procedure with estimator (10) and when using the fully Bayesian procedure proposed by Carvalho, Polson and Scott (2010) with a half-Cauchy prior on $\tau$. For each data point, 100 replicates of an $n$-dimensional vector sampled from a $\mathcal{N}(\theta_{0},I_{n})$ distribution were created, where $\theta_{0}$ had either 20, 40 or 200 (5%, 10% or 50%) entries equal to an integer $A$ ranging from 1 to 10, and all the other entries equal to zero. The full Bayesian version was implemented using the code provided in (Scott, 2010), and the coordinatewise posterior mean was used as the estimator of $\theta_{0}$. For the empirical Bayes procedure, the estimator (10) was used with $c_{1}=2$ and $c_{2}=1$. Performance was measured by squared error loss, which was averaged across replicates to create Figure 2.

With $\tau=\frac{1}{n}$, this leads to the ‘universal threshold’ of $\sqrt{2\sigma^{2}\log n}$, or with $\tau=\left(\frac{p_{n}}{n}\right)^{\alpha}$, a ‘threshold’ at $\sqrt{2\alpha\sigma^{2}\log(n/p_{n})}$. Based on this property and the proofs of the main results, we can divide the underlying parameters into three cases:

Those that are exactly or close to zero, where the observations are shrunk close to zero;

Those that are larger than the threshold, where the horseshoe estimator essentially behaves like the identity;

Those that are close to the ‘threshold’, where the horseshoe estimator is most likely to shrink the observations too much.

The horseshoe estimator performs well in cases (i) and (ii) due to its pole at zero and its heavy tails respectively. The hardest parameters to recover from the noise are those that are close to the threshold, and these are the ones that affect the estimation risk the most. This phenomenon explains the peaks in the graphs of Figure 2 around $A=3.5$.

Concluding remarks

The choice of the global shrinkage parameter $\tau$ is critical towards ensuring the right amount of shrinkage of the observations to recover the underlying mean vector. The value of $\tau=\frac{p_{n}}{n}\sqrt{\log(n/p_{n})}$ was found to be optimal. Theorem 4.1 indicates that quite a wide range of estimators for $\tau$ will work well, especially in cases where the underlying mean vector is sparse. Of course, it should not come as a surprise that an estimator designed to recover sparse vectors will work especially well if the truth is indeed sparse. An interesting extension to this work would be to investigate whether the posterior concentration properties of the horseshoe prior still remain when a hyperprior is placed on $\tau$. The result that $\tau=\frac{p_{n}}{n}$ (up to a log factor) yields optimal rates, together with the simulation results, suggests that in a fully Bayesian approach, a prior on $\tau$ which is restricted to $$ may perform better than the suggested half-Cauchy prior.

The simulation results also indicate that mean vectors with the nonzero means close to the universal threshold are the hardest to recover. In future simulations involving shrinkage rules, it would therefore be interesting to study the challenging case where all the nonzero parameters are at this threshold. The performance of the empirical Bayes estimator (10) leaves something to be desired around the threshold. In additional numerical experiments (not shown), we tried two other estimators of $\tau$. The first was the ‘oracle estimator’ $\widehat{\tau}=\frac{p_{n}}{n}$. For values of the nonzero means well past the ‘threshold’, the behaviour of this estimator was very similar to that of (10). However, before the threshold, the squared error loss of the empirical procedure with the oracle estimator was between that of the full Bayes estimator and empirical Bayes with estimator (10). The second estimator was the mean of the samples of $\tau$ from the full Bayes estimator. The resulting squared error loss was remarkably close to that of the full Bayes estimator, for all values of the nonzero means. Neither of these two estimators is of much practical use. However, their range of behaviours suggests room for improvement over the estimator (10), and it would be worthwhile to study more refined estimators for $\tau$.

An interesting question is what aspects of the horseshoe prior are truly essential towards optimal posterior contraction properties. Our proofs do not elucidate whether the pole at zero of the horseshoe prior is required, or if a prior with heavy tails, and in a sense ‘sufficient’ mass at zero would work as well. The failure of the Lasso to concentrate around the true mean vector at the minimax rate does indicate that heavy tails in itself may not be sufficient, and adding mass at zero solves this problem (Castillo, Schmidt-Hieber and Van der Vaart, 2014; Castillo and Van der Vaart, 2012). It is possible that the pole at zero is inessential, in particular if the global tuning parameter is chosen carefully, for instance by empirical Bayes. If the tuning parameter is chosen by a full Bayes method, the peak may be more essential, depending on its prior.

Acknowledgements

The authors would like to thank two anonymous referees for their helpful suggestions, as well as James Scott for his advice on implementing the full Bayesian version of the horseshoe estimator.

Appendix A Proofs

This section begins with Lemma A.1, providing bounds on some of the degenerate hypergeometric functions appearing in the posterior mean and posterior variance. This is followed by two lemmas that are needed for the proofs of Theorems 3.1 and 3.2: Lemma A.2 provides two upper bounds on the horseshoe estimator and Lemma A.3 gives a bound on the absolute value of the difference between the horseshoe estimator and an observation. We then proceed to the proof of Theorem 3.1, after which Lemma A.4 provides upper bounds on the posterior variance. These upper bounds are then used in the proof of Theorem 3.2. The proof of Theorem 3.4 is given next, followed by Lemmas A.5 and A.6 supporting the proof of Theorem 3.5. This section concludes with the proofs of Theorem 4.1 and Lemma A.7, which both concern the empirical Bayes procedure discussed in section 4.

where (11) and (13) hold for $\tau<1/\sqrt{a}$, (12) holds for $\tau<1$, and (14) and (15) hold for $\tau<1/a$.

Write $\xi=y^{2}/(2\sigma^{2})$. We first note that for $z\geq\tau^{2}$, we have $z\leq\tau^{2}+(1-\tau^{2})z\leq 2z$, while for $z\leq\tau^{2}$, we have $\tau^{2}\leq\tau^{2}+(1-\tau^{2})z\leq 2\tau^{2}$. Hence, we can bound $I_{k}$ from above by:

and from below by half of that quantity. We bound the integral over $[0,\tau^{2}]$ in all cases by bounding the factor $e^{\xi z}$ by 1 or $e^{\tau^{2}\xi}$. For the integral over $[\tau^{2},1]$, we first substitute $u=\xi z$, yielding: $\int_{\tau^{2}}^{1}z^{k-1}e^{\xi z}dz=\xi^{-k}\int_{\tau^{2}\xi}^{\xi}u^{k-1}e^{u}du$. For (11) and (13), we split the domain of integration into $\left[\tau^{2}\xi,\frac{\xi}{a}\right]$ and $\left[\frac{\xi}{a},\xi\right]$. For $I_{\frac{3}{2}}$, we bound by:

yielding (11). Similarly, for $I_{\frac{1}{2}}$:

resulting in (13). The bound (12) is obtained similarly, but without splitting up $[\tau^{2}\xi,\xi]$ further, by the inequality

For the bounds on $I_{-\frac{1}{2}}$, we split up the domain of integration $\left[\tau^{2}\xi,\xi\right]$ into $\left[\tau^{2}\xi,\tau\xi\right],\left[\tau\xi,\frac{\xi}{a}\right]$ and $\left[\frac{\xi}{a},\xi\right]$, and then bound by:

If $\tau^{2}<1$, the posterior mean of the horseshoe prior can be bounded above by:

$T_{\tau}(y)\leq ye^{\frac{y^{2}}{2\sigma^{2}}}f(\tau)$, where $f$ is such that $f(\tau)\leq\frac{2}{3}\tau$;

$T_{\tau}(y)\leq y\frac{\frac{2}{3}e^{\tau^{2}\frac{y^{2}}{2\sigma^{2}}}\tau+2e^{\frac{y^{2}}{2a\sigma^{2}}}\left(\frac{1}{\sqrt{a}}-\tau\right)+\frac{2\sqrt{a}\sigma^{2}}{y^{2}}\left(e^{\frac{y^{2}}{2\sigma^{2}}}-e^{\frac{y^{2}}{2a\sigma^{2}}}\right)}{\frac{1}{\tau}+e^{\tau^{2}\frac{y^{2}}{2\sigma^{2}}}\left(\frac{1}{\tau}-\frac{1}{\sqrt{\tau}}\right)+\frac{a\sigma^{2}\sqrt{a}}{y^{2}}\left(e^{\frac{y^{2}}{2a\sigma^{2}}}-e^{\tau\frac{y^{2}}{2\sigma^{2}}}\right)+\frac{\sigma^{2}}{y^{2}}\left(e^{\frac{y^{2}}{2\sigma^{2}}}-e^{\frac{y^{2}}{2a\sigma^{2}}}\right)}$ , for any $a>1$ and $\tau<\tfrac{1}{a}$.

We bound the integrals in the numerator and denominator of expression (1). For the first upper bound, we will use the fact that for $0\leq z\leq 1$, $e^{\frac{y^{2}}{2\sigma^{2}}z}$ is bounded below by 1 and above by $e^{\frac{y^{2}}{2\sigma^{2}}}$. The posterior mean can therefore be bounded by:

By Shafer’s inequality for the arctangent (Shafer, 1966):

which completes the proof for the first upper bound. For the second inequality, we note that, in the notation of Lemma A.1, $T_{\tau}(y)=y\frac{I_{\frac{1}{2}}(y)}{I_{-\frac{1}{2}}(y)}$. The bounds in Lemma A.1 yield the stated inequality. ∎

For $\tau^{2}<1$, the absolute value of the difference between the horseshoe estimator and an observation $y$ can be bounded by a function $h(y,\tau)$ such that for any $c>2$:

We assume $y>0$ without loss of generality. By a change of variables of $x=1-z$:

By following the proof of Watson’s lemma provided in Miller (2006), we can find bounds on the numerator and denominator of the above expression. First define $g(x)=(1-x)^{-\frac{1}{2}}\frac{1}{\textbf{1}-(1-\tau^{2})x}$ and note that by Taylor’s theorem, $g(x)=g(0)+xg^{\prime}(\xi_{x}),$ where $\xi_{x}$ is between 0 and $x$. Let $s$ be any number between 0 and 1. Because $g^{\prime\prime}(x)$ is not negative for $x\in[0,1)$, we have that for $x\in[0,s]$, $s\in(0,1)$: $g^{\prime}(0)\leq g^{\prime}(x)\leq g^{\prime}(s)$. The numerator can then be bounded by:

where $h_{1}(y,\sigma,s)=4\sigma^{4}-2\sigma^{2}(sy^{2}+2\sigma^{2})e^{-\frac{sy^{2}}{2\sigma^{2}}}$, $h_{2}(y,\sigma,s)=16\sigma^{6}-2\sigma^{2}(s^{2}y^{4}+4s\sigma^{2}y^{2}+8\sigma^{4})e^{-\frac{sy^{2}}{2\sigma^{2}}}$ and $h_{3}(\tau)=\arctan\left(\frac{\sqrt{1-\tau^{2}}}{\tau}\right)\tau^{-1}(1-\tau^{2})^{-\frac{3}{2}}-(1-\tau^{2})^{-1}$. The denominator can similarly be bounded by:

where $h_{4}(y,\sigma,s)=2\sigma^{2}-2\sigma^{2}e^{-\frac{sy^{2}}{2\sigma^{2}}}$ and $h_{5}(y,\sigma,s)=4\sigma^{4}-2\sigma^{2}e^{-\frac{sy^{2}}{2\sigma^{2}}}(sy^{2}+2\sigma^{2})$. Hence:

For any fixed $\tau$, this bound tends to zero as $y$ tends to infinity. If $\tau\to 0$, the term containing $h_{3}(\tau)$ could potentially diverge. For $s=\frac{2}{3}$ and $y=\sqrt{c\sigma^{2}\log(1/\tau)}$, where $c$ is a positive constant, this term displays the following limiting behaviour as $\tau\to 0$:

because $\lim_{\tau\downarrow 0}\arctan\left(\frac{\sqrt{1-\tau^{2}}}{\tau}\right)(1-\tau^{2})^{-\frac{3}{2}}=\frac{\pi}{2}$, $\lim_{\tau\downarrow 0}\frac{\tau}{1-\tau^{2}}=0$ and the factor $\left(c\sigma^{2}\log(1/\tau)\right)^{\frac{3}{2}}\tau^{\frac{c}{3}-1}$ tends to zero as $\tau\downarrow 0$ if $\frac{c}{3}-1>0$ and infinity otherwise. The condition $c>3$ is related to the choice of $s=\frac{2}{3}$ and can be improved to any constant strictly greater than $2$ by choosing $s$ appropriately close to one. Hence, we find that the absolute value of the difference between the posterior mean and an observation can be bounded by a function $h(y,\tau)$ with the desired property. ∎

Nonzero parameters Denote $\zeta_{\tau}=\sqrt{2\sigma^{2}\log(1/\tau)}$. We will show

for all nonzero $\theta_{i}$, which can be done by bounding $\sup_{y}|T_{\tau}(y)-y|$:

Lemma A.3 yields the following bound on the difference between the observation and the horseshoe estimator: $|T_{\tau}(y)-y|\leq h(y,\tau)$, where $h(y,\tau)$ is such that $\lim_{\tau\downarrow 0}\sup_{|y|>c\zeta_{\tau}}h(y,\tau)=0$ for any $c>1$. Combining this with the inequality $|T_{\tau}(y)-y|\leq|y|$, we have as $\tau\to 0$:

which implies (16), as $|T_{\tau}(y)|\leq|y|$:

Parameters equal to zero We split up the term for the zero means into two parts:

where $\zeta_{\tau}=\sqrt{2\sigma^{2}\log(1/\tau)}$. For the first term, we have, by the first bound in Lemma A.2:

where the identity $\frac{d}{dy}ye^{\frac{y^{2}}{2\sigma^{2}}}=\frac{y^{2}}{\sigma^{2}}e^{\frac{y^{2}}{2\sigma^{2}}}+e^{\frac{y^{2}}{2\sigma^{2}}}$ was used to bound $\int_{-\zeta_{\tau}}^{\zeta_{\tau}}y^{2}e^{\frac{y^{2}}{2\sigma^{2}}}dy$. For the second term, because $|T_{\tau}(y)|\leq|y|$ for all $y$, we have by the identity $y^{2}\phi(y)=\phi(y)-\frac{d}{dy}[y\phi(y)]$, and by Mills’ ratio:

where the last inequality holds for $\zeta_{\tau}>\sigma^{2}$. If we apply this inequality and combine this upper bound with the upper bound on the first term, we find, for $\zeta_{\tau}>\sigma^{2}$ (corresponding to $\tau<e^{-\frac{\sigma^{2}}{2}}$):

Hence, for $\tau<e^{-\frac{\sigma^{2}}{2}}$:

Conclusion By (16) and (19), we find for $\tau<e^{-\frac{\sigma^{2}}{2}}$:

The posterior variance when using the horseshoe prior can be expressed as:

${\rm Var}(\theta\mid y)\leq\sigma^{2}+y^{2}$;

${\rm Var}(\theta\mid y)\leq\left(\frac{\sigma^{2}}{y}+y\right)T_{\tau}(y)-T_{\tau}(y)^{2}$.

where $m(y)$ is the density of the marginal distribution of $y$. Equality (20) can be found by combining the expressions

with the equality $T_{\tau}(y)=y+\sigma^{2}\frac{m^{\prime}(y)}{m(y)}$. The first upper bound is implied by the property $|T_{\tau}(y)|<|y|$ and the fact that $(1-z)^{2}\leq 1$ for $z\in$. The second upper bound can be demonstrated by noting that $(1-z)^{2}\leq 1-z$ for $z\in$ and hence:

Hence, as $\tau\to 0$: ${\rm Var}(\theta|y)\lesssim\sigma^{2}$, for any $|y|$ that increases as least as fast as $\zeta_{\tau}$ when $\tau$ decreases. Now suppose $|y|\leq\zeta_{\tau}$. Then, by the bound ${\rm Var}(\theta\mid y)\leq\sigma^{2}+y^{2}$ from Lemma A.4, we find:

Zero means By the bound ${\rm Var}(\theta\mid y)\leq\sigma^{2}+y^{2}$, we find for $c\geq 1$:

For $|y|<c\zeta_{\tau}$, we consider the upper bound ${\rm Var}(\theta\mid y)\leq\left(\frac{\sigma^{2}}{y}+y\right)T_{\tau}(y)-T_{\tau}(y)^{2}$ from Lemma A.4. From this bound, we get ${\rm Var}(\theta\mid y)\leq\frac{\sigma^{2}}{y}T_{\tau}(y)+yT_{\tau}(y)$. Hence:

We bound the first integral from (A) by applying the first bound on $T_{\tau}(y)$ from Lemma A.2:

because $f(\tau)\leq\frac{2}{3}\tau$. For the second term in (A), we first note that the second bound from Lemma A.2 can be relaxed to:

By expanding $(1-z)^{2}z^{-\frac{1}{2}}=z^{-\frac{1}{2}}-2z^{\frac{1}{2}}+z^{\frac{3}{2}}$, we see that the final term in (20) is equal to:

As $\frac{T_{\tau}(y)}{y}$ is non-negative, we can bound the posterior variance from below by the final two terms in (20). By the above equality, this yields the following lower bound:

where $I_{k}$ is as in Lemma A.1. We now use the bounds from Lemma A.1 with $a=2$ and take $\xi$ equal to $c\log(1/\tau)$ for some nonnegative constant $c$. Then $e^{\xi}=\frac{1}{\tau^{c}}$ and $e^{\frac{\xi}{2}}=\frac{1}{\tau^{\frac{c}{2}}}$. Taking for each bound on $I_{k}$, $k\in\left\{\tfrac{3}{2},\tfrac{1}{2},-\tfrac{1}{2}\right\}$, the term that diverges fastest as $\tau$ approaches zero, we find that the lower bound is asymptotically of the order:

By the substitution $x=y^{2}/\zeta_{\tau}^{2},dy=\frac{\sigma\sqrt{\log(1/\tau)}}{\sqrt{2x}}dx$, we find:

where in the last step, we used $\tau^{\tau x}\geq\tau^{\tau}\geq e^{-\frac{1}{e}}$ for $x\in,\tau\in(0,1]$. By plugging this into (25), we find that as $\tau\to 0$:

finishing the proof for the first statement of the theorem.

We first consider $|y_{i}|>\zeta_{\tau}$. By the first bound of Lemma A.4:

The first integral from (A) can be split into two parts by splitting up the factor $\sigma^{2}+y^{2}$, the first of which can be bounded, by substituting $x=(y-a_{n})/\sigma$ and applying Mills’ ratio:

The second of these integrals is, by $y^{2}=(y-a_{n})^{2}-a_{n}^{2}+2a_{n}y$, equal to:

By substituting $x=-y$ in the second integral from (A) and then applying the same inequalities to it as to the first integral, the following bound is obtained:

If $a_{n}\lesssim 1/\zeta_{\tau}$, then $e^{-\frac{a_{n}^{2}}{2\sigma^{2}}}e^{\frac{\zeta_{\tau}a_{n}}{\sigma}}=\mathcal{O}(1)$ and $\zeta_{\tau}\pm a_{n}=\mathcal{O}(\zeta_{\tau})$, yielding an upper bound on (33) of order $\tau\zeta_{\tau}$.

We now consider $|y_{i}|\leq\zeta_{\tau}$. We use the second bound of Lemma A.4:

Applying inequality $\frac{1}{y}T_{\tau}(y)\leq\frac{2}{3}\tau e^{\frac{y^{2}}{2\sigma^{2}}}$ from Lemma A.2 to the first integral yields the bound:

If $a_{n}\lesssim 1/\zeta_{\tau}$, we have $a_{n}\zeta_{\tau}=\mathcal{O}(1)$ and thus this term will be of order $\tau\zeta_{\tau}$. For the second integral from (A), we use bound (23). This leads to three integrals to be bounded, $I_{1},I_{2}$ en $I_{3}$.

$I_{1},I_{2}$ and $I_{3}$ will all be of no larger order than $\tau\zeta_{\tau}$ if $a_{n}\lesssim 1/\zeta_{\tau}$. ∎

For $k=0$, the statement is immediate. By integration by parts the integral is seen to be equal to $y^{k}e^{y}-e-\int_{1}^{y}ku^{k-1}e^{u}du$. For $k\neq 0$, the latter integral is bounded above by

This is further bounded above by a multiple of $(1\vee y^{k-1})e^{y/2}+y^{k-1}e^{y}$. ∎

Let $I_{k}$ be as in Lemma A.1. There exist functions $R_{k}$ with $\sup_{\zeta_{\tau}/4\leq y\leq 4\zeta_{\tau}}|R_{k}(y)|\to 0$ for $k>0$ and $k=-\tfrac{1}{2}$, such that,

We split the integral in the definition of $I_{k}$ over the intervals $[0,\tau^{2}]$ and $[\tau^{2},1]$. The first interval contributes, uniformly in $y\tau\to 0$,

by the substitution $u=z/\tau^{2}$. The integral tends to $\int_{0}^{1}\frac{u^{k}}{1+u}\,du$, by the dominated convergence theorem, for any $k>-1$. The second interval contributes, with the substitution $u=(y^{2}/2\sigma^{2})z$:

In the second integral the argument satisfies $u\geq 1$, and hence $u/((y^{2}\tau^{2}/(2\sigma^{2})+(1-\tau^{2}))\to 1$, uniformly in $u$ and $y\tau\to 0$. Hence

as $y\to\infty$, by Lemma A.5. For the first integral we separately consider the cases $k>0$ and $k=-1/2$. If $k>0$, then $\int_{0}^{1}u^{k-1}e^{u}du$ converges, and hence, by the dominated convergence theorem, uniformly in $y\tau\to 0$,

If $k=-1/2$ , then we substitute $v=2\sigma^{2}u/(\tau^{2}y^{2})$ and rewrite the integral as

Denote $\zeta_{\tau}=\sqrt{2\sigma^{2}\log(1/\tau)}$ and assume that $\theta_{i}=\gamma\zeta_{\tau}$ for $i=1,\ldots,p_{n}$ and $\theta_{i}=0$ for $i=p_{n}+1,\ldots,n$. We prove (7) by proving that there exists a positive constant $c_{1}(\gamma)$ such that

If (36) holds, we have, by Jensen’s inequality:

as $\tau\to 0$. In addition, we have $T_{\tau}(y)=yI_{\frac{1}{2}}(y)/I_{-\frac{1}{2}}(y)$. For $|y|=\sqrt{2\sigma^{2}c\log(1/\tau)}$, with $c>1$, the lower bound (12) on $I_{\frac{1}{2}}(y)$ behaves as $(\sigma^{2}/y^{2})e^{\frac{y^{2}}{2\sigma^{2}}}$, while the upper bound (15) on $I_{-\frac{1}{2}}(y)$ behaves as $(2a\sqrt{a}\sigma^{2}/y^{2})e^{\frac{y^{2}}{2\sigma^{2}}}$, as $\tau\to 0$. Therefore, for $|y|>\zeta_{\tau}$, we have $T_{\tau}(y)\gtrsim y$. Thus, we can bound by:

By combining the lower bounds (37) and (38) with the upper bound (2), we arrive at (7). For the posterior variance, we already have $\sum_{i=p_{n}+1}^{n}{\rm Var}(\theta_{i}\mid Y_{i})\asymp(n-p_{n})\tau\zeta_{\tau}$ by (24) and (26). Expression (8) can therefore be proven by showing that there exists a positive constant $c_{2}(\gamma)$ such that:

We shall show that the first and third integrals are negligible, while the second gives the approximation in (36). On the domain of the second integral, we have $\zeta_{\tau}/4\leq\zeta_{\tau}+y\leq 4\zeta_{\tau}$, so we can apply Lemma A.6 to see that this integral is asymptotic to

where $c_{1}=\int_{0}^{\infty}z^{-1/2}(1-z)^{-1}dz$ and $c_{2}=\int_{0}^{1}z^{1/2}(1-z)^{-1}dz$. On $[-\zeta_{\tau}/2,3\zeta_{\tau}]$:

By the substitution $u=\zeta_{\tau}y-2\sigma^{2}\log\zeta_{\tau}$, the remaining integral is equal to, with $a_{\tau}=-\frac{\zeta_{\tau}^{2}}{2}-2\sigma^{2}\log\zeta_{\tau}$ and $b_{\tau}=3\zeta_{\tau}^{2}-2\sigma^{2}\log\zeta_{\tau}$:

by the dominated convergence theorem. This yields the approximation in (36), with $c_{1}(\gamma)=(2\sigma/\sqrt{2\pi})\int_{-\infty}^{\infty}e^{\frac{\gamma u}{\sigma^{2}}}/(c_{1}+2\sigma^{2}e^{\frac{u}{\sigma^{2}}})du$.

For the first integral in (40), we use bound 1 from Lemma A.2, and obtain a bound on its absolute value equal to

where the last inequality follows by integration by parts. This is of much smaller order than the second integral from (40). In the third integral of (40), we bound $I_{\frac{1}{2}}(\zeta_{\tau}+y)/I_{-\frac{1}{2}}(\zeta_{\tau}+y)$ by 1, giving the upper bound

by Mills’ ratio. This is also of much smaller order than the second integral from (40), thus concluding the proof of (36).

Proof of (39) By expanding the term $(1-z)^{2}$ in the numerator of the final term of (20), the posterior variance can be seen to be equal to:

Because $I_{\frac{1}{2}}(y)/I_{-\frac{1}{2}}(y)$ can be interpreted as the mean of the density proportional to $z\to z^{-1/2}e^{y^{2}z/(2\sigma^{2})}/(\tau^{2}+(1-\tau^{2})z)$, and $I_{\frac{3}{2}}(y)/I_{-\frac{1}{2}}(y)$ as the second moment, it follows that the term in square brackets in (43) is nonnegative. By (43), we write:

The first term in (A) is as (40), except without the factor $(\zeta_{\tau}+y)$. Following the same steps as the proof of (36), we see that it is smaller than a multiple of $\zeta_{\tau}^{-1}$ times the bound on (40), so it is of the order $\zeta_{\tau}^{2\gamma-3}\tau^{(1-\gamma)^{2}}$. The first and third integrals of the second term of (A) are also negligible. For the first, we use that the expression in square brackets is nonnegative and bounded above by $I_{\frac{3}{2}}(y)/I_{-\frac{1}{2}}(y)$, which in turn is bounded above by $I_{\frac{1}{2}}(y)/I_{-\frac{1}{2}}(y)$. We bound as in (42), with the difference that the leading factor is $(\zeta_{\tau}+y)^{2}$ instead of $(\zeta_{\tau}+y)$. This leads to the order $\zeta_{\tau}\tau^{(1-\gamma)^{2}+\gamma}$, much smaller than the claimed rate. For the third integral, we can bound the term in square brackets by 1 and use Mills’ ratio to see that it is of the order $\zeta_{\tau}\tau^{(4-\gamma)^{2}}$.

We are left with the middle integral of the second term of (A). On the domain of this integral, by Lemma A.6:

where $c_{3}=\int_{0}^{1}z^{3/2}(1+z)^{-1}dz$, and $c_{1}$ is as in (41). We see that $I_{\frac{3}{2}}(y)/I_{-\frac{1}{2}}(y)$ and $I_{\frac{1}{2}}(y)/I_{-\frac{1}{2}}(y)$ are asymptotic to the same function on this domain. Since $A/(A+B)-A^{2}/(A+B)^{2}=AB/(A+B)^{2}$, it follows that up to $\mathcal{O}(\tau)$, the middle integral is asymptotic to

We substitute $u=\zeta_{\tau}y-2\sigma^{2}\log\zeta_{\tau}$ to reduce this to

This is asymptotic to expression (39), with $c_{2}(\gamma)=(2\sigma c_{1}/\sqrt{2\pi})\int_{-\infty}^{\infty}e^{\frac{\gamma u}{\sigma^{2}}}/(c_{1}+2\sigma^{2}e^{\frac{u}{\sigma^{2}}})^{2}du$.

where (17) was used in the second line, and $Z$ follows a standard normal distribution. If $A$ is such that $\widehat{\tau}>\tau$ holds with probability one on $A$, we can use the inequality $\zeta_{\widehat{\tau}}<\zeta_{\tau}$ if $\widehat{\tau}>\tau$ to find:

and then bound each of terms on the right hand side. For the nonzero means, we take $c=1$, while for the zero means, we consider $c\geq 1$. Note that for $\zeta_{\widehat{\tau}}$ to be well-defined, we need $\widehat{\tau}\leq 1$ and consequently, when we consider $\widehat{\tau}>c\tau$, we must have $c\tau<1$.

If $\widehat{\tau}\leq\tau$, the inequality $\zeta_{\widehat{\tau}}^{2}\leq\zeta_{\tau}^{2}$ needed for (46) does not hold. For this case, we assume that $\widehat{\tau}\geq g(n,p_{n})$ with probability one, for some function $g(n,p_{n})$, corresponding to $\zeta_{\widehat{\tau}}\leq\sqrt{-2\sigma^{2}\log g(n,p_{n})}$. Then we find by (45):

By (47) and (48), we have for $\theta_{i}\neq 0$:

By this inequality combined with the inequality $\log(y+1)<y$, we have:

By (46) and (A), we get for any $c\geq 1$ such that $c\tau<1$:

Now suppose $\widehat{\tau}\leq c\tau$ for some $c\geq 1$ such that $c\tau<1$. First note that $|T_{\tau}(y)|$ increases monotonically in $\tau$, as is clear from

Because sign$(T_{\widehat{\tau}}(y_{i}))$ = sign$(T_{c\tau}(y_{i}))$ and $0\leq|T_{\widehat{\tau}}(y_{i})|\leq|T_{c\tau}(y_{i})|$, we have:

And thus, by (18), we have for $\theta_{i}=0$:

The function $x\log\left(1+\frac{1}{x}\right)$ is monotonically increasing in $x$ for $x\in$. Hence, with the choice $\tau=\frac{p_{n}}{n}$ or $\tau=\frac{p_{n}}{n}\sqrt{\log(n/p_{n})}$, the conditions stated in the theorem are sufficient for (54) to be bounded by the minimax squared error rate in the worst case.

Suppose $Y_{i}\sim\mathcal{N}(\theta_{i},\sigma^{2}),i=1,\ldots,p_{n}$ and $Y_{i}\sim\mathcal{N}(0,\sigma^{2}),i=p_{n}+1,\ldots,n$ and define

For $c_{1}=2$, this quantity will exceed $\log\frac{n}{p_{n}}$ if $p_{n}\gtrsim\log n$. If $c_{1}>2$, we require $((c_{2}-1)p_{n}+1)(\tfrac{c_{1}}{2}-1)\geq 1$, which is certainly satisfied if $p_{n}\to\infty$.

References