Hypothesis Testing in High-Dimensional Regression under the Gaussian Random Design Model: Asymptotic Theory

Introduction

In matrix form, letting $\bm{y}=(y_{1},\dots,y_{n})^{\sf T}$ and denoting by ${\bm{X}}$ the matrix with rows $\bm{x}_{1}^{\sf T}$,$\cdots$, $\bm{x}_{n}^{\sf T}$ we have

We are interested in high-dimensional settings where the number of parameters exceeds the sample size, i.e., $p>n$, but the number of non-zero entries of $\bm{\theta}_{0}$ (to be denoted by $s_{0}$) is smaller than $p$. In this situation, a recurring problem is to select the non-zero entries of $\bm{\theta}_{0}$ that hence can provide a succinct explanation of the data. The vast literature on this topic is briefly overviewed in Section 1.1.

The Gaussian design assumption arises naturally in some important applications. Consider for instance the problem of learning a high-dimensional Gaussian graphical model from data. In this case we are given i.i.d. samples $\bm{z}_{1},\bm{z}_{2},\dots,\bm{z}_{n}\sim{\sf N}(0,\bm{K}^{-1})$, with $\bm{K}$ a sparse positive definite matrix whose non-zero entries encode the underlying graph structure. As first shown by Meinshausen and Bühlmann , the $i$-th row of $\bm{K}$ can be estimated by performing linear regression of the $i$-th entry of the samples ${\bm{z}}_{1},{\bm{z}}_{2},\dots,{\bm{z}}_{n}$ onto the other entries . This reduces the problem to a high-dimensional regression model under Gaussian designs. Standard Gaussian designs were also shown to provide useful insights for compressed sensing applications .

In statistics and signal processing applications, it is unrealistic to assume that the set of nonzero entries of $\bm{\theta}_{0}$ can be determined with absolute certainty. The present paper focuses on the problem of quantifying the uncertainty associated to the entries of $\bm{\theta}_{0}$. More specifically, we are interested in testing null-hypotheses of the form:

for $i\in[p]\equiv\{1,2,\dots,p\}$ and assigning p-values for these tests. Rejecting $H_{0,i}$ is equivalent to stating that $\theta_{0,i}\neq 0$.

Any hypothesis testing procedure faces two types of errors: false positives or type I errors (incorrectly rejecting $H_{0,i}$, while $\theta_{0,i}=0$), and false negatives or type II errors (failing to reject $H_{0,i}$, while $\theta_{0,i}\neq 0$). The probabilities of these two types of errors will be denoted, respectively, by $\alpha$ and $\beta$ (see Section 2.1 for a more precise definition). The quantity $1-\beta$ is also referred to as the power of the test, and $\alpha$ as its significance level. It is trivial to achieve $\alpha$ arbitrarily small if we allow for $\beta=1$ (never reject $H_{0,i}$) or $\beta$ arbitrarily small if we allow for $\alpha=1$ (always reject $H_{0,i}$). This paper aims at optimizing the trade-off between power $1-\beta$ and significance $\alpha$.

Without further assumptions on the problem structure, the trade-off is trivial and no non-trivial lower bound on $1-\beta$ can be established. Indeed we can take $\theta_{0,i}\neq 0$ arbitrarily close to $0$, thus making $H_{0,i}$ in practice indistinguishable from its complement. We will therefore assume that, whenever $\theta_{0,i}\neq 0$, we have $|\theta_{0,i}|>\mu$ as well. The smallest value of $\mu$ such that the power and significance reach some fixed non-trivial value (e.g., $\alpha=0.05$ and $1-\beta\geq 0.9$) has a particularly compelling interpretation, and provides an answer to the following question: What is the minimum magnitude of $\theta_{0,i}$ to be able to distinguish it from the noise level, with a given degree of confidence?

In the case of orthogonal designs we have $n=p$ and $\bm{X}^{{\sf T}}\bm{X}=n{\rm I}_{n\times n}$. By an orthogonal transformation, we can limit ourselves to ${\bm{X}}=\sqrt{n}\,{\rm I}_{n\times n}$, i.e., $y_{i}=\sqrt{n}\,\theta_{0,i}+w_{i}$. Hence testing hypothesis $H_{0,i}$ reduces to testing for the mean of a univariate Gaussian.

It is easy to see that we can distinguish the $i$-th entry from noise only if its size is at least of order $\sigma/\sqrt{n}$. More precisely, for any $\alpha\in(0,1)$, $\beta\in(0,\alpha)$, we can achieve significance $\alpha$ and power $1-\beta$ if and only if $|\theta_{0,i}|\geq c(\alpha,\beta)\,\sigma/\sqrt{n}$ for some constant $c(\alpha,\beta)$ [7, Section 3.9].

To move away from the orthogonal case, consider standard Gaussian designs. Several papers studied the estimation problem in this setting . The conclusion is that there exist computationally efficient estimators $\bm{\widehat{\theta}}$ that are consistent (in high-dimensional sense) for $n\geq c_{1}s_{0}\log(p/s_{0})$, with $c_{1}$ a numerical constant. By far the most popular such estimator is the Lasso or Basis Pursuit Denoiser .

On the other hand, no practical estimator is known that is consistent under a significantly smaller sample size (impossibility results have been proven in this direction, see e.g. ). We expect hypothesis testing to require at least as large sample size as point estimation, i.e. $n\geq c_{0}s_{0}\log(p/s_{0})$ for some $c_{0}=c_{0}(\alpha,\beta)$.

These simple remarks motivate the following seemingly simple question:

Assume standard Gaussian design $\bm{X}$, and fix $\alpha,\beta\in(0,1)$. Are there constants $c=c(\alpha,\beta)$, $c_{1}=c_{1}(\alpha,\beta)$ and a hypothesis testing procedure achieving the desired significance and power for all $\mu\geq c\sigma/\sqrt{n}$, $n\geq c_{1}s_{0}\log(p/s_{0})$?

Despite the seemingly idealized setting, the answer to this question is highly non-trivial. To document this point, we consider in Appendix C two hypothesis testing methods that were recently proposed by Zhang and Zhang , and by Bühlmann . These approaches apply to a broader class of design matrices ${\bm{X}}$ that satisfy the restricted eigenvalue property . We show that, when specialized to the case of standard Gaussian designs $\bm{x}_{i}\sim{\sf N}(0,{\rm I}_{p\times p})$, these methods require $|\theta_{0,i}|\geq\mu=c\,\max\{\sigma s_{0}\log p/\,n,\sigma/\sqrt{n}\}$ to reject hypothesis $H_{0,i}$ with a given degree of confidence (with $c$ being a constant independent of the problem dimensions). In other words, these methods are guaranteed to succeed only if the coefficient to be tested is larger than the ideal scale $\sigma/\sqrt{n}$, by a diverging factor of order $s_{0}\log p/\sqrt{n}$. In particular, the results of do not allow to answer the above question.

In this paper, we answer positively to this question. As in , our approach is based on the Lasso estimator

We use the solution to this problem to construct a debiased estimator of the form

A similar approach was developed independently in and (after a a preprint version of the present paper became available online) in . Apart from differences in the construction of $\bm{M}$, the three papers differ crucially in the assumptions and the regime analyzed, and establish results that are not directly comparable. In the present paper we assume a specific (random) model for the design matrix $\bm{X}$. In contrast and assume deterministic designs, or random designs with general unknown covariance.

On the other hand, we are able to analyze a regime that is significantly beyond reach of the mathematical techniques of , even for the very special case of standard Gaussian designs. Namely, for standard designs, we consider $\mu$ of order $\sigma/\sqrt{n}$, and $n$ of order $s_{0}\log(p/s_{0})$.

This regime is both challenging and interesting because $\theta_{0,i}$ (when non-vanishing) is of the same order as the noise level. Indeed our analysis requires an exact asymptotic distributional characterization of the problem (4).

The contributions of this paper are organized as follows:

We state the problem formally, by taking a minimax point of view. Based on this formulation, we prove a general upper bound on the minimax power of tests with a given significance level $\alpha$. We then specialize this bound to the case of standard Gaussian design matrices, showing formally that no test can achieve non-trivial significance $\alpha$, and power $1-\beta$, unless $|\theta_{0,i}|\geq\mu_{{\tiny\rm UB}}=c\sigma/\sqrt{n}$, with $c$ a dimension-independent constant.

We define a hypothesis testing procedure that is well-suited for the case of standard Gaussian designs, $\bm{\Sigma}={\rm I}_{p\times p}$. We prove that this test achieves a ‘nearly-optimal’ power-significance trade-off in a properly defined asymptotic sense. Here ‘nearly optimal’ means that the trade-off has the same form as the previous upper bound, except that $\mu_{{\tiny\rm UB}}$ is replaced by $\mu=C\mu_{{\tiny\rm UB}}$ with $C$ a universal constant. In particular, we provide a positive answer to the open question discussed above.

Our analysis builds on an exact asymptotic characterization of the Lasso estimator, first developed in .

We introduce a generalization of the previous hypothesis testing method to Gaussian designs with general covariance matrix $\bm{\Sigma}$. In this case we cannot establish validity in the regime $n\geq c_{1}s_{0}\log(p/s_{0})$, since a rigorous generalization of the distributional result of is not available.

However: $(1)$ We prove that such a generalized distributional limit holds under the stronger assumption that $n$ is much larger than $s_{0}(\log p)^{2}$ (see Theorem 4.5). $(2)$ We show that this distributional limit can be derived from the powerful replica heuristics in statistical physics for the regime $n\geq c_{1}s_{0}\log(p/s_{0})$. (See Section 4 for further discussion of the validity of this heuristics.)

Conditional on this standard distributional limit holding, we prove that the proposed procedure is nearly optimal in this case as well.

We validate our approach on both synthetic and real data in Sections 3.4, 4.6 and Section 6, comparing it with the methods of . Simulations suggest that the latter are indeed overly conservative in the present setting, resulting in suboptimal statistical power. (As emphasized above, the methods of apply to a broader class of design matrices ${\bm{X}}$.)

Let us stress that the present treatment has two important limitations. First, it is asymptotic: it would be important to develop non-asymptotic bounds. Second, for the case of general designs, it requires to know or estimate the design covariance $\bm{\Sigma}$. In Section 4.5 we discuss a simple approach to this problem for sparse $\bm{\Sigma}$. A full study of this issue is however beyond the scope of the present paper.

After a a preprint version of the present paper became available online, several papers appeared that partially address these limitations. In particular make use of debiased estimators of the form (5), and have much weaker assumptions on the design ${\bm{X}}$. Note however that these papers require a significantly larger sample size, namely $n\geq(s_{0}\log p)^{2}$. Hence, even limiting ourselves to standard designs, the results presented here are not comparable to the ones of , and instead complement them. We refer to Section 5 for further discussion of the relation.

In contrast we work within the Gaussian random design model, and focus on the asymptotics $s_{0},p,n\to\infty$ with $s_{0}/p\to{\varepsilon}\in(0,1)$ and $n/p\to\delta\in(0,1)$. The study of this type of high-dimensional asymptotics was pioneered by Donoho and Tanner , who assumed standard Gaussian designs and focused on exact recovery in absence of noise. The estimation error in presence of noise was characterized in . Further work in the same or related setting includes .

Wainwright also considered the Gaussian design model and established upper and lower thresholds $n_{{\tiny\rm UB}}(p,s_{0};\bm{\Sigma})$, $n_{{\tiny\rm LB}}(p,s_{0};\bm{\Sigma})$ for correct recovery of ${\rm supp}(\bm{\theta}_{0})$ in noise $\sigma>0$, under an additional condition on $\mu\equiv\min_{i\in{\rm supp}(\bm{\theta}_{0})}|\theta_{0,i}|$. The thresholds $n_{{\tiny\rm UB}}(p,s_{0};\bm{\Sigma})$, $n_{{\tiny\rm LB}}(p,s_{0};\bm{\Sigma})$ are of order $s_{0}\log p$ for many covariance structures $\bm{\Sigma}$, provided $\mu\geq C\sqrt{(\log p)/n}$ for some constant $C>0$. Correct support recovery depends, in a crucial way, on the irrepresentability condition of .

Let us stress that the results on support recovery offer limited insight into optimal hypothesis testing procedures. Under the conditions that guarantee exact support recovery, both type I and type II error rates tend to $0$ rapidly as $n,p,s_{0}\to\infty$, thus making it difficult to study the trade-off between statistical significance and power. Here we are interested in triples $n,p,s_{0}$ for which $\alpha$ and $\beta$ stay bounded. As discussed in the previous section, the regime of interest (for standard Gaussian designs) is $c_{1}s_{0}\log(p/s_{0})\leq n\leq c_{2}s_{0}\log(p)$. At the lower end the number of observations $n$ is so small that essentially nothing can be inferred about ${\rm supp}(\bm{\theta}_{0})$ using optimally tuned Lasso estimator, and therefore a nontrivial power $1-\beta>\alpha$ cannot be achieved. At the upper end, the number of samples is sufficient enough to recover ${\rm supp}(\bm{\theta}_{0})$ with high probability, leading to arbitrary small errors $\alpha,\beta$

Let us finally mention that resampling methods provide an alternative path to assess statistical significance. A general framework to implement this idea is provided by the stability selection method of . However, specializing the approach and analysis of to the present context does not provide guarantees superior to , that are more directly comparable to the present work.

2 Notations

Throughout, $\phi(x)=e^{-x^{2}/2}/\sqrt{2\pi}$ is the Gaussian density and $\Phi(x)\equiv\int_{-\infty}^{x}\phi(u){\rm d}u$ is the Gaussian distribution. For two functions $f(n)$ and $g(n)$, with $g(n)\geq 0$, the notation $f(n)=\Omega(g(n))$ means that $f$ is bounded below by $g$ asymptotically, namely, there exists constant $C>0$ and integer $n_{0}>0$, such that $f(n)\geq Cg(n)$ for $n>n_{0}$. Further, $f(n)=O(g(n))$ means that $f$ is bounded above by $g$ asymptotically, namely, for some constants $C<\infty$ and integer $n_{0}>0$, $f(n)\leq C|g(n)|$ for all $n>n_{0}$. Finally $f(n)=\Theta(g(n))$ if both $f(n)=\Omega(g(n))$ and $f(n)=O(g(n))$.

Minimax formulation

In this section we define the hypothesis testing problem, and introduce a minimax criterion for evaluating hypothesis testing procedures. In subsection 2.2 we state our upper bound on the minimax power and, in subsection 2.3, we outline the prof argument, that is based on a reduction to binary hypothesis testing.

We consider the minimax criterion to measure the quality of a testing procedure. In order to define it formally, we first need to establish some notations.

A testing procedure for the family of hypotheses $H_{0,i}$, cf. Eq. (3), is given by a family of measurable functions

As mentioned above, we will measure the quality of a test $T$ in terms of its significance level $\alpha$ (probability of type I errors) and power $1-\beta$ ($\beta$ is the probability of type II errors). A type I error (false rejection of the null) leads one to conclude that a relationship between the response vector ${\bm{y}}$ and a column of the design matrix ${\bm{X}}$ exists when in reality it does not. On the other hand, a type II error (the failure to reject a false null hypothesis) leads one to miss an existing relationship.

Adopting a minimax point of view, we require that these metrics are achieved uniformly over $s_{0}$-sparse vectors. Formally, for $\mu>0$, we let

The minimax power for testing hypothesis $H_{0,i}$ against the alternative $|\theta_{i}|\geq\mu$ is given by the function $1-\beta^{\rm opt}_{i}(\,\cdot\,;\mu):\to$ where, for $\alpha\in$

The following are straightforward yet useful properties.

The optimal power $\alpha\mapsto 1-\beta^{\rm opt}_{i}(\alpha;\mu)$ is non-decreasing. Further, by using a test such that $T_{i,{\bm{X}}}({\bm{y}})=1$ with probability $\alpha$ independently of ${\bm{y}}$, ${\bm{X}}$, we conclude that $1-\beta^{\rm opt}_{i}(\alpha;\mu)\geq\alpha$.

To prove the first property, notice that, for any $\alpha\leq\alpha^{\prime}$ we have $1-\beta_{i}(\alpha;\mu)\leq 1-\beta_{i}(\alpha^{\prime};\mu)$. Indeed $1-\beta_{i}(\alpha^{\prime};\mu)$ is obtained by taking the supremum in Eq. (9) over a family of tests that includes those over which the supremum is taken for $1-\beta_{i}(\alpha;\mu)$.

2 Upper bound on the minimax power

It is easy to check that, for any $\alpha>0$, $u\mapsto G(\alpha,u)$ is continuous and monotone increasing. For $u$ fixed $\alpha\mapsto G(\alpha,u)$ is continuous and monotone increasing. Finally $G(\alpha,0)=\alpha$ and $\lim_{u\to\infty}G(\alpha,u)=1$.

We then have the following upper bound on the optimal power of random Gaussian designs. (We refer to Section 7.3 for the proof.)

The next corollary specializes the above result to the case of standard Gaussian designs. (The proof is immediate and hence we omit it.)

For $i\in[p]$, let $1-\beta^{\rm opt}_{i}(\alpha;\mu)$ be the minimax power of a standard Gaussian design ${\bm{X}}$ with covariance matrix $\bm{\Sigma}={\rm I}_{p\times p}$, cf. Definition 2.1. Then, for any $\xi\in[0,(3/2)\sqrt{n-s_{0}+1}]$ we have

It is instructive to look at the last result from a slightly different point of view. Given $\alpha\in(0,1)$ and $1-\beta\in(\alpha,1)$, how big does the entry $\mu$ need to be so that $1-\beta^{\rm opt}_{i}(\alpha;\mu)\geq 1-\beta$? It follows from Corollary 2.4 that to achieve a pair $(\alpha,\beta)$ as above we require $\mu\geq\mu_{{\tiny\rm UB}}=c\sigma/\sqrt{n}$ for some $c=c(\alpha,\beta)$.

Previous work requires $\mu\geq c\,\max\{\sigma s_{0}\log p/\,n,\sigma/\sqrt{n}\}$ to achieve the same goal although for deterministic designs ${\bm{X}}$ (see Appendix C). This motivates the central question of the present paper (already stated in the introduction): Can hypothesis testing be performed in the ideal regime $\mu\geq c\sigma/\sqrt{n}$?

As further clarified in the next section and in Section 7.1, Theorem 2.3 by an oracle-based argument. Namely, we upper bound the power of any hypothesis testing method, by the power of an oracle that knows, for each coordinates $j\in[p]\setminus j$, whether $\theta_{0,j}\in{\rm supp}(\bm{\theta}_{0})$ or not. In other words the procedure has access to ${\rm supp}(\bm{\theta}_{0})\backslash\{i\}$. At first sight, this oracle appears exceedingly powerful, and hence the bound might be loose. Surprisingly, the bound turns out to be tight, at least in an asymptotic sense, as demonstrated in Section 3.

Let us finally mention that a bound similar to the present one was announced independently –and from a different viewpoint– in .

3 Proof outline

The proof of Theorem 2.3 is based on a simple reduction to the binary hypothesis testing problem. We first introduce the binary testing problem, in which the vector of coefficients $\bm{\theta}$ is chosen randomly according to one of two distributions.

We denote by $1-\beta^{\rm bin}_{i,{\bm{X}}}(\,\cdot\,;Q)$ the optimal power for the binary hypothesis testing problem $\bm{\theta}_{0}\sim Q_{0}$ versus $\bm{\theta}_{0}\sim Q_{1}$, namely:

The reduction is stated in the next lemma.

Let $Q_{0}$, $Q_{1}$ be any two probability measures supported, respectively, on ${\cal R}_{0}$ and ${\cal R}_{1}$ as per Definition 2.5. Then, the minimax power for testing hypothesis $H_{0,i}$ under the random design model, cf. Definition 2.1, is bounded as

Here expectation is taken with respect to the law of ${\bm{X}}$ and the $\inf$ is over all measurable functions ${\bm{X}}\mapsto\alpha_{{\bm{X}}}$.

The binary hypothesis testing problem is characterized in the next lemma by reducing it to a simple regression problem. For $S\subseteq[p]$, we denote by ${\rm{\bf P}}_{S}$ the orthogonal projector on the linear space spanned by the columns $\{\bm{\widetilde{x}}_{i}\}_{i\in S}$. We also let ${\rm{\bf P}}^{\perp}_{S}={\rm I}_{n\times n}-{\rm{\bf P}}_{S}$ be the projector on the orthogonal subspace.

If $|S|<s_{0}$ then for any $\xi>0$ there exists distributions $Q_{0}$, $Q_{1}$ as per Definition 2.5, depending on $i$, $S$, $\mu$ but not on ${\bm{X}}$, such that $\beta^{\rm bin}_{i,{\bm{X}}}(\alpha;Q)\geq\beta^{\rm oracle}_{i,{\bm{X}}}(\alpha;S,\mu)-\xi$.

The proof of this Lemma is presented in Section 7.2.

The proof of Theorem 2.3 follows from Lemmas 2.6 and 2.7, cf. Section 7.3.

Hypothesis testing for standard Gaussian designs

In this section we describe our hypothesis testing procedure (that we refer to as SDL-test) in the case of standard Gaussian designs, see subsection 3.1. In subsection 3.2, we develop asymptotic bounds on the probability of type I and type II errors. The test is shown to nearly achieve the ideal tradeoff between significance level $\alpha$ and power $1-\beta$, using the upper bound stated in the previous section.

Our results are based on a characterization of the high-dimensional behavior of the Lasso estimator, developed in . For the reader’s convenience, and to provide further context, we recall this result in subsection 3.3. Finally, subsection 3.4 discusses some numerical experiments.

Our SDL-test procedure for standard Gaussian designs is described in Table 1.

The key is the construction of the unbiased estimator $\bm{\widehat{\theta}}^{u}$ in step 3. The asymptotic analysis developed in and in the next section establishes that $\bm{\widehat{\theta}}^{u}$ is an asymptotically unbiased estimator of $\bm{\theta}_{0}$, and the empirical distribution of $\{\widehat{\theta}^{u}_{i}-\theta_{0,i}\}_{i=1}^{p}$ is asymptotically normal with variance $\tau^{2}$. Further, the variance $\tau^{2}$ can be consistently estimated using the residual vector $\bm{r}$. These results establish that (in a sense that will be made precise next) the regression model (2) is asymptotically equivalent to a simpler sequence model

with noise having zero mean. In particular, under the null hypothesis $H_{0,i}$, $\widehat{\theta}^{u}_{i}$ is asymptotically gaussian with mean $0$ and variance $\tau^{2}$. This motivates rejecting the null if $|\widehat{\theta}^{u}_{i}|\geq\tau\Phi^{-1}(1-\alpha/2)$.

2 Asymptotic analysis

Note that this definition assumes the coefficients $\theta_{0,i}$ are of order one, while the noise is scaled as $\sigma(p)^{2}=\Theta(n)$. Equivalently, we could have assumed $\theta_{0,i}=\Theta(1/\sqrt{n})$ and $\sigma^{2}(p)=\Theta(1)$: the two settings only differ by a scaling of ${\bm{y}}$. We favor the first scaling as it simplifies somewhat the notation in the following.

As before, we will measure the quality of the proposed test in terms of its significance level (size) $\alpha$ and power $1-\beta$. Recall that $\alpha$ and $\beta$ respectively indicate the type I error (false positive) and type II error (false negative) rates. The following theorem establishes that the $P_{i}$’s are indeed valid p-values, i.e., allow to control type I errors. Throughout $S_{0}(p)=\{i\in[p]:\theta_{0,i}(p)\neq 0\}$ is the support of $\bm{\theta}_{0}(p)$.

A more general form of Theorem 3.2 (cf. Theorem 4.3) is proved in Section 7. We indeed prove the stronger claim that the following holds true almost surely

The result of Theorem 3.2 follows then by taking the expectation of both sides of Eq. (20) and using bounded convergence theorem and exchangeability of the columns of ${\bm{X}}$.

Our next theorem proves a lower bound for the power of the proposed test. In order to obtain a non-trivial result, we need to make suitable assumption on the parameter vectors $\bm{\theta}_{0}=\bm{\theta}_{0}(p)$. In particular, we need to assume that the non-zero entries of $\bm{\theta}_{0}$ are lower bounded in magnitude. If this were not the case, it would be impossible to distinguish arbitrarily small parameters $\theta_{0,i}$ from $\theta_{0,i}=0$. (In Appendix B, we also provide an explicit formula for the regularization parameter $\lambda=\lambda(p_{\Theta_{0}},\sigma,{\varepsilon},\delta)$ that achieves this power.)

There exists a (deterministic) choice of $\lambda=\lambda(\sigma,{\varepsilon})$ such that the following happens.

where $\tau_{*}=\tau_{*}(\sigma_{0},{\varepsilon},\delta)$ is defined as follows

Here, $M({\varepsilon})$ is given by the following parametric expression in terms of the parameter $\kappa\in(0,\infty)$:

Theorem 3.3 is proved in Section 7. We indeed prove the stronger claim that the following holds true almost surely:

The result of Theorem 3.3 follows then by taking the expectation of both sides of Eq. (24) and using exchangeability of the columns of ${\bm{X}}$.

Again, it is convenient to rephrase Theorem 3.3 in terms of the minimum value of $\mu$ for which we can achieve statistical power $1-\beta\in(\alpha,1)$ at significance level $\alpha$. It is known that $M({\varepsilon})=2{\varepsilon}\log(1/{\varepsilon})\,(1+O({\varepsilon}))$ . Hence, for $n\geq 2\,s_{0}\log(p/s_{0})\,(1+O(s_{0}/p))$, we have $\tau_{*}^{2}=O(1)$. Since $\lim_{u\to\infty}G(\alpha,u)=1$, any pre-assigned statistical power can be achieved by taking $\mu\geq C({\varepsilon},\delta)\sigma/\sqrt{n}$ which matches the fundamental limit established in the previous section.

Let us finally comment on the choice of the regularization parameter $\lambda$. Theorem 3.2 holds irrespective of $\lambda$, as long as it is kept fixed in the asymptotic limit. In other words, control of type I errors is fairly insensitive to the regularization parameters. On the other hand, to achieve optimal minimax power, it is necessary to tune $\lambda$ to the correct value. The tuned value of $\lambda=\lambda(p_{\Theta_{0}},\sigma,{\varepsilon},\delta)$ for the standard Gaussian sequence model is provided in Appendix A. Further, the factor $\sigma$ (and hence the need to estimate the noise level) can be omitted if –instead of the Lasso– we use the scaled Lasso . In subsection 3.4, we discuss another way of choosing $\lambda$ that also avoid estimating the noise level.

3 Gaussian limit

Theorems 3.2 and 3.3 are based on an asymptotic distributional characterization of the Lasso estimator developed in . We restate it here for the reader’s convenience.

with ${\sf d}=(1-\|\bm{\widehat{\theta}}\|_{0}/n)^{-1}$.

In particular, this result implies that the empirical distribution of $\{\widehat{\theta}^{u}_{i}-\theta_{0,i}\}_{i=1}^{p}$ is asymptotically normal with variance $\tau_{0}^{2}$. This naturally motivates the use of $|\widehat{\theta}^{u}_{i}|/\tau_{0}$ as a test statistics for hypothesis $H_{0,i}:\,\theta_{0,i}=0$.

The definitions of ${\sf d}$ and $\tau$ in step 2 are also motivated by Theorem 3.4. In particular, ${\sf d}(\bm{y}-{\bm{X}}\bm{\widehat{\theta}})/\sqrt{n}$ is asymptotically normal with variance $\tau_{0}^{2}$. This is used in step 2, where $\tau$ is just the robust median absolute deviation (MAD) estimator (we choose this estimator since it is more resilient to outliers than the sample variance ).

4 Numerical experiments

As an illustration, we generated synthetic data from the linear model (1) with $\bm{w}\sim{\sf N}(0,{\rm I}_{p\times p})$ and the following configurations.

Design matrix: For pairs of values $(n,p)=\{(300,1000),(600,1000),(600,2000)\}$, the design matrix is generated from a realization of $n$ i.i.d. rows $\bm{x}_{i}\sim{\sf N}(0,{\rm I}_{p\times p})$.

Regression parameters: We consider active sets $S_{0}$ with $|S_{0}|=s_{0}\in\{10,20,25,50,100\}$, chosen uniformly at random from the index set $\{1,\cdots,p\}$. We also consider two different strengths of active parameters $\theta_{0,i}=\mu$, for $i\in S_{0}$, with $\mu\in\{0.1,0.15\}$.

We examine the performance of SDL-test (cf. Table 1) at significance levels $\alpha=0.025,0.05$. The experiments are done using glmnet-package in R that fits the entire Lasso path for linear regression models. Let ${\varepsilon}=s_{0}/p$ and $\delta=n/p$. We do not assume ${\varepsilon}$ is known, but rather estimate it as $\bar{{\varepsilon}}=0.25\,\delta/\log(2/\delta)$. The value of $\bar{{\varepsilon}}$ is half the maximum sparsity level ${\varepsilon}$ for the given $\delta$ such that the Lasso estimator can correctly recover the parameter vector if the measurements were noiseless . Provided it makes sense to use Lasso at all, $\bar{{\varepsilon}}$ is thus a reasonable ballpark estimate.

The regularization parameter $\lambda$ is chosen as to satisfy

where $\tau$ and ${\sf d}$ are determined in step 2 of the procedure. Here $\kappa_{*}=\kappa_{*}(\bar{{\varepsilon}})$ is the minimax threshold value for estimation using soft thresholding in the Gaussian sequence model, see and Remark B.1. Note that $\tau$ and ${\sf d}$ in the equation above depend implicitly upon $\lambda$. Since glmnet returns the entire Lasso path, the value of $\lambda$ solving the above equation can be computed by the bisection method.

As mentioned above, the control of type I error is fairly robust for a wide range of values of $\lambda$. However, the above is an educated guess based on the analysis of . We also tried the values of $\lambda$ proposed for instance in on the basis of oracle inequalities.

Figure 2 shows the results of SDL-test and the method of for parameter values $p=1000,n=600,s_{0}=25,\mu=0.15$, and significance levels $\alpha\in\{0.025,0.05\}$. Each point in the plot corresponds to one realization of this configuration (there are a total of $10$ realizations). We also depict the theoretical curve $(\alpha,G(\alpha,\mu_{0}/\tau_{*}))$, predicted by Theorem 3.3. The empirical results are in good agreement with the asymptotic prediction.

We compare SDL-test with the ridge-based regression method and the low dimensional projection estimator (LDPE ) . Table 2 summarizes the results for a few configurations $(p,n,s_{0},\mu)$, and $\alpha=0.05$. Simulation results for a larger number of configurations and $\alpha=0.05,0.025$ are reported in Tables 8 and 9 in Appendix E.

As demonstrated by these results, LDPE and the ridge-based regression are both overly conservative. Namely, they achieve smaller type I error than the prescribed level $\alpha$ and this comes at the cost of a smaller statistical power than our testing procedure. This is to be expected since the approach of and cover a broader class of design matrices ${\bm{X}}$, and are not tailored to random designs.

Note that being overly conservative is a drawback, when this comes at the expense of statistical power. The data analysts should be able to decide the level of statistical significance $\alpha$, and obtain optimal statistical power at that level.

The reader might wonder whether the loss in statistical power of methods in and is entirely due to the fact that these methods achieve a smaller number of false positives than requested. In Fig. 3, we run SDL-test , ridge-based regression , and LDPE for $\alpha\in\{0.01,0.02,\cdots,0.1\}$ and for $10$ realizations of the problem per each value of $\alpha$. We plot the average type I error and the average power of each method versus $\alpha$. As we see even for the same empirical fraction of type I errors, SDL-test results in a higher statistical power.

Hypothesis testing for nonstandard Gaussian designs

In this section, we generalize our testing procedure to nonstandard Gaussian design models where the rows of the design matrix ${\bm{X}}$ are drawn independently from distribution ${\sf N}(0,\bm{\Sigma})$.

We first describe the generalized SDL-test procedure in subsection 4.1 under the assumption that $\bm{\Sigma}$ is known. In subsection 4.2, we show that this generalization can be justified from a certain generalization of the Gaussian limit theorem 3.4 to nonstandard Gaussian designs.

Establishing such a generalization of Theorem 3.4 appears extremely challenging. We nevertheless show that such a limit theorem follows from the replica method of statistical physics in section 4.4. We also show that a version of this limit theorem is relatively straightforward in the regime $s_{0}=o(n/(\log p)^{2})$.

Finally, in Section 4.5 we discuss a procedure for estimating the covariance $\bm{\Sigma}$ (cf. Subroutine in Table 4). Appendix F proposes an alternative implementation that does not estimate $\bm{\Sigma}$ but instead bounds the effect of unknown $\bm{\Sigma}$.

The hypothesis testing procedure SDL-testfor general Gaussian designs is defined in Table 3.

The basic intuition of this generalization is that $(\widehat{\theta}^{u}_{i}-\widehat{\theta}_{0,i})/(\tau[(\bm{\Sigma}^{-1})_{ii}]^{1/2})$ is expected to be asymptotically ${\sf N}(0,1)$, whence the definition of (two-sided) p-values $P_{i}$ follows as in step 4. Parameters ${\sf d}$ and $\tau$ in step 2 are defined in the same manner to the standard Gaussian designs.

2 Asymptotic analysis

Let $v_{i}=(\theta_{0,i},(\widehat{\theta}^{u}_{i}-\theta_{0,i})/\tau,(\bm{\Sigma}^{-1})_{ii})$, for $1\leq i\leq p$, and $\nu^{(p)}$ be the empirical distribution of $\{v_{i}\}_{i=1}^{p}$ defined as

We will next show that the SDL-test procedure is appropriate for any random design model for which the standard distributional limit holds. Our first theorem is a generalization of Theorem 3.2 to this setting.

The proof of Theorem 4.3 is deferred to Section 7. In the proof, we show the stronger result that the following holds true almost surely

The result of Theorem 4.3 follows then by taking the expectation of both sides of Eq. (31) and using bounded convergence theorem.

The following theorem characterizes the power of SDL-test for general $\bm{\Sigma}$, and under the assumption that a standard distributional limit holds .

Theorem 4.4 is proved in Section 7. We indeed prove the stronger result that the following holds true almost surely

We also notice that in contrast to Theorem 3.3, where $\tau_{*}$ has an explicit formula that leads to an analytical lower bound for the power (for a suitable choice of $\lambda$), in Theorem 4.4, $\tau$ depends upon $\lambda$ implicitly and can be estimated from the data as in step 3 of SDL-test procedure. The result of Theorem 4.4 holds for any value of $\lambda$.

3 Gaussian limit for n≫s0​(log⁡p)2n\gg s_{0}(\log p)^{2}

In the following theorem we show that if sample size $n$ asymptotically dominates $s_{0}(\log p)^{2}$, then the standard distributional limit can be established rigorously.

$n(p)\leq p$, and $s_{0}(\log p)^{2}/n(p)\to 0$;

There exist constants $c_{\rm min},c_{\rm max}>0$ such that the eigenvalues of $\bm{\Sigma}$ lie in the interval $[c_{\rm min},c_{\rm max}]$: $c_{\rm min}\leq\lambda_{\rm min}(\bm{\Sigma})\leq\lambda_{\rm max}(\bm{\Sigma})\leq c_{\rm max}$;

The empirical distribution of $\{(\bm{\Sigma}^{-1})_{ii})\}_{1\leq i\leq p}$ converges weakly to the probability distribution of the random variable $\Upsilon$;

The regularization parameter is $\lambda=C_{*}\sigma\sqrt{(\log p)/n}$ for $C_{*}=C_{*}(c_{\rm min},c_{\rm max})$ a sufficiently large constant.

Then the sequence has a standard distributional limit with ${\sf d}=(1-\|\bm{\widehat{\theta}}(\lambda)\|_{0}/n)^{-1}$ and $\tau=\sigma_{0}$. Alternatively, $\tau$ can be taken to be a solution of Eq. (37) below.

Theorem 4.5 is proved in Section 7.7. The proof uses techniques from our conference paper .

Notice that this result does allow to control type I errors using Theorem 4.3, but does not allow to lower bound the power, using Theorem 4.4, since $|S_{0}(p)|/p\to 0$. A lower bound on the power under the same assumptions presented in this section can be found in . In the present paper we focus instead on the case $|S_{0}(p)|/p$ bounded away from $0$.

4 Gaussian limit via the replica heuristics for smaller sample size nn

As mentioned above, the standard distributional limit follows from Theorem 3.4 for $\Sigma={\rm I}_{p\times p}$. Even in this simple case, the proof is rather challenging . Partial generalization to non-gaussian designs and other convex problems appeared recently in and , each requiring over 50 pages of proofs.

where the the limit exists by the above assumptions on the convergence of ${\mathfrak{E}}^{(p)}(a,b)$. Then, the parameters $\tau$ and ${\sf d}$ of the standard distributional limit are obtained by setting ${\sf d}=(1-\bm{\widehat{\theta}}/n)^{-1}$ and solving the following with respect to $\tau^{2}$:

In other words, the replica method indicates that the standard distributional limit holds for a large class of non-diagonal covariance structures $\bm{\Sigma}$. It is worth stressing that convergence assumption for the sequence ${\mathfrak{E}}^{(p)}(a,b)$ is quite mild, and is satisfied by a large family of covariance matrices. For instance, it can be proved that it holds for block-diagonal matrices $\bm{\Sigma}$ as long as the blocks have bounded length and the blocks empirical distribution converges.

The replica method is a non-rigorous but highly sophisticated calculation procedure that has proved successful in a number of very difficult problems in probability theory and probabilistic combinatorics. Attempts to make the replica method rigorous have been pursued over the last 30 years by some world-leading mathematicians . This effort achieved spectacular successes, but so far does not provide tools to prove the above replica claim. In particular, the rigorous work mainly focuses on ‘i.i.d. randomness’, corresponding to the case covered by Theorem 3.4.

Over the last ten years, the replica method has been used to derive a number of fascinating results in information theory and communications theory, see e.g. . More recently, several groups used it successfully in the analysis of high-dimensional sparse regression under standard Gaussian designs . The rigorous analysis of ours and other groups subsequently confirmed these heuristic calculations in several cases.

There is a fundamental reason that makes establishing the standard distributional limit a challenging task. This requires in fact to characterize the distribution of the estimator (4) in a regime where the standard deviation of $\widehat{\theta}_{i}$ is of the same order as its mean. Further, $\widehat{\theta}_{i}$ does not converge to the true value $\theta_{0,i}$, hence making perturbative arguments ineffective.

The analysis becomes easier for a larger number of samples. In Theorem 4.5 below we will show that (a suitable version of) the standard distributional holds for $n$ asymptotically larger than $s_{0}(\log p)^{2}$. This uses methods from our companion paper .

5 Covariance estimation

So far we assumed that the design covariance $\bm{\Sigma}$ is known. This setting is relevant for semi-supervised learning applications, where the data analyst has access to a large number $N\gg p$ of ‘unlabeled examples’. These are i.i.d. feature vectors $\bm{u}_{1}$, $\bm{u}_{2}$,…$\bm{u}_{N}$ with $\bm{u}_{1}\sim{\sf N}(0,\bm{\Sigma})$ distributed as $\bm{x}_{1}$, for which the response variable $y_{i}$ is not available. In this case $\bm{\Sigma}$ can be estimated accurately by $N^{-1}\sum_{i=1}^{n}\bm{u}_{i}\bm{u}_{i}^{{\sf T}}$. We refer to for further background on such applications.

In other applications, $\bm{\Sigma}$ is unknown and no additional data is available. In this case we proceed as follows:

We estimate $\Sigma$ from the design matrix $\bm{X}$ (equivalently, from the feature vectors $\bm{x}_{1}$, $\bm{x}_{2}$, …$\bm{x}_{n}$). We let $\bm{\widehat{\Sigma}}$ denote the resulting estimate.

We use $\bm{\widehat{\Sigma}}$ instead of $\bm{\Sigma}$ in step 3 of our hypothesis testing procedure.

The problem of estimating covariance matrices in high-dimensional setting has attracted considerable attention in the past. Several estimation methods provide a consistent estimate $\bm{\widehat{\Sigma}}$, under suitable structural assumptions on $\bm{\Sigma}$. For instance if $\bm{\Sigma}^{-1}$ is sparse, one can apply the graphical model method of , the regression approach of , or CLIME estimator , to name a few.

Since the covariance estimation problem is not the focus of our paper, we will test the above approach using a very simple covariance estimation method. Namely, we assume that $\bm{\Sigma}$ is sparse and estimate it by thresholding the empirical covariance. A detailed description of this estimator is given in Table 4. We refer to for a theoretical analysis of this type of methods. Note that the Lasso is unlikely to perform well if the columns of ${\bm{X}}$ are highly correlated and hence the assumption of sparse $\bm{\Sigma}$ is very natural. On the other hand, we would like to emphasize that this covariance thresholding estimation is only one among many possible approaches.

As an additional contribution, in Appendix F we describe an alternative covariance-free procedure that only uses bounds on $\bm{\Sigma}$ where the bounds are estimated from the data.

In our numerical experiments, we use the estimated covariance returned by Subroutine. As shown in the next section, computed p-values appear to be fairly robust with respect to errors in the estimation of $\bm{\Sigma}$. It would be interesting to develop a rigorous analysis of SDL-test that accounts for the covariance estimation error.

6 Numerical experiments

Elements below the diagonal are given by the symmetry condition $\Sigma_{kj}=\Sigma_{jk}$. (Notice that this is a circulant matrix.)

In Fig. 4(a), we compare SDL-test with the ridge-based regression method proposed in . While the type I errors of SDL-test are in good match with the chosen significance level $\alpha$, the method of is conservative. As in the case of standard Gaussian designs, this results in significantly smaller type I errors than $\alpha$ and smaller average power in return. Also, in Fig. 5, we run SDL-test , ridge-based regression , and LDPE for $\alpha\in\{0.01,0.02,\cdots,0.1\}$ and for $10$ realizations of the problem per each value of $\alpha$. We plot the average type I error and the average power of each method versus $\alpha$. As we see, similar to the case of standard Gaussian designs, even for the same empirical fraction of type I errors, SDL-test results in a higher statistical power.

Table 5 summarizes the performances of the these methods for a few configurations $(p,n,s_{0},\mu)$, and $\alpha=0.05$. Simulation results for a larger number of configurations and $\alpha=0.05,0.025$ are reported in Tables 10 and 11 in Appendix E.

Let $\bm{z}=(z_{i})_{i=1}^{p}$ denote the vector with entries $z_{i}\equiv\widehat{\theta}^{u}_{i}/(\tau[(\bm{\Sigma}^{-1})_{ii}]^{1/2})$. In Fig. 4(b) we plot the normalized histograms of $\bm{z}_{S_{0}}$ (in red) and $\bm{z}_{S_{0}^{c}}$ (in white), where $\bm{z}_{S_{0}}$ and $\bm{z}_{S_{0}^{c}}$ respectively denote the restrictions of $\bm{z}$ to the active set $S_{0}$ and the inactive set $S_{0}^{c}$. The plot clearly exhibits the fact that $\bm{z}_{S^{c}_{0}}$ has (asymptotically) standard normal distribution and the histogram of $\bm{z}_{S_{0}}$ appears as a distinguishable bump. This is the core intuition in defining SDL-test.

Discussion

In this section we compare our contribution with related work in order to put it in proper perspective. We first compare it with other recent debiasing methods in subsection 5.1. In subsection 5.2 we then discuss the role of of the factor ${\sf d}$ in our definition of $\bm{\widehat{\theta}}^{u}$: this is an important difference with respect to the methods of . We finally contrast the Gaussian limit in Theorem 3.4 and Le Cam’s local asymptotic normality theory, that plays a pivotal role in classical statistics.

As explained several times in the previous sections, the key step in our procedure is to correct the Lasso estimator through a debiasing procedure. For the reader’s convenience, we copy here the definition of the latter:

The approach of is similar in that it is based on debiased estimator of the form

where $\bm{M}$ is computed from the design matrix ${\bm{X}}$. The authors of propose to compute $\bm{M}$ by doing sparse regression of each column of ${\bm{X}}$ onto the others.

After a first version of the present paper became available as an online preprint, de Geer, Bühlmann and Ritov studied an approach similar to (and to ours) in a random design setting. They provide guarantees under the assumptions that $\bm{\Sigma}^{-1}$ is sparse and that the sample size $n$ asymptotically dominates $(s_{0}\log p)^{2}$. The authors also establish asymptotic optimality of their method in terms of semiparametric efficiency. The semiparametric setting is also at the center of .

A further development over the approaches of was proposed by the present authors in . This paper constructs the matrix $\bm{M}$ by solving an optimization problem that controls the bias of $\bm{\widehat{\theta}^{*}}$ and minimize its variance meanwhile. This method does not require any sparsity assumption on $\bm{\Sigma}$ or $\bm{\Sigma}^{-1}$, but still requires sample size $n$ to asymptotically dominate $(s_{0}\log p)^{2}$.

It is interesting to compare and contrast the results of , with the contribution of the present paper. (Let us emphasize that appeared after submission of the present work.)

The approach of guarantees control of type I error, and optimality for non-Gaussian designs. (Both of require however sparsity of $\bm{\Sigma}^{-1}$.)

In contrast, our results are fully rigorous only in the special case $\bm{\Sigma}={\rm I}$.

Neither of the papers requires knowledge of covariance $\bm{\Sigma}$. The method in estimates $\bm{\Sigma}^{-1}$ assuming that it is sparse, however the method does not require such estimation.

In contrast, our generalization to arbitrary Gaussian designs postulates knowledge of $\bm{\Sigma}$. (Further this generalization relies on the standard distributional limit assumption.)

The work of focuses on random designs, but requires $n$ much larger than $(s_{0}\log p)^{2}$. This is roughly the square of the number of samples needed for consistent estimation.

In contrast, we achieve similar power, and confidence intervals with optimal sample size $n=O(s_{0}\log(p/s_{0}))$.

In summary, the present work is complementary to the one in in that it provides a sharper characterization, within a more restrictive setting. Together, these papers provide support for the use of debiasing methods of the form (42).

2 Role of the factor 𝖽{\sf d}

It is worth stressing one subtle, yet interesting, difference between the methods of of and the one of the present paper. In both cases, a debiased estimator is constructed using Eq. (42). However:

The approach of sets $\bm{M}$ to be an estimate of $(\bm{\Sigma}^{-1})$. In the idealized situation where $\bm{\Sigma}$ is known, this construction reduces to setting $\bm{M}=\bm{\Sigma}^{-1}$.

In contrast, our prescription (41) amounts to setting $\bm{M}={\sf d}\,\bm{\Sigma}^{-1}$, with ${\sf d}=(1-\|\bm{\widehat{\theta}}\|_{0}/n)^{-1}$. In other words, we choose $\bm{M}$ as a scaled version of the inverse covariance.

The mathematical reason for the specific scaling factor is elucidated by the proof of Theorem 3.4 in . Here we limit ourselves to illustrating through numerical simulations that this factor is indeed crucial to ensure the normality of $(\widehat{\theta}_{i}^{u}-\theta_{0,i})$ in the regime $n=\Theta(s_{0}\log(p/s_{0}))$.

We consider the same setup as in Section 4.6 where the rows of the design matrix are generated independently from ${\sf N}(0,\bm{\Sigma})$ with $\bm{\Sigma}_{jk}$ given by (40) for $j\leq k$. We fix undersampling ratio $\delta=n/p$ and sparsity level ${\varepsilon}=s_{0}/p$ and consider values $p\in\{250,500,750,\cdots,3500\}$. We also take active sets $S_{0}$ with $|S_{0}|=s_{0}$ chosen uniformly at random from the index set $\{1,\cdots,p\}$ and set $\theta_{0,i}=0.15$ for $i\in S$.

The goal is to illustrate the effect of the scaling factor ${\sf d}$ on the empirical distribution of $(\widehat{\theta}^{u}_{i}-\theta_{0,i})$, for large $n,p,s_{0}$. As we will see, the effect becomes more pronounced as the ratio $n/s_{0}=\delta/{\varepsilon}$ (i.e. the number of samples per non-zero coefficient) becomes smaller. As above, we use $\bm{\widehat{\theta}}^{u}$ for the unbiased estimator developed in this paper (which amounts to Eq. (42) with $\bm{M}={\sf d}\bm{\Sigma}^{-1}$). We will use $\bm{\widehat{\theta}}^{{\sf d}=1}$ for the ‘ideal’ unbiased estimator corresponding to the proposal of (which amounts to Eq. (42) with $\bm{M}=\bm{\Sigma}^{-1}$).

In Fig. 7, we plot the histogram of $\bm{v}$ for $p=3000$ and using both $\bm{\widehat{\theta}}^{u}$ and $\bm{\widehat{\theta}}^{{\sf d}=1}$. Again, the plots clearly demonstrate importance of ${\sf d}$ in obtaining a Gaussian behavior.

(${\varepsilon}=0.02,\delta=0.6$). Figures 6(b) and 8 show similar plots for this case. As we see, the effect of ${\sf d}$ becomes less noticeable here. The reason is that we expect $\|\bm{\widehat{\theta}}\|_{0}/n=O(s_{0}/n)$, and ${\sf d}=(1-\|\bm{\widehat{\theta}}\|_{0}/n)^{-1}=1+O(s_{0}/n)\approx 1$ for $s_{0}$ much smaller than $n$.

3 Comparison with Local Asymptotic Normality

Our approach is based on an asymptotic distributional characterization of the Lasso estimator, cf. Theorem 3.4. Simplifying, the Lasso estimator is in correspondence with a debiased estimator $\bm{\widehat{\theta}}^{u}$ that is asymptotically normal in the sense of finite-dimensional distributions. This is analogous to what happens in classical statistics, where local asymptotic normality (LAN) can be used to characterize an estimator distribution, and hence derive test statistics .

This analogy is only superficial, and the mathematical phenomenon underlying Theorem 3.4 is altogether different from the one in local asymptotic normality. We refer to for a more complete understanding, and only mention a few points:

LAN theory holds in the low-dimensional limit, where the number of parameters $p$ is much smaller than the number of samples $n$. Even more, the focus is on $p$ fixed, and $n\to\infty$.

In contrast, the Gaussian limit in Theorem 3.4 holds with $p$ proportional to $n$.

The starting point of LAN theory is low-dimensional consistency, namely $\bm{\widehat{\theta}}\to\bm{\theta}_{0}$ as $n\to\infty$. As a consequence, the distribution of $\bm{\widehat{\theta}}$ can be characterized by a local approximation around $\bm{\theta}_{0}$.

In contrast, in the high-dimensional asymptotic regime of Theorem 3.4, the mean square error per coordinate $\|\bm{\widehat{\theta}}-\bm{\theta}_{0}\|_{2}^{2}/p$ remains bounded away from zero . As a consequence, normality does not follow from local approximation.

Indeed, in the present case, the Lasso estimator (which is of course a special case of M-estimator) $\bm{\widehat{\theta}}$ is not normal. Only the debiased estimator $\bm{\widehat{\theta}}^{u}$ is asymptotically normal. Further, while LAN theory holds quite generally in the classical asymptotics, the present theory is more sensitive to the properties of the design matrix ${\bm{X}}$.

Real data application

We tested our method on the UCI communities and crimes dataset . This concerns the prediction of the rate of violent crime in different communities within US, based on other demographic attributes of the communities. The dataset consists of a response variable along with 122 predictive attributes for 1994 communities. Covariates are quantitative, including e.g., the fraction of urban population or the median family income. We consider a linear model as in (2) and hypotheses $H_{0,i}$. Rejection of $H_{0,i}$ indicates that the $i$-th attribute is significant in predicting the response variable.

In order to evaluate various hypothesis testing procedures, we need to know the true significant variables. To this end, we let $\bm{\theta}_{0}=({\bm{X}}_{\rm tot}^{{\sf T}}{\bm{X}}_{\rm tot})^{-1}{\bm{X}}_{\rm tot}^{{\sf T}}{\bm{y}}$ be the least-square estimator, using the whole data set. Figure 9 shows the the entries of $\bm{\theta}_{0}$. Clearly, only a few entries have non negligible values which correspond to the significant attributes. In computing type I errors and powers, we take the elements in $\bm{\theta}_{0}$ with magnitude larger than $0.04$ as active and the others as inactive.

In order to validate our approach in the high-dimensional $p>n$ regime, we take random subsamples of the communities (hence subsamples of the rows of ${\bm{X}}_{\rm tot}$) of size $n=84$. We compare SDL-test with the method of , over $20$ realizations and significance levels $\alpha=0.01,0.025,0.05$. The fraction of type I errors and statistical power is computed by comparing to $\bm{\theta}_{0}$. Table 6 summarizes the results. As the reader can see, Buhlmann’s method is very conservative yielding to no type-I errors and but much smaller power than SDL-test.

In table 7, we report the relevant features obtained from the whole dataset as described above, corresponding to the nonzero entries in $\bm{\theta}_{0}$. We also report the features identified as relevant by SDL-test and those identified as relevant by Ridge-based regression method, from one random subsample of communities of size $n=84$. Features description is available in .

Finally, in Fig. 10 we plot the normalized histograms of $\bm{v}_{S_{0}}$ (in red) and $\bm{V}_{S_{0}^{c}}$ (in white). Recall that $\bm{v}=(v_{i})_{i=1}^{p}$ denotes the vector with $v_{i}\equiv\widehat{\theta}^{u}_{i}/(\tau[(\bm{\Sigma}^{-1})_{ii}]^{1/2})$. Further, $\bm{v}_{S_{0}}$ and $\bm{v}_{S_{0}^{c}}$ respectively denote the restrictions of $\bm{v}$ to the active set $S_{0}$ and the inactive set $S_{0}^{c}$. This plot demonstrates that $\bm{v}_{S^{c}_{0}}$ has roughly standard normal distribution as predicted by the theory.

Proofs

We now take expectation of these inequalities with respect to $\bm{\theta}\sim Q_{0}$ (in the first case) and $\bm{\theta}\sim Q_{1}$ (in the second case) and we get, with the notation introduced in the Definition 2.5,

The thesis follows since $\xi>0$ is arbitrary.

2 Proof of Lemma 2.7

it is a straightforward calculation to drive the power of this test as

where the function $G(\alpha,u)$ is defined as per Eq. (10). Next we show that the power of this test converges to $1-\beta^{\rm oracle}_{i,{\bm{X}}}(\alpha;S,\mu)$ as $a\to\infty$. Hence the claim is proved by taking $a\geq a(\xi)$ for some $a(\xi)$ large enough.

where the second step follows from matrix inversion lemma. Clearly, as $a\to\infty$, the right hand side of the above equation converges to $(1/\sigma)\,{\rm{\bf P}}^{\perp}_{S}$. Therefore, the power converges to $1-\beta^{\rm oracle}_{i,{\bm{X}}}(\alpha;S,\mu)=G(\alpha,\mu\sigma^{-1}\|{\rm{\bf P}}^{\perp}_{S}\bm{\widetilde{x}}_{i}\|)$.

3 Proof of Theorem 2.3

Let $u_{{\bm{X}}}\equiv\mu\|{\rm{\bf P}}^{\perp}_{S}\bm{\widetilde{x}}_{i}\|_{2}/\sigma$. By Lemma 2.6 and 2.7, we have,

with the $\sup$ taken over measurable functions ${\bm{X}}\mapsto\alpha_{{\bm{X}}}$, and $G(\alpha,u)$ defined as per Eq. (10).

Since $\bm{\widetilde{x}}_{i}$ and ${\bm{X}}_{S}$ are jointly Gaussian, we have

with $\bm{z}_{i}\sim{\sf N}(0,{\rm I}_{n\times n})$ independent of ${\bm{X}}_{S}$. It follows that

4 Proof of Theorem 3.3

In addition, since $\mu_{0}\sigma_{0}/2$ is a continuity point of the distribution of $\Theta_{0}$, we have

Now, we take the expectation of both sides of Eq. (56) with respect to the law of random design ${\bm{X}}$ and random noise $\bm{w}$. Changing the order of limit and expectation by applying dominated convergence theorem and using linearity of expectation, we obtain

5 Proof of Theorem 4.3

converges weakly to ${\sf N}(0,1)$. Hence,

Applying the same argument as in the proof of Theorem 3.3, we obtain the following by taking the expectation of both sides of the above equation

6 Proof of Theorem 4.4

Similar to the proof of Theorem 3.3, by taking the expectation of both sides of the above inequality we get

7 Proof of Theorem 4.5

In order to prove the claim, we will establish the following (corresponding to the the case $\Theta_{0}=0$ of Definition 4.1):

If $\tau$ solves Eq. (37), then $\tau^{2}\to\sigma^{2}_{0}$ as $p\to\infty$.

Recalling $\bm{r}\equiv{\sf d}(\bm{y}-\bm{X}\bm{\widehat{\theta}})/\sqrt{n}$, the empirical distribution of $\{r_{i}\}_{1\leq i\leq n}$ converges weakly to ${\sf N}(0,\sigma_{0}^{2})$.

We will prove these three claims after some preliminary remarks. First notice that, by [59, Theorem 6] (and using assumptions $(i)$ and $(iii)$) ${\bm{X}}$ satisfies the restricted eigenvalue property RE$(s_{0},3s_{0},3)$ of with a $p$-independent constant $\kappa=\kappa(c_{\rm min},c_{\rm max})>0$, almost surely for all $p$ large enough. (Indeed Theorem 6 of ensures that this holds with probability at least $1-e^{-\Omega(n(p))}$, and hence almost surely for all $p$ large enough by Borel-Cantelli lemma.)

We can therefore apply [18, Theorem 7.2] to conclude that there exists a constant $C_{0}$ such that, almost surely for all $p$ large enough, we have

(Here we used $\sigma^{2}\leq 2n\sigma_{0}^{2}$ for all $p$ large enough.) In particular, from Eq. (67) and assumption $(i)$, it follows that $\lim_{p\to\infty}\|\widehat{\theta}\|_{0}/n=0$ and hence, almost surely,

By Eq. (68), we can assume ${\sf d}\in(1/2,2)$ for all $p$ large enough. By Eq. (37) it is sufficient to show that ${\sf E}_{1}(\tau^{2},b)\to 0$ uniformly for $b\in[1/2,2]$, $\tau\in[0,M\sigma_{0}]$, for some $M\geq 2$. Since $\|\bm{\theta}\|_{0}/p\to 0$, and by dominated convergence, we have

Therefore, substituting ${\bm{y}}=\tau\bm{\Sigma}^{-1/2}\bm{z}$, we have

The random variables $(\bm{\Sigma}^{1/2}{\bm{z}})_{i}$ are ${\sf N}(0,\Sigma_{ii})$. Therefore by union bound, since $\Sigma_{ii}\leq c_{\rm max}$, for $Z\sim{\sf N}(0,1)$, we have

7.2 Claim 2

Let $\bm{z}=\bm{\Sigma}^{-1}{\bm{X}}^{\sf T}\bm{w}/n$. Conditional on ${\bm{X}}$, we have

Using the assumption $\sigma^{2}/n\to\sigma_{0}^{2}$ and employing [33, Lemma 7.2], we have, almost surely,

Consequently, we have, for almost every sequence of matrices ${\bm{X}}$, letting $Z\sim{\sf N}(0,1)$ independent of ${\bm{X}}$

(Here, the first identity follows from Eq. (74), the second from Eq. (75) and the Lipschitz continuity of $\psi$, and the last from assumption $(iv)$, together with the fact that $\psi$ is bounded Lipschitz.)

Next, applying Gaussian isoperimetry to the conditional measure of $\bm{z}$ given ${\bm{X}}$ (noting that $\|\bm{C}\|_{2}\leq C_{1}$ almost surely for all $n$ large enough and some constant $C_{1}<\infty$), and to the Lipschitz function $\bm{z}\mapsto\Psi({\bm{z}})\equiv p^{-1}\sum_{i=1}^{p}\psi(0,z_{i},(\bm{\Sigma}^{-1})_{ii})$, we have

almost surely for all $n$ large enough. Using Borel-Cantelli lemma, we conclude that, almost surely

Substituting ${\bm{y}}={\bm{X}}\theta_{0}+\bm{w}$ in definition of $\bm{\widehat{\theta}}^{u}$, we get

where we recall that $\bm{\widehat{\Sigma}}\equiv({\bm{X}}^{{\sf T}}{\bm{X}})/n$ and we defined

The proof is therefore concluded if we can show that, almost surely,

In order to simplify the notation, and since the last argument plays no role, we let $\psi_{i}(x,y)\equiv\psi(x,y,(\bm{\Sigma}^{-1})_{ii})$. Without loss of generality we will assume that $\|\psi_{i}\|_{\infty}\leq 1$, and that the Lipschitz modulus of $\psi_{i}$ is at most one.