Li 2019 - Kullback-Leibler Risk Bound and Bias Theory

Summary

Thm 3.2: for sparse ${\Omega }_0$ the GHS Cesaro-average KL risk bound has a log factor of order $\log \{\log (\cdot )/(\dots )\}$ on each zero entry (prior mass $\sim \log (n^{1/2}p)/(n^{1/2}p)$) rather than $1/(n^{1/2}p)$, so its bound is smaller than for any prior bounded near the true value (lasso, SCAD). Thm 4.1: the GHS posterior mean of an entry is $(1-E[Z])\hat{\omega }$ with $E[Z]=O(1/{\hat{\omega }}^2)$, hence within $O(1/\hat{\omega })$ of least squares, while the double-exponential posterior mean stays biased by a constant.

Overview

All methods’ bounds diverge when $p^2/n\to \infty$; GHS diverges more slowly.

Main Content

Lemma 3.1 (Barron 1988)

With KL neighbourhood $A_{ϵ}=\{\Omega :D(p_{{\Omega }_0}\Vert p_{\Omega })\le ϵ\}$ and prior $\nu$ with $\nu (A_{ϵ})>0$ for all $ϵ>0$, the Cesaro-average risk of the Bayes posterior-mean density estimate ${\hat{p}}_n$ satisfies $R_n=\frac{1}{n}{\sum }_{j=1}^n\mathbb{E}D(p_{{\Omega }_0}\Vert {\hat{p}}_j)\le ϵ-\frac{1}{n}\log \nu (A_{ϵ})$ for all $ϵ>0$.

Theorem 3.2 (KL risk bound under GHS)

Let $y\sim N(0,{\Omega }_0)$, ${\sum }_{i,j}{\sigma }_{ij}^0=Mp$ for a constant $M$, and assume a cube of half-width $2/(M{n}^{1/2}p)$ around ${\Omega }_0$ lies in ${\mathcal{S}}_p$. Then $R_n\le \frac{1}{n}-\frac{1}{n}\log \nu (A_{1/n})$ and: (1) under GHS, $\log \nu (A_{1/n})>p_0\log \left\{\frac{1}{M{n}^{1/2}p}{C}_1\log \right(2M{n}^{1/2}p\left)\right\}+p_1\log \frac{C_2}{n^{1/2}p}$ with $p_0$ zero and $p_1$ non-zero entries of ${\Omega }_0$; (2) for any other prior density continuous, bounded and strictly positive near ${\omega }_{ij}^0$: $\log \nu (A_{1/n})>p^2\log \frac{K_1}{n^{1/2}p}$.

Reading: densities bounded near zero (Bayesian lasso; the SCAD-penalty-equivalent density by Polson & Scott 2012) get prior measure of order $1/(n^{1/2}p)$ per dimension, whereas the horseshoe’s unbounded spike gives order $\log (n^{1/2}p)/(n^{1/2}p)$ on each zero entry and $1/(n^{1/2}p)$ on non-zero ones. Rissanen (1986) lower bound: no estimator avoids divergence in high dimension; recommendation: screen variables first.

Theorem 4.1 (element-wise bias)

Assume the diagonal is known; by partial regression $y_p\mid y_{-p}\sim N(-{\omega }_{pp}^{-1}{\omega }_{p(-p)}{y}_{-p},{\omega }_{pp}^{-1})$, and the least-squares estimate ${\hat{\omega }}_{p(-p)}\mid Y_{(-p)}\sim N({\omega }_{p(-p)}^0,{\omega }_{pp}^0(Y_{(-p)}^{\top }{Y}_{(-p)}{)}^{-1})$. Scale ${\omega }_{pj}$ and ${\hat{\omega }}_{pj}$ by $\{{\omega }_{pp}^0(Y_{(-p)}^{\top }{Y}_{(-p)}{)}_{jj}^{-1}{\}}^{-1/2}$ to get ${\omega }_{pj}^{\prime },{\hat{\omega }}_{pj}^{\prime }$. (1) GHS: $E({\omega }_{pj}^{\prime }\mid Y,\tau )=(1-E(Z_{pj})){\hat{\omega }}_{pj}^{\prime }$, $Z_{pj}\sim$ Compound Confluent Hypergeometric $(1,\frac{1}{2},1,{\hat{\omega }}_{pj}^{\prime 2}/2,1,{\theta }_{pj}^{-1})$ on $(0,1)$, with ${\theta }_{pj}={\omega }_{pp}^0(Y_{(-p)}^{\top }{Y}_{(-p)}{)}_{jj}^{-1}{\tau }^{-2}$. If ${\hat{\omega }}_{pj}^{\prime 2}/2>1$: $E(Z_{pj})<4(C_1+C_2){\theta }_{pj}(1+{\hat{\omega }}_{pj}^{\prime 2}/2)/{\hat{\omega }}_{pj}^{\prime 4}$, with $C_1=1-2/e\approx 0.26$, $C_2=\Gamma (1/2)\Gamma (2)/\Gamma (2.5)=0.75$; so $E(Z_{pj})=O(1/{\hat{\omega }}_{pj}^{\prime 2})$ as ${\hat{\omega }}_{pj}^{\prime }\to \infty$. (2) Double-exponential prior: $E({\omega }_{pj}^{\prime }\mid Y)={\hat{\omega }}_{pj}^{\prime }+\frac{d}{d{\hat{\omega }}_{pj}^{\prime }}\log m_{\text{lasso}}({\hat{\omega }}_{pj}^{\prime })$, with limit $\pm a$, $a=2^{1/2}/(nv)$, $v$ the prior variance: a non-vanishing shrinkage bias. (3) ${\hat{\omega }}_{pj}^{\prime 2}\mid Y_{(-p)}\sim \text{noncentral }{\chi }_1^2({\omega }_{pj}^{\prime 2})$; ${\omega }_{pj}^{\prime 2}={\omega }_{pj}^{02}{\omega }_{pp}^{0-1}\{(Y_{(-p)}^{\top }{Y}_{(-p)}{)}_{jj}^{-1}{\}}^{-1}$ with the last factor Gamma-distributed (mean prop. to $n-p+2$, mode prop. to $n-p$) when $n>p-1$.

Interpretation: with $n\gg p$ and a non-zero true entry, GHS is close to unbiased with high probability; BGL’s bias does not vanish and is larger when sparse (large $a$). Bias is smaller when ${\omega }_{pj0}^2{\omega }_{pp0}^{-1}$ is large even for $n<p$. Proofs: Supplement S.1 (Thm 3.2: cube measure via the horseshoe’s Thm 1 of Carvalho et al. 2010) and S.2 (Thm 4.1, via Bhadra et al. 2016 and Carvalho et al. 2010 Thm 2).

Examples

Shrinkage factor: for the scalar horseshoe, $E[\kappa \mid y]$ with $\kappa =1/(1+{\lambda }^2{\tau }^2)$ has a U-shaped (horseshoe) prior; large $|\hat{\omega }|$ pushes $\kappa \to 0$ (no shrinkage). Numerical illustration with $\tau =0.1$: posterior mean $\approx \hat{\omega }(1-E[\kappa ])$; quick Monte Carlo is left to the reader.

import numpy as np
rng = np.random.default_rng(0)
def hs_post_mean(y, tau=0.1, s=1.0, M=400000):      # normal-means horseshoe, importance sampling over lambda
    lam = np.abs(rng.standard_cauchy(M)); v = (lam*tau)**2
    w = np.exp(-0.5*y**2/(s**2 + v))/np.sqrt(s**2 + v)           # marginal likelihood of y given lambda
    shrink = v/(s**2 + v)                                       # E[theta | y, lambda] = shrink * y
    return y*np.sum(w*shrink)/np.sum(w)
print([round(hs_post_mean(y), 2) for y in (0.5, 1.5, 3., 6., 12.)])   # approaches y for large y

Connections

Limits of shrinkage vs selection: Meinshausen Buhlmann 2010 - Stability Selection - Overview controls false positives differently.

See Also

Li 2019 - Graphical Horseshoe Estimator - Overview