Honest confidence sets for high-dimensional regression by projection and shrinkage

TL;DR

This paper introduces a novel projection and shrinkage method to construct honest confidence sets in high-dimensional linear regression, effectively handling unknown sparsity and achieving optimal adaptation.

Contribution

It develops a new approach that separates signals into strong and weak groups, ensuring honest confidence sets without sparsity constraints and adapting to the true sparsity.

Findings

Outperforms existing methods in finite sample scenarios.

Achieves honest confidence sets over the full parameter space.

Adapts to the optimal rate when the true parameter is sparse.

Abstract

The issue of honesty in constructing confidence sets arises in nonparametric regression. While optimal rate in nonparametric estimation can be achieved and utilized to construct sharp confidence sets, severe degradation of confidence level often happens after estimating the degree of smoothness. Similarly, for high-dimensional regression, oracle inequalities for sparse estimators could be utilized to construct sharp confidence sets. Yet the degree of sparsity itself is unknown and needs to be estimated, causing the honesty problem. To resolve this issue, we develop a novel method to construct honest confidence sets for sparse high-dimensional linear regression. The key idea in our construction is to separate signals into a strong and a weak group, and then construct confidence sets for each group separately. This is achieved by a projection and shrinkage approach, the latter implemented via Stein estimation and the associated Stein unbiased risk estimate. Our confidence set is honest over the full parameter space without any sparsity constraints, while its diameter adapts to the optimal rate of $n^{-1/4}$ when the true parameter is indeed sparse. Through extensive numerical comparisons, we demonstrate that our method outperforms other competitors with big margins for finite samples, including oracle methods built upon the true sparsity of the underlying model.

Figures

Equations

$$y=X\beta +\epsilon ,$$

$$\underset{n\to \infty }{\mathrm{liminf}}\underset{\beta \in {\mathbb{R}}^p}{\inf }{\mathbb{P}}_{\beta }\left\{X\beta \in \hat{C}\right\}\ge 1-\alpha ,$$

$$\frac{1}{n}\Vert X\hat{\beta }-X\beta {\Vert }^2\le c\frac{s\log p}{n}$$

$$|{\hat{C}}_a|=O_p\left(n^{-1/4}+\sqrt{s\log p/n}\right)$$

$$\hat{\mu }(y;\tilde{\mu })$$

$$\hat{L}(y;\tilde{\mu })$$

$$\underset{n\to \infty }{\mathrm{liminf}}\underset{\mu \in {\mathbb{R}}^n}{\inf }{\mathbb{P}}_{\mu }\left\{\left|{\sigma }^2\hat{L}-n^{-1}\Vert \hat{\mu }-\mu {\Vert }^2\right|\le c_s(\alpha ){\sigma }^2{n}^{-1/2}\right\}\ge 1-\alpha ,$$

$${\mu }_A=P_A\mu ,{\mu }_{\perp }=P_A^{\perp }\mu =({\mathbf{I}}_n-P_A)\mu ,$$

$$\hat{C}=$$

$$\Vert {\hat{\mu }}_A-{\mu }_A{\Vert }^2=\Vert P_A\epsilon {\Vert }^2\mid A\sim {\sigma }^2{\chi }_k^2.$$

$$r_A^2=c_1{\tilde{r}}_A^2=c_1{\sigma }^2{\chi }_{k,\alpha /2}^2/n,$$

$$\mathbb{P}\left\{\frac{\Vert P_A\mu -{\hat{\mu }}_A{\Vert }^2}{n{r}_A^2}\le 1/c_1\right\}=1-\alpha /2.$$

$${\hat{\mu }}_{\perp }$$

$$\hat{L}$$

$$B=(n-k){\sigma }^2/\Vert y_{\perp }{\Vert }^2.$$

$$\underset{(n-k)\to \infty }{\mathrm{liminf}}\underset{\beta \in {\mathbb{R}}^p}{\inf }\mathbb{P}\left\{\left|{\sigma }^2\hat{L}-(n-k{)}^{-1}\Vert {\hat{\mu }}_{\perp }-{\mu }_{\perp }{\Vert }^2\right|\le c_s(\alpha ){\sigma }^2(n-k{)}^{-1/2}\right\}\ge 1-\alpha ,$$

$$r_{\perp }^2=c_2{\tilde{r}}_{\perp }^2=c_2\frac{n-k}{n}{\sigma }^2\left\{\hat{L}+c_s(\alpha /2)(n-k{)}^{-1/2}\right\},$$

$$\underset{(n-k)\to \infty }{\mathrm{liminf}}\underset{\beta \in {\mathbb{R}}^p}{\inf }\mathbb{P}\left\{\frac{\Vert {\mu }_{\perp }-{\hat{\mu }}_{\perp }{\Vert }^2}{n{r}_{\perp }^2}\le 1/c_2\right\}\ge 1-\alpha /2.$$

$$\mathbb{E}|\hat{C}{|}^2=$$

$$\mathbb{E}|\hat{C}{|}^2\le 4{\sigma }^2\max \left\{c_1\frac{{\chi }_{k,\alpha /2}^2}{n},c_2\frac{n-k}{n}\left(1+c_s(\alpha /2)(n-k{)}^{-1/2}\right)\right\}$$

$$r_A^2{\asymp }_{p}k/n,r_{\perp }^2=O_p\left(\frac{\sqrt{n-k}}{n}+\frac{\Vert {\mu }_{\perp }{\Vert }^2}{n}\right).$$

$$|\hat{C}|=2(r_A\vee r_{\perp }){\asymp }_p{n}^{-1/4}.$$

$$c_{\ast }\le \frac{\Vert X_{A}v{\Vert }^2}{n\Vert v{\Vert }^2}\le c^{\ast },\text{for all }A\text{ with }|A|=s^{\ast }\text{ and all nonzero }v\in {\mathbb{R}}^{s^{\ast }}.$$

$$|\hat{C}|=O_p\left\{(n^{-1/4}+\Vert {\beta }_{A^c}\Vert )\vee \sqrt{k/n}\right\}$$

$$\hat{\beta }=\hat{\beta }(y^{\prime },X^{\prime };\lambda ):=\underset{\beta \in {\mathbb{R}}^p}{\mathrm{argmin}}\left[\frac{1}{2n}\Vert y^{\prime }-X^{\prime }\beta {\Vert }^2+\lambda \Vert \beta {\Vert }_1\right],$$

$$A=\{j:{\hat{\beta }}_j\ne 0\},$$

$$|\hat{C}|=O_p\left(n^{-1/4}+\sqrt{s\log p/n}\right).$$

$$|\hat{C}|=O_p\left(n^{-1/4}\vee \sqrt{s/n}\right).$$

$$\rho (t;\lambda ,\gamma )={\int }_0^{|t|}{\left(1-\frac{u}{\gamma \lambda }\right)}_{+}du=\{\begin{array}{ll}|t|-t^2/(2\gamma \lambda )& \text{if}|t|\le \gamma \lambda \\ \gamma \lambda /2& \text{if}|t|>\gamma \lambda \end{array},$$

$${\hat{\beta }}_{\lambda ,\gamma }^{\text{mcp}}={\hat{\beta }}_{\lambda ,\gamma }^{\text{mcp}}(y^{\prime },X^{\prime }):=\underset{\beta \in {\mathbb{R}}^p}{\mathrm{argmin}}\left[\frac{1}{2n}\Vert y^{\prime }-X^{\prime }\beta {\Vert }^2+\lambda \sum _{j=1}^p\rho (|{\beta }_j|;\lambda ,\gamma )\right].$$

$$|\hat{C}|=O_p\left(n^{-1/4}\vee \sqrt{s/n}\right).$$

Full text

Keywords: adaptive confidence set, high-dimensional inference, sparse linear regression, Stein estimate.

1 Introduction

Consider high-dimensional linear regression

[TABLE]

where $y\in\mathbb{R}^{n}$, $X=[X_{1}|\cdots|X_{p}]\in\mathbb{R}^{n\times p}$, $\beta\in\mathbb{R}^{p}$, $\varepsilon\sim\mathcal{N}_{n}(0,\sigma^{2}\mathbf{I}_{n})$ and $p>n$. While there is a rich body of research on parameter estimation under this model concerning signal sparsity (e.g. Bickel et al. (2009); Zhang and Huang (2008); Negahban et al. (2012)), how to construct confidence sets remains elusive. In this work, we focus on confidence sets for the mean $\mu=X\beta$ with the following two properties: First, the confidence set $\widehat{C}$ is (asymptotically) honest over all possible parameters. That is, for a given confidence level $1-\alpha$,

[TABLE]

where $\mathbb{P}_{\beta}$ is taken with respect to the distribution of $y\sim\mathcal{N}_{n}(X\beta,\sigma^{2}\mathbf{I}_{n})$, regarding $X$ as fixed. Second, the diameter of $\widehat{C}$ is able to adapt to the sparsity and the strength of $\beta$. In practical applications, sparsity assumptions are very hard to verify, and for many data sets they are at most a good approximation. The first property guarantees that our confidence sets reach the nominal coverage without imposing any sparsity assumption, while the second property allows us to leverage sparse estimation when $\beta$ is indeed sparse.

Our problem is related to the construction of confidence sets in nonparametric regression, for which a line of work has laid down important theoretic foundations and provided methods of construction (Li 1989; Beran and Dümbgen 1998; Hoffman and Lepski 2002; Juditsky and Lambert-Lacroix 2003; Baraud 2004; Genovese and Wasserman 2005; Robins and van der Vaart 2006; Cai and Low 2006; Bull and Nickl 2013). Despite such notable advances, lack of numerical support casts doubt on the merit of borrowing these nonparametric regression methods directly for sparse regression. Taking the adaptive method based on sample splitting in Robins and van der Vaart (2006) as an example, an honest confidence set for $\mu$ can be constructed as $\widehat{C}_{a}=\{\mu\in\mathbb{R}^{n}:n^{-1/2}\|\mu-X\hat{\beta}\|\leq r_{n}\}$, where $X\hat{\beta}$ is an initial estimate independent of $y$, and its (normalized) diameter $|\widehat{C}_{a}|\mathop{:=}2r_{n}=O_{p}(n^{-1/4}+n^{-1/2}\|X\hat{\beta}-X\beta\|)$. A common choice for $\hat{\beta}$ under model (1.1) for $p>n$ is a sparse estimator, such as the lasso (Tibshirani 1996) or $\ell_{0}$-penalized least-squares estimator. With high probability, the prediction loss of the lasso estimator typically satisfies

[TABLE]

for some $c>0$, uniformly for all $\beta\in\mathscr{B}(s):=\{v\in\mathbb{R}^{p}:\|v\|_{0}\leq s\}$; see for example Bickel et al. (2009). Under this choice, the diameter $|\widehat{C}_{a}|$ is of the order

[TABLE]

for all $\beta\in\mathscr{B}(s)$. For a precise statement, see Theorem 6 below. This method has nice theoretical properties when $s=o(n/\log p)$. But even for moderately sparse signals with $s/n\to\delta\in(0,1)$, the bound on the right side of (1.4) approaches $\infty$ as $p>n\to\infty$ and thus offers little insight into the performance of the confidence set. The upper bound (1.3) also critically depends on the regularization parameter used for the initial estimate $\hat{\beta}$. In fact, our numerical results show that, for finite samples with $(s,n,p)=(10,200,800)$, this confidence set can be worse than a naive $\chi^{2}$ region $\{\mu:\|y-\mu\|^{2}\leq\sigma^{2}\chi^{2}_{n,\alpha}\}$, where $\chi^{2}_{n,\alpha}$ denotes the $1-\alpha$ quantile of the $\chi^{2}$ distribution with $n$ degrees of freedom. A similar issue occurs in the related but different problem of constructing confidence sets for $\beta$. Nickl and van de Geer (2013) have shown that one can construct a confidence set for $\beta$ that is honest over $\mathscr{B}(k_{1})$ for $k_{1}=o(n/\log p)$, and for any $s\leq k_{1}$, the diameter is on the same order as that in (1.4) for any $\beta\in\mathscr{B}(s)$. Compared to the unrestricted honesty in (1.2) over the entire space $\mathbb{R}^{p}$, the restriction on the honesty region to $\mathscr{B}(k_{1})$ also reflects the challenge faced in the construction of confidence sets when $p>n$. Recently, Ewald and Schneider (2018) provide an exact formula to compute a lower bound of the coverage rate of a confidence set centered at the lasso, over the entire parameter space for any significance level $\alpha\in(0,1)$, and vice versa; however, low dimension $(p<n)$ is a vital condition in their proof, making it impossible to generalize their idea to the high-dimensional problem that we are studying.

The construction of confidence sets is fundamentally different from the problem of inferring error bounds for a sparse estimator (Nickl and van de Geer 2013). It is seen from (1.4) that no matter how sparse the true $\beta$ is, the diameter of $\widehat{C}_{a}$ cannot converge at a rate faster than $n^{-1/4}$. Indeed, results in Li (1989) imply that, for the linear model (1.1) with $p\geq n$, the diameter of an honest confidence set for $\mu$, in the sense of (1.2), cannot adapt at any rate $o(n^{-1/4})$. This is in sharp contrast to error bounds for a sparse estimator, such as that in (1.3), which can decay at a much faster rate when $\beta$ is sufficiently sparse. It is not desired to construct confidence sets directly from error bounds like (1.3) even we only require honesty for $\beta\in\mathscr{B}(k_{1})$ with a given $k_{1}=o(n/\log p)$, because its diameter, on the order of $\sqrt{k_{1}\log p/n}$, cannot adapt to any sparser $\beta\in\mathscr{B}(s)$ for $s<k_{1}$.

Motivated by these challenges, we propose a new two-step method to construct a confidence set for $\mu=X\beta$, allowing the dimension $p\gg n$ in (1.1). The basic idea of our method is to estimate the radius of the confidence set separately for strong and weak signals defined by the magnitude of $|\beta_{j}|$. Using a sparse estimate, such as the lasso, one can recover the set $A$ of large $|\beta_{j}|$ accurately and expect a small radius for a confidence ball for $\mu_{A}$, the projection of $\mu$ onto the subspace spanned by $X_{j},j\in A$. By construction, $(\mu-\mu_{A})$ is composed of weak signals. Thus, in the second step, we shrink our estimate of this part towards zero by Stein’s method and construct a confidence set with Stein’s unbiased risk estimate (Stein 1981). Combining the inferential advantages of sparse estimators and Stein estimators, our method overcomes many of the aforementioned difficulties. First, our confidence set is honest for all $\beta\in\mathbb{R}^{p}$, and its diameter is well under control for all possible values of $\beta$ including the dense case. Second, by using elastic radii our confidence set, an ellipsoid in general, can adapt to signal strength and sparsity. The radius for strong signals adapts to the sparsity of the underlying model via sparse estimation or model selection, while the radius for weak signals adapts according to the degree of shrinkage of the Stein estimate. Without any signal strength assumption, the diameter of our confidence set is $O_{p}(n^{-1/4}+\sqrt{s\log p/n})$, the same as (1.4), for $\beta\in\mathscr{B}(s)$. It may further reduce to $O_{p}(n^{-1/4}+\sqrt{s/n})$ under an assumption on the separability between the strong and the weak signals, which shrinks to the optimal rate $n^{-1/4}$ when the signal sparsity $s=O(\sqrt{n})$, as opposed to $s=O(\sqrt{n}/\log p)$ in (1.4). Third, we provide a data-driven selection of the set $A$ from multiple candidates, which protects our method from a bad choice and thus makes it very robust. We demonstrate with extensive numerical results that our method can construct much smaller confidence sets than other competing methods, including the adaptive method (Robins and van der Vaart 2006) discussed above and oracle approaches making use of the true sparsity of $\beta$ (the oracle). These results highlight the practical usefulness of our method.

Note that the construction of confidence sets for $\mu=X\beta$ is different in nature from the construction of confidence intervals for an individual $\beta_{j}$ or a low-dimensional projection of $\beta$. For the latter, the optimal rate of an interval length can be $n^{-1/2}$ when $\beta$ is sufficiently sparse (Schneider 2016; Cai and Guo 2017), such as the intervals constructed by de-biased lasso methods (Zhang and Zhang 2014; van de Geer et al. 2014; Javanmard and Montanari 2014). Although simultaneous inference methods have been proposed based on bootstrapping de-biased lasso estimates (Zhang and Cheng 2017; Dezeure et al. 2017), these methods are shown to achieve the desired coverage only for extremely sparse $\beta$ such that $\|\beta\|_{0}=o(\sqrt{n/(\log p)^{3}})$, which severely limits their practical application.

The remainder of this paper is organized as follows: Section 2 develops our two-step Stein method in details, including its theoretical properties and algorithmic implementation. To demonstrate the advantage of our method, we develop in Section 3 a few competing methods making use of the lasso prediction or the oracle of the true sparsity. Extensive numerical comparisons are provided in Section 4 to show the superior performance of our two-step Stein method, relative to the competitors, in a variety of sparsity settings, including when $\beta$ is quite dense. The paper is concluded in Section 5 with further discussions. Proofs of all theoretical results are deferred to the Supplementary Material.

Throughout the paper, we always assume model (1.1) with $\varepsilon\sim\mathcal{N}_{n}(0,\sigma^{2}\mathbf{I}_{n})$ unless otherwise noted. We denote by $\mathbb{P}_{\beta}$ the distribution of $[y\mid X]$ and $\mathbb{E}_{\beta}$ the corresponding expectations, where the subscript $\beta$ may be dropped when its meaning is clear from the context. Denote by $[p]$ the index set $\{1,\ldots,p\}$ and by $|A|$ the size of a set $A\subseteq[p]$. Write $a_{n}=\Omega(b_{n})$ if $b_{n}=O(a_{n})$ and $a_{n}\asymp b_{n}$ if $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$. We use $\Omega_{p}(.)$ and $\asymp_{p}$ if the above statements hold in probability. For a vector $v=(v_{j})_{1:m}$, let $v_{A}=(v_{j})_{j\in A}$ be the restriction of $v$ to the components in $A$. For a matrix $M=[M_{1}\mid\ldots\mid M_{m}]$, where $M_{j}$ is the $j$th column, denote by $M_{A}=(M_{j})_{j\in A}$ the submatrix consisting of columns in $A$. For $a,b\in\mathbb{R}^{n}$, $\langle a,b\rangle\mathop{:=}a^{\mathsf{T}}b$ is the inner product. Define $a\vee b\mathop{:=}\max\{a,b\}$ and $a\wedge b\mathop{:=}\min\{a,b\}$ for $a,b\in\mathbb{R}$.

2 Two-step Stein method

Dividing $\beta$ into strong and weak signals, our method constructs a confidence set $\widehat{C}(y)$ with an ellipsoid shape for $X\beta$ that is honest as defined in (1.2). Note that under a high-dimensional asymptotic framework, all variables $X=X(n)$, $y=y(n)$, $\beta=\beta(n)$ and $s=s_{n}$ depend on $n$ as $p=p_{n}\gg n\to\infty$, while $X(n)$ is regarded as a fixed design matrix for each $n$. We often suppress the dependence on $n$ to simplify the notation.

2.1 Preliminaries on Stein estimation

We will use a simplified Stein estimate (Li 1989) to construct the confidence set for weak signals. For a linear estimate $\tilde{\mu}=T_{n}y$, where $y\sim\mathcal{N}_{n}(\mu,\sigma^{2}\mathbf{I}_{n})$ and $T_{n}\in\mathbb{R}^{n\times n}$, let $R_{n}=\mathbf{I}_{n}-T_{n},$ and define

[TABLE]

where $\hat{\mu}(y;\tilde{\mu})$ is the Stein estimate associated with the initial estimate $\tilde{\mu}$ and $\sigma^{2}\hat{L}(y;\tilde{\mu})$ is the Stein unbiased risk estimate (SURE). Li (1989) proved the uniform consistency of $\hat{L}$.

Lemma 1** (Theorem 3.1 in Li (1989)).**

Assume that $y\sim\mathcal{N}_{n}(\mu,\sigma^{2}\mathbf{I}_{n})$. For any $\alpha\in(0,1)$, there exists a constant $c_{s}(\alpha)>0$ such that

[TABLE]

where $\hat{\mu}$ and $\hat{L}$ are defined in (2.1) and (2.2).

2.2 Method of construction

Now, consider the linear model (1.1) and let $\mu=X\beta$. Given a pre-constructed candidate set $A=A_{n}\subseteq[p]$, independent of $(X,y)$, define

[TABLE]

where $P_{A}$ is the orthogonal projection from $\mathbb{R}^{n}$ onto $\operatorname{span}(X_{A})$ and $P_{A}^{\perp}$ is the projection to the orthogonal complement. A good candidate set $A$ is supposed to include all strong signals, say $A=\{j:|\beta_{j}|>\tau\}$. With such a choice, $\|\mu_{\perp}\|$ will be small. Typically, we split our data set into two halves, $(X,y)$ and $(X^{\prime},y^{\prime})$, and apply a model selection method on $(X^{\prime},y^{\prime})$ to construct the set $A$. See Section 2.3 for more detailed discussion.

We estimate $\mu_{A}$ and $\mu_{\perp}$, respectively, by $\hat{\mu}_{A}$ and $\hat{\mu}_{\perp}$, compute radii $r_{A}$ and $r_{\perp}$, and construct a $(1-\alpha)$ confidence set $\widehat{C}$ for $\mu$ in the form of

[TABLE]

Note that $\widehat{C}$ is an ellipsoid in $\mathbb{R}^{n}$, where $r_{A}=r_{A}(\alpha)$ and $r_{\perp}=r_{\perp}(\alpha)$ correspond to the major and minor axes, respectively. Our method consists of a projection and a shrinkage step:

Step 1: Projection. Let $\hat{\mu}_{A}=P_{A}y$ and $k=\operatorname{rank}(X_{A})\leq|A|$. Since $A$ is independent of $(y,X)$, we have

[TABLE]

Thus, we choose

[TABLE]

where $\chi^{2}_{k,\alpha/2}$ is the $(1-\alpha/2)$ quantile of the $\chi^{2}_{k}$ distribution and $c_{1}>1$ is a constant, so that

[TABLE]

Step 2: Shrinkage. Let $y_{\perp}=P_{A}^{\perp}y$. As mentioned above, under a good choice of $A$ that contains strong signals, $\|\mu_{\perp}\|$ is expected to be small. Therefore, we shrink $y_{\perp}$ towards zero via Stein estimation to construct $\hat{\mu}_{\perp}$. Note that $y_{\perp}$ is in an $(n-k)$-dimensional subspace of $\mathbb{R}^{n}$. Letting $\tilde{\mu}=0$ and $R_{n}=P_{A}^{\perp}$ in (2.1) and (2.2), we obtain

[TABLE]

where the shrinkage factor

[TABLE]

It then follows from Lemma 1 that

[TABLE]

for any sequence of $A=A_{n}$ as long as $(n-k)\to\infty$. Therefore, if we choose

[TABLE]

where $c_{2}>1$ is a constant, we have

[TABLE]

In practical implementation, we estimate the constant $c_{s}(\alpha)$ in (2.11) by simulation, which will be discussed in Section 2.5.

If $1/c_{1}+1/c_{2}=1$, confidence set (2.4) made up from (2.7) and (2.13) is honest and the expectation of its (normalized) diameter $|\widehat{C}|\mathop{:=}2(r_{A}\vee r_{\perp})$ can be calculated explicitly for all $\beta\in\mathbb{R}^{p}$:

Theorem 1**.**

Assume $1/c_{1}+1/c_{2}=1$, $A$ is independent of $(y,X)$ with $\operatorname{rank}(X_{A})=k$, and $(n-k)\to\infty$ as $n\to\infty$. Then the confidence set $\widehat{C}$ (2.4) constructed by the two-step Stein method is honest in the sense of (1.2). Furthermore, the squared diameter of $\widehat{C}$ has expectation

[TABLE]

where $\chi^{2}_{n-k}(\rho)$ follows a noncentral $\chi^{2}$ distribution with $n-k$ degrees of freedom and non-centrality parameter $\rho=\|\mu_{\perp}\|^{2}/\sigma^{2}$.

In the above result, we did not impose any assumptions on $A$ except $(n-k)\to\infty$, which allows many choices of $A$. Our confidence set $\widehat{C}$ is honest as in (1.2) and its diameter is under control for all $\beta\in\mathbb{R}^{p}$. Since $\mathbb{E}[1/\chi^{2}_{n-k}(\rho)]>0$, a uniform but very loose upper bound

[TABLE]

holds for all $\beta\in\mathbb{R}^{p}$. In particular, when $\beta$ is dense, the diameter will be comparable to that of the naive $\chi^{2}$ region. As corroborated with the numerical results in Section 4.4, this protects our method from inferior performance when sparsity assumptions are violated, making it robust to different data sets. Next, we will show that our confidence set is adaptive: When $\beta$ is indeed sparse with separable strong and weak signals, the radii $r_{A}$ and $r_{\perp}$ will adapt to the optimal rate with a proper choice of $A$ that contains strong signals.

2.3 Adaptation of the diameter

To simplify our analysis, we set $c_{1}=c_{2}=2$ in this section so that they can be ignored when calculating the convergence rates of $r_{A}$ and $r_{\perp}$. These rates do not change as long as $c_{1}$ and $c_{2}$ stay as constants when $n\to\infty$. Lemma 2 specifies conditions for the diameter of $\widehat{C}$ to converge at the optimal rate $n^{-1/4}$.

Lemma 2**.**

Suppose that $k=\operatorname{rank}(X_{A})$ and $\lVert\mu_{\perp}\rVert=o(\sqrt{n-k})$. Then

[TABLE]

Therefore, if $k=O(\sqrt{n})$ and $\|\mu_{\perp}\|=O(n^{1/4})$, then the diameter of $\widehat{C}$

[TABLE]

The $\ell_{2}$ norm of the weak signals $\lVert\mu_{\perp}\rVert$ can be bounded by $\|\beta_{A^{c}}\|$ under the sparse Riesz condition on $X$ and a sparsity assumption on $\beta$. A design matrix $X$ satisfies the sparse Riesz condition (Zhang and Huang 2008) with rank $s^{*}$ and spectrum bounds $0<c_{*}<c^{*}<\infty$, denoted by $\textup{SRC}(s^{*},c_{*},c^{*})$, if

[TABLE]

Under our asymptotic framework, $s^{*}$, $c^{*}$ and $c_{*}$ are allowed to depend on $n$.

Theorem 2**.**

Suppose $X$ satisfies $\textup{SRC}(s^{*},c_{*},c^{*})$ with $s^{*}\geq|\operatorname{supp}(\beta)\cap A^{c}|$, and let $k=\operatorname{rank}(X_{A})$. If $\limsup_{n}c^{*}<\infty$, $k=o(n)$ and $\lVert\beta_{A^{c}}\rVert=o(1)$, then

[TABLE]

for the two-step Stein method. In particular, $|\widehat{C}|\asymp_{p}n^{-1/4}$ if $k=O(\sqrt{n})$ and $\lVert\beta_{A^{c}}\rVert=O(n^{-1/4})$.

Remark 1*.*

Let us take a closer look at the conditions in this theorem for $|\widehat{C}|\asymp_{p}n^{-1/4}$. Suppose that $\beta$ has $O(\sqrt{n})$ strong coefficients that can be reliably detected by a model selection method, while all other signals are weak such that $\lVert\beta_{A^{c}}\rVert=O(n^{-1/4})$. Then we can have $k\leq|A|=O(\sqrt{n})$ with high probability. This shows that the sparsity $s=\|\beta\|_{0}$ is allowed to be $O(\sqrt{n})$. The only additional constraint on $s$ comes from the assumption $\textup{SRC}(s^{*},c_{*},c^{*})$ with $s^{*}\geq s$, which holds for Gaussian designs if $s\log p=o(n)$ (Zhang and Huang 2008). Compared to (1.4) which requires $s\log p=O(\sqrt{n})$, we have relaxed the sparsity assumption on $\beta$ to attain the optimal rate $n^{-1/4}$ by imposing a mild condition on the decay rate of the weak signals $\|\beta_{A^{c}}\|$.

Now we discuss a few methods to find $A$ so that our confidence sets can adapt to the sparsity and signal strength of $\beta$. We split the whole data set into $(X,y)$ and $(X^{\prime},y^{\prime})$, with respective sample sizes $n$ and $n^{\prime}$, so that they are independent. Henceforth, we assume an even partition with $n^{\prime}=n$, which simplifies the notation and is commonly used in practice, unless otherwise noted. The first method is to apply lasso on $(X^{\prime},y^{\prime})$:

[TABLE]

where $\lambda$ is a tuning parameter. Then choose

[TABLE]

that is, we define strong signals by the support of the lasso. This choice of $A$ is justified by the following corollary. Let $A_{0}=\operatorname{supp}(\beta)$ and $S_{0}=\{j\in A_{0}:|\beta_{j}|\geq K\sqrt{s\log p/n}\}$ for a sufficiently large $K$.

Corollary 3**.**

Suppose that $X$ and $X^{\prime}$ satisfy $\textup{SRC}(s^{*},c_{*},c^{*})$, where $0<c_{*}<c^{*}$ are constants. Let the confidence set $\widehat{C}$ (2.4) be constructed by the two-step Stein method with $A$ chosen by (2.18) and $\lambda=c_{0}\sigma\sqrt{c^{*}\log p/n},c_{0}>2\sqrt{2}$. Assume $s\leq(s^{*}-1)/(2+4c^{*}/c_{*})$ and $s\log p=o(n)$. Then for any $\beta\in\mathscr{B}(s)$ we have

[TABLE]

If in addition $\lVert\beta_{A_{0}\setminus S_{0}}\rVert=O(n^{-1/4})$, then

[TABLE]

The rate of $|\widehat{C}|$ in (2.19) does not depend on any assumption on signal strength, and it is identical to (1.4). However, our method can achieve a faster rate (2.20) if $\lVert\beta_{A_{0}\setminus S_{0}}\rVert=O(n^{-1/4})$. Together with the definition of $S_{0}$, this essentially imposes a separability assumption between the strong and the weak signals when $s\log p\gg\sqrt{n}$.

To weaken the beta-min condition on strong signals in $S_{0}$, we may apply a better model selection method to define $A$, such as using the minimax concave penalty (MCP) (Zhang 2010):

[TABLE]

for $\gamma>1$. Accordingly, a regularized least-squares estimate is defined by

[TABLE]

Suppose we choose $A=\operatorname{supp}(\hat{\beta}^{\textup{mcp}}_{\lambda,\gamma})$ in our two-step Stein method. The model selection consistency of $\hat{\beta}^{\textup{mcp}}_{\lambda,\gamma}$ makes it possible for $|\widehat{C}|$ to adapt at the rate (2.20) under the same SRC assumption but a weaker beta-min condition than Corollary 3.

Corollary 4**.**

Suppose that $X$ and $X^{\prime}$ satisfy $\textup{SRC}(s^{*},c_{*},c^{*})$, where $0<c_{*}<c^{*}$ are constants, $s^{*}\geq(c^{*}/c_{*}+1/2)s$, and $s\log p=o(n)$. Choose a sequence of $(\lambda_{n},\gamma_{n})$ satisfying $\lambda_{n}\gg\sqrt{\log p/n}$ and $\gamma_{n}\geq c_{*}^{-1}\sqrt{4+c_{*}/c^{*}}$. If $\beta\in\mathscr{B}(s)$ and $\inf_{A_{0}}|\beta_{j}|\geq(\gamma_{n}+1)\lambda_{n}$, then $\mathbb{P}\{\operatorname{supp}(\hat{\beta}^{\textup{mcp}}_{\lambda_{n},\gamma_{n}})=A_{0}\}\to 1$, and consequently the $\widehat{C}$ constructed by the two-step Stein method with $A=\operatorname{supp}(\hat{\beta}^{\textup{mcp}}_{\lambda_{n},\gamma_{n}})$ has diameter

[TABLE]

Remark 2*.*

Compared to (1.4) for confidence sets centering at a sparse estimator, the diameter of our method in (2.20) and (2.23) converges faster by a factor of $(\log p)^{1/2}$ when $s=\Omega(\sqrt{n})$. Accordingly, our method achieves the optimal rate when $s=O(\sqrt{n})$ instead of $s=O(\sqrt{n}/\log p)$ as for (1.4). Under a high-dimensional setting with $p\gg n$, say $p=\exp(n^{a})$ for $a\in(0,1/2)$, this improvement in rate can be very substantial, which is supported by our numerical results. The faster rate of our method is made possible by its adaption to both signal strength and sparsity, while the rate of (1.4) is obtained by adaption to sparsity only (cf. Theorem 6). We emphasize that our method achieves the adaptive rates in the above results, while being uniformly honest over the entire $\mathbb{R}^{p}$ (Theorem 1). One could construct a confidence set with diameter $O_{p}(\sqrt{s/n})$ using only the covariates selected by a consistent model selection method, which would be faster than the rate (2.23). However, such a confidence set is not honest over $\mathbb{R}^{p}$, because it cannot reach the nominal coverage rate for those $\beta$ that do not satisfy the required beta-min condition for model selection consistency. Our method overcomes this difficulty with the shrinkage step, based on the uniform consistency of the SURE (Lemma 1).

Remark 3*.*

For an uneven partition of the whole data set, the conclusions of Corollaries 3 and 4 still hold as long as both $n^{\prime}\asymp n\to\infty$. However, it is a common and reasonable choice to have $n=n^{\prime}$, since $(X^{\prime},y^{\prime})$ and $(X,y)$ can be swapped to construct a confidence set for $X^{\prime}\beta$, making full use of the whole data set.

2.4 Multiple candidate sets

It is common to have multiple choices for the candidate set $A$ in our two-step Stein method. Let

[TABLE]

be a collection of candidate sets. We can apply the two-step Stein method to construct $M=M_{n}$ confidence sets for $\mu$, denoted by $\widehat{C}_{m}$, and then choose an optimal set $\widehat{C}_{m^{*}}$ by certain criterion such as minimizing the volume or the diameter. Furthermore, the cardinality of $\mathcal{H}$ may be unbounded as $n$ increases, i.e., $M_{n}\to\infty$. In what follows, we show that under mild conditions, (2.7) and (2.13) hold uniformly for all $A\in\mathcal{H}$ after modifying $r_{A}$ and $r_{\perp}$ accordingly, which implies $\widehat{C}_{m^{*}}$ is asymptotically honest.

Put $k=\operatorname{rank}(X_{A})$ for $A\in\mathcal{H}$ and $k_{\max}=\max_{A\in\mathcal{H}}k$. Intuitively, the cardinality of $\mathcal{H}$ (i.e. $M$) and the maximum size of $A$ in $\mathcal{H}$ (i.e. $k_{\max}$) determine the radii and the coverage probability of $\widehat{C}_{m}$.

For strong signals, we apply the following concentration inequality to show (2.7) holds uniformly:

Lemma 3**.**

Suppose $\chi^{2}_{n}$ follows a $\chi^{2}$ distribution with $n$ degrees of freedom. Then for any $\delta>0$,

[TABLE]

This lemma with a union bound implies

[TABLE]

Then choosing

[TABLE]

as the radius for strong signals, we have

[TABLE]

For weak signals, we establish (2.13) uniformly over $\mathcal{H}$ via the following result:

Lemma 4**.**

Suppose all components of $\varepsilon$ in (1.1), $\varepsilon_{i},i=1,\ldots,n$, have mean [math], common second, forth and sixth moments and their eighth moments are bounded by some constant $d$. For any $\delta>0$ there exists a positive number $D$ depending on $d$ such that

[TABLE]

The proof of Lemma 4 mainly follows the ideas in Li (1985). In our model with $\varepsilon\sim\mathcal{N}_{n}(0,\sigma^{2}\mathbf{I}_{n})$, the first term on the right hand side of (2.26) simplifies to

[TABLE]

via Lemma 3. Assume that the cardinality of $\mathcal{H}$ and the maximum size of $A\in\mathcal{H}$ satisfy $M\ll(n-k_{\max})^{2}$. To achieve the desired coverage for weak signals, it is sufficient to pick $\delta$ such that $\delta^{2}=\Omega(\log M)$ and $\delta^{4}=\Omega(M)$. Therefore, we can set

[TABLE]

for some constant $c_{m}(\alpha/2)>0$, and the corresponding radius

[TABLE]

for any $A\in\mathcal{H}$, so that the upper bound in (2.26) is $\leq\alpha/2$. Now we generalize Theorem 1 to establish asymptotic honesty uniformly over $\mathcal{H}$:

Theorem 5**.**

Given $\mathcal{H}$, construct confidence sets $\widehat{C}_{m},m=1,\ldots,M$, with $r_{A}$ and $r_{\perp}$ as in (2.25) and (2.27), respectively, for $A=A_{m}$. Suppose $\lim_{n\to\infty}M/(n-k_{\max})^{2}=0$, $1/c_{1}+1/c_{2}=1$, and each $A_{m}$ is independent of $(X,y)$. Then the confidence sets $\widehat{C}_{m}$ are uniformly honest over $\mathcal{H}$, i.e.,

[TABLE]

Consequently, $\widehat{C}_{m^{*}}$ chosen by any criterion is asymptotically honest.

Remark 4*.*

The increment of $r_{A}^{2}$ in (2.25), $2\sqrt{k\log(4M/\alpha)}/n$, reflects the cost for achieving uniform honesty over $\mathcal{H}$. But this factor will not cause a slower rate for $r_{A}$ if $\log M=O_{p}(k)$, where the $k$ here is the size of the selected candidate set $A_{m^{*}}$. Compared with (2.12), the factor $M^{1/4}/\sqrt{n-k}$ in (2.27), also the cost for uniform honesty, will in general lead to slower convergence of $r_{\perp}$. However, this is a worthwhile price to protect our method from an improper candidate set $A$ that does not satisfy the assumptions in Theorem 2. For example, if the candidate set $A$ misses some strong signals, we may end up with $\hat{L}\asymp_{p}1$ and the radius of weak signals $r_{\perp}$ will not converge to [math] at all. Such bad choices of $A$ will be excluded if $\widehat{C}_{m^{*}}$ is chosen by minimizing its volume over $\mathcal{H}$. In this sense, our method provides a data-driven selection of an optimal candidate set.

To construct $\mathcal{H}$, we threshold the lasso $\hat{\beta}$ in (2.17) calculated from $(X^{\prime},y^{\prime})$ to obtain

[TABLE]

for a sequence of threshold values $\tau_{m}=a_{m}\lambda$, e.g. $a_{m}\in[0,4]$. It is possible for two different $\tau_{m}$ to define the same $A$, which will be counted once in $\mathcal{H}$. By setting $\tau_{m}=0$ for some $m$, $A=\operatorname{supp}(\hat{\beta})$ will be included in $\mathcal{H}$, though it may not be selected as the optimal $\widehat{C}_{m^{*}}$. In the proof of Corollary 3, we have shown $\lVert\hat{\beta}\rVert_{0}=O_{p}(\sqrt{n})$, and therefore both $M$ and $k_{\max}$ are $O_{p}(\sqrt{n})$, which means $M\ll(n-k_{\max})^{2}$ with high probability. As a result, we can guarantee uniform honesty over all $\widehat{C}_{m}$. Other choices of $\mathcal{H}$ are possible, such as stepwise variable selection with BIC. It is possible that $A=\varnothing$ for a large value of $\tau_{m}$. In this special case, $r_{A}=0$, so the confidence set reduces to a ball, i.e., $\left\{\mu\in\mathbb{R}^{n}:{\|\mu-\hat{\mu}_{\perp}\|^{2}}\leq nr_{\perp}^{2}\right\}$.

2.5 Algorithm and implementation

We implement our method with a sequence of candidate sets $A_{m}$ defined by (2.28). Given the data set, $\sigma^{2}$, $\lambda$ in (2.17) and threshold values $\{a_{m}\lambda\}_{1\leq m\leq M}$, this section describes some technique details in our algorithm to construct the confidence set (2.4) by the two-step Stein method.

Data splitting. We split the original data set into $(X^{\prime},y^{\prime})$ and $(X,y)$. Apply lasso on $(X^{\prime},y^{\prime})$ to get $\hat{\beta}$ in (2.17) with the tuning parameter $\lambda$. Threshold $\hat{\beta}$ by $\tau_{m}=a_{m}\lambda$ for $m=1,\ldots,M$ in (2.28) to define candidate sets $A_{m}$. Note that $A_{m}$, $m=1,\ldots,M$, are independent of $(X,y)$.

Choice of $c_{1}$ and $c_{2}$. When $A\neq\varnothing$, we consider two criteria to choose the constants $c_{1}$ in (2.6) and $c_{2}$ in (2.12). The first criterion is to minimize the log-volume of $\widehat{C}$, namely,

[TABLE]

up to an additive constant, which becomes a constrained optimization problem

[TABLE]

where $\tilde{r}_{A}$ and $\tilde{r}_{\perp}$ are defined in (2.6) and (2.12) and $E>2$ is a pre-determined upper bound. It is easy to obtain the solution

[TABLE]

For all numerical results in this paper, we use $E=10$. Without the constraint $c_{1},c_{2}\leq E$, the minimizer would be $(c_{1},c_{2})=(n/k,n/(n-k))$ so that under the conditions of Corollary 3, $r_{A}=\sqrt{n/k}\tilde{r}_{A}\asymp_{p}1$ and thus the diameter $|\widehat{C}|$ would not converge to [math]. Therefore, a finite upper bound $E$ must be imposed.

The second criterion is to minimize the diameter $|\widehat{C}|$

[TABLE]

which yields the solution

[TABLE]

As a result, we have $r_{A}=r_{\perp}=(\tilde{r}_{A}^{2}+\tilde{r}_{\perp}^{2})^{1/2}$ and the confidence set reduces to a ball.

Computation of $c_{s}(\alpha)$. For any candidate set $A$, the radius $r_{\perp}$ (2.12) depends on the constant $c_{s}(\alpha)$, which is essentially the quantile of the deviation between $\sigma^{2}\hat{L}$ and the loss of the Stein estimator $\hat{\mu}_{\perp}$. We use the following simulation procedure to estimate $c_{s}(\alpha)$: First draw $\check{Y}_{j}\sim\mathcal{N}_{n}(0,\sigma^{2}\mathbf{I}_{n})$ for $j=1,2,\ldots,N$. For each $j$, compute

[TABLE]

Then the $(1-\alpha)$ quantile of the empirical distribution of

[TABLE]

is a consistent estimator of $c_{s}(\alpha)$ as long as $\|\mu_{\perp}\|=o(\sqrt{n})$, which is the case under the assumptions of Corollary 3. Expression (2.34) can be written as a function of a $\chi^{2}_{n}$ random variable, which simplifies its simulation.

Clearly, the estimate of $c_{s}(\alpha)$ does not depend on $A$ and is used for any candidate set $A\in\mathcal{H}$ in our implementation. Moreover, we find the multiple set adjustments on the radii, i.e., the factors of $(\log M)^{1/2}$ and $M^{1/4}$, are usually negligible given a reasonable sample size, say $n\geq 100$. Therefore, we simply use the radii $r_{A}$ and $r_{\perp}$ in (2.6) and (2.12) for each $A\in\mathcal{H}$.

Algorithm 1 summarizes the two-step Stein method with multiple candidate sets $A_{m}$.

Remark 5*.*

In the calculation of $r_{\perp}$ and $c_{s}(\alpha)$, we use truncated SURE for $\hat{L}=(1-B)_{+}$ in (2.9) and similarly for $\check{L}_{j}$ in (2.33). Such a truncated rule has been used for the James-Stein estimator (Efron and Morris 1973) and does not affect the asymptotic validity of our method.

For all numerical results in this paper, we assume the noise variance $\sigma^{2}$ is known. In real applications, one may use sample splitting to estimate $\hat{\sigma}=\hat{\sigma}(y^{\prime},X^{\prime})$ from $(X^{\prime},y^{\prime})$ and then plug $\sigma=\hat{\sigma}$ into the construction of confidence sets. As long as $\hat{\sigma}$ is consistent, all the asymptotic results in this work still hold. For high-dimensional linear models, the scaled lasso provides a consistent $\hat{\sigma}$ (Sun and Zhang 2012).

3 Competing methods

To illustrate the effectiveness of our two-step Stein method, we first present three alternative procedures that can be derived by extending ideas from construction of nonparametric regression confidence sets in conjunction with lasso estimation. Since all of them make use of oracle properties, we review an error bound for lasso prediction due to Bickel et al. (2009).

3.1 Lasso prediction error

Given $X$, $y$ and $\lambda>0$, consider the lasso estimator $\hat{\beta}=\hat{\beta}(y,X;\lambda)$ defined as in (2.17). Let $\omega(X)=\max_{j}(\|X_{j}\|^{2}/n)$. Error bounds of lasso prediction have been established under the restricted eigenvalue assumption (Bickel et al. 2009). For $S\subseteq[p]$ and $c_{0}>0$, define the cone

[TABLE]

We say the design matrix $X$ satisfies $\textup{RE}(s,c_{0})$, for $s\in[p]$ and $c_{0}>0$, if

[TABLE]

Lemma 5** (Theorem 7.2 in Bickel et al. (2009)).**

Let $n\geq 1$ and $p\geq 2$. Suppose that $\|\beta\|_{0}\leq s$ and $X$ satisfies Assumption $\textup{RE}(s,3)$. Choose $\lambda=K\sigma\sqrt{\log(p)/n}$ for $K>2\sqrt{2}$. Then we have

[TABLE]

Remark 6*.*

The original theorem in Bickel et al. (2009) assumes that all the diagonal elements of the Gram matrix $X^{\mathsf{T}}X/n$ are $1$ for simplicity, while we remove this assumption by including the term $\omega(X)$.

3.2 Another adaptive method

Here we develop another adaptive method following the procedure in Section 3 of Robins and van der Vaart (2006), which constructs a confidence set for $\mu$ from $y\sim\mathcal{N}_{n}(\mu,\sigma^{2}\mathbf{I}_{n})$ via sample splitting. Applied to the linear model (1.1), the method can be described as follows. Split the original data set into $(X^{\prime},y^{\prime})$ and $(X,y)$, of which the former is used to obtain an initial lasso estimate $\hat{\beta}=\hat{\beta}(y^{\prime},X^{\prime};\lambda)$ (2.17), and the latter is used to compute two quantities

[TABLE]

where $R_{n}$ is an estimate of the loss $\|X\beta-X\hat{\beta}\|^{2}/n$. Then, a confidence ball for $\mu=X\beta$ is constructed in the form of

[TABLE]

where $z_{\alpha}$ is the $(1-\alpha)$ quantile of the standard normal distribution. Note that $\hat{\tau}_{n}$ in (3.5) contains the term $\|\mu-X\hat{\beta}\|$ as well so an explicit form of the confidence ball is

[TABLE]

where $r_{a}$ is the radius.

To establish the convergence rate of the diameter of $\widehat{C}_{a}$, we need an assumption, similar to RE$(s,c_{0})$, on the restricted maximum eigenvalue of $X^{\mathsf{T}}X/n$ over the cone $\mathscr{C}(S,c_{0})$ (3.1). For $s\in[p]$ and $c_{0}>0$, let

[TABLE]

Theorem 6**.**

The $(1-\alpha)$ confidence set $\widehat{C}_{a}$ (3.5) is honest for all $\beta\in\mathbb{R}^{p}$. Suppose $s\log p=o(n)$, the sequence $X=X(n)$ satisfies

[TABLE]

and so does the sequence $X^{\prime}=X^{\prime}(n)$. Then with a proper choice of $\lambda\asymp\sqrt{\log p/n}$, for any $\beta\in\mathscr{B}(s)$ the diameter

[TABLE]

These properties have been informally discussed in the introduction (Section 1). Although $\widehat{C}_{a}$ is also honest over the entire parameter space, the upper bound on its diameter critically depends on the sparsity of $\beta$. The scaling $s\log p=o(n)$ is the minimum requirement for the lasso to be consistent in estimating $\mu$ or $\beta$. In general, this scaling is also needed for the RE assumption to hold with $\liminf_{n}\kappa(2s,3;X)>0$ (Negahban et al. 2012) and for the upper bound on $|\widehat{C}_{a}|$ to be informative. This is different from the universal bound (2.15) on $\mathbb{E}|\widehat{C}|^{2}$ for the two-step method. The diameter $|\widehat{C}_{a}|$ adapts to the optimal rate for sufficiently sparse $\beta$ as $s\log p=O(\sqrt{n})$; see Remark 2 for related discussion. Our numerical results in Section 4.4 demonstrate that $|\widehat{C}_{a}|$ can be 10 times larger than the diameter of our two-step Stein method when $\beta$ is not sparse.

3.3 An oracle lasso method

We calculate the lasso $\hat{\beta}=\hat{\beta}(y,X;\lambda)$ from the whole data set without sample splitting, which we denote by $(X,y)$ in this subsection.

Assuming the true sparsity $s_{\beta}=\|\beta\|_{0}$ is known (the oracle), a $(1-\alpha)$ confidence ball for $X\beta$ is constructed as

[TABLE]

where $c_{o}(\alpha)$ is a constant depending on the design matrix $X$ and the tuning parameter $\lambda$. We estimate $c_{o}(\alpha)$ by a similar procedure to be described in Section 3.4 for a two-step lasso method. Although there are sharper upper bounds, e.g. $O(s_{\beta}\log(p/s_{\beta})/n)$, for lasso prediction error (e.g. Chapter 11 in Hastie et al. (2015)), our choice of $\lambda$ is tuned to achieve the desired coverage rate in our numerical results and thus the corresponding $r_{o}$ is already optimized in this sense.

It should be pointed out that the oracle lasso is not implementable in practice since the true sparsity $s_{\beta}$ is unknown. In theory, it can build a confidence set with a diameter on the order of $(s_{\beta}\log p/n)^{1/2}$, potentially faster than the rate $n^{-1/4}$, however, the constant $c_{o}(\alpha)$ can be large and difficult to approximate. Indeed, in comparison with the oracle lasso, our method often constructs confidence sets with a smaller volume even under highly sparse settings, which highlights the practical usefulness of our two-step method.

3.4 A two-step lasso method

To appreciate the advantage of using Stein estimates in the shrinkage step of our construction, we compare our method with a two-step lasso method, in which we replace the Stein estimate by the lasso to build a confidence set for $\mu_{\perp}$, the mean for weak signals. Consider the two-step method in Section 2.2 with a given candidate set $A$. Let $k=\operatorname{rank}(X_{A})$ and further assume $A$ contains strong signals only, that is, $A\subseteq\operatorname{supp}(\beta)$. We use the same method to find $\hat{\mu}_{A}$ and $r_{A}$ (2.6) in the projection step. Like the oracle lasso, we assume the true sparsity $s_{\beta}=\|\beta\|_{0}$ is given and construct a confidence set for $\mu_{\perp}$ based on the error bound for lasso prediction.

Apply lasso on $(P_{A}^{\perp}X,y_{\perp})=(P_{A}^{\perp}X,P_{A}^{\perp}y)$ with a tuning parameter

[TABLE]

to find the estimate

[TABLE]

It is natural to estimate the center $\mu_{\perp}=P_{A}^{\perp}\mu$ by the lasso prediction $\hat{\mu}_{\perp}=P_{A}^{\perp}X\tilde{\beta}$. As a corollary of Lemma 5, we find an error bound for $\|\hat{\mu}_{\perp}-\mu_{\perp}\|^{2}$:

Corollary 7**.**

Let $n\geq 1$ and $p\geq 2$. Suppose that $\|\beta\|_{0}\leq s$ and Assumption $\textup{RE}(s,3)$ holds for $X$. Choose $\lambda_{2}$ as in (3.7). Then for any fixed $A\subseteq\operatorname{supp}(\beta)$ with $k=\operatorname{rank}(X_{A})<s$, we have

[TABLE]

Accordingly, the radius for weak signals is chosen as

[TABLE]

where $c_{l}(\alpha/2)=c_{l}(\alpha/2;P_{A}^{\perp}X)$ is a constant. Lastly, we combine $(\hat{\mu}_{\perp},r_{\perp})$ with $(\hat{\mu}_{A},r_{A})$ as in (2.4) to define the confidence set $\widehat{C}$.

Again we use sample splitting to define the candidate set ${A}$ by thresholding the lasso estimate $\hat{\beta}(y^{\prime},X^{\prime};\lambda)$ in (2.17) with a threshold value $\tau=\Omega_{p}(\|\hat{\beta}-\beta\|_{\infty})$ so that $\mathbb{P}\left({A}\subseteq\operatorname{supp}(\beta)\right)\to 1$, satisfying the assumption in Corollary 7. Upper bounds on $\|\hat{\beta}-\beta\|_{\infty}$ are available under certain conditions; see, for example, Theorem 11.3 in Hastie et al. (2015).

Remark 7*.*

Suppose $\beta$ is sufficiently sparse so that $s_{\beta}\log p\ll\sqrt{n}$. Then, it follows that both $r_{A}$ and $r_{\perp}$ of the two-step lasso converge faster than the rate of $n^{-1/4}$. This is not surprising and shows the advantage of the oracle knowledge of the true sparsity $s_{\beta}$. Of course, in practice we do not know $s_{\beta}$ and therefore, this two-step lasso method, like the oracle lasso, is not implementable for real problems. The numerical comparisons in the next section will show that our two-step Stein method, which does not use the true sparsity in its construction, is more appealing than the two-step lasso: Its adaptation to the underlying sparsity is comparable to the two-step lasso, while its coverage turns out to be much more robust.

We follow the same procedure as the two-step Stein method to implement the two-step lasso method with multiple candidate sets $A_{m},m=1,\ldots,M$ — threshold $\hat{\beta}(y^{\prime},X^{\prime};\lambda)$ with a sequence of threshold values to construct $A_{m}$ (2.28) and then choose the confidence set with the minimum volume or diameter. The main difference is how to approximate $c_{l}(\alpha)$ in (3.10), which is done by the following approach.

We first use $b=\max_{i\in[p]}(X_{i}^{\prime\mathsf{T}}y^{\prime})/\|X_{i}^{\prime}\|^{2}$ as a rough upper bound for $\|\beta\|_{\infty}$. For $j=1,2,\ldots,N$, we draw an $s_{\beta}$-sparse vector, $\gamma_{j}\in\mathbb{R}^{p}$, of which the nonzero components follow $\mathcal{U}(-b,b)$. Then we sample $Y_{j}^{*}\sim\mathcal{N}_{n}(X\gamma_{j},\sigma^{2}\mathbf{I}_{n})$ and calculate lasso estimate $\hat{\gamma}_{j}(\lambda)=\hat{\beta}(Y_{j}^{*},X;\lambda)$ as in (2.17) with the tuning parameter $\lambda$ for all $j$. Let $c_{j}={\|X(\hat{\gamma}_{j}(\lambda)-\gamma_{j})\|^{2}}/{(\sigma^{2}s_{\beta}\log p)}$. For a large $N$, $c_{l}(\alpha)$ can be approximated by the $(1-\alpha)$ quantile of $\{c_{j}\}$. Here, $\lambda=\nu\cdot K\sigma^{2}\sqrt{\log p/n}$, where $\nu\leq 1$ is a pre-determined constant. This choice is slightly smaller than the theoretical value in Lemma 5, but gives a stable estimate of $c_{l}(\alpha)$ with the desired coverage. As we calculate $b$ with $(X^{\prime},y^{\prime})$ in the above, our estimate of $c_{l}(\alpha)$ is independent of the response $y$. It is possible that a candidate set $A_{m}$ defined by (2.28) may contain $s$ or more predictors. In this case, we will only include the largest $s-1$ predictors in terms of their absolute lasso coefficients, as Corollary 7 requires $|A_{m}|<s$.

4 Numerical results

We will first compare our method with the above competing methods when $\beta$ is sparse relative to the sample size, i.e., $s/n$ is small, and then consider the more challenging settings in which the sparsity $s$ is comparable to $n$.

4.1 Simulation setup

The rows of $X$ and $X^{\prime}$, both of size $n\times p$, are independently drawn from $\mathcal{N}_{p}(0,\Sigma)$ and the columns are normalized to have an identical $\ell_{2}$ norm. We use three designs for $\Sigma$ as in Dezeure et al. (2015):

[TABLE]

The support of $\beta$ is randomly chosen and its $s$ nonzero components are generated in two ways:

They are drawn independently from a uniform distribution $\mathcal{U}(-b,b)$. 2. 2.

Half of the nonzero components follow $\mathcal{U}(-b,b)$ while the other half following $\mathcal{U}(-0.2,0.2)$, so there are two signal strengths under this setting.

Lastly, $y$ and $y^{\prime}$ are drawn from $\mathcal{N}_{n}(X\beta,\sigma^{2}\mathbf{I}_{n})$ and $\mathcal{N}_{n}(X^{\prime}\beta,\sigma^{2}\mathbf{I}_{n})$, respectively. In our results, we chose $n=n^{\prime}=200$, $p=800$, $\sigma^{2}=1$ and $s=10$, and $b$ took 10 values evenly spaced between $(0,1)$ and $(1,5)$. In total, we had $60$ simulation settings, each including one design for $\Sigma$, one way of generating $\beta$, and one value for $b$. Under each setting, $100$ data sets were generated independently, so that the total number of data sets used in this simulation study was 6,000.

The confidence level $1-\alpha$ was set to $0.95$. The threshold values $\{a_{m}\}$ in (2.28) were evenly spaced from [math] to $4$ with a step of $0.05$. All the competing methods use lasso in some of the steps, and the tuning parameter $\lambda$ was chosen by three approaches:

1. the minimum theoretical value in Bickel et al. (2009), $\lambda_{val}=2\sqrt{2}\sigma\sqrt{\log p/n}$,
2. cross validation $\lambda_{cv}$, and 3) one standard error rule $\lambda_{1se}$. For the one standard error rule, we choose the largest $\lambda$ whose test error in cross validation is within one standard error of the error for $\lambda_{cv}$. Since it is time-consuming to approximate $c_{o}(\alpha)=c_{o}(\alpha;X,\lambda)$ for the oracle lasso when $\lambda$ is chosen by a data-dependent way, we set $c_{o}(\alpha;X,\lambda_{cv})=\eta_{1}c_{o}(\alpha;X,\lambda_{val})$ and $c_{o}(\alpha;X,\lambda_{1se})=\eta_{2}c_{o}(\alpha;X,\lambda_{val})$, where the factors $\eta_{k}$ were chosen such that the overall coverage rate across data sets simulated with $b>0.3$ was around the desired level.

Unlike the adaptive method in Section 3.2 and our two-step methods, the oracle lasso method does not require sample splitting. Consequently, a confidence set is constructed based on the whole data set including both $(X,Y)$ and $(X^{\prime},Y^{\prime})$ for a fair comparison. We compare the geometric average radius $\bar{r}=(r_{A}^{|A|}r_{\perp}^{n-|A|})^{1/n}$ of our two-step methods with $r_{a}$ of the adaptive method and $r_{o}$ of the oracle lasso. This is equivalent to comparing the volumes of the confidence sets.

4.2 Results on the two-step Stein method

In this subsection we compare the two-step Stein method with the adaptive method and the oracle lasso. The constants $c_{1}$ and $c_{2}$ of our method were chosen by minimizing the volume in (2.29) with upper bound $E=10$.

Figure 1 compares the geometric average radius $\bar{r}$ among the three methods against the signal strength $b$ under the first way of drawing $\beta$. Every point in a panel was computed by averaging $\bar{r}$ from $100$ data sets under a particular simulation setting. It is seen from the figure that $\bar{r}$ by our method was dramatically smaller than the other two methods for almost every setting. This suggests that the volumes of our confidence sets were orders of magnitude smaller than the other two methods, as the ratio of the radii will be raised to the power of $n=200$ for comparing volumes. When $X$ was drawn from the equal correlation (Equi.corr) design, $\bar{r}$ of the oracle lasso and the adaptive methods kept increasing as $b$ increased, while $\bar{r}$ by our method became stable after $b>2$. Overall, the equal correlation design was more challenging than the other two designs, for which our method outperformed the other two methods with the largest margin. Unlike the other two methods, our method was less sensitive to the choices of $\lambda$ and the designs of $X$. Essentially, $r_{A}$ and $r_{\perp}$ by our method are determined by the candidate set $A$. Even if a different $\lambda$ is used, our method can choose adaptively an optimal $A$ close to $\operatorname{supp}(\beta)$, showing the advantage of using multiple candidate sets.

In a similar way, Figure 2 plots $\bar{r}$ against $b$ in the second scenario of drawing $\beta$. When $b$ is large (e.g, $b\geq 1$), the $\beta$ contains a mixture of weak and strong signals. Again, we see that $\bar{r}$ of our method was smaller than the other two competitors for most settings. The average radius by our method often decreased as $b>1$, which shows that our method can properly distinguish strong signals and weak signals.

The coverage rates, each computed from $100$ data sets, for each of the three ways of choosing $\lambda$ are summarized in Figure 3. We pooled the results from three types of design matrices together in the figure, because the coverage rates distributed similarly across them. The coverage rates of our method matched the desired $95\%$ confidence level very well, with coverage rate $>0.9$ for $96\%$ of the cases. This result is particularly satisfactory for a quite small sample size of $n=200$. The adaptive method also showed a good coverage, but slightly more conservative than the desired level. The oracle lasso had the most variable coverage rate across different settings when $\lambda$ was selected in a data-dependent way ($\lambda_{cv}$ or $\lambda_{1se}$). In fact, its coverage could drop below $0.5$ for these two cases (not shown in the figure). This shows the difficulty in practice to construct stable confidence sets using error bounds like (3.3) even with a known sparsity. Together with the results in Figures 1 and 2, this comparison demonstrates the advantage of the proposed two-step Stein method: It builds much smaller confidence sets, while closely matching the desired confidence level. In particular, our confidence sets were uniformly smaller than those by the adaptive method (Section 3.2) for all simulation settings and all choices of $\lambda$.

4.3 Comparison with the two-step lasso method

We discussed in Section 2.5 two ways to choose $c_{1}$ and $c_{2}$, that is, by minimizing the volume or by minimizing the diameter of the confidence set for our proposed two-step framework. Here we compare the two-step Stein method and the two-step lasso, each with the two ways to choose the constants. The two-step Stein method by minimizing the volume (abbreviated as TSV) is the same method used in the previous comparison. Similarly, we use the short-hand TSD, TLV, and TLD for the two-step Stein method by minimizing diameter, the two-step lasso method by minimizing volume and by minimizing diameter, respectively. The true sparsity $s=10$ was given to the two-step lasso methods. Only the first scenario of generating $\beta$ was considered in this comparison, since most results in the second scenario were similar. Figure 4 shows the plots of radius against $b$ by the four methods under different settings, while Figure 5 reports the distribution of the coverage rates. The two-step lasso methods apply the lasso twice, one to generate candidate sets $A_{m}$ and the other to compute $\hat{\mu}_{\perp}$ and $r_{\perp}$ for weak signals. To clarify, the three ways of choosing $\lambda$ in these figures refer to the step to generate candidate sets $A_{m}$, while $\lambda_{2}$ in (3.8) was set to $\nu K\sigma^{2}\sqrt{\log(p-|A|)/(n-|A|)}$, where $\nu=0.5$ in our simulation.

We make the following observations from the two figures. First, the two-step Stein methods showed a substantially more satisfactory coverage than the two-step lasso methods. The coverage was close to 0.95 for both TSV and TSD, while the coverage rates of TLV and TLD had a much larger variance and were especially poor when $\lambda$ was chosen via cross validation. The confidence sets by the two-step lasso methods had a slightly smaller average radius than the two-step Stein methods for the Toeplitz and the exponential decay designs. However, given their low and unstable coverage rates, this does not imply the two-step lasso methods constructed better confidence sets. Recall that $|\widehat{C}|=O_{p}(n^{-1/4}\vee\sqrt{s/n})$ for the two-step Stein methods and $|\widehat{C}|=O_{p}(\sqrt{s\log p/n})$ for the two-step lasso methods. The signals were very sparse in our simulation, with $s=10$ much smaller than $p$, favorable for the two-step lasso methods. Even so, we find the two-step Stein methods very competitive, noting that the radii of both TSV and TSD were actually comparable or slightly smaller than the two-step lasso methods for the equal correlation designs, in which the predictors were highly correlated. This comparison demonstrates that the two-step Stein method is more appealing in practice, as it does not require any prior knowledge about the underlying sparsity but gives a better and more stable coverage. Second, both ways of choosing the constants $c_{1}$ and $c_{2}$ worked well for the two-step Stein method. On the contrary, it is seen from Figure 5 that the coverage rate of TLV was significantly lower than that of TLD in the bottom two panels. Lastly, between using $\lambda_{cv}$ and $\lambda_{1se}$ in the lasso for defining candidate sets $A_{m}$, we recommend the latter, as it tends to give comparable radii but a better coverage, especially for the two-step lasso.

We also compared the performance between the oracle lasso method and TLD, both constructing confidence sets based on the lasso prediction (3.3) with a known sparsity. The coverage rates of the two methods were quite comparable as reported in Figures 3 and 5. The geometric average radius of the oracle lasso method (Figure 1) was $2$ to $5$ times that of TLD (Figure 4). The difference was especially significant when the signal strength was high (large $b$). This comparison confirms that, by separating strong and weak signals, our two-step framework can greatly improve the efficiency of the constructed confidence sets.

4.4 Dense signal settings

We have shown the advantages of our two-step Stein method in the last two subsections under sparse settings. Recall that the dimension of our data was $(n,p)=(200,800)$ with sparsity $s=10$ for $\beta$ in the previous comparisons. The goal of this subsection is to illustrate the stable performance of our method when the true signal is dense. As such, we changed the sparsity to $s=100$ for the first way of generating $\beta$ and $s=200$ for the second way of generating $\beta$. We focused on the equal correlation design, which was the most difficult one among the three designs. With the same set of values for the signal strength $b$, we had 20 distinct parameter settings for data generation in this comparison, and again we simulated 100 data sets under each setting. The tuning parameter $\lambda$ was selected as $\lambda_{1se}$ for all the results here.

Figure 6 compares the geometric average $\bar{r}$ against $b$ and the coverage among the adaptive method, the oracle lasso and our two-step Stein method. In all the scenarios reported in panels (a) and (b), our method outperformed the other two methods with very big margins in terms of the volume of a confidence set. For $b>1$, the radius of our method approached the naive radius ${(\chi^{2}_{n,\alpha}/n)^{1/2}}$ as suggested by Theorem 1, while the radii of the oracle lasso and the adaptive methods kept increasing to much greater than the naive $\chi^{2}$ radius. This shows that the two competing methods failed to construct acceptable confidence sets when the signal was dense. Since the sparsity level $s$ is comparable to $n$ for the data sets here, the upper bounds for the diameters of these two methods, $|\widehat{C}_{o}|=O_{p}(\sqrt{s\log p/n})$ and $|\widehat{C}_{a}|=O_{p}(n^{-1/4}+\sqrt{s\log p/n})$, are no longer useful or even valid. It is seen from Figure 6(c) and (d) that the coverage rates of the two-step Stein method were much better than the oracle lasso, but slightly lower than the adaptive method. Nevertheless, our confidence sets still maintained a minimum coverage of $0.9$ in most cases, which is quite satisfactory given the way smaller diameters than the adaptive method.

To understand the behavior of our method in this dense signal setting, we examined the number of variables selected as strong signals in the set $A$, i.e., $k=|A|$. Figure 7 displays the box plot of $k$ across 100 data sets for each value of $b$ under the first way to generate $\beta$. When $b\leq 1$, our two-step method still chose a nonempty candidate set, but $k$ dropped to [math] for $b\geq 2$, i.e., $A=\varnothing$. Note that the radius of our method will be close to the naive $\chi^{2}$ radius when $k=n$ or $k=0$; see (2.14) in Theorem 1. When the signal strength $b\leq 1$, some small nonzero coefficients are close to zero so $\beta$ is effectively quite sparse, in which case the lasso can select a good subset $A$ of strong signals. On the contrary, when $b$ is large, the lasso will not be able to select a majority of the strong signals, leaving $\|\mu_{\perp}\|=\|P_{A}^{\perp}\mu\|$ too big. In this setting, our method automatically adjusts its “optimal” choice to $A=\varnothing$, constructing a confidence set centered at the Stein estimate $\hat{\mu}(y;0)$ (2.8) with radius estimated via the SURE.

5 Discussion

For high-dimensional regression, oracle inequalities for sparse estimators cannot be directly utilized to construct honest and adaptive confidence sets due to the unknown signal sparsity. To overcome this difficulty, we have developed a two-step Stein method, via projection and shrinkage, to construct confidence sets for $\mu=X\beta$ in (1.1) by separating signals into a strong group and a weak group. Not only is honesty achieved over the full parameter space $\mathbb{R}^{p}$, but also our confidence sets can adapt to the sparsity and strength of $\beta$. We also implemented an adaptive way to choose a proper subspace for the projection step among multiple candidate sets, which protects our method from a poor separation between strong and weak signals. Our two-step Stein method showed very satisfactory performance in extensive numeric comparisons, outperforming other competing methods under various parameter settings.

The focus of this work is on the confidence set for $\mu=X\beta$. Although related, it is different from the problem of inference on $\beta$. In general, it is difficult to infer a confidence set for $\beta$ from the confidence set for $X\beta$ without any constraint on $X$ and $\beta$, because $X$ does not have a full column rank under the high-dimensional setting. However, if we know that $\|\beta\|_{0}\leq s$, then a confidence set $\widehat{C}$ for $\mu$ can be converted into a confidence set for $\beta$ as $\widehat{B}:=\{\beta\in\mathscr{B}(s):X\beta\in\widehat{C}\}$, which is the union of $s$-dimensional subspaces intersecting $\widehat{C}$. It is interesting future work to study the convergence rate of $\widehat{B}$ and related computational issues, such as how to draw $\beta$ from $\widehat{B}$. On the other hand, if $X$ satisfies $\textup{SRC}(s,c_{*},c^{*})$, then

[TABLE]

A hypothesis test about the mean $X\beta$ can be carried out by using the confidence set $\widehat{C}$ to obtain a lower bound on $\|X\beta\|$, which carries over to a lower bound on $\|\beta\|$ with the above inequality and thus can be used to perform a test about $\beta$. See Nickl and van de Geer (2013) for a related discussion. We have also demonstrated that our method works well even when the underlying $\beta$ is dense, e.g. $\|\beta\|_{0}\asymp n$, which is important for practical applications. See Bradic et al. (2018) for recent theoretical results on high-dimensional inference for non-sparse $\beta$.

Another future direction is to incorporate the confidence set $\widehat{C}$ with the method of estimator augmentation (Zhou 2014; Zhou and Min 2017) for lasso-based inference. Estimator augmentation can be used to simulate from the sampling distribution of the lasso without solving the lasso problem repeatedly, provided a point estimate of $\mu=X\beta$. Given $\widehat{C}$, one may randomize the point estimate of $\mu$ by sampling from the confidence set, which has been shown to improve the inferential performance of estimator augmentation (Min and Zhou 2019).

Appendix A Proofs

Proof of Lemma 2.

By the law of large number, we have

[TABLE]

where $\Phi^{-1}$ is the inverse of the cumulative distribution function of $\mathcal{N}(0,1)$. It follows from (2.6) and (A.1) that

[TABLE]

Let $\varepsilon_{\perp}=P_{A}^{\perp}\varepsilon$. Under the normality assumption of $\varepsilon$, we have

[TABLE]

It follows, by noting $\lVert\mu_{\perp}\rVert=o(\sqrt{n-k})$, that

[TABLE]

By plugging (A.3) in (2.12), we obtain

[TABLE]

If $k=O_{p}(\sqrt{n})$ and $\|\mu_{\perp}\|=O(n^{1/4})$, it follows from (A.2) and (A) that $|\widehat{C}|\asymp_{p}n^{-1/4}$. ∎

Proof of Theorem 2.

Under sparse Riesz condition, letting $G=A^{c}\cap\operatorname{supp}(\beta)$, we have

[TABLE]

which, together with $k=o(n)$ and $\|\beta_{A^{c}}\|=o(1)$, implies $\|\mu_{\perp}\|=o(\sqrt{n})=o(\sqrt{n-k})$. Thus, by Lemma 2, $r_{\perp}^{2}=O_{p}(n^{-1/2}+\lVert\beta_{A^{c}}\rVert^{2})$ and the rest of the proof is straightforward. ∎

Proof of Corollary 3.

Under the choice of $\lambda$ in this corollary and the assumption that $s\leq(s^{*}-1)/(2+4c^{*}/c_{*})$, Theorem 1 and Theorem 3 in Zhang and Huang (2008) imply that, for any $\epsilon>0$, there exists $N$ such that when $n>N$,

[TABLE]

where $M_{1}^{*}$ and $M_{2}^{*}$ are two constants depending on $c_{0}$, $c_{*}$ and $c^{*}$. It follows from (A.5) that

[TABLE]

Thus, we have

[TABLE]

Now, all the conditions in Theorem 2 are satisfied, leading to (2.19). Further, (A.6) implies that $S_{0}\subset A$ and thus $\|\beta_{A^{c}}\|=\|\beta_{A^{c}\cap A_{0}}\|\leq\|\beta_{A_{0}\backslash S_{0}}\|=O(n^{-1/4})$ with probability at least $1-\epsilon$. Consequently, (2.20) follows from (2.16). ∎

Proof of Corollary 4.

If $\mathbb{P}(A=A_{0})\to 1$, then the rate of $|\widehat{C}|$ in (2.23) follows immediately from (2.16) in Theorem 2. Thus, it remains to show that $\hat{\beta}^{\textup{mcp}}_{\lambda_{n},\gamma_{n}}=\hat{\beta}^{\textup{mcp}}_{\lambda,\gamma}(y^{\prime},X^{\prime})$ (2.22) is model selection consistent by verifying the conditions of the following corollary, which is a simplified version of Corollary 4.2 in Huang et al. (2012).

Corollary 8**.**

Let $\lambda_{\min}$ be the smallest eigenvalue of $(X^{\prime}_{A_{0}})^{\mathsf{T}}X^{\prime}_{A_{0}}/n$, $\tau_{n}=\sigma\sqrt{2\log s/(n\lambda_{\min})}$ and $\lambda^{*}=2\sigma\sqrt{2c^{*}\log(p-s)/n}$. Suppose that $X^{\prime}$ satisfies $\textup{SRC}(s^{*},c_{*},c^{*})$, where $0<c_{*}<c^{*}$ are constants and $s^{*}\geq(c^{*}/c_{*}+1/2)s$. If a sequence of $(\lambda_{n},\gamma_{n})$ satisfies $\inf_{A_{0}}|\beta_{j}|\geq\gamma_{n}\lambda_{n}+a_{n}\tau_{n}$ with $a_{n}\to\infty$, $\lambda_{n}\geq a_{n}\lambda^{*}$, $n\lambda_{n}^{2}/(4c^{*})>\sigma^{2}$ and $\gamma_{n}\geq c_{*}^{-1}\sqrt{4+c_{*}/c^{*}}$, then $\mathbb{P}\{\operatorname{supp}(\hat{\beta}^{\textup{mcp}}_{\lambda_{n},\gamma_{n}})=A_{0}\}\to 1$.

Under the SRC assumption $\lambda_{\min}$ is bounded from below by $c_{*}>0$. It follows from $\tau_{n}=O(\sqrt{\log s/n})$, $\lambda^{*}=O(\sqrt{\log p/n})$ and $\lambda_{n}\gg\sqrt{\log p/n}$ that there exists $a_{n}\to\infty$ such that $\lambda_{n}\geq a_{n}(\lambda^{*}\vee\tau_{n})$. Then we have the following: $\inf_{A_{0}}|\beta_{j}|\geq(\gamma_{n}+1)\lambda_{n}\geq\gamma_{n}\lambda_{n}+a_{n}\tau_{n}$, $\lambda_{n}\geq a_{n}\lambda^{*}$, and $n\lambda_{n}^{2}/(4c^{*})\gg\log p>\sigma^{2}$ when $n$ is sufficiently large. Thus all the conditions in Corollary 8 are satisfied under the assumptions of Corollary 4. This completes the proof.

Technically, we did not invoke the assumption $s\log p=o(n)$ in the proof. But it is required for the sparse Riesz condition to hold (e.g. for Gaussian designs). ∎

Proof of Lemma 3.

We have the following inequalities for any positive $x$ and degree of freedom of $n$ from Lemma 1 in Laurent and Massart (2000):

[TABLE]

The solutions of $2\sqrt{n}\sqrt{x_{1}}+2x_{1}=\sqrt{n}\delta$ and $2\sqrt{n}\sqrt{x_{2}}=\sqrt{n}\delta$ are plugged in (A.7) and (A.8) to obtain

[TABLE]

so that

[TABLE]

To finish the proof, we will show that

[TABLE]

is bounded by $\delta^{2}$ for any $n$. Replacing $\sqrt{1+2\delta/\sqrt{n}}$ with its Taylor expansion $1+\delta/\sqrt{n}+O(\delta^{2}/n)$ in (A.9), we get $f(n)=\delta^{2}+O(n^{-1/2})\to\delta^{2}$, as $n\to\infty$. If $f(n)$ is monotonically increasing in $n$, then $\delta^{2}$ is a tight upper bound of $f(n)$ for all $n$. Lastly, to prove the monotonicity, it suffices to show the derivative

[TABLE]

which can be verified easily. Now the proof is completed. ∎

Proof of Lemma 4.

Let

[TABLE]

A few steps of derivation shows that

[TABLE]

It follows from (A.10) that

[TABLE]

where the second probability on the right hand side is bounded by

[TABLE]

To prove the theorem, it suffices to show that all three probabilities in (A.11) can be bounded by either $D/(n-k)^{2}$ or $D/\delta^{4}$ for some constant $D>0$. Before that, we introduce the following three inequalities derived from Theorem 2 in Whittle (1960):

[TABLE]

for some constant $D_{1}$ depending on the moments of $\varepsilon_{i}$. In our case, $D_{1}$ only depends on the upper bound $d$ of the eighth moment. The first term of (A.11) can be bounded by

[TABLE]

Similarly, using (A.13) and (A.14), we can also show that both the second and the third terms are bounded by $D_{2}/(\sigma^{2}\delta^{4})$ for some $D_{2}>0$ depending only on $d$. Lastly, the proof is finished by letting $D=(16D_{1})\vee(D_{2}/\sigma^{2})$. ∎

Proof of Theorem 6.

The honesty of $\widehat{C}_{a}$ in (3.5) is guaranteed by Theorem 3.1 and Proposition 2.1 in Robins and van der Vaart (2006) with the only assumption $y/\sqrt{n}\sim\mathcal{N}_{n}(\mu/\sqrt{n},\sigma^{2}\mathbf{I}_{n}/n)$. It is not difficult to verify that $(X^{\prime},y^{\prime})$ satisfies all the conditions in Corollary B.2 and Theorem 7.2 of Bickel et al. (2009). Thus, with probability approaching one, we have $\|\hat{\beta}-\beta\|^{2}=O(s\log p/n)$ and $(\hat{\beta}-\beta)\in\mathscr{C}(A_{0},3)$, as defined in (3.1), with $A_{0}=\operatorname{supp}(\beta)$. By the definition of $\zeta(s,3;X)$, this implies that

[TABLE]

Again, by Theorem 3.1 in Robins and van der Vaart (2006), we have

[TABLE]

which completes the proof. ∎

Proof of Corollary 7.

Rewrite orthogonal matrix $P_{A}^{\perp}=VV^{\mathsf{T}}$, where $V\in\mathbb{R}^{n\times(n-k)}$ consists of orthogonal unit column vectors. Write the lasso estimate in (3.8) as $\tilde{\beta}=F(y_{\perp},P_{A}^{\perp}X;n\lambda_{2})$, where $F$ is understood as a mapping with a parameter $n\lambda_{2}>0$. Since $P_{A}^{\perp}X_{A}=0$, the loss in (3.8) becomes

[TABLE]

which demonstrates that $\tilde{\beta}_{A}=0$ and $\tilde{\beta}_{A^{c}}=F(V^{\mathsf{T}}y,V^{\mathsf{T}}X_{A^{c}};n\lambda_{2})$. Moreover, we have

[TABLE]

We will verify that the lasso problem, $\tilde{\beta}_{A^{c}}=F(V^{\mathsf{T}}y,V^{\mathsf{T}}X_{A^{c}};n\lambda_{2})$, satisfies all the assumptions in Lemma 5 so that we can apply (3.3) to bound the prediction error on the left side of (A.16). Since $A\subseteq\operatorname{supp}(\beta)$, we have $\|\beta_{A^{c}}\|_{0}\leq s-k$. Next, we show $V^{\mathsf{T}}X_{A^{c}}\in\mathbb{R}^{(n-k)\times(p-k)}$ satisfies $\textup{RE}(s-k,3)$. Let $D$ be any subset of $[p-k]$ such that $|D|\leq(s-k)$. For any nonzero $\gamma\in\mathbb{R}^{p-k}$ in the cone $\mathscr{C}(D,3)$, a vector $\delta=(\eta,\gamma)\in\mathbb{R}^{p}$ can always be constructed satisfying

[TABLE]

since $P_{A}X_{A^{c}}\gamma\in\operatorname{span}(X_{A})$. Define a mapping $g:i\mapsto i+|A|$ for $i\in[p]$ and let $B=[|A|]\cup g(D)\subset[p]$. Then $|B|=|A|+|D|\leq s$, and $\delta\in\mathscr{C}(B,3)$ because

[TABLE]

where the second step is due to $\gamma\in\mathscr{C}(D,3)$. Based on that $X$ satisfies $\textup{RE}(s,3)$, we arrive at the following inequality:

[TABLE]

which shows that $\textup{RE}(s-k,3)$ holds for $V^{\mathsf{T}}X_{A^{c}}$ and $\kappa(s-k,3;V^{\mathsf{T}}X_{A^{c}})\geq\sqrt{n/(n-k)}\kappa(s,3;X)$. Lastly, $n\lambda_{2}=K\sigma\sqrt{n\log(p-k)}\geq K\sigma\sqrt{(n-k)\log(p-k)}$, as required in Lemma 5.

So far, we have shown that $(V^{\mathsf{T}}X_{A^{c}},V^{\mathsf{T}}y)$ and $\lambda_{2}$ satisfy all the conditions in Lemma 5, which with (A.16) implies that

[TABLE]

for any $A\subseteq\operatorname{supp}(\beta)$. Then inequality (3.9) immediately follows by noting that $\omega(V^{\mathsf{T}}X_{A^{c}})\leq\omega(X)$ and substituting $\kappa(s-k,3;V^{\mathsf{T}}X_{A^{c}})$ with $\sqrt{n/(n-k)}\kappa(s,3;X)$. ∎