Unified Bayesian theory of sparse linear regression with nuisance parameters^††thanks: Research is partially supported by a Faculty Research and Professional Development Grant from College of Sciences of North Carolina State University.

As briefly discussed above, the form in (1) is general and includes many interesting statistical models. Here we provide specific examples belonging to (1) in detail. In Section 5, these examples will be used to apply the main results developed in this study.

The rest of this paper is organized as follows. In Section 2, some notations are introduced and a prior distribution on sparse regression coefficients is specified. Sections 3–4 provide our main results on the posterior contraction, the Bernstein-von Mises phenomenon, and selection consistency of the posterior distribution. In Section 5, our general theorems are applied to the examples considered above to derive the posterior asymptotic properties in each specific example. All technical proofs are provided in Appendix.

Here we describe the notations we use throughout this paper. For a vector $\theta=(\theta_{j})\in\mathbb{R}^{p}$ and a set $S\subset\{1,\dots,p\}$ of indices, we write $S_{\theta}=\{j:\theta_{j}\neq 0\}$ to denote the support of $\theta$, $s\coloneqq|S|$ (or $s_{\theta}\coloneqq|S_{\theta}|$) to denote the cardinality of $S$ (or $S_{\theta}$), and $\theta_{S}=\{\theta_{j}:j\in S\}$ and $\theta_{S^{c}}=\{\theta_{j}:j\notin S\}$ to separate components of $\theta$ using $S$. In particular, the support of the true parameter $\theta_{0}$ and its cardinality are written as $S_{0}$ and $s_{0}\coloneqq|S_{0}|$, respectively. The notation $\lVert\theta\rVert_{q}=(\sum_{j}|\theta_{j}|^{q})^{1/q}$, $1\leq q<\infty$, stands for the $\ell_{q}$-norm and $\lVert\theta\rVert_{\infty}=\max_{j}|\theta_{j}|$ denotes the maximum norm. We write $\rho_{\min}(A)$ and $\rho_{\max}(A)$ for the minimum and maximum eigenvalues of a square matrix $A$, respectively. For a matrix $X=(\!(x_{ij})\!)$, let $\lVert X\rVert_{\rm sp}=\rho_{\max}^{1/2}(X^{T}X)$ stand for the spectral norm and $\lVert X\rVert_{\rm F}=(\sum_{i,j}x_{ij}^{2})^{1/2}$ stand for the Frobenius norm of $X$. We also define a matrix norm $\lVert X\rVert_{\ast}=\max_{j}\lVert X_{\cdot j}\rVert_{2}$ for $X_{\cdot j}$ the $j$th column of $X$, which is used for compatibility conditions. The column space of $X$ is denoted by ${\rm span}(X)$. For further convenience, we write $\varsigma_{\min}(X)=\rho_{\min}^{1/2}(X^{T}X)$ for the minimum singular value of $X$. The notation $X_{S}$ means the submatrix of $X$ with columns chosen by $S$. For sequences $a_{n}$ and $b_{n}$, $a_{n}\lesssim b_{n}$ (or $b_{n}\gtrsim a_{n}$) stands for $a_{n}\leq Cb_{n}$ for some constant $C>0$ independent of $n$, and $a_{n}\asymp b_{n}$ means $a_{n}\lesssim b_{n}\lesssim a_{n}$. These inequalities are also used for relations involving constant sequences.

For given parameters $\theta$ and $\eta$, we write the joint density as $p_{\theta,\eta}=\prod_{i=1}^{n}p_{\theta,\eta,i}$ for $p_{\theta,\eta,i}$ the density of the $i$th observation vector $Y_{i}$. In particular, the true joint density is expressed as $p_{0}=\prod_{i=1}^{n}p_{0,i}$ for $p_{0,i}\coloneqq p_{\theta_{0},\eta_{0},i}$ with the true parameters $\theta_{0}$ and $\eta_{0}$. The notation $\mathbb{E}_{0}$ denotes the expectation operator with the true density $p_{0}$. For two probability measures $P$ and $Q$, let $\lVert P-Q\rVert_{\rm TV}$ denote the total variation between $P$ and $Q$. For two $n$-variate densities $f\coloneqq\prod_{i=1}^{n}f_{i}$ and $g\coloneqq\prod_{i=1}^{n}g_{i}$ of independent variables, denote the average Rényi divergence (of order $1/2$) by $R_{n}(f,g)=-n^{-1}\sum_{i=1}^{n}\log\int\sqrt{f_{i}g_{i}}$.

For any $\eta_{1},\eta_{2}\in\mathbb{H}$, we define $d_{n}^{2}(\eta_{1},\eta_{2})=d_{A,n}^{2}(\eta_{1},\eta_{2})+d_{B,n}^{2}(\eta_{1},\eta_{2})$ for the two squared pseudo-metrics:

For compatibility conditions, the uniform compatibility number $\phi_{1}$ and the smallest scaled singular value $\phi_{2}$ are defined as

We write $Y^{(n)}=(Y_{1}^{T},\dots,Y_{n}^{T})^{T}$ for the observation vector, $n_{\ast}=\sum_{i=1}^{n}m_{i}$ for the dimension of $Y^{(n)}$, and $\Theta=\mathbb{R}^{p}$ for the parameter space of $\theta$. Lastly, for a (pseudo-)metric space $({\cal F},d)$, let $N(\epsilon,{\cal F},d)$ denote the $\epsilon$-covering number, the minimal number of $\epsilon$-balls that cover $\cal F$.

In this subsection, we specify a prior distribution for the high-dimensional regression coefficients $\theta$. A prior for $\eta$ should satisfy the conditions required for the main results, so its specific characterization is deferred to Section 3. On the other hand, the prior for $\theta$ specified here is always good for our purposes and satisfies all requirements.

We first select a dimension $s$ from a prior $\pi_{p}$, and then randomly choose $S\subset\{1,\dots,p\}$ for given $s$. A nonzero part $\theta_{S}$ of $\theta$ is then selected from a prior $g_{S}$ on $\mathbb{R}^{s}$ while $\theta_{S^{c}}$ is fixed to zero. The resulting prior specification for $(S,\theta)$ is formulated as

where $\delta_{0}$ is the Dirac measure at zero on $\mathbb{R}^{p-s}$ with suppressed dimensionality. For the prior $\pi_{p}$ on the model dimensions, we consider a prior satisfying the following: for some constants $A_{1},A_{2},A_{3},A_{4}>0$,

Examples of priors satisfying (3) can be found in Castillo and van der Vaart, [9] and Castillo et al., [8]. For the prior $g_{S}$, the $s$-fold product of the exponential power density is considered, where the regularization parameter is allowed to vary with $p$ and $\lVert X\rVert_{\ast}$, i.e.,

for some constants $L_{1},L_{2},L_{3}>0$. The order of $\lambda$ is important in that it determines the boundedness requirement of the true signal $\theta_{0}$ (see condition (C3) below). A particularly interesting case is obtained when $\lambda$ is set to the lower bound $\lVert X\rVert_{\ast}/(L_{1}p^{L_{2}})$. Then the boundedness condition becomes very mild by choosing $L_{2}$ sufficiently large. When $\lambda$ is set to the upper bound, the boundedness condition is still reasonably mild. However, it can actually be relaxed if the true signal is known to be small enough, though we do not pursue this generalization in this study. In Section 4, we shall see that values of $\lambda$ that do not increase too fast are in fact necessary for a distributional approximation and selection consistency.

The goal of this subsection is to study posterior contraction of $\theta$ relative to the $\ell_{1}$- and $\ell_{2}$-metrics. To do so, we derive the posterior contraction rate with respect to the average Rényi divergence $R_{n}(f,g)$, and then the rates for $\theta$ relative to more concrete metrics will be recovered from the Rényi contraction.

To proceed, we first need to examine a dimensionality property of the support of $\theta$. The following theorem shows that the posterior distribution is concentrated on models of relatively small sizes.

Compared to the literature [e.g., 8, 23, 3], the rate in Theorem 1 is floored by the extra term $n\bar{\epsilon}_{n}^{2}/\log p$. This arises from the presence of a nuisance parameter in the model formulation. To minimize its impact, a prior on $\eta$ should be chosen such that (C2) holds for as small $\bar{\epsilon}_{n}$ as possible; a suitable choice induces the (nearly) optimal contraction rate.

Using the basic results in Theorem 1, the next theorem obtains the rate at which the posterior distribution contracts at the truth with respect to the average Rényi divergence. The theorem requires additional assumptions on a prior.

1. (C5)

   For $s_{\star}\coloneqq s_{0}\vee(n\bar{\epsilon}_{n}^{2}/\log p)$ with $\bar{\epsilon}_{n}$ satisfying (C2), a sufficiently large $B>0$, and some sequences $\gamma_{n}$ and $\epsilon_{n}\geq\sqrt{s_{\star}\log(p\vee\overline{m}\vee\gamma_{n})/n}$ satisfying $\epsilon_{n}^{2}/\overline{m}\rightarrow 0$, there exists a subset ${\cal H}_{n}\subset{\cal H}$ such that

   	$\displaystyle\min_{1\leq i\leq n}\inf_{\eta\in{\cal H}_{n}}\rho_{\min}(\Delta_{\eta,i})$	$\displaystyle\geq\frac{1}{\gamma_{n}},$		(6)
   	$\displaystyle\log N\left(\frac{1}{6\overline{m}\gamma_{n}n^{3/2}},{\cal H}_{n},d_{n}\right)$	$\displaystyle\lesssim n\epsilon_{n}^{2},$		(7)
   	$\displaystyle e^{Bs_{\star}\log p}\Pi({\cal H}\setminus{\cal H}_{n})$	$\displaystyle\rightarrow 0.$		(8)

The above conditions are related to the classical ones in the literature (e.g., see Theorem 2.1. of Ghosal et al., [15]). Condition (6) requires that for every $i\leq n$, the minimum eigenvalue of $\Delta_{\eta,i}$ is not too small on a sieve ${\cal H}_{n}$. Although $\gamma_{n}$ can be any positive sequence, a sequence increasing exponentially fast makes the entropy in (7) too large, resulting in a suboptimal rate $\epsilon_{n}$. If $\gamma_{n}$ can be chosen to be smaller than $p$ and $\overline{m}$, then this does not lead to any deterioration of the rate in $\epsilon_{n}$. The entropy condition (7) is actually stronger than needed. Scrutinizing the proof of the theorem, one can see that the entropy appearing in the theorem is obtained using pieces that are smaller than those giving the exponentially powerful test in Lemma 2 in Appendix. However, the covering number with those pieces looks more complicated and the form in (7) suffices for all examples in the present paper. Lastly, condition (8) implies that the outside of a sieve ${\cal H}_{n}$ should possess sufficiently small prior mass to kill the factor $s_{\star}\log p$ arising from the lower bound of the denominator of the posterior distribution. In fact, conditions similar to (C2), (7) and (8) are also required for the prior of $\theta$. By reading the proof, it is easy to see that the prior (2) explicitly satisfies the analogous conditions on an appropriately chosen sieve.

We want to sharpen the rate $\epsilon_{n}\geq\sqrt{s_{\star}\log(p\vee\overline{m}\vee\gamma_{n})/n}$ as much as possible. In most instances, $\gamma_{n}$ can be chosen such that $\log\gamma_{n}\lesssim\log p$. This is trivially satisfied if $\gamma_{n}$ is some polynomial in $n$ as in the examples in this paper. If $p$ is known to increase much faster than $n$, e.g., $\log p\asymp n^{c}$ for some $c\in(0,1)$, then $\gamma_{n}$ need not be a polynomial in $n$ and the condition can be met more easily with a sequence that grows even faster. Note also that we typically have $\log\overline{m}\lesssim\log p$ in most cases. These postulates lead to $\epsilon_{n}\geq\sqrt{(s_{\star}\log p)/n}$. Indeed, it is often possible to choose $\epsilon_{n}=\sqrt{(s_{\star}\log p)/n}$, which is commonly guaranteed by choosing an appropriate sieve ${\cal H}_{n}$ and a prior. The condition will be made precise in (C5^∗) below for recovery and we only consider the situation that $\epsilon_{n}=\sqrt{(s_{\star}\log p)/n}$ in what follows.

Although Theorem 2 provides the basic results for posterior contraction, it does not give precise interpretations for the parameters $\theta$ and $\eta$ themselves, because of the abstruse expression of the average Rényi divergence. The contraction rates with respect to more concrete metrics are recovered under some additional conditions. Under the additional assumption $a_{n}\epsilon_{n}^{2}\rightarrow 0$, it can be shown that Theorem 1 and Theorem 2 explicitly imply that for the set

with a sufficiently large constant $M_{1}$, the posterior mass of ${\cal A}_{n}$ goes to one in probability (see the proof of Theorem 3). To complete the recovery, we need to separate the sum of squares of the mean into $\lVert X(\theta-\theta_{0})\rVert_{2}$ and $nd_{A,n}^{2}(\eta,\eta_{0})$, which requires an additional condition. The conditions required for the recovery are clarified as follows.

1. (C5^∗)

   While $\log\overline{m}\lesssim\log p$, (C5) holds for $\gamma_{n}$ and $\epsilon_{n}=\sqrt{(s_{\star}\log p)/n}$ such that $\log\gamma_{n}\lesssim\log p$ and $a_{n}\epsilon_{n}^{2}\rightarrow 0$ with $a_{n}$ satisfying (C1).

2. (C6)

   For $s_{\star}$ satisfying (C5^∗), there exists $\eta_{\ast}\in\mathbb{H}$ such that

   	$\displaystyle\liminf_{n\geq 1}\inf_{(\theta,\eta)\in{\cal A}_{n}}\frac{\sum_{i=1}^{n}(\theta-\theta_{0})^{T}X_{i}^{T}(\xi_{\eta,i}-\xi_{\eta_{\ast},i})}{\lVert X(\theta-\theta_{0})\rVert_{2}^{2}+nd_{A,n}^{2}(\eta,\eta_{\ast})}$	$\displaystyle>-\frac{1}{2},$
   	$\displaystyle d_{A,n}(\eta_{\ast},\eta_{0})$	$\displaystyle\lesssim\sqrt{\frac{s_{\star}\log p}{n}},$

   where $\epsilon_{n}$ in $\mathcal{A}_{n}$ satisfies $\epsilon_{n}=\sqrt{(s_{\star}\log p)/n}$.

By expanding the quadratic term for the mean in ${\cal A}_{n}$, one can see that the separation is possible if (C6) is satisfied. Clearly, (C6) is trivially satisfied if the model has only $X\theta$ for its mean, in which we take $\xi_{\eta,i}-\xi_{\eta_{\ast},i}=\xi_{\eta_{\ast},i}-\xi_{\eta_{0},i}=0$ for every $i\leq n$. In many cases where there exists $\eta^{\prime}\in{\cal H}$ such that $d_{A,n}(\eta^{\prime},\eta_{0})=0$, we can often take $\eta_{\ast}=\eta^{\prime}$ for the second inequality of (C6) to hold automatically.

The following theorem shows that the posterior distribution of $\theta$ and $\eta$ contracts at their respective true values at some rates, relative to more easily comprehensible metrics than the average Rényi divergence. In the expressions, if $K_{1}s_{\star}+s_{0}<1$, the compatibility numbers should be understood be equal to 1 for interpretation.

The thresholds for contraction depend upon the compatibility conditions, which make their implication somewhat vague. As $K_{1}s_{\star}+s_{0}$ is much smaller than $n_{\ast}$, it is not unreasonable to assume that $\phi_{1}(K_{1}s_{\star}+s_{0})$ and $\phi_{2}(K_{1}s_{\star}+s_{0})$ are bounded away from zero, whence the compatibility number is removed from the rates. We refer to Example 7 of Castillo et al., [8] for more discussion. In the next subsection, we will see that one of these restrictions is actually necessary for shape approximation or selection consistency.

Notwithstanding the lack of formal study of minimax rates with additional complications, we still want to match our rates for $\theta$ with those in simple linear regression, which we call the “optimal” rates. In this sense, Theorem 3 only provides the suboptimal rates for $\theta$ if $s_{0}=o(s_{\star})$. Although the theorem gives the optimal results if $s_{0}\log p\gtrsim n\bar{\epsilon}_{n}^{2}$, it is practically hard to check this condition as $s_{0}$ is unknown. If $s_{0}$ is known to be nonzero, the desired conclusion is trivially achieved as soon as $n\bar{\epsilon}_{n}^{2}/\log p\lesssim 1$. The following corollary, however, shows that the optimal rates are still available even if $s_{0}=0$, with restrictions on $\bar{\epsilon}_{n}$ and the prior.

The corollary is useful in limited situations, especially when a parametric rate is available for a nuisance parameter. Even if $n\bar{\epsilon}_{n}^{2}=\log n$, we need to further assume that $\log n=o(\log p)$, i.e., the ultra high-dimensional setup, to conclude that (a) holds, while we can always apply (b) because $\log n\lesssim\log p$. Although assertion (b) holds for any $s_{0}\geq 0$ if $A_{4}$ is chosen sufficiently large, its specific threshold is not directly available. Indeed, by carefully reading the proof of Theorem 1 together with Lemma 1 in Appendix, one can see that the threshold depends on unknown constant bounds for the eigenvalues of the true covariance matrix in (C4). Still, (b) holds for any $A_{4}>0$ if $s_{0}>0$. We believe that the assumption $s_{0}>0$ is very mild, and hence simply apply (b) with this assumption to conclude the optimal contraction for models with finite dimensional nuisance parameters. The optimal rates can still be achieved for any $s_{0}\geq 0$ by verifying the conditions in the following subsection. With finite dimensional nuisance parameters, we do not pursue this direction as it seems an overkill considering the mildness of the assumption $s_{0}>0$, though those conditions are actually required for the Bernstein-von Mises theorem and selection consistency in Section 4.

In semiparametric situations with high- or infinite-dimensional nuisance parameters, none of (a) and (b) generally works unless $p$ increases sufficiently fast. Still, the optimal rates can be achieved under stronger conditions using the semiparametric theory, as the following subsection provides.

Recall that only suboptimal rates may be available from Theorem 3 if $s_{0}\log p\lesssim n\bar{\epsilon}_{n}^{2}$. In many semiparametric situations, however, it is often possible to obtain parametric rates for finite dimensional parameters under stronger conditions, even when there are infinite-dimensional nuisance parameters in a model [4, 7]. It has also been shown that a similar argument holds in some high-dimensional semiparametric regression models [10]. Therefore, it is naturally of interest to examine under what conditions we can replace $s_{\star}$ by $s_{0}$ in the rates for $\theta$, even if $s_{0}\log p\lesssim n\bar{\epsilon}_{n}^{2}$. Similar to other semiparametric settings [4, 10], this can be established by the semiparametric theory, but requires stronger conditions than those in traditional fixed dimensional parametric cases because of the high-dimensions of the parameters in our setup.

To proceed, some additional conditions are required for technical reasons, which are made for the size of $\bar{\epsilon}_{n}$ as the optimal rates are automatically attained if $s_{0}\log p\gtrsim n\bar{\epsilon}_{n}^{2}$. Still, in a practical sense, the conditions almost always need to be verified to reach the optimal rates, since only oracle rates are generally available and we do not know which term is greater.

In what follows, we write $\bar{s}_{\star}\coloneqq n\bar{\epsilon}_{n}^{2}/\log p$ for $\bar{\epsilon}_{n}$ satisfying the conditions of Theorem 3 through the definition of $\epsilon_{n}$. We first assume the following condition on the uniform compatibility number.

1. (C3)

   For a sufficiently large $M$, the uniform compatibility number $\phi_{1}(M\bar{s}_{\star}+s_{0})$ is bounded away from zero.

This condition is weaker than assuming that the smallest scaled singular value $\phi_{2}(M\bar{s}_{\star}+s_{0})$ is bounded away from zero, as we have $\phi_{1}(s)\geq\phi_{2}(s)$ for any $s>0$ by the Cauchy-Schwarz inequality. We will also resort on a slightly stronger condition with respect to $\phi_{1}$ for a distributional approximation in the following section. In this sense, our condition is weaker than those for Theorem 4 of Castillo et al., [8]. Condition (C3) is not restrictive as (C5^∗) requires $s_{\star}=o(n)$; we again refer to Example 7 of Castillo et al., [8].