A Flexible Framework for Hypothesis Testing in High-dimensions

TL;DR

This paper introduces a versatile framework for hypothesis testing in high-dimensional linear regression, enabling valid inference for complex hypotheses when the number of parameters exceeds samples, with proven error control and high power.

Contribution

It develops a general testing framework for high-dimensional models that handles various hypotheses and provides minimax optimal confidence intervals for linear functionals.

Findings

Controls false positive rate near nominal level

Demonstrates high statistical power in experiments

Provides minimax rate optimal confidence intervals

Abstract

Hypothesis testing in the linear regression model is a fundamental statistical problem. We consider linear regression in the high-dimensional regime where the number of parameters exceeds the number of samples ($p> n$). In order to make informative inference, we assume that the model is approximately sparse, that is the effect of covariates on the response can be well approximated by conditioning on a relatively small number of covariates whose identities are unknown. We develop a framework for testing very general hypotheses regarding the model parameters. Our framework encompasses testing whether the parameter lies in a convex cone, testing the signal strength, and testing arbitrary functionals of the parameter. We show that the proposed procedure controls the type I error, and also analyze the power of the procedure. Our numerical experiments confirm our theoretical findings and demonstrate that we control false positive rate (type I error) near the nominal level, and have high power. By duality between hypotheses testing and confidence intervals, the proposed framework can be used to obtain valid confidence intervals for various functionals of the model parameters. For linear functionals, the length of confidence intervals is shown to be minimax rate optimal.

Tables

Table 1. (a) Type I error ( % percent \% )

$c\backslash \rho$	$0.2$	$0.4$	$0.6$	$0.8$
$0.02$	0.00	0.004	1.33	2.33
$0.04$	0.33	1.66	2.33	3.00
$0.06$	1.66	2.00	3.00	3.66
$0.08$	3.33	4.33	3.66	4.66
$0.1$	3.00	4.00	4.66	4.33

Table 2. (a) Type I error ( % percent \% )

$c\backslash \rho$	$0.2$	$0.4$	$0.6$	$0.8$
$0.02$	0.00	0.004	1.33	2.33
$0.04$	0.33	1.66	2.33	3.00
$0.06$	1.66	2.00	3.00	3.66
$0.08$	3.33	4.33	3.66	4.66
$0.1$	3.00	4.00	4.66	4.33

Table 3. (b) Statistical power ( % percent \% )

$c\backslash \rho$	$0.2$	$0.4$	$0.6$	$0.8$
$0.2$	8.00	10.66	18.66	14.33
$0.3$	17.33	24.66	28.66	35.33
$0.4$	86.00	93.33	92.66	84.66
$0.5$	90.00	88.00	97.33	86.66
$0.6$	100.00	88.33	100.00	100.00

Table 4. (a) Type I error ( % percent \% )

$b\backslash \rho$	$0.2$	$0.4$	$0.6$	$0.8$
$1$	2.00	2.00	2.00	3.33
$0.8$	0.66	2.33	2.33	2.66
$0.6$	3.00	3.66	1.00	2.66
$0.4$	2.66	2.33	1.33	2.00
$0.2$	2.33	1.66	2.33	3.66

Table 5. (a) Type I error ( % percent \% )

$b\backslash \rho$	$0.2$	$0.4$	$0.6$	$0.8$
$1$	2.00	2.00	2.00	3.33
$0.8$	0.66	2.33	2.33	2.66
$0.6$	3.00	3.66	1.00	2.66
$0.4$	2.66	2.33	1.33	2.00
$0.2$	2.33	1.66	2.33	3.66

Table 6. (b) Statistical power ( % percent \% )

$b\backslash \rho$	$0.2$	$0.4$	$0.6$	$0.8$
$-0.2$	35.33	68.00	78.00	80.00
$-0.4$	99.33	100.00	100.00	100.00
$-0.6$	100.00	100.00	100.00	100.00
$-0.8$	100.00	100.00	100.00	100.00
$-1$	100.00	100.00	100.00	100.00

Table 7. Table 3: Coverage rate of the confidence intervals for ξ 𝖳 ​ θ 0 superscript 𝜉 𝖳 subscript 𝜃 0 \xi^{\sf T}\theta_{0} and ‖ θ 0 ‖ 2 2 superscript subscript norm subscript 𝜃 0 2 2 \|\theta_{0}\|_{2}^{2} computed as in ( 60 ) for the real data experiment and at various noise levels σ 𝜎 \sigma .

$\sigma$	$1$	$5$	$10$
${\xi }^{𝖳}{\theta }_0$	0.96	0.94	0.93
${\Vert {\theta }_0\Vert }_2^2$	0.95	0.93	0.94

Equations

$$y_i={\theta }_0^{\mathsf{T}}{x}_i+w_i,w_i\sim \mathsf{N}(0,{\sigma }^2).$$

$$y=X{\theta }_0+w,w\sim \mathsf{N}(0,{\sigma }^2{\mathrm{I}}_{n\times n}).$$

$$H_0:{\theta }_0\in {\Omega }_0\text{ versus }{H}_A:{\theta }_0\notin {\Omega }_0,,$$

$${\hat{\theta }}^n(y,X;\lambda )\equiv \arg \underset{\theta \in {\mathbb{R}}^p}{\max }\left\{\frac{1}{2n}\Vert y-X\theta {\Vert }_2^2+\lambda \Vert \theta {\Vert }_1\right\}.$$

$$H_0:\Vert Q{\theta }_0{\Vert }_2\in {\Omega }_0,$$

$$\eta ({\theta }_0,\sigma )=\frac{\mathbb{E}((x_i^{\mathsf{T}}{\theta }_0{)}^2)}{\mathrm{Var}(y_i)}=\frac{\mu }{1+\mu },$$

$$\Vert X{\Vert }_{{\psi }_2}=\underset{q\ge 1}{\sup }{q}^{-1/2}(\mathbb{E}|X{|}^q{)}^{1/q}.$$

$$\Vert X{\Vert }_{{\psi }_2}=\underset{\Vert x\Vert \le 1}{\sup }\Vert ⟨X,x⟩{\Vert }_{{\psi }_2}.$$

$$\displaystyle\phi^{2}(J,S)\equiv\min_{{\theta}\in\mathbb{R}^{p}}\Big{\{}\frac{|S|\langle{\theta},J{\theta}\rangle}{\|{\theta}_{S}\|_{1}^{2}}:\quad{\theta}\in\mathbb{R}^{p},\,\|{\theta}_{S^{c}}\|_{1}\leq 3\|{\theta}_{S}\|_{1}\Big{\}}\,.$$

$${\mathcal{P}}_U(\theta )=\sum _{i=1}^k⟨\theta ,u_i⟩u_i=U{U}^{\mathsf{T}}\theta ,$$

$${\tilde{H}}_0:{\mathcal{P}}_U({\theta }_0)\in {\mathcal{P}}_U({\Omega }_0).$$

$$\underset{\theta \in {\Omega }_0}{\sup }{\mathbb{P}}_{\theta }(R=1)\le \underset{{\mathcal{P}}_U(\theta )\in {\mathcal{P}}_U({\Omega }_0)}{\sup }{\mathbb{P}}_{\theta }(R=1)\le \alpha .$$

$$\{{\hat{\theta }}^n(\lambda ),\hat{\sigma }(\lambda )\}=\underset{\theta \in {\mathbb{R}}^p,\sigma >0}{\mathrm{argmin}}\left\{\frac{1}{2\sigma n}\Vert y-X\theta {\Vert }_2^2+\frac{\sigma }{2}+\lambda \Vert \theta {\Vert }_1\right\}.$$

$${\hat{\gamma }}^{\mathrm{d}}=U^{\mathsf{T}}\hat{\theta }+\frac{1}{n}{G}^{\mathsf{T}}{X}^{\mathsf{T}}(y-X\hat{\theta }),$$

$$\begin{array}{rl}& \text{minimize}{g}^{\mathsf{T}}\hat{\Sigma }g\\ & \text{subject to}\Vert \hat{\Sigma }g-u_i{\Vert }_{\infty }\le \mu \end{array}$$

$$\sqrt{n}({\hat{\gamma }}^{\mathrm{d}}-U^{\mathsf{T}}{\theta }_0)=Z+\Delta ,Z\sim \mathsf{N}(0,{\sigma }^2{G}^{\mathsf{T}}\hat{\Sigma }G),\Delta =\sqrt{n}(G^{\mathsf{T}}\hat{\Sigma }-U^{\mathsf{T}})({\theta }_0-\hat{\theta }).$$

$$\mathbb{P}\left(\Vert \Delta {\Vert }_{\infty }\ge \frac{c\mu \sigma s_0}{{\varphi }_0^2}\sqrt{\log p}\right)\le 2p^{-c_0}+2e^{-n/16},c_0=\frac{c^2}{32K}-1.$$

$$\begin{array}{rl}& \text{minimize}{g}^{\mathsf{T}}\hat{\Sigma }g\\ & \text{subject to}\Vert \hat{\Sigma }g-u_i{\Vert }_{\infty }\le \mu ,\\ & |u_i^{\mathsf{T}}\hat{\Sigma }g-1|\le \mu .\end{array}$$

$$Q^{(n)}\equiv \frac{{\hat{\sigma }}^2}{n}(G^{\mathsf{T}}\hat{\Sigma }G),D^{(n)}\equiv \mathrm{diag}(\{Q_{ii}^{(n)}{\}}^{-1/2}).$$

$$\begin{array}{rl}{\theta }^{\mathrm{p}}=& \underset{\theta \in {\mathbb{R}}^p}{\mathrm{argmin}}\Vert D({\hat{\gamma }}^{\mathrm{d}}-U^{\mathsf{T}}\theta ){\Vert }_{\infty }\\ & \text{subject to}\theta \in {\Omega }_0.\end{array}$$

$$T_n=\Vert D({\hat{\gamma }}^{\mathrm{d}}-U^{\mathsf{T}}{\theta }^{\mathrm{p}}){\Vert }_{\infty }$$

$$R_X(y)=\{\begin{array}{ll}1& \text{ if }{T}_n\ge z_{\alpha /(2k)}(\text{reject }{\tilde{H}}_0)\\ & \\ 0& \text{ otherwise }(\text{fail to reject }{\tilde{H}}_0).\end{array}$$

$$\mathbb{P}\left(\Vert D({\hat{\gamma }}^{\mathrm{d}}-U^{\mathsf{T}}{\theta }_0){\Vert }_{\infty }\le t\right)-\mathbb{P}\left(\Vert DZ{\Vert }_{\infty }\le t\right)=o_P(1).$$

$${\mu }_{\ast }(X;U;G)\equiv |\hat{\Sigma }G-U{|}_{\infty },$$

$$\displaystyle\alpha_{n}(R)\equiv\sup\Big{\{}{\mathbb{P}}_{{\theta}_{0}}(R=1):\,\,{\theta}_{0}\in\Omega_{0},\;\|{\theta}_{0}\|_{0}\leq s_{0}(n)\Big{\}}\,.$$

$$\underset{n\to \infty }{\mathrm{limsup}}{\alpha }_n(R_X)\le \alpha .$$

$$\mathrm{d}({\theta }_0,{\Omega }_0;U)=\underset{\theta \in {\Omega }_0}{\inf }\Vert U^{\mathsf{T}}(\theta -{\theta }_0){\Vert }_{\infty }.$$

$$\displaystyle\beta_{n}(R)\equiv\sup\Big{\{}{\mathbb{P}}_{{\theta}_{0}}(R=0):\,\,\|{\theta}_{0}\|_{0}\leq s_{0}(n),\;{\rm d}({\theta}_{0},\Omega_{0};U)\geq\eta\Big{\}}$$

$$\underset{n\to \infty }{\mathrm{liminf}}\frac{1-{\beta }_n(R_X)}{1-{\beta }_n^{\ast }(\eta )}\ge 1,1-{\beta }_n^{\ast }(\eta )\equiv F{\left(\alpha ,\frac{\sqrt{n}\eta }{\hat{\sigma }{m}_0},k\right)}_{+}$$

Full text

1 Introduction

Consider the high-dimensional regression model where we are given $n$ i.i.d. pairs $(y_{1},x_{1})$, $(y_{2},x_{2})$, $\cdots$, $(y_{n},x_{n})$ with $y_{i}\in\mathbb{R}$, and $x_{i}\in\mathbb{R}^{p}$, denoting the response values and the feature vectors, respectively. The linear regression model posits that response values are generated as

[TABLE]

Here ${\theta}_{0}\in\mathbb{R}^{p}$ is a vector of parameters to be estimated. In matrix form, letting $y=(y_{1},\dots,y_{n})^{\sf T}$ and denoting by $X$ the matrix with rows $x_{1}^{\sf T}$,$\cdots$, $x_{n}^{\sf T}$ we have

[TABLE]

We are interested in high-dimensional models where the number of parameters $p$ may far exceed the sample size $n$. To make informative inference feasible in this setting, we assume sparsity structure for the model, that is $\theta_{0}$ has only a few ($s_{0}<n$) number of nonzero entries, whose identities are unknown.

Our goal in this paper is to understand various parameter structures of the high-dimensional model. Specifically, we develop a flexible framework for testing null hypothesis of the form

[TABLE]

for a general set $\Omega_{0}\subset\mathbb{R}^{p}$. Remarkably, we make no additional assumptions (such as convexity or connectedness) on $\Omega_{0}$.

In Section 5, we will relax the sparsity assumption on the model parameters to the approximate sparsity. Consider the linear model $y=X\theta_{*}+w$, where $\theta_{*}\in\mathbb{R}^{p}$ is not necessarily sparse. The approximate sparsity posits that even if the true signal $X\theta_{*}$ cannot be written as a sparse linear combination of the covariates, there exists at least one sparse linear combination of the covariates that gets close to the true signal. Formally, we assume that there exists a vector $\theta_{0}\in\mathbb{R}^{p}$ such that $\|\theta_{0}\|_{0}=s_{0}$, and $\|X\theta_{*}-X\theta_{0}\|=o_{P}(1)$. Note that this notion of approximate sparsity is similar to but stronger than the one introduced in [BTW07, BCCH12].333In [BCCH12], the approximate sparsity assumption allows $\|X\theta_{*}-X\theta_{0}\|=O_{P}(\sqrt{s_{0}})$, while here we are imposing stronger requirement $\|X\theta_{*}-X\theta_{0}\|=o_{P}(1)$.

In addition, in Section 6 we extend our analysis to non-gaussian noise.

1.1 Motivation

High-dimensional models are ubiquitous in many areas of applications. Examples range from signal processing (e.g. compressed sensing), to recommender systems (collaborative filtering), to statistical network analysis, to predictive analytics, etc. The widespread interest in these applications has spurred remarkable progress in the area of high-dimensional data analysis [CT07, BRT09, BvdG11]. Given that the number of parameters goes beyond the sample size, there is no hope to design reasonable estimators without making further assumption on the structure of model parameters. A natural such assumption is sparsity, which posits that only $s_{0}$ of the parameters ${\theta}_{0,i}$ are nonzero, and $s_{0}\leq n$. A prominent approach in this setting for estimating the model parameters is via the Lasso estimator [Tib96, CD95] defined by

[TABLE]

(We will omit the arguments of $\widehat{\theta}^{n}(y,X;\lambda)$ whenever clear from the context.)

To date, the majority of work on high-dimensional parametric models has focused on point estimation such as consistency for prediction [GR04], oracle inequalities and estimation of parameter vector [CT07, BRT09, RWY09], model selection [MB06, ZY06, Wai09], and variable screening [FL08]. The work [BTW07] extended the oracle inequalities for the lasso to the setting of weak sparsity and weak approximation, where the effect of covariates on the response can be controlled up to a small approximation error by conditioning on a relatively small number of covariates, whose identities are unknown. The minimax rate for estimating the parameters in the high-dimensional linear model was studied in [YZ10, RWY11], assuming that the true model parameters belong to some $\ell_{q}$ ball.

Despite this remarkable progress, the fundamental problem of statistical significance is far less understood in the high-dimensional setting. Uncertainty assessment is particularly important when one seeks subtle statistical patterns about the model parameters ${\theta}_{0}$.

Below, we discuss some important examples of high-dimensional inference that can be performed when provided a methodology for testing hypothesis of the form (3).

Example 1 (Testing ${\theta}_{\min}$ condition) Support recovery in high-dimension concerns the problem of finding a set $\widehat{S}\subseteq\{1,2,\dots,p\}$, such that ${\mathbb{P}}(\widehat{S}=S)$ is large, where $S\equiv\{i:\,{\theta}_{0,i}\neq 0,\,1\leq i\leq p\}$. Work on support recovery requires the nonzero parameters be large enough to be detected. Specifically, for exact support recovery meaning that ${\mathbb{P}}(\hat{S}\neq S)\to 1$, it is assumed that $\min_{i\in S}|{\theta}_{0,i}|=\Omega(\sqrt{(\log p)/n})$. This assumption is often referred to as ${\theta}_{\min}$ condition and is shown to be necessary for exact support recovery [MY09, ZY06, FL01, ZY06, Wai09, MB06].

Relaxing the task of exact support recovery, let $\alpha$ and $\beta$ be the type I and type II error rates in detecting nonzero (active) parameters of the model. In [JM14b], it is shown that even for gaussian design matrices, any hypothesis testing rule with nontrivial power $1-\beta>\alpha$ requires $\min_{i\in S}|{\theta}_{0,i}|=\Omega(1/\sqrt{n})$. Despite ${\theta}_{\min}$ assumption is commonplace, it is not verifiable in practice and hence it calls for developing methodologies that can test whether such condition holds true.

For a vector ${\theta}\in\mathbb{R}^{p}$, define support of ${\theta}$ as ${\rm supp}({\theta})=\{1\leq i\leq p:\,{\theta}_{i}\neq 0\}$. In (3), letting $\Omega_{0}=\{{\theta}\in\mathbb{R}^{p}:\,\min_{i\in{\rm supp}({\theta})}|{\theta}_{i}|\geq c\}$, we can test ${\theta}_{\min}$ condition for any given $c\geq 0$ and at a pre-assigned significance level $\alpha$.

Example 2 (Confidence intervals for quadratic forms) We can apply our method to test hypothesis of form

[TABLE]

for some given set $\Omega_{0}\subseteq[0,\infty)$ and a given matrix $Q\in\mathbb{R}^{m\times p}$. By duality between hypothesis testing and confidence interval, we can also construct confidence intervals for quadratic forms $\|Q{\theta}_{0}\|$.

In the case of $Q={\rm I}$, this yields inference on the signal strength $\|{\theta}\|_{2}^{2}$. As noted in [JBC17], armed with such testing method one can also provide confidence intervals for the estimation error, namely $\|{\widehat{\theta}}-{\theta}_{0}\|_{2}^{2}$. Specifically, we split the collected samples into two independent groups $(y^{(0)},X^{(0)})$ and $(y^{(1)},X^{(1)})$, and construct an estimate ${\widehat{\theta}}$ just by using the first group. Letting $\tilde{y}\equiv y^{(1)}-X^{(1)}{\widehat{\theta}}$, we obtain a linear regression model $\tilde{y}=X^{(1)}({\theta}_{0}-{\widehat{\theta}})+w$. Further, if ${\widehat{\theta}}$ is a sparse estimate, then ${\theta}_{0}-{\widehat{\theta}}$ is also sparse. Therefore, inference on the signal strength on the obtained model is similar to inference on the error size $\|{\theta}_{0}-{\widehat{\theta}}\|_{2}^{2}$.

Inference on quadratic forms turns out to be closely related to a number of well-studied problems, such as estimate of the noise level $\sigma^{2}$ and the proportion of explained variation [FGH12, BEM13, Dic14, JBC17, VG18, GWCL19]. To expand on this point, suppose that attributes $x_{i}$ are drawn i.i.d. from a gaussian distribution with covariance $\Sigma$, and the noise level $\sigma^{2}$ is unknown. Then, ${\rm Var}(y_{i})=\sigma^{2}+\|\Sigma^{1/2}{\theta}_{0}\|_{2}^{2}$. Since $\|y\|_{2}^{2}/{\rm Var}(y_{i})$ follows a $\chi^{2}$ distribution with $n$ degrees of freedom, we have $\|y\|_{2}^{2}/n={\rm Var}(y_{i})[1+O_{P}(n^{-1/2})]$. Hence, task of inference on the quadratic form $\|\Sigma^{1/2}{\theta}_{0}\|_{2}^{2}$ and the noise level $\sigma^{2}$ are intimately related. This is also related to the proportion of explained variation defined as

[TABLE]

with $\mu=(1/\sigma^{2})\|\Sigma^{1/2}{\theta}_{0}\|_{2}^{2}$ the signal-to-noise ratio. This quantity is of crucial importance in genetic variability [VHW08] as it somewhat quantifies the proportion of variance in a trait (response) that is explained by genes (design matrix) rather than environment (noise part).

Example 3 (Testing individual parameters ${\theta}_{0,i}$) Recently, there has been a significant interest in testing individual hypothesis $H_{0,i}:{\theta}_{i}=0$, in the high-dimensional regime. This is a challenging problem because obtaining an exact characterization of the probability distribution of the parameter estimates in the high-dimensional regime is notoriously hard.

A successful approach is based on debiasing the regularized estimators. The resulting debiased estimator is amenable to distributional characterization which can be used for inference on individual parameters [JM14a, JM14b, ZZ14, VdGBRD14, JM13]. Our methodology for testing hypothesis of form (3) is built upon the debiasing idea. It also recovers the debiasing approach for $\Omega_{0}=\{{\theta}\in\mathbb{R}^{p}:\,{\theta}_{i}=0\}$.

Example 4 (Confidence intervals for predictions) For a new sample $\xi$, we can perform inference on the response value $\xi^{\sf T}{\theta}_{0}$ by letting $\Omega_{0}=\{{\theta}:\,\xi^{\sf T}{\theta}_{0}=c\}$ for a given value $c$. Further, by duality between hypothesis testing and confidence intervals, we can construct confidence interval for $\xi^{\sf T}{\theta}_{0}$. We refer to Section 7 for further details.

Example 5 (Confidence intervals for $f({\theta}_{0})$) Let $f:\mathbb{R}^{p}\to\mathbb{R}$ be an arbitrary function. By letting $\Omega_{0}=\{{\theta}:\,f({\theta}_{0})=c\}$ we can test different values of $f({\theta}_{0})$. Further, by employing the duality relationship between hypothesis testing and confidence intervals, we can construct confidence intervals for $f({\theta}_{0})$. Note that Examples 3, 4 are special cases of $f({\theta}_{0})=e_{i}^{\sf T}{\theta}_{0}$ and $f({\theta}_{0})=\xi^{\sf T}{\theta}_{0}$. Here $e_{i}$ is the $i$-th standard basis element with one at the $i$-th entry and zero everywhere else.

Example 6 (Testing over convex cones) For a given cone ${\cal C}$, our framework allows us to test whether ${\theta}_{0}$ belongs to ${\cal C}$. Some examples that naturally arise in studying treatment effects are nonnegative cone ${\cal C}_{\geq 0}=\{{\theta}\in\mathbb{R}^{p}:\,{\theta}_{i}\geq 0\,\,\text{for all }1\leq i\leq p\}$, and monotone cone ${\cal C}_{M}=\{{\theta}\in\mathbb{R}^{p}:\,{\theta}_{1}\leq{\theta}_{2}\leq\dotsc\leq{\theta}_{p}\}$. Letting ${\theta}_{i}$ denote the mean of treatment $i$, by testing ${\theta}_{0}\in{\cal C}_{\geq 0}$, one can test whether all the treatments in the study are harmless. Another case is when treatments correspond to an ordered set of dosages of the same drug. Then, one might reason that if the drug is of any effect, its effect should follow a monotone relationship with its dosage. This hypothesis can be cast as ${\theta}_{0}\in{\cal C}_{M}$. Such testing problems over cones have been studied for gaussian sequence models by [Kud63, RW78, RCLN86], and very recently by [WWG19].

1.2 Other Related work

Testing in the high-dimensional linear model has experienced a resurgence in the past few years. Most closely related to us is the line of work on debiasing/desparsifying pioneered by [ZZ14, VdGBRD14, JM14a]. These papers propose a debiased estimator $\widehat{\theta}^{{\rm d}}$ such that every coordinate $\widehat{\theta}^{{\rm d}}_{i}$ is approximately gaussian under the condition that ${s_{0}^{2}(\log p)}/{n}\to 0$, which is in turn used to test single coordinates of ${\theta}_{0}$, $H_{0}:{\theta}_{0,i}=0$, and construct confidence intervals for ${\theta}_{0,i}$. In a parallel line of work, [BCH11, BCFVH17, BCH13, BCH14] have also designed an asymptotically gaussian pivot via the post-double-selection lasso, under the same sample size condition of ${s_{0}^{2}(\log p)}/{n}\to 0$. [CG17] established that the sample size conditions required by debiasing and post-double-selection are minimax optimal meaning to construct a confidence interval of length $O({1}/{\sqrt{n}})$ for a coordinate of ${\theta}_{0}$ requires ${s_{0}^{2}(\log p)}/{n}\to 0$.

The debiasing and post-double-selection approaches have also been applied to a wide variety of other models for testing ${\theta}_{0,i}$ including missing data linear regression [WWBS19], quantile regression [ZKL14], and graphical models [RSZZ15, CRZZ16, WK16, BK18].

In the multiple testing realm, the debiasing approach has been used to control directional FDR [JJ19]. Other methods such as FDR-thresholding and SLOPE procedures controls the false discovery rate (FDR) when the design matrix $X$ is orthogonal [SC16, BvdBS*+*15, ABDJ06]. In the non-orthogonal setting, the knockoff procedure [BC15] controls FDR whenever $n\geq 2p$, and the noise is isotropic; In [JS16], knockoff was generalized to also control for the family-wise error rate. More recently, [CFJL18] developed the model-free knockoff which allows for $p>n$ when the distribution of $X$ is known.

In parallel, there have been developments in selective inference, namely inference for the variables that the lasso selects. [LSST16, TTLT16] developed exact tests for the regression coefficients corresponding to variables that lasso selects. This was further generalized to a wide variety of polyhedral model selection procedures including marginal screening and orthogonal matching pursuit in [LT14]. [TT18, FST14, HPM*+*16] developed more powerful and general selective inference procedures by introducing noise in the selection procedure. To allow for selective inference in the high-dimensional setting, [LSST16] combined the polyhedral selection procedure with the debiased lasso to construct selectively valid confidence intervals for ${\theta}_{0,i}$ when ${s_{0}(\log p)}/{\sqrt{n}}\to 0$.

Much of the previous work has focused on testing coordinates or one-dimensional projections of ${\theta}_{0}$. An exception is the work [NvdG13] which studies the problem of constructing confidence sets for the high dimensional linear models, so that the confidence sets are honest over the family of sparse parameters, under i.i.d gaussian designs. Our work increases the applicability of the debiasing approach by allowing for general hypothesis, ${\theta}_{0}\in\Omega_{0}$. The set $\Omega_{0}$ can be non-convex or even disconnected. Our setup encompasses a broad range of testing problems and it is shown to be minimax optimal for special cases such as $\Omega=\{{\theta}:{\theta}_{i}=0\}$ and $\Omega_{0}=\{{\theta}:\xi^{\sf T}{\theta}=c\}$.

The authors in [ZB17] have studied the problem (3) independently and indeed [ZB17] was posted online around the same time that the first draft of our paper was released. This work also leverages the idea of debiasing but greatly differs from this work, both in methodology and theory, which we now discuss. In [ZB17], the debiased estimator is constructed in the standard basis (as compared to ours which is done in a lower dimensional subspace) and is followed by an $\ell_{1}$ projection to construct the test statistic. The test statistic involves a data dependent vector and the method uses bootstrap to approximate the distribution of the test statistic and set the critical values. In terms of theory, [ZB17] shows that the proposed method controls the type I error at the desired level assuming that $\log p=o(n^{1/8})$ and $s_{0}=o(n^{1/4}/\sqrt{\log p})$ (See Theorem 1 therein), while we prove such result for our test under $s_{0}=o(\sqrt{n}/\log p)$. It is shown in [ZB17] that the rule achieves asymptotic power one provided that the signal strength (measured in term of the $\ell_{\infty}$ distance of ${\theta}_{0}$ from $\Omega_{0}$) asymptotically dominates $n^{-1/4}$. In comparison, in Theorem 3.4 we establish a lower bound of the power for all values of the signal strength and as a corollary of that we show the method achieves power one if the signal strength dominates $n^{-1/2}$ asymptotically.

1.3 Organization of the paper

In the remaining part of the introduction, we present the notations and a few preliminary definitions. The rest of the paper presents the following contributions:

• •

Section 2. We explain our testing methodology. It consists of constructing a debiased estimator for the projections of the model parameters in a lower dimensional subspace. It is then followed by an $\ell_{\infty}$ projection to form the test statistic.

• •

Section 3. We present our main results. Specifically, we show that our method controls false positive rate under a pre-assigned $\alpha$ level. We also derive an analytical lower bound for the statistical power of our test. In case of $\Omega_{0}=\{{\theta}\in\mathbb{R}^{d}:\,{\theta}_{i}=0\}$ (Example 3), it matches the bound proposed in [JM14a, Theorem 3.5], which is also shown to be minimax optimal.

• •

Section 5. We explain the notion of approximate sparsity and discuss how our results can be extended to allow for approximately sparse models.

• •

Section 6. We relax the gaussianity assumption on the noise component and discuss how to address possibly non-gaussian noise under proper moment conditions.

• •

Section 7. We provide applications of our framework for some special cases: Inference on linear predictions, quadratic forms of the parameters and testing the $\theta_{\min}$ condition. In Section 7.1, we discuss the existing literature for these subproblems and compare it to our proposed methodology.

• •

Section 8. We provide numerical experiments to corroborate our findings and evaluate type I error and statistical power of our test under various settings.

• •

Section 9. Proof of Theorems are given in this section, while the proof of technical lemmas are deferred to appendices.

1.4 Notations

We start by adapting some simple notations that will be used throughout the paper, along with some basic definitions from the literature on high-dimensional regression.

We use $e_{i}$ to refer to the $i$-th standard basis element, e.g., $e_{1}=(1,0,\dotsc,0)$. For a vector $v$, ${\rm supp}(v)$ represents the positions of nonzero entries of $v$. For a vector $\theta$ and a subset $S$, $\theta_{S}$ is the restriction of $\theta$ to indices in $S$. For an integer $p\geq 1$, we use the notation $[p]=\{1,\cdots,p\}$. We write $\|v\|_{p}$ for the standard $\ell_{p}$ norm of a vector $v$, i.e., $\|v\|_{p}=(\sum_{i}|v_{i}|^{p})^{1/p}$ and $\|v\|_{0}$ for the number of nonzero entries of $v$. Whenever the subscript $p$ is not mentioned it should be read as $\ell_{2}$ norm. For a matrix $A$, we denote by $|A|_{\infty}\equiv\max_{i\leq m,j\leq n}|A_{ij}|$, the maximum absolute value of entries of $A$. Further, its maximum and minimum singular values are respectively indicated by by $\sigma_{\max}(A)$ and $\sigma_{\min}(A)$. Throughout, $\Phi(x)\equiv\int_{-\infty}^{x}e^{-t^{2}/2}{\rm d}t/\sqrt{2\pi}$ denotes the CDF of the standard normal distribution. We also denote the $z$-values $z_{\alpha}=\Phi^{-1}(1-\alpha)$.

The term “with high probability” means with probability converging to one as $n\to\infty$ and for two functions $f(n)$ and $g(n)$, the notation $f(n)=o(g(n))$ means that $g$ ‘dominates’ $f$ asymptotically, namely, for every fixed positive $C$, there exists $n(C)$ such that $f(n)\leq Cg(n)$ for $n>n(C)$. Likewise, $f(n)=O(g(n))$ indicates that $f$ is ‘bounded’ above by $g$ asymptotically, i.e., $f(n)\leq Cg(n)$ for some positive constant $C$. Analogously, we use he notations $o_{P}(\cdot)$ and $O_{P}(\cdot)$ to indicate asymptotic behavior is probability as the sample size $n$grows to infinity.

Let $\widehat{\Sigma}=(X^{\sf T}X)/n\in\mathbb{R}^{p\times p}$ be the sample covariance of the design $X\in\mathbb{R}^{n\times p}$. In the high-dimensional setting, where $p$ exceeds $n$, $\widehat{\Sigma}$ is singular. As common in high-dimensional statistics, we assume compatibility condition which requires $\widehat{\Sigma}$ to be nonsingular in a restricted set of directions.

We use the notation $\|\cdot\|_{\psi_{2}}$ to refer to the sub-gaussian norm. Specifically, for a random variable $X$, we let

[TABLE]

For a random vector $X\in\mathbb{R}^{m}$, its sub-gaussian norm is defined as

[TABLE]

Definition 1.1**.**

For a symmetric matrix $J\in\mathbb{R}^{p\times p}$ and a set $S\subseteq[p]$, the compatibility condition is defined as

[TABLE]

Matrix $J$ is said to satisfy compatibility condition for a set $S\subseteq[p]$, with constant $\phi_{0}$ if $\phi(J,S)\geq\phi_{0}$.

2 Projection statistic

Depending on the structure of $\Omega_{0}$ it may be useful to instead of testing the null hypothesis $H_{0}:\theta_{0}\in\Omega_{0}$, we test it in a lower dimensional space. Consider an $k$-dimensional subspace represented by an orthonormal basis $\{u_{1},\dotsc,u_{k}\}$, with $u_{i}\in\mathbb{R}^{p}$. For this section, we assume that the basis $\{u_{1},\ldots,u_{k}\}$ is predetermined and fixed. In Section 4, we discuss how to choose the subspace depending on $\Omega_{0}$ to maximize the power of the test. The projection onto this subspace is given by

[TABLE]

where $U=[u_{1},\dotsc,u_{k}]\in\mathbb{R}^{p\times k}$. We also use the notation $\mathcal{P}_{U}(\Omega_{0})=\{\mathcal{P}_{U}(\theta):\,\theta\in\Omega_{0}\}$ to denote the projection of $\Omega_{0}$ onto the subspace $U$. Define the hypothesis

[TABLE]

Under the null $H_{0}$, $\tilde{H}_{0}$ also holds, so controlling the type-I error of $\tilde{H}_{0}$ also controls the type-I error of $H_{0}$. In the following we propose a testing rule $R\in\{0,1\}$ for the null hypothesis $\tilde{H}_{0}$ and show that it controls type-I error below a pre-assigned level $\alpha$. Consequently,

[TABLE]

For now, we consider an arbitrary fixed subspace $U$, and then after we analyze the statistical power of our test we provide guidelines on how to choose $U$ to increase the power.

In order to test $\tilde{H}_{0}$ we construct a test statistic based on the debiasing approach.

We first let $\{{\widehat{\theta}},\widehat{\sigma}\}$ be the scaled Lasso estimator [SZ12] given by

[TABLE]

This optimization simultaneously gives an estimate of ${\theta}_{0}$ and $\sigma$. We use regularization parameter $\lambda=\sqrt{2.05(\log p)/n}$. Due to the $\ell_{1}$ penalization, the lasso estimator ${\widehat{\theta}}$ is biased towards small $\ell_{1}$ norm, and so is the projection $\mathcal{P}_{U}(\theta_{0})$. We view $\mathcal{P}_{U}(\theta_{0})$ in the basis $U$, namely $\gamma_{0}=U^{\sf T}\theta_{0}$ and construct a debiased estimator for it in the following way:

[TABLE]

with the decorrelating matrix ${G}=[g_{1}|\dotsc|g_{k}]\in\mathbb{R}^{p\times k}$, where each $g_{i}$ is obtained by solving the optimization problems for each $1\leq i\leq k$:

[TABLE]

Note that the decorrelating matrix $G\in\mathbb{R}^{p\times p}$ is a function of $X$, but not of $y$. We next state a lemma that provides a a bias-variance decomposition for $\widehat{\gamma}^{{\rm d}}$ and brings insight about the form of debiasing given by (11).

Lemma 2.1**.**

Let $X\in\mathbb{R}^{n\times p}$ be any (deterministic) design matrix. Assuming that optimization problem (12) is feasible for $i\in[k]$, let $\widehat{\gamma}^{{\rm d}}=\widehat{\gamma}^{{\rm d}}(\lambda)$ be a general debiased estimator as per Eq (11). Then, setting $Z=G^{\sf T}X^{\sf T}w/\sqrt{n}$, with $w$ the noise vector in the regression (2), we have

[TABLE]

Further, assume that $X$ satisfies the compatibly condition for the set $S={\rm supp}(\theta_{0})$, $|S|\leq s_{0}$, with constant $\phi_{0}$, and let $K\equiv\max_{i\in[p]}(X^{\sf T}X/n)_{ii}$. Then, choosing $\lambda=c\sqrt{(\log p)/n}$, we have

[TABLE]

Lemma 2.1 can be proved in a similar way to Theorem 2.3 of [JM14a] and its proof is omitted here. The decomposition (13) explains the rationale behind optimization (12). Indeed the convex program (12) aims at optimizing two objectives. On one hand, the constraint controls the term $|G^{\sf T}\widehat{\Sigma}-U^{\sf T}|_{\infty}$, which by Lemma 2.1 controls the bias term $\|\Delta\|_{\infty}$. On the other hand, it minimizes the objective function $g^{\sf T}\widehat{\Sigma}g$, which controls the variance of $\widehat{\gamma}^{{\rm d}}_{i}$. Therefore, the parameter $\mu$ in optimization (12) controls the bias-variance tradeoff and should be chosen large enough to ensure that (12) is feasible. (See Section 3.1 for further discussion.)

Remark 2.2**.**

In the special case of $k=1$ and $u=e_{i}$, the debiased estimator (11) reduces to the one introduced in [JM14a]. For the special case of $k=1$, it becomes similar to the estimator proposed by [CG17] that is used to construct confidence intervals for linear functionals of $\theta_{0}$. Note that the proposed debiasing procedure incurs small bias in the infinity norm with respect to the rotated basis, $\left\|\widehat{\gamma}^{{\rm d}}-U^{\sf T}\theta_{0}\right\|_{\infty}$, as opposed to the standard debiasing procedure  [JM14a, JM14b, ZZ14, VdGBRD14, JM13] which incurs small bias, in the infinity norm, with respect to the original basis, and not necessarily in the rotated basis.

The following assumption ensures that the entries of noise $Z$ have non-vanishing variances.

Assumption 2.3**.**

We have $\lim\inf_{n\to\infty}\min_{i\in[k]}(G^{\sf T}\widehat{\Sigma}G)_{i,i}\geq c_{0}>0$, for some positive constant $c_{0}$.

The above assumption entails the decorrelating matrix $G$, where our proposal constructs via optimization (12). In the following lemma, we provide a sufficient condition for the above assumption to hold.

Lemma 2.4**.**

Suppose that $\lim\sup_{n\to\infty}\mu(\max_{i\in[k]}\,\|u_{i}\|_{1})\leq c<1$ and $\lim\sup_{n\to\infty}\max_{i\in[k]}(u_{i}^{\sf T}\widehat{\Sigma}u_{i})<C<\infty$, for some constant $c,C$. Then, Assumption 2.3 holds.

We refer to Appendix A.1 for the proof of Lemma 2.4.

Remark 2.5**.**

The very recent work [CCG19] uses the debiasing approach for inference on individualized treatment effect (and for general linear function $u^{T}\theta_{0}$ ). The proposed mechanism slightly differs from (12) in that it includes an extra constraint. By this trick, the proposed mechanism of [CCG19] can be used for inference on a broad family of loading vector $u$. We can follow the same idea and replace optimization (12) by the following optimization

[TABLE]

This way Assumption 2.3 is automatically satisfied (See [CCG19, Lemma 1] for the details).

Define the shorthand

[TABLE]

To ease the notation, we hereafter drop the superscript $(n)$. We next construct a test statistic $T_{n}$ so that the large values of $T_{n}$ provide evidence against the null hypothesis. For this, consider the $\ell_{\infty}$ projection estimator given by

[TABLE]

We then define the test statistic to be the optimal value of (17), i.e.,

[TABLE]

The reason for using $\ell_{\infty}$ norm in the projection is that the bias term of $\widehat{\gamma}^{{\rm d}}$ is controlled in $\ell_{\infty}$ norm (See Lemma 2.1.) The decision rule is then based on the test statistic:

[TABLE]

The above procedure generalizes the debiasing approach of [JM14a]. Specifically, for $\Omega_{0}=\{{\theta}:\theta_{1}=0\}=\{0\}\times\mathbb{R}^{p-1}$ and $U=e_{1}e_{1}^{\sf T}$, the test rule becomes the one proposed by [JM14a] for testing hypothesis of the form $H_{0}:\theta_{0,1}=0$ versus its alternative.

Remark 2.6**.**

Using Lemma 2.1, under the null hypothesis $H_{0}:\,\theta_{0}\in\Omega_{0}$, we have that $D(\widehat{\gamma}^{{\rm d}}-U^{\sf T}{\theta^{{\rm p}}})$ is (asymptotically) stochastically dominated by $DZ$, whose entries are dependent and are distributed as standard normal. The choice of threshold $z_{\alpha/(2k)}$ in (19) comes from using this observation and union bounding to control the (two-sided) tail of $\|DZ\|_{\infty}$. Given that Lemma 2.1 also characterizes the dependency structure of the entries of $DZ$, we can pursue another (less conservative) approach to choose the rejection threshold. As an implication of Lemma 2.1, and since $k$ (dimension of $Z$) is fixed, we have that for all $t\in\mathbb{R}$,

[TABLE]

Under the null hypothesis $H_{0}$, we have $\|D(\widehat{\gamma}^{{\rm d}}-U^{\sf T}{\theta^{{\rm p}}})\|_{\infty}\leq\|D(\widehat{\gamma}^{{\rm d}}-U^{\sf T}\theta_{0})\|_{\infty}$, and by (20), the distribution of $\|D(\widehat{\gamma}^{{\rm d}}-U^{\sf T}\theta_{0})\|_{\infty}$ is asymptotically equal to the maximum of dependent standard normal variables $\|DZ\|_{\infty}$, whose distribution can be easily simulated since the covariance of the multivariate gaussian vector $DZ$ is known.

In the next section, we prove that decision rule (19) controls type-I error below the target level $\alpha$ provided the basis $U$ is independent of the samples $(y_{i},x_{i})$, $1\leq i\leq n$. We also develop a lower bound on the statistical power of the testing rule and use that to choose the basis $U$.

3 Main results

3.1 Controlling false positive rate

Definition 3.1**.**

Consider a given triple $(X;U;{G})$ where $X\in\mathbb{R}^{n\times p}$, $U\in\mathbb{R}^{p\times k}$ with $U^{\sf T}U=I$ and ${G}\in\mathbb{R}^{p\times k}$. The generalized coherence parameter of $(X;U;{G})$ denoted by $\mu_{*}(X;U;{G})$ is given by

[TABLE]

where $\widehat{\Sigma}=(X^{\sf T}X)/n$ is the sample covariance of $X$. The minimum generalized coherence of $(X;U)$ is $\mu_{\min}(X;U)=\min_{{G}\in\mathbb{R}^{p\times k}}\mu_{*}(X;U;{G})$.

Note that choosing $\mu\geq\mu_{\min}(X;U)$, the optimization (12) becomes feasible.

We take a minimax perspective and require that the probability of type I error (false positive) to be controlled uniformly over $s_{0}$-sparse vectors.

For a testing rule $R\in\{0,1\}$ and a set $\Omega_{0}$, we define

[TABLE]

Our first result shows validity of our test for general set $\Omega_{0}$ under deterministic designs.

Theorem 3.2**.**

Consider a sequence of design matrices $X\in\mathbb{R}^{n\times p}$, with dimensions $n\to\infty$, $p=p(n)\to\infty$ satisfying the following assumptions. For each $n$, the sample covariance $\widehat{\Sigma}=(X^{\sf T}X)/n$ satisfies compatibility condition for the set $S_{0}={\rm supp}({\theta}_{0})$, with a constant $\phi_{0}>0$. Also, assume that $K\geq\max_{i\in[p]}\widehat{\Sigma}_{ii}$ for some constant $K>0$. Also consider a sequence of matrices $U\in\mathbb{R}^{p\times k}$, with fixed $k$ and $p=p(n)\to\infty$, such that $U^{\sf T}U=I_{k}$.

Consider the linear regression (2) and let $\widehat{\theta}^{n}$ and $\widehat{\sigma}$ be obtained by scaled Lasso, given by (10), with $\lambda=c\sqrt{(\log p)/n}$. Construct a debiased estimator $\widehat{\gamma}^{{\rm d}}$ as in (11) using $\mu\geq\mu_{\min}(X;U)$, where $\mu_{\min}(X;U)$ is the minimum generalized coherence parameter as per Definition 3.1, and suppose that Assumption 2.3 holds. Choose $c^{2}>32K$ and suppose that $s_{0}=o(\min\{{1}/(\mu\sqrt{\log p}),{n}/{\log p}\})$. For the test $R_{X}$ defined in Equation (19), and for any $\alpha\in[0,1]$, we have

[TABLE]

We next prove validity of our test for general set $\Omega_{0}$ under random designs.

Theorem 3.3**.**

Let $\Sigma\in\mathbb{R}^{p\times p}$ such that $\sigma_{\min}(\Sigma)\geq C_{\min}>0$ and $\sigma_{\max}(\Sigma)\leq C_{\max}<\infty$ and $\max_{i\in[p]}\Sigma_{ii}\leq 1$. Suppose that $X\Sigma^{-1/2}$ has independent sub-gaussian rows, with mean zero and sub-gaussian norm $\|\Sigma^{-1/2}x_{1}\|_{\psi_{2}}=\kappa$, for some constant $\kappa>0$.

Let $\widehat{\theta}^{n}$ and $\widehat{\sigma}$ be obtained by scaled Lasso, given by (10), with $\lambda=c\sqrt{(\log p)/n}$, and $c^{2}>48$. Consider an arbitrary $U\in\mathbb{R}^{p\times k}$, with $U^{\sf T}U=I$, that is independent of the samples $\{(x_{i},y_{i})\}_{i=1}^{n}$. Construct a debiased estimator $\widehat{\gamma}^{{\rm d}}$ as in (11) with $\mu=a\sqrt{(\log p)/n}$ and $a^{2}>48e^{2}\kappa^{4}C_{\max}/C_{\min}$. In addition, suppose that $\lim\sup_{n\to\infty}\mu(\max_{i\in[k]}\,\|u_{i}\|_{1})\leq c^{\prime}$, for some constant $0<c^{\prime}<1$ and $s_{0}=o(\sqrt{n}/{\log p})$.

For the test $R_{X}$ defined in Equation (19), and for any $\alpha\in[0,1]$, we have

[TABLE]

We refer to Section 9 for the proof of Theorem 3.2 and 3.3.

3.2 Statistical power

We next analyze the statistical power of our test. Before proceeding, note that without further assumption, we cannot achieve any non-trivial power, namely, power of $\alpha$ which is obtained by a rule that randomly rejects null hypothesis with probability $\alpha$. Indeed, by choosing ${\theta}_{0}\notin\Omega_{0}$ but arbitrarily close to $\Omega_{0}$, once can make $H_{0}$ essentially indistinguishable from $H_{A}$. Taking this point into account, for a set $\Omega_{0}\subseteq\mathbb{R}^{p}$ and ${\theta}_{0}\in\mathbb{R}^{p}$, we define the distance ${\rm d}({\theta}_{0},\Omega_{0})$ as

[TABLE]

We will assume that, under alternative hypothesis, ${\rm d}({\theta}_{0},\Omega_{0};U)\geq\eta$ as well. Define

[TABLE]

Quantity $\beta_{n}$ is the probability of type II error (false negative) and $1-\beta_{n}$ is the statistical power of the test.

Theorem 3.4**.**

Let $R_{X}$ be the test defined in Equation (19). Under the conditions of Theorem 3.3, for all $\alpha\in[0,1]$:

[TABLE]

where we define $m_{0}$ as

[TABLE]

Further, for $\alpha\in[0,1]$, $x\in\mathbb{R}_{+}$, and integer $k\geq 1$, the function $F(\alpha,x,k)$ is defined as follows:

[TABLE]

The proof of Theorem 3.4 is given in Section 9.3.

Note that for any fixed $k\geq 1$ and $\alpha>0$, the function $x\mapsto{F}(\alpha,x,k)$ is continuous and monotone increasing, i.e., the larger ${\rm d}({\theta}_{0},\Omega_{0})$ the higher power is achieved. Also, in order to achieve a specific power $\beta>\alpha$, our scheme requires $\eta>c_{\beta}m_{0}(\sigma/\sqrt{n})$, for some constant $c_{\beta}$ that depends on the desired power $\beta$. In addition, if $\eta\sqrt{n}\to\infty$, the rule achieves asymptotic power one.

It is worth noting that in case of testing individual parameters $H_{0,i}:{\theta}_{0,i}=0$ (corresponding to $\Omega_{0}=\{{\theta}\in\mathbb{R}^{p}:\,{\theta}_{i}=0\}$ and $k=1$), we recover the power lower bound established in [JM14a], which by comparing to the minimax trade-off studied in [JM14b], is optimal up to a constant.

4 Choice of subspace $U$

Before we start this section, let us stress again that by Theorems 3.2 and 3.3, the proposed testing rule controls type-I error below the desired level $\alpha$, for any choice of $U\in\mathbb{R}^{p\times k}$, with $1\leq k\leq p$ and $U^{\sf T}U={\rm I}$ that is independent of $X$. Here, we provide guidelines for choosing $U$ that yields high power. To this end we use the result of Theorem 3.4.

Note that

[TABLE]

where we recall that $\sigma_{\min}(\Sigma)>C_{\min}>0$ and $\|u_{i}\|=1$, for $i\in[k]$. Hence,

[TABLE]

We propose to choose $U$ by maximizing the right-hand side of (30), which by Theorem 3.4 serves as a lower bound for the power of the test. Nevertheless, the above optimization involves $\theta_{0}$ which is unknown. To cope with this issue, we use the Lasso estimate $\widehat{\theta}$ via the following procedure:

We randomly split the data $(y,X)$ into two subsamples $(y^{(1)},X^{(1)})$ and $(y^{(2)},X^{(2)})$ each with sample size $n_{0}=n/2$. We let ${\widehat{\theta}}^{(1)}$ be the optimizer of the scaled Lasso applied to $(y^{(1)},X^{(1)})$. 2. 2.

We choose $U\in\mathbb{R}^{p\times k}$ by solving the following optimization:

[TABLE] 3. 3.

We construct the debiased estimator using the data $(y^{(2)},X^{(2)})$. Specifically, set $\widehat{\Sigma}^{(2)}\equiv(1/n_{0})(X^{(2)})^{\sf T}(X^{(2)})$ and let $g_{i}$ be the solution of the following optimization problems for each $1\leq i\leq k$:

[TABLE]

Define the decorrelating matrix ${G}=[g_{1}|\dotsc|g_{k}]\in\mathbb{R}^{p\times k}$ and let ${\widehat{\theta}}^{(2)}$ be the optimizer of the scaled Lasso applied to $(y^{(2)},X^{(2)})$. Let

[TABLE] 4. 4.

Set $Q\equiv(\widehat{\sigma}^{2}/n)(G^{\sf T}\widehat{\Sigma}^{(2)}G)$ and $D\equiv{\rm diag}(\{Q_{ii}\}^{-1/2})$. Find the $\ell_{\infty}$ projection as

[TABLE] 5. 5.

Define the test statistics $T_{n}=\|D(\widehat{\gamma}^{{\rm d}}-U^{\sf T}\theta^{\rm p})\|_{\infty}$. The testing rule is given by

[TABLE]

Note that the data splitting above ensures that $U$ is independent of $(y^{(2)},X^{(2)})$, which is required for our analysis (See Theorems 3.2, 3.3 and 3.4.)

4.1 Convex sets $\Omega_{0}$

When the set $\Omega_{0}$ is convex, step (2) in the above procedure can be greatly simplified. Indeed, we can only focus on $k=1$ in this case.

Lemma 4.1**.**

Define the set $\mathcal{J}$ of matrices as

[TABLE]

If $\Omega_{0}$ is convex then there exists a unit norm $u^{*}\in\mathbb{R}^{p\times 1}$ such that $u^{*}\in\mathcal{J}$.

Proof of Lemma 4.1 is given in Appendix A.3.

Focusing on $k=1$, optimization (31) reduces to the following optimization over $u\in\mathbb{R}^{p\times 1}$:

[TABLE]

The function $x\mapsto{F}(\alpha,x,k)$ is monotone increasing in $x$ and by substituting for ${\rm d}(\theta_{0},\Omega_{0};u)$, this becomes equivalent to the following problem:

[TABLE]

Given that the objective is linear in $u$ and $\theta$, and the set $\Omega_{0}$ is convex we can apply the Von Neumann’s minimax theorem and change the order of $\max$ and $\min$:

[TABLE]

Denote the orthogonal projection of ${\widehat{\theta}}^{(1)}$ onto $\Omega_{0}$ by $\mathcal{P}_{\Omega_{0}}({\widehat{\theta}}^{(1)})=\arg\min_{{\theta}\in\Omega_{0}}\|{\theta}-{\widehat{\theta}}^{(1)}\|_{2}$. Then it is straightforward to see that the optimal $u$ is given by

[TABLE]

with $\mathcal{P}^{\perp}_{\Omega_{0}}({\widehat{\theta}}^{(1)})={\widehat{\theta}}^{(1)}-\mathcal{P}_{\Omega_{0}}({\widehat{\theta}}^{(1)})$.

We remind again that the type I error is controlled at the desired level for any $U\in\mathbb{R}^{p\times k}$ with $U^{\sf T}U={\rm I}$ that is independent of $(y,X)$. The choice of $u$ in (40) is a guideline for increasing power in case of convex sets $\Omega_{0}$.

Remark 4.2**.**

Let us stress again that convexity assumption of set $\Omega_{0}$ is crucial in deriving the recipe (40). To build further insight, we provide a concert example of a non-convex $\Omega_{0}$ and argue that $k=1$ is not the right choice. Let $\Omega_{0}=\Omega_{1}\cup\Omega_{2}$, where $\Omega_{i}=\{x\in\mathbb{R}^{p}:\,|x_{i}|\leq a,\,|x_{j}|\leq 3a,\,\text{ for }j\neq i\}$, for $i=1,2$ and a fixed constant $a>0$. Let $\theta_{0}=(2a,2a,0,\dotsc,0)\in\mathbb{R}^{p}$. Observe that $\Omega_{0}$ is not convex and $\theta_{0}\notin\Omega_{0}$. By choosing $k=p$ and $U={\rm I}_{p\times p}$, we have ${\rm d}(\theta_{0},\Omega_{0},U)=a$ and hence our method achieves non-trivial power. However, we argue that setting $k=1$, our method cannot do better than random guessing. Specifically, we show that for any vector $u\in\mathbb{R}^{p}$, we ${\rm d}(\theta_{0},\Omega_{0},u)=0$. By symmetry, assume that $|u_{1}|\leq|u_{2}|$. Note that the point $z_{0}=\pm(a\operatorname{sign}(u_{1}),3a\operatorname{sign}(u_{2}),\dotsc,3a\operatorname{sign}(u_{p}))\in\Omega_{1}\subset\Omega_{0}$. Further, $u^{\sf T}z_{0}=\pm(a|u_{1}|+3a|u_{2}|+\dotsc+3a|u_{p}|)$. By convexity of $\Omega_{1}$, we have that $\mathcal{P}_{u}(\Omega_{1})\supseteq A$ where $A=\{\alpha u:|\alpha|\leq|u^{\sf T}z_{0}|\}$. In addition, we have $u^{\sf T}\theta_{0}=2a(u_{1}+u_{2})$ and using the assumption $|u_{1}|\leq|u_{2}|$, we get $|u^{\sf T}\theta_{0}|\leq|u^{\sf T}z_{0}|$. Therefore, $\mathcal{P}_{u}(\theta_{0})\in A\subseteq\mathcal{P}_{u}(\Omega_{1})\subset\mathcal{P}_{u}(\Omega_{0})$. This implies that ${\rm d}(\theta_{0},\Omega_{0},u)=0$, meaning that we cannot do better than random guessing if the inference is done in the on-dimensional projected space. We refer to Figure 1 for a schematic illustration of this example in $p=2$.

5 Approximate sparsity

With the aim of broadening the application of our proposed method, we relax the sparsity assumption of the model to a so-called approximate sparsity structure. Consider the linear model

[TABLE]

with $w\sim{\sf N}(0,\sigma^{2}{\rm I}_{n\times n})$, and $\theta_{*}\in\mathbb{R}^{p}$ the unknown model parameters that is not necessarily sparse. However, we assume that there exists at least one sparse linear combination of the covariates that gets close to the true signal. This is formally stated as the approximate sparsity stated below, which is similar to the one introduced by [BCCH12].

Assumption 5.1**.**

(Approximately Sparse Model).* The signal $X\theta_{*}$ is well approximated by a linear combination of unknown $s_{0}\geq 1$ covariates:*

[TABLE]

The approximate sparsity assumption in [BCCH12] is weaker than the one we are imposing here, as the former allows for $\|r\|=O_{P}(\sqrt{s_{0}})$.

The next assumption is also introduced by [BCCH12], under the name of “RF condition”. This is basically an assumption on the moments of covariates and the noise component. In stating that we borrow the following empirical process notation from [BCCH12]: ${\mathbb{E}}_{n}[f]\equiv{\mathbb{E}}_{n}[f(z_{i})]\equiv\sum_{i=1}^{n}f(z_{i})/n$ and $\bar{{\mathbb{E}}}[f]\equiv{\mathbb{E}}{\mathbb{E}}_{n}[f]={\mathbb{E}}{\mathbb{E}}_{n}[f(z_{i})]=\sum_{i=1}^{n}{\mathbb{E}}[f(z_{i})]/n$.

Assumption 5.2**.**

(Moment Condition).* Suppose that the following moment conditions holds:*

• $(i)$

For a constant $C_{2}>0$, $\bar{{\mathbb{E}}}[y_{i}^{2}]+\bar{{\mathbb{E}}}[X_{ij}^{2}y_{i}^{2}]+1/\bar{{\mathbb{E}}}[X_{ij}^{2}w_{i}^{2}]\leq C_{2}$.

• $(ii)$

We have $\max_{j\in[p]}\bar{{\mathbb{E}}}[|X_{ij}^{3}w_{i}^{3}|]\leq o(\sqrt{n/(\log p)^{3}})$, and also $s_{0}\log p=o(n)$.

• $(iii)$

$\underset{i\in[n],j\in[p]}{\max}X_{ij}^{2}(s_{0}\log p)/n\to 0$, in probability and $\underset{j\in[p]}{\max}|({\mathbb{E}}_{n}-\bar{{\mathbb{E}}})[X_{ij}^{2}w_{i}^{2}]|+|({\mathbb{E}}_{n}-\bar{{\mathbb{E}}})[X_{ij}^{2}y_{i}^{2}]|\to 0$, in probability.

The above moment condition was proposed in [BCCH12] where they bound the estimate error of selection methods such as Lasso under approximate sparsity condition. Our lemma below provides a set of alternative conditions that, for sub-gaussian designs, imply the Moment condition 5.2.

Lemma 5.3**.**

Suppose that the design $X$ has independent sub-gaussian centered rows with uniformly bounded sub-gaussian norm ($\|x_{i}\|_{\psi_{2}}\leq C$). Assume that $y_{i}$ and $w_{i}$ have uniformly bounded conditional moments of order $4$, that is ${\mathbb{E}}(y_{i}^{4}|x_{i})\leq C^{\prime}$ and ${\mathbb{E}}(w_{i}^{4}|x_{i})\leq C^{\prime\prime}$, for $i\in[n]$. In addition, suppose that $s_{0}=o(n/\log^{2}(p))$ and $\log p=o(n^{1/3})$. Then the Moment Condition 5.2 holds.

We refer to Appendix A.11 for the proof of Lemma 5.3.

Iterated Lasso. Following [BCCH12], we consider a weighed Lasso estimator of $\theta_{0}$. Formally, let ${\widehat{\theta}}$ be given by

[TABLE]

where the regularization $\lambda$ is chosen as

[TABLE]

The weights $\gamma_{i}$, $j\in[p]$ are ideally chosen as $\gamma_{j}=\sqrt{{\mathbb{E}}_{n}[X_{ij}^{2}w_{i}^{2}]}$. But since the noise terms $w_{i}$ are unobserved this ideal option is not realizable. Hence, we use an iterative method proposed in [BCCH12, BCH14] to set the weights $\gamma_{i}$. (The resulting Lasso estimator ${\widehat{\theta}}$ is referred to as ‘iterated Lasso’ in [BCCH12, BCH14].) The details of the procedure is described in Algorithm 1.

Our next theorem is analogous to Theorem 3.3 and shows our procedure controls the type-I error for random designs under approximately sparse models.

Theorem 5.4**.**

Let $\Sigma\in\mathbb{R}^{p\times p}$ such that $\sigma_{\min}(\Sigma)\geq C_{\min}>0$ and $\sigma_{\max}(\Sigma)\leq C_{\max}<\infty$ and $\max_{i\in[p]}\Sigma_{ii}\leq 1$. Suppose that the regression model (41) is approximately sparse (Assumption 5.1), and assume that the responses $y_{i}$ have uniformly bounded conditional moment of order 4, that is ${\mathbb{E}}(y_{i}^{4}|x_{i})\leq C^{\prime}$ for $i\in[n]$ and a constant $C^{\prime}>0$ independent of $n$.

Let ${\widehat{\theta}}$ be the iterated Lasso estimator using data $(y,X)$, given by (43). Consider an arbitrary $U\in\mathbb{R}^{p\times k}$, with $U^{\sf T}U=I$, that is independent of the samples $\{(x_{i},y_{i})\}_{i=1}^{n}$. Construct a debiased estimator $\widehat{\gamma}^{{\rm d}}$ as in (11) with $\mu=a\sqrt{(\log p)/n}$, and $a^{2}>48e^{2}\kappa^{4}C_{\max}/C_{\min}$. In addition, suppose that $\lim\sup_{n\to\infty}\mu(\max_{i\in[k]}\,\|u_{i}\|_{1})\leq c^{\prime}$, for some constant $0<c^{\prime}<1$, $s_{0}=o(\sqrt{n}/{\log p})$ and $\log p=o(n^{1/3})$.

For the test $R_{X}$ defined in Equation (19), and for any $\alpha\in[0,1]$, we have

[TABLE]

We refer to Section 9.4 for the proof of Theorem 5.4.

6 Extension to Non-Gaussian Noise

Our analysis can be extended to the case of non-gaussian noise measurements. Specifically, suppose that the noise term $w_{i}$ satisfies

[TABLE]

for some constants $a,B>0$, and $1\leq i\leq n$.

Recall that our analysis is based on a bias-variance decomposition of the estimate $\widehat{\gamma}^{{\rm d}}$ as in Lemma 2.1. The bias term $\|\Delta\|_{\infty}$ can be bounded as

[TABLE]

The first term does not involve the noise term $w$ and can be treated as before. For bounding $\|\theta_{0}-{\widehat{\theta}}\|_{1}$, we used the result of [BCCH12, Theorem 1] (See Proposition 9.7 in the Appendix) that also applies to non-gaussian noise as long as the moment conditions (Assumption 5.2) hold, which by Lemma 5.3, for sub-gaussian designs it reduces to requiring the noise variables $w_{i}$ have bounded conditional moment of order 4.

So the remaining part is characterizing the limiting distribution of $Z$. To this end, we will show that the Lindeberg condition holds and hence $Z$ admits an asymptotically normal distribution by virtue of central limit theorem.

Similar to the approach taken in [JM14a], we slightly modify our construction of the decorrelating matrix $G$ to ensure the Lindeberg condition holds. Let ${G}=[g_{1}|\dotsc|g_{k}]\in\mathbb{R}^{p\times k}$, where each $g_{i}$ is obtained by solving the following optimization problems for each $1\leq i\leq k$:

[TABLE]

Our following proposition shows that $Z$ admits an asymptotically normal distribution in the non-gaussian setting.

Proposition 6.1**.**

Suppose that the noise variables $w_{i}$ are independent with ${\mathbb{E}}(w_{i}|X)=0$, ${\mathbb{E}}(w_{i}^{2}|X)=\sigma^{2}$ and ${\mathbb{E}}(|w_{i}|^{4+a}|X)\leq B$ for some $a>4\beta/(1-2\beta)$. Let $G=[g_{1}|\dotsc|g_{k}]\in\mathbb{R}^{p\times k}$ be the matrix constructed by solving optimization problem (47). For $i\in[p]$, define

[TABLE]

Suppose that the assumptions of Theorem 5.4 hold. Then, for any sequence $i=i(n)\in[p]$, and any $x\in\mathbb{R}$, we have

[TABLE]

with $\Phi(x)$ indicating the cdf of standard normal variable.

We refer to Appendix A.4 for the proof of Proposition 6.1.

7 Discussion

It is useful to study the proposed methodology for some specific choices of $\Omega_{0}$ and discuss its optimality.

Example 1 (Predictions). Fix an arbitrary $c\in\mathbb{R}$ and consider the set $\Omega_{0}=\{\theta:\xi^{\sf T}\theta=c\}$. This corresponds to the set where the (noiseless) unobserved response on the new feature vector $\xi$ is $c$. We can use our methodology to test $H_{0}:\theta_{0}\in\Omega_{0}$ versus its alternative. Further, by duality of hypothesis testing and confidence intervals, our methodology provides confidence intervals for a linear functional of the form $\xi^{\sf T}\theta_{0}$.

Computing $u$ from (40) in this case gives $u=\xi/\|\xi\|$. Since $\xi$ is independent of $(y,X)$, the data splitting step in the procedure becomes superfluous. By duality, we construct $(1-\alpha)$ confidence interval for $\xi^{\sf T}{\theta}_{0}$ by finding the range of values $c$ such that the rule fails to reject $H_{0}$ at level $\alpha$. This is formalized in the next lemma.

Lemma 7.1**.**

Consider a sequence of design matrices $X\in\mathbb{R}^{n\times p}$, with dimensions $n,p\to\infty$, $p=p(n)\to\infty$ satisfying the assumptions of Theorem 3.2. For given $\alpha\in(0,1)$, define $C(\alpha)=[c_{\min},c_{\max}]$ with

[TABLE]

where $\widehat{\gamma}^{{\rm d}}$ is the debiased estimator given by (33) with $u=\xi/\|\xi\|$. Then,

[TABLE]

We refer to Appendix A.5 for the proof of Lemma 7.1. The constructed confidence interval has length of rate $\|\xi\|/\sqrt{n}$. In [CG17], it is shown that the minimax expected length of confidence intervals for $\xi^{\sf T}\theta_{0}$, with a sparse vector $\xi$ (i.e., $\|\xi\|_{0}=O(s_{0})$) is $\|\xi\|(1/\sqrt{n}+s_{0}(\log p)/n)$. Therefore, in the regime $s_{0}=o(\sqrt{n}/\log p)$, which is the focus of the current paper, the constructed confidence intervals are minimax rate optimal. It is worth noting that the confidence interval defined in Lemma 7.1 is similar to the one proposed by [CG17]. For the case of non-sparse $\xi$, [CG17] establishes the minimax rate $\|\xi\|_{\infty}s_{0}\sqrt{(\log p)/n}$ for the expected length of confidence interval for $\xi^{\sf T}\theta_{0}$, and hence our construction (49) has an optimality gap in this case.

Example 2 (Quadratic forms). As another example we apply our framework to testing squared-$\ell_{2}$ norm of $\theta_{0}$. Consider the set $\Omega_{0}(c)=\{\theta:\|\theta\|_{2}^{2}=c\}$, where $c\geq 0$ is a fixed arbitrary constant. We use the proposed framework to test the null hypothesis $H_{0}:\theta_{0}\in\Omega_{0}(c)$. Computing $u$ from (40) in this case gives $u={\widehat{\theta}}^{(1)}/\|{\widehat{\theta}}^{(1)}\|$. We next use the duality between hypothesis testing and confidence intervals to construct confidence intervals for $\|\theta_{0}\|_{2}^{2}$.

Lemma 7.2**.**

Consider a sequence of design matrices $X\in\mathbb{R}^{n\times p}$, with dimensions $n,p\to\infty$, $p=p(n)\to\infty$ satisfying the assumptions of Theorem 3.3). For given $\alpha\in(0,1)$, define $C(\alpha)=[c_{\min},c_{\max}]$ with

[TABLE]

where $a_{+}=\max(a,0)$ and $\widehat{\gamma}^{{\rm d}}$ is the debiased estimator given by (33) with $u={\widehat{\theta}}^{(1)}/\|{\widehat{\theta}}^{(1)}\|$. Then,

[TABLE]

We give the proof of Lemma 7.2 in Appendix A.6.

Example 3 (Testing $\theta_{\min}$ condition). For a given $c>0$, define the set $\Omega_{0}=\{{\theta}\in\mathbb{R}^{p}:\;\min_{j\in{\rm supp}({\theta})}|{\theta}_{j}|\geq c\}$. Apart from the importance of this example as discussed in the introduction, it differs from previous example in that the set $\Omega_{0}$ is non-convex and disconnected. Recall that the guideline (40) was provided for convex sets $\Omega_{0}$, which is not true in this example.

Before proposing a choice of $U$ for this example, we state a lemma.

Lemma 7.3**.**

Let $v\in\mathbb{R}^{p}$ and define ${\theta}\in\mathbb{R}^{p}$ with ${\theta}_{i}={\cal S}(v_{i},c)$, where

[TABLE]

Then ${\theta}$ is a solution to $\min_{{\theta}\in\mathbb{R}^{p}}\|D(v-{\theta})\|_{\infty}$, subject to ${\theta}\in\Omega_{0}$, for any diagonal matrix $D$.

Proof of Lemma 7.3 is straightforward and is omitted.

In the numerical experiments, we apply our framework for this example with $k=1$ and $U=u\in\mathbb{R}^{p}$ given by:

[TABLE]

We refer to Appendix A.7 for a justification for this choice. By using Lemma 7.3, the test statistic in this case amounts to $T_{n}=|d(\widehat{\gamma}^{{\rm d}}-{\cal S}(\widehat{\gamma}^{{\rm d}},c))|$ (See step 5 of the algorithm presented in Section 4).

7.1 Prior art

The inference problem (3) studied in this paper is very general and encompasses several important problems such as the examples discussed in Section 1.1. For specific choices of set $\Omega_{0}$, one may use the structure of the set $\Omega_{0}$ to come up with methods with higher statistical power. However, in the sequel we argue that for three classes of inferential problems, our proposed framework either recovers the previously proposed methods for that specific problem, or have comparable performance. We also contrast the underlying assumptions of our framework and those of other methods designed for these specialized problems.

1. Inference on prediction: As discussed in Section 7, for inference on linear functions $\gamma_{0}=\xi^{\sf T}\theta_{0}$ (predictions), our framework proposes $u=\xi/\|\xi\|$ and construct a debiased estimator of $\gamma_{0}$ taking the following form

[TABLE]

with $g$ is obtained by solving optimization (12). As argued for the case of random designs with population covariance $\Sigma$, this implies $g\approx\Sigma^{-1}\xi/\|\xi\|$. As also discussed earlier in the introduction and previous section, a similar approach has been used by [CG17] and they prove that the resulting confidence interval would be minimax rate optimal. It is indeed an appealing property of our method that, despite its generality, it recovers the method of [CG17] for this specific case and enjoys minimax optimality.

$\bullet$ Assumptions: In terms of assumptions, [CG17] focuses on high-dimensional linear models with gaussian designs (rows of design matrix are drawn i.i.d from a multivariate normal distribution), sparse parameter vector and gaussian measurement noise. Our analysis in Section 3 considers sub-gaussian random designs (Theorem 3.3) and coherent fixed design (Theorem 3.2). We also extended our analysis to approximately sparse models (Section 5) and non-gaussian noise (Section 6).

$\bullet$ Least-favorable one-dimensional sub-model: It is worth noting that the form of debiasing (57) for linear functionals of $\theta$ can also be derived from the perspective of least-favorable scores discussed in an earlier work [ZZ14]. Akin to the semi-parametric models, consider the one-dimensional sub-model $\{\theta_{0}+u\phi,\,|\phi|<{\varepsilon}_{*}\}$ with ${\varepsilon}_{*}\to 0$, $\phi$ scalar and $u\in\mathbb{R}^{p}$. By imposing the constraint $\xi^{\sf T}u=1$, we have $\xi^{\sf T}(\theta_{0}+u\phi)-\xi^{\sf T}\theta_{0}=\phi$. The idea of [ZZ14] is to look for the least favorable submodels at $\theta_{0}$, given by $\theta_{0}+u\phi$ with $u_{0}$ the direction that minimizes Fishers information. For the log-likelihood $\ell_{i}(\theta_{0})=\ell(\theta_{0}|y_{i},x_{i})$, recall that the Fisher information operator at $\theta$ is defined as $F=-{\mathbb{E}}(\ddot{\ell_{i}}(\theta))$ and for linear regression with gaussian errors, we have $F=\tfrac{1}{\sigma}^{2}{\mathbb{E}}(x_{i}x_{i}^{\sf T})=\tfrac{1}{\sigma}^{2}\Sigma$. The least-favorable direction in the sub-model is then given by

[TABLE]

Following [ZZ14], one can construct a low-dimensional projection estimator (LDPE) as a one-step maximum likelihood correction of $\widehat{\theta}$ in the direction of the least favorable sub-model $u$ as follows

[TABLE]

where in the last step we replaced the denominator by its expectation. Comparing (58) with (57) we see that (up to a normalization by $\|\xi\|$) they are the same if $g=\Sigma^{-1}\xi$. However, $\Sigma$ is unknown in general and optimization (12) try to find $g\approx\Sigma^{-1}\xi$ that also minimizes the variance of the obtained debiased estimator.

$\bullet$ Choice of $k$ and effect of sample splitting: Our procedure uses sample splitting to find the best subspace $U$ for the sake of statistical power. On one side, the sample splitting incurs loss in power as we are using only half of data points. On the other side, the purpose of sample splitting was to choose $U$ so as to increase the power. To understand this trade-off we consider the following inference problem. Consider a function $h:\mathbb{R}^{p}\mapsto\mathbb{R}^{q}$ defined as $h(\theta)=(\xi_{1}^{\sf T}\theta,\dotsc,\xi_{q}^{\sf T}\theta)$, for a linearly independent set $\{\xi_{1},\dotsc,\xi_{q}\}$. The goal is to do inference on the value of $h(\theta_{0})$. We consider the following two methods of choosing $U$ in constructing the debiased estimator:

Method 1: We let $k=q$ and $U$ be a basis for the space spanned by $\{\xi_{1},\dotsc,\xi_{q}\}$. This method does not require any sample splitting. 2. 2.

Method 2: Define $\Omega_{0}=\{\theta:h(\theta)=c\}$, for a given $c>0$. Since $\Omega_{0}(c)$ is convex, our methodology sets $k=1$ and chooses $u$ as in (40). Here we require sample splitting for $q\geq 2$. (cf. Section 4.1)

Note that the two methods become identical for $q=1$. We next compare (the analytical lower bound on) the statistical power of these two methods for choosing $U$. Let $\eta_{u}={\rm d}(\widehat{\theta},\Omega_{0};u)$ and $\eta_{u}={\rm d}(\widehat{\theta},\Omega_{0};U)$, with $u$ given by (40) and $U$ a basis for the space $\{\xi_{1},\dotsc,\xi_{q}\}$. Using Theorem 3.4 and Equation (30), the lower bound for the power of method 1 and method 2 are respectively given by $F(\alpha,\tfrac{1}{\widehat{\sigma}}\sqrt{nC_{\min}}\eta_{U},q)$ and $F(\alpha,\frac{1}{\sqrt{2}\widehat{\sigma}}\sqrt{nC_{\min}}\eta_{u},1)$. Furthermore, by Equation (90) we have $\eta_{u}\geq\eta_{U}$ and since $F(\alpha,x,k)$ is increasing in $x$, we get $F(\alpha,\frac{1}{\sqrt{2}\widehat{\sigma}}\sqrt{nC_{\min}}\eta_{u},1)\geq F(\alpha,\frac{1}{\sqrt{2}\widehat{\sigma}}\sqrt{nC_{\min}}\eta_{U},1)$. In summary, we have

[TABLE]

The above lower bounds nicely capture the tradeoff between the choice of $k$ and the sample splitting. The function $F(\alpha,x,k)$ is decreasing in $k$ which supports the use of $k=1$, but the function is increasing in the $x$ and hence decreases under sample splitting. To understand this tradeoff we basically need to compare $F(\alpha,x,1)$ and $F(\alpha,\sqrt{2}x,q)$, with $x=\frac{1}{\sqrt{2}\widehat{\sigma}}\sqrt{nC_{\min}}\eta_{U}$. In Figure 2, we have plotted these curves for $\alpha=0.05$ and several values of $q$. As we see for small values of signal strength $x$, method 2 ($k=1$ and sample splitting) outperforms, while for larger signal strength $x$, method 1 ($k>1$ and no sample splitting) prevails.

2. Inference on quadratic forms of parameters: The work [JBC17] proposed EigenPrism, a procedure to construct two-sided confidence interval for the signal squared magnitude $\|\theta_{0}\|^{2}$. An appealing property of this procedure is that, albeit its applicability to the high-dimensional setting ($p>n$), it does not make any assumption on the coefficient sparsity. However, it is theoretically justified only for standard gaussian designs where $X_{ij}\sim{\sf N}(0,1)$, independently. As explained in [JBC17], this assumption is crucial because it ensures that in the SVD of $X=UDV^{\sf T}$, the columns of $V$ are uniformly distributed on the unit sphere, and hence allows for computing the expectation and variance of inner products of columns of $V$ with $\theta_{0}$, which constitutes a main building component of EigenPrism. By contrast, our procedure (when specialized to inference on quadratic forms of parameters as discussed in Section 7, Example 2) applies to a much broader family of sub-gaussian random designs, but assumes the coefficient sparsity $s_{0}=o(\sqrt{n}/\log p)$.

In the limit $n,p\to\infty$ and $n/p\to\gamma\in(0,1)$, the length of confidence intervals constructed by EigenPrism for $\|\theta_{0}\|^{2}$ works out at $C_{\gamma}(\|\theta_{0}\|^{2}+\sigma^{2})\frac{z_{\alpha/2}}{\sqrt{n}}$, with $C_{\gamma}$ a numerical constant defined based on Marcenko-Pastur distribution with parameter $\gamma$. By comparison, using Lemma 7.2, the confidence interval obtained by our method is of length $2L<\frac{2z_{\alpha/2}}{\sqrt{C_{\min}}}\|\widehat{\theta}^{(1)}\|\tfrac{\sigma}{\sqrt{n}}$. As we see the length of confidence intervals for $\|\theta_{0}\|^{2}$ from both methods scale at rate $1/\sqrt{n}$.

3. Inference on individual parameters: As discussed in Section 1.1, for the special case of inference on an individual model parameter, our approach recovers the debiasing method of [JM14a]. Similar debiasing approach (with different construction of the the decorrelating matrix, using node-wise regression) was proposed in [ZZ14, VdGBRD14] and its validity is proved under the assumption that the precision matrix $\Sigma^{-1}$ is sparse. The work [BCH14] has proposed a significantly different approach for doing inference on an an individual parameter, called “post-double selection”. Suppose that we are interested in parameter $\theta_{i}$. This method consists of two selection steps: 1) Let $I_{1}$ be the covariates selected by Lasso in regressing columns $i$ of the design matrix on the other columns; 2) Let $I_{2}$ denote the covariates selected by Lasso in regressing $y$ on the design $X$. The estimation of parameter $\theta_{i}$ is then defined as the least squares estimator obtained by regression $y$ on $x_{i}$ and the selected features $I_{1}\cup I_{2}$ (we may expand this set to also include other features that the statistician thinks are relevant). It is shown that the post-double estimator obeys an asymptotically normal distribution.

The limiting distribution of the post-double estimator is characterized under approximate sparsity structure and also applies to non-gaussian noise as well, as far as some moment conditions (similar to Assumption 5.2) hold. Let us stress that the approximate sparsity assumption in [BCH14] is much weaker than ours in that it allows for $\|r\|=O_{P}(\sqrt{s_{0}})$, while we require $\|r\|=o_{P}(1)$. In addition, the analysis of the post-double estimator extends to possibly heteroscedastic noise distributions.

8 Numerical illustration

In this section, we examine the performance of our inference framework in terms of coverage rate and length of confidence intervals, type I error and statistical power under different setups. We consider linear model (2) where the design matrix $X\in\mathbb{R}^{n\times p}$ has i.i.d rows generated from ${\sf N}(0,\Sigma)$, with $\Sigma\in\mathbb{R}^{p\times p}$ being the toeplitz matrix $\Sigma_{i,j}=\rho^{|i-j|}$. For coefficient parameter ${\theta}_{0}$, we consider a uniformly random support (set of nonzero parameters) $S\subseteq[p]$, with $|S|=s_{0}$. The measurement errors are $w_{i}\sim{\sf N}(0,1)$.

8.1 Testing $\theta_{\min}$ condition

We consider the set $\Omega_{0}=\{{\theta}:\min_{j\in{\rm supp}({\theta}_{0})}|{\theta}_{0,j}|\geq c\}$ and the null hypothesis $H_{0}:\,\theta_{0}\in\Omega_{0}$. As explained in Section 7 (Example 3), the set $\Omega_{0}$ is non-convex (indeed disconnected) and we consider one-dimensional projection of the problems along the direction $u$ given by (56) for this example. For the scaled Lasso estimator $\widehat{\theta}^{n}$, given by (10), we set the regularization parameter $\lambda=\sqrt{2.05(\log p)/n}$. Further, the parameter $\mu$ in constructing the debiased estimator (see optimization problem (12)) is set to $\mu=2\sqrt{(\log p)/n}$. We set $p=1000$, $n=600$, $s_{0}=10$. The nonzero parameters $\theta_{0,i}$, $i\in S$, are chosen as $0.1,0.2,\dotsc,1$. We set $\alpha=0.05$ and vary the values of $c$ and $\rho$. The rejection probabilities are computed based on 300 random samples for each value of pair $(c,\rho)$. When $c\leq 0.1$, $H_{0}$ holds and thus the rejection probability corresponds to the type I error. When $c>0.1$, the rejection probability corresponds to the power of the test. The results are reported in Table 1. As we see in Table 1(a), type I error is controlled below the desired level $\alpha=0.05$. Also, as evident in Table 1(b), the power of our test increases at a very fast rate as $c$ increases.

8.2 Confidence intervals for linear functions

We use our methodology to construct $95\%$ confidence intervals for functions of the form $\xi^{\sf T}{\theta}_{0}$. We set $p=3000$, $s_{0}=30$ and choose the correlation parameter $\rho=0.5$. The model parameters are set as follows. We set $\theta_{0,j}=0.5$ for $j=1,\dotsc,s_{0}$, and $\theta_{0,j}=0.5/(j-s_{0}+1)$, for $j=s_{0}+1,\dotsc,p$.

We construct confidence intervals according to Lemma (7.1). We choose fives vectors $\xi_{1},\xi_{2},\dotsc,\xi_{5}$ as eigenvectors of $\Sigma$ with well-separated eigenvalues. Specifically, sorting the eigenvalues of $\Sigma$ as $\sigma_{1}\geq\sigma_{2}\geq\dotsc\geq\sigma_{3000}$, we choose the eigenvectors corresponding to $\sigma_{1}$, $\sigma_{750}$, $\sigma_{1500}$, $\sigma_{2250}$, $\sigma_{3000}$. For each $\xi_{i}$, we vary $n$ in $\{1000,1200,1400,\dotsc,2600\}$. For each configuration $(\xi_{i},n)$, we consider $300$ independent realizations of measurement noise and on each realization, we construct $95\%$ confidence interval for $\xi_{i}^{\sf T}{\theta}_{0}$ based on Lemma (7.1).

In Figure 3(a), we plot the average coverage probability of constructed confidence intervals for each configuration. Each curve corresponds to one of the vectors $\xi_{i}$. As we see, the coverage probability for all of them and across different values of $n$ is close to the nominal value.

In Figure 3(b), we plot the average length of confidence intervals as we vary the sample size $n$ in the log-log scale. As evident from the figure, the length of confidence intervals scales as $1/\sqrt{n}$.

8.3 Testing for the non-negative cone

Define \Omega_{0}=\{\theta:\theta_{i}\geq 0\text{ for all i}\} as the non-negative cone. In this section, we test whether $\theta_{0}\in\Omega_{0}$ versus $\theta_{0}\notin\Omega_{0}$. The null model is generated as follows. The nonzero entries in support $S$ are chosen as $b,b/2,b/3,\dotsc,b/s_{0}$, where $s_{0}=|S|$ and $b>0$. The entries outside $S$ are set to zero. The alternative model is generated similar where $b$ is replaced by $-b$. As in the previous sections, the design matrix $X\in\mathbb{R}^{n\times p}$ has i.i.d rows generated from ${\sf N}(0,\Sigma)$, with $\Sigma\in\mathbb{R}^{p\times p}$ being the toeplitz matrix $\Sigma_{i,j}=\rho^{|i-j|}$, and measurement errors $w_{i}\sim{\sf N}(0,1)$, with parameters $(n,s_{0},p)=(600,10,1000)$. We set $\alpha=0.05$ and vary the values of $b$ and $\rho$. The rejection probabilities are computed based on 300 random samples for each value of pair $(b,\rho)$.

The simulation report in Table 2 shows that the type I error is controlled below the target level $\alpha=0.05$. Per statistical power, the method achieves power at least $99\%$ for $|b|\geq 0.4$. Note that we have a very difficult alternative in the sense that only a small fraction of the coordinates $({s_{0}}/{d})$ is negative with small magnitudes ranging in $[b/10,b]$, so it is a very mild violation of the null, yet our algorithm still has high power.

8.4 Real data experiment

We measure the performance of our testing procedure on a riboflavin data set, which is publicly available by [BKM14] and can be downloaded via the ‘${{\sf hdi}}$’ R-package. The data set includes $p=4088$ predictors corresponding to the genes and $n=71$ samples. The response variable indicates the logarithm of the riboflavin production rate and the covariates are the logarithm of the expression levels of the genes. We model the riboflavin production rate by a linear model. We first fit the Lasso solution ${\widehat{\theta}}$ using the glmnet package [FHT10] and then generate $N=100$ instances of the problem as $y^{(i)}=X{\widehat{\theta}}+w^{(i)}$, where $w^{(i)}\sim{\sf N}(0,\sigma^{2}{\rm I}_{n})$. In other words, we treat ${\widehat{\theta}}$ as the true parameter $\theta_{0}$ and generate new data by resampling the noise.

We run two sets of experiments on this data.

CI for predictions. We fix a vector $\xi\in\mathbb{R}^{p}$ that is generated as $\xi_{i}\sim{\sf N}(0,1/\sqrt{p})$, independently for $i\in[p]$. On each problem instance $(i)$, we construct confidence interval ${{\rm CI}}^{(i)}$ for $\xi^{T}\theta_{0}$, using Lemma 7.1. We compute the coverage rate as

[TABLE]

CI for squared norm. On each problem instance $(i)$, we construct confidence interval for $\|\theta_{0}\|_{2}^{2}$, using Lemma 7.2 and compute the coverage rate given by (60).

The results are reported in Table 3. As we see for various values of noise standard deviation $\sigma$, the coverage rates of the constructed intervals remain close to the nominal value. In Figure 4, we depict the constructed confidence intervals for $40$ random problem instances, in each experiment.

9 Proof of Theorems

9.1 Proof of Theorem 3.2

We first prove a lemma to bound the estimation error of $\widehat{\sigma}$ returned by the scaled Lasso. The following lemma uses the analysis of [SZ12] and its proof is given in Appendix A.8 for reader’s convenience.

Lemma 9.1**.**

Under the assumptions of Theorem 3.2, let $\widehat{\sigma}=\widehat{\sigma}(\lambda)$ be the scaled Lasso estimator of the noise level, with $\lambda=c\sqrt{(\log p)/n}$ and define $\sigma_{*}=\|w\|/\sqrt{n}$. Then, $\widehat{\sigma}$ satisfies

[TABLE]

Armed with Lemmas 9.1 and 2.1 we are ready to prove Theorem 3.2. Under $H_{0}$, we have ${\theta}_{0}\in\Omega_{0}$ and hence by invoking Lemma 2.1, we have

[TABLE]

Note that for $\tilde{Z}\equiv\widehat{\sigma}DZ/(\sigma\sqrt{n})\in\mathbb{R}^{k}$, we have $\tilde{Z}_{i}\sim{\sf N}(0,1)$. The entries of $\tilde{Z}$ are correlated though.

Fix $\epsilon>0$ and apply Equation (62) to write

[TABLE]

For the second term, we proceed as follows

[TABLE]

Now, note that $\sigma^{*}\to\sigma$, in probability, as $n$ tends to infinity. Therefore, by applying Lemma (9.1) and using the assumption $s_{0}=o(n/\log p)$, we get

[TABLE]

Using this in (9.1), we have

[TABLE]

We next note that by definition (16), and using the assumption $\lim\inf_{n\to\infty}\min_{i\in[k]}(G^{\sf T}\widehat{\Sigma}G)_{ii}\geq c_{0}>0$, we have from which we obtain

[TABLE]

By Equation (65), we have ${\mathbb{P}}((\sigma/\widehat{\sigma})\geq 2)\to 0$. In addition, since $s_{0}=o(1/(\mu\sqrt{\log p}))$, for $n$ and $p$ large enough, we have $c\mu s_{0}\sqrt{\log p}/\phi_{0}^{2}\leq\epsilon\sqrt{c_{0}}/2$. Hence by (14),

[TABLE]

By substituting (68) in (9.1), we get

[TABLE]

By union bounding over the entries of $\tilde{Z}$, we get

[TABLE]

Observe that the above holds for any $\epsilon>0$, and that the right-hand side is bounded pointwise for all $\epsilon$. Therefore, by applying the dominated convergence theorem, we get

[TABLE]

The result follows by choosing $x=\Phi^{-1}(1-\alpha/(2k))$.

9.2 Proof of Theorem 3.3

For $\phi_{0},s_{0},K\geq 0$, let $\mathcal{E}_{n}=\mathcal{E}_{n}(\phi_{0},s_{0},K)$ be the event that the compatibility condition holds for $\widehat{\Sigma}=(X^{{\sf T}}X/n)$, for all sets $S\subseteq[p]$, $|S|\leq s_{0}$ with constant $\phi_{0}>0$, and that $\max_{i\in[p]}\,\widehat{\Sigma}_{i,i}\leq K$. Explicitly

[TABLE]

Then, by result of [RZ13, Theorem 6] (see also [JM14a, Theorem 2.4(a)]), random designs satisfy the compatibility condition with constant $\phi_{0}=\sqrt{C_{\min}}/2$, provided that $n\geq\nu s_{0}\log(p/s_{0})$, where $\nu=c\kappa^{4}(C_{\max}/C_{\min})$, for a constant $c>0$. More precisely,

[TABLE]

where $c_{1}=c_{1}(c)>0$ is a constant.

We next provide an explicit upper bound for the minimum generalized coherence $\mu_{\min}(X;U)$ (cf. Definition 3.1) for random designs.

Proposition 9.2** ( [JM14a]).**

Let $\Sigma\in\mathbb{R}^{p\times p}$ be such that $\sigma_{\min}(\Sigma)\geq C_{\min}>0$ and $\sigma_{\max}(\Sigma)\leq C_{\max}<\infty$ and $\max_{i\in[p]}\Sigma_{ii}\leq 1$. Suppose that $X\Sigma^{-1/2}$ has independent sub-gaussian rows, with mean zero and sub-gaussian norm $\|\Sigma^{-1/2}x_{1}\|_{\psi_{2}}=\kappa$, for some constant $\kappa>0$. For $U\in\mathbb{R}^{p\times k}$ independent of $X$ satisfying $U^{\sf T}U=I$, and for fixed constant $a>0$, define

[TABLE]

In other words, $\mathcal{G}_{n}(a)$ is the event that problem (12) is feasible for $\mu=a\sqrt{(\log p)/n}$. Then, for $n\geq a^{2}C_{\min}\log p/(4e^{2}C_{\max}\kappa^{4})$, the following holds true with high probability

[TABLE]

We refer to Appendix A.9 for the proof of Proposition 9.2.

The last step is to prove that Assumption 2.3 holds. In doing that, we use Lemma 2.4. Note that the first condition of this lemma holds by assumption of the theorem. To prove the second condition, we use the following result.

Lemma 9.3**.**

Let $\Sigma\in\mathbb{R}^{p\times p}$ such that $\sigma_{\min}(\Sigma)\geq C_{\min}>0$ and $\sigma_{\max}(\Sigma)\leq C_{\max}<\infty$ and $\max_{i\in[p]}\Sigma_{ii}\leq 1$. Suppose that $X\Sigma^{-1/2}$ has independent sub-gaussian rows, with mean zero and sub-gaussian norm $\|\Sigma^{-1/2}x_{1}\|_{\psi_{2}}=\kappa$, for some constant $\kappa>0$. Let $\widehat{\Sigma}\equiv(X^{\sf T}X)/n$. For $u_{i}\in\mathbb{R}^{p}$ independent of $X$, we have

[TABLE]

for a constant $C>0$ depending on $\kappa$, $C_{\max}$, and $c>2$ depending on $C$.

We refer to Appendix A.2 for the proof of Lemma 9.3. The second condition of Lemma 2.4 follows from $u_{i}^{\sf T}\Sigma u_{i}\leq C_{\max}\|u_{i}\|^{2}=C_{\max}$, union bounding over $i\in[k]$ and Lemma 9.3 (along with Borel-Cantelli Lemma).

Putting the three probabilistic bounds (72), (74) and (75) together in Theorem 3.2, we obtain that for random designs with $s_{0}=o(\sqrt{n}/(\log p))$, we have $\underset{n\to\infty}{\lim\sup}\;\;\alpha_{n}(R_{X})\leq\alpha$.

9.3 Proof of Theorem 3.4

We start by stating a lemma that will be used later in the proof.

Lemma 9.4**.**

Under the assumptions of Theorem 3.3, for any $i\in[k]$ we have

[TABLE]

where $c$ is a constant depending on $a,C$ and by a suitable choice of them, we have $c\geq 2$.

We refer to Appendix A.10 for the proof of Lemma 9.4.

Corollary 9.5**.**

Assuming the setting of Theorem 3.3, by an application of Borel-Cantelli lemma and using Lemma 9.4, of any $i\in[k]$ we have almost surely

[TABLE]

Recalling the definition of $m_{0}$, given by (28), we have the following corollary.

Corollary 9.6**.**

Recalling the definition of $m_{0}$ given by (28), for any $i\in[k]$, we have almost surely

[TABLE]

Let $z_{*}\equiv\Phi^{-1}(1-\alpha/(2k))$ and write

[TABLE]

We define the shorthands $v\equiv DU^{\sf T}(\theta^{\rm p}-{\theta}_{0})$ and $\tilde{v}\equiv D(\widehat{\gamma}^{{\rm d}}-U^{\sf T}{\theta}_{0})$. Note that $v,\tilde{v}\in\mathbb{R}^{k}$. We further let $i^{\star}\equiv\arg\max_{i\in[k]}|v_{i}|$. Then, we can write

[TABLE]

By a very similar argument we used to derive Equation (69), we can show that for any fixed $i\in[k]$ and all $x\in\mathbb{R}$, we have

[TABLE]

In words, each coordinate of $\tilde{v}$ asymptotically admits a standard normal distribution.

The other remark we want to make is about the quantity $\|{v}\|_{\infty}$, which will be a key factor in determining the power of the test. Because $\theta^{\rm p}\in\Omega_{0}$, we have

[TABLE]

Continuing with (78), we write

[TABLE]

where $(a)$ follows from Equations (79) and (81); $(b)$ holds because of Corollary 9.6 and Equation (80). Here $Z$ is a standard normal variable; $(c)$ follows by substituting for $z_{*}$.

9.4 Proof of Theorem 5.4

The proof goes along the same lines of the proof of Theorem 3.2 and 3.3.

Defining $r=X\theta_{*}-X\theta_{0}$ and by plugging in for $y=X\theta_{*}+w=X\theta_{0}+r+w$ in the definition (11), we get

[TABLE]

with

[TABLE]

Sine $w\sim{\sf N}(0,\sigma^{2}{\rm I}_{n\times n})$, we have $Z|X\sim{\sf N}(0,\sigma^{2}G^{\sf T}\widehat{\Sigma}G)$. We next bound $\|\Delta\|_{\infty}$.

It is straightforward to see that the assumptions of Theorem 5.4 implies the assumption of Lemma 5.3 and hence by the result of the lemma, the moment conditions (Assumption 5.2) hold. To deal with $\Delta_{1}$, we use the following result from [BCCH12] that bounds the $\ell_{1}$ error of the iterated Lasso estimator under the Assumptions 5.1 and 5.2.

Proposition 9.7**.**

([BCCH12, Theorem 1]) Suppose that in the regression model (41), Assumption 5.1 (approximate sparsity) and Assumption 5.2 (Moment Conditions) hold. Let ${\widehat{\theta}}$ be the iterated lasso estimator (43) with weights $\gamma_{j}$ specified by Algorithm 44. Then, ${\widehat{\theta}}$ satisfies

[TABLE]

with high probability, for some finite constant $C>0$.

Now let ${\cal E}_{n}$ be the probability event that $\|{\widehat{\theta}}-\theta_{0}\|_{1}\leq CC_{\min}^{-1}s_{0}\sqrt{(\log p)/n}$. Recall the event $\mathcal{G}_{n}(a)$ from (73) and define ${\cal F}_{n}\equiv\mathcal{G}_{n}(a)\cap{\cal E}_{n}$. Then, by using Propositions 9.2 and 9.7, we have the ${\cal F}_{n}$ happens with high probability. Further, on the event ${\cal F}_{n}$ we have

[TABLE]

We next bound $\Delta_{2}$. Write

[TABLE]

Using lemma 9.4, we have

[TABLE]

with $C^{\prime}=1/C_{\min}+C$, and with probability at least $1-2p^{-c}$, for $c\geq 2$. By union bounding over $i\in[k]$, we get

[TABLE]

with probability at least $1-2kp^{-c}\geq 1-2p^{-c+1}$. Using Assumption 5.1, $\|r\|=o_{P}(1)$, which gives us

[TABLE]

Combining (85) and (86), we have

[TABLE]

Hence $\|\Delta\|_{\infty}=o_{p}(1)$ and $Z|X$ is asymptotically normally distributed. Having this result, we can then follows the lines of the proof of Theorem 3.3 to show that our procedure controls the type I error, that is $\lim\sup_{n\to\infty}\alpha_{n}(R_{X})\leq\alpha$.

Acknowledgements

A. Javanmard was partially supported by an Outlier Research in Business (iORB) grant from the USC Marshall School of Business, a Google Faculty Research award and the NSF CAREER award DMS-1844481.

Appendix A Proof of Technical Lemmas

A.1 Proof of Lemma 2.4

We start by providing a non-asymptotic lower bound on $(G^{\sf T}\widehat{\Sigma}G)_{i,i}$.

Lemma A.1**.**

Let $G$ be the matrix with rows $g_{i}^{\sf T}$ obtained by solving optimization (12). Then, for all $i\in[p]$,

[TABLE]

Using this lemma we write

[TABLE]

which completes the proof.

A.1.1 Proof of Lemma A.1

The proof proceeds as the proof of [JM14a, Lemma 3.1]. Let $C_{i}(\mu)$ be the solution of optimization (12). We write

[TABLE]

Hence, for feasible $\tilde{g}$ and any $c\geq 0$, and by using that $\|u_{i}\|=1$ for $i\in[k]$,

[TABLE]

Then by minimizing over all feasible $\tilde{g}$,

[TABLE]

The minimum of the right hand side is achieved for $g=cu_{i}/2$ which implies that

[TABLE]

The claim follows by optimizing over $c\geq 0$.

A.2 Proof of Lemma 9.3

Fix $i\in[k]$ and write $u_{i}^{\sf T}\widehat{\Sigma}u_{i}=\frac{1}{n}\sum_{\ell=1}^{n}(u_{i}^{\sf T}x_{\ell})^{2}$. Let $V_{\ell}=u_{i}^{\sf T}x_{\ell}$ then ${\mathbb{E}}(V_{\ell}^{2})=u_{i}^{\sf T}\Sigma u_{i}$. Further, using the sub-gaussian assumption on the covariates $x_{i}$, we have

[TABLE]

Let $S_{\ell}=V_{\ell}^{2}-u_{i}^{\sf T}\Sigma u_{i}$. Then $S_{\ell}$ is zero mean and its sub-exponential norm can be bounded as $\|S_{\ell}\|_{\psi_{1}}\leq 2\|V_{\ell}^{2}\|_{\psi_{1}}\leq 2\|V_{\ell}^{2}\|_{\psi_{1}}\leq 4\|V_{\ell}\|_{\psi_{2}}^{2}\leq 4\kappa^{2}C_{\max}\equiv C^{\prime}$. Therefore, by an application of Bernstein inequality for centered sub-exponential random variables [Ver12] (similar to the proof of Lemma A.3), we have that for ${\varepsilon}\leq eC^{\prime}$ ,

[TABLE]

For ${\varepsilon}=C\sqrt{(\log p)/n}$ and assuming $n\geq[C/(eC^{\prime})]^{2}\log p$, we obtain

[TABLE]

The result follows.

A.3 Proof of Lemma 4.1

Consider the following two optimization problems:

[TABLE]

Let ${{\sf OPT}}_{1}$ and ${{\sf OPT}}_{2}$ respectively denote the optimal value of problems $({{\rm P}}_{1})$ and $({{\rm P}}_{2})$. Clearly ${{\sf OPT}}_{1}\geq{{\sf OPT}}_{2}$. We next show the reverse side.

First note that

[TABLE]

Since the right-hand side is linear in $v$ and $\theta$, and $\Omega_{0}$ is convex, by Von Neumann’s minimax theorem, we have

[TABLE]

Let $\tilde{v}=Uv$. Since $U$ has orthonormal columns we have $\|\tilde{v}\|_{2}=\|v\|_{2}\leq\|v\|_{1}\leq 1$. Using this observation along with Equations (87) and (88), we get

[TABLE]

Therefore, for any $U\in\mathcal{J}$, there exists unit norm vector $u\in\mathbb{R}^{p}$, such that

[TABLE]

Before we proceed with the rest of the proof we state a lemma about the function $G$.

Lemma A.2**.**

The function $k\mapsto{F}(\alpha,x,k)$ is strictly decreasing in $k$.

Now choose any $U\in\mathcal{J}$ and choose unit norm $u$ that satisfies (90). Then,

[TABLE]

where the first inequality follows from monotonicity of ${F}(\alpha,x,k)$ in $x$ and the second inequality follow from Lemma A.2. This implies that ${{\sf OPT}}_{1}\leq{{\sf OPT}}_{2}$.

Therefore ${{\sf OPT}}_{1}={{\sf OPT}}_{2}$ which completes the proof. Indeed, we have proved a stronger claim that $\mathcal{J}$ only includes one-dimensional subspaces ($k=1$). This follows readily from the above proof and the fact that ${F}(\alpha,x,k)$ is strictly decreasing in $k$ as per Lemma A.2.

A.3.1 Proof of Lemma A.2

Recall the definition of ${F}$ given by

[TABLE]

Let $z=\Phi^{-1}(1-\alpha/(2y))$. We then have

[TABLE]

where $\varphi(t)\equiv e^{-t^{2}/2}{\rm d}t/\sqrt{2\pi}$ is the standard normal density function. Since $z>0$ and $\Phi$ is monotone increasing, it is easy to see that $(\partial/\partial y){F}(\alpha,x,y)<0$ for $y>0$.

A.4 Proof of Proposition 6.1

Write

[TABLE]

Note that conditional on $X$, the random variables $\zeta_{\ell}$ are zero mean and independent. In addition, $\sum_{\ell=1}^{n}{\mathbb{E}}(\zeta_{\ell}^{2}|X)=n$. Let $c_{n}=(g_{i}^{\sf T}\widehat{\Sigma}g_{i})^{1/2}$. Similar to the proof of Theorem 3.3, by using Lemma 9.3 and 2.4, Assumption 2.3 holds and hence

[TABLE]

for some positive constant $c_{0}$. We are now ready to prove that the Lindeberg condition holds. If optimization (47) is feasible for $i\in[k]$, then $|\zeta_{\ell}|\leq(\sigma c_{n})^{-1}\|Xg_{i}\|_{\infty}\|w\|_{\infty}\leq(\sigma c_{n})^{-1}n^{\beta}\|w\|_{\infty}$. Therefore,

[TABLE]

where the last limit follows since $a>4\beta/(1-2\beta)$ and ${\mathbb{E}}(|w_{1}|^{4+a})$ is finite.

What we are left to prove is that optimization (47) is feasible for all $i\in[k]$, with high probability. This follows by showing that $\Sigma^{-1}u_{i}$ is a feasible solution to (47) using the tail bound inequality for sub-gaussian variables. This is very similar to the argument presented in the proof of [JM14a, Lemma 6.3] and is omitted here.

A.5 Proof of Lemma 7.1

By computing $u$ from (40) in case of $\Omega_{0}=\{{\theta}:\,\langle\xi,{\theta}\rangle=c\}$, we have $u=\xi/\|\xi\|$. Let $q=(\widehat{\sigma}^{2}/n)(g^{\sf T}\widehat{\Sigma}^{(2)}g+10^{-4})$ and $d=q^{-1/2}$. Then, the test statistics (18) becomes

[TABLE]

because ${\theta^{{\rm p}}}\in\Omega_{0}$.

By duality of hypothesis testing and confidence intervals, the $(1-\alpha)$ confidence interval of $\langle\xi,{\theta}_{0}\rangle$, denoted by $C(\alpha)$, consists of all values $c$ such that we fail to reject $H_{0}$ at level $\alpha$. Namely, $C(\alpha)=[c_{\min},c_{\max}]$ such that $c\in C(\alpha)$ if and only if $T_{n}<z_{\alpha/2}$. Plugging for $d$ this yields

[TABLE]

The proof is complete.

A.6 Proof of Lemma 7.2

For $\phi_{0},s_{0},K\in\mathbb{R}_{\geq 0}$, define the set $\mathcal{E}_{n}(\phi_{0},s_{0},K)$ as follows:

[TABLE]

By using [JM14a, Theorem 2.4 (a)], there exists constant $c_{*}\leq 2000$, such that for $n\geq Cs_{0}\log(p/s_{0})$, $C=4c_{*}(C_{\max}\kappa^{4}/C_{\min})$ and $\phi_{0}=C_{\min}^{1/2}$, $K\geq 1+20\kappa^{2}\sqrt{(\log p)/n}$, we have

[TABLE]

Define the set $\Omega_{1}\equiv\{\theta\in\mathbb{R}^{p}:\,\|{\widehat{\theta}}^{(1)}-\theta\|_{2}\leq\frac{\sqrt{20s_{0}}}{\phi_{0}^{2}}\lambda\}$. Using the result of [BvdG11, Lemma 6.10], we have that for $\lambda\geq 8\sigma\sqrt{K(1+c_{0})(\log p)/n}$,

[TABLE]

(Note that ${\widehat{\theta}}^{(1)}$ is computed based on $n/2$ samples.)

Let $\mathcal{G}$ be the probability event that both of the above high probability events hold, that is $\mathcal{G}\equiv\{X\in\mathcal{E}_{n}(\phi_{0},s_{0},K)\}\cap\{\theta_{0}\in\Omega_{1}\}$. Therefore, ${\mathbb{P}}(\mathcal{G})\geq 1-4e^{-c_{1}n}-2p^{-c_{0}}$.

Assuming $\mathcal{G}$, we rewrite the $\ell_{\infty}$ projection in (34) for this case with $u={\widehat{\theta}}^{(1)}/\|{\widehat{\theta}}^{(1)}\|$ and $k=1$, as follows:

[TABLE]

and the test statistics is given by $T_{n}=|d(\widehat{\gamma}^{{\rm d}}-({\widehat{\theta}}^{(1)})^{\sf T}\theta^{\rm p}/\|{\widehat{\theta}}^{(1)}\|)|$. By duality of hypothesis testing and confidence intervals, we need to find the range of values of $c$, such that $T_{n}\leq z_{\alpha/2}$ (i.e, the test rule fails to reject the null hypothesis). Note that $T_{n}\leq z_{\alpha/2}$ if and only if

[TABLE]

Writing $({\widehat{\theta}}^{(1)})^{\sf T}{\theta^{{\rm p}}}=\tfrac{1}{2}(\|{\theta^{{\rm p}}}\|^{2}+\|{\widehat{\theta}}^{(1)}\|^{2}-\|{\theta^{{\rm p}}}-{\widehat{\theta}}^{(1)}\|^{2})$ and using that fact that ${\theta^{{\rm p}}}\in\Omega_{0}(c)\cap\Omega_{1}$, the above inequality yields

[TABLE]

with $C\equiv\frac{6}{4}0\sigma^{2}K(1+c_{0})\phi_{0}^{-4}$. Rearranging the terms and substituting for $d$, we get $c\in C(\alpha)=[c_{\min},c_{\max}]$, where $c_{\min}$ and $c_{\max}$ are given by

[TABLE]

with $L$ given by (53). As shown in the proof of Theorem 3.2, Assumption 2.3 holds which implies that $L\gtrsim 1/\sqrt{n}$. In addition, by our assumption $s_{0}=o(\sqrt{n}/\log p)$, which results in $\frac{s_{0}\log p}{n}=o(L)$.

Regarding the coverage probability, we use the duality of hypothesis testing and confidence intervals to obtain

[TABLE]

Hence,

[TABLE]

A.7 Choice of $U$ for testing $\theta_{\min}$ condition

Here we provide a justification for the choice of $U$, given by (56), for testing $\theta_{\min}$ condition. Recall that in this case $\Omega_{0}=\{\theta\in\mathbb{R}^{p}:\min_{j\in{\rm supp}(\theta)}|{\theta}_{j}|\geq c\}$. Instead of directly solving optimization (31), which is hard due to non-convexity of $\Omega_{0}$, we first develop a lower bound and find $U$ that maximizes the lower bound.

The lower bound is obtained by fixing $k=1$ in the optimization (31). The problem then amounts to

[TABLE]

which by plugging in for ${\rm d}({\widehat{\theta}}^{(1)},\Omega_{0};u)$ is equivalent to

[TABLE]

We claim that the optimal $u$ should be one of the standard basis element. To see this, consider $u\neq e_{i}$, for $i\in[p]$. Then, there exists a vector $v\in\mathbb{R}^{p}$ such that $v_{j}\neq 0$ for all $j\in[p]$ and $v^{\sf T}u=0$. Choose $\lambda\in\mathbb{R}$ large enough such that all the coordinates of ${\theta}={\widehat{\theta}}^{(1)}+\lambda v$ have magnitude larger than $c$. Therefore, ${\theta}\in\Omega_{0}$ and $u^{\sf T}(\theta-{\widehat{\theta}}^{(1)})=0$.

Setting $u=e_{i}$, the objective becomes $\inf_{\theta\in\Omega_{0}}|\theta_{i}-{\widehat{\theta}}^{(1)}_{i}|=|{\cal S}({\widehat{\theta}}_{i}^{(1)},c)-{\widehat{\theta}}_{i}^{(1)}|$, by Lemma 7.3. Therefore, the optimal value of objective is achieved for $i=i^{\star}$ given by (56).

A.8 Proof of Lemma 9.1

We apply [SZ12, Theorem 1], where using their notation with their $\lambda_{0}$ replaced by $\lambda$, $\xi=3$, $T={\rm supp}({\theta}_{0})$, $\kappa(\xi,T)\geq\phi_{0}$, $\eta_{*}(\sigma^{*}\lambda,\xi)\leq 4s_{0}\lambda^{2}/\phi_{0}^{2}$. By a straightforward manipulation of Eq. (13) in [SZ12], we have for $\|X^{\sf T}w/(n\sigma^{*})\|_{\infty}\leq\lambda/2$,

[TABLE]

Note that

[TABLE]

We define $v_{j}=w^{\sf T}Xe_{j}/(\sqrt{n}\sigma)$. Since $v_{j}\sim{\sf N}(0,\widehat{\Sigma}_{jj})$ by applying a standard tail bound on the supremum of $p$ gaussian random variables, we get

[TABLE]

For the second term, note that

[TABLE]

with $Z_{j}\sim{\sf N}(0,1)$ independent. By a standard tail bound for $\chi^{2}$ random variables we have

[TABLE]

Combining (99), (100) in (98), we get that

[TABLE]

which yields the desired result.

A.9 Proof of Proposition 9.2

Note that by Definition 3.1, clearly

[TABLE]

Therefore the statement follows readily from the following lemma.

Lemma A.3**.**

Consider a random design matrix $X\in\mathbb{R}^{n\times p}$, with i.i.d. rows having mean zero and population covariance $\Sigma$. Assume that

• $(i)$

We have $\sigma_{\min}(\Sigma)\geq C_{\min}>0$, and $\sigma_{\max}(\Sigma)\leq C_{\max}<\infty$.

• $(ii)$

The rows of $X\Sigma^{-1/2}$ are sub-gaussian with $\kappa=\|\Sigma^{-1/2}x_{1}\|_{\psi_{2}}$.

Let $\widehat{\Sigma}=(X^{\sf T}X)/n$ be the empirical covariance. Then, for any fixed $U\in\mathbb{R}^{p\times k}$ independent of $X$ satisfying $U^{\sf T}U=I$, and for any fixed constant $a>0$, the following holds true

[TABLE]

with $c_{2}=(a^{2}C_{\min})/(24e^{2}\kappa^{4}C_{\max})-2$.

Proof of Lemma A.3.

The proof is an application of the Bernstein-type inequality for sub-exponential random variables [Ver12]. Define $\tilde{x}_{\ell}=\Sigma^{-1/2}x_{\ell}$, for $\ell\in[n]$, and write

[TABLE]

Fix $i,j\in[p]$, and for $\ell\in[n]$, let $v^{(ij)}_{\ell}=(e_{i}^{\sf T}\Sigma^{1/2}\tilde{x}_{\ell})(\tilde{x}_{\ell}^{\sf T}\Sigma^{-1/2}u_{j})-u_{j,i}$, where $u_{j,i}$ denotes the $i$-th component of $u_{j}$. Notice that ${\mathbb{E}}(v^{(ij)}_{\ell})=0$, and the $v^{(ij)}_{\ell}$ are independent for $\ell\in[n]$, since $U$ is independent of $X$. In addition, $H_{i,j}=(1/n)\sum_{\ell=1}^{n}v^{(ij)}_{\ell}$. By [Ver12, Remark 5.18], we have

[TABLE]

Moreover, for any two random variables $X$ and $Y$, we have

[TABLE]

Hence, by assumption $(ii)$, we obtain

[TABLE]

Define $\kappa^{\prime}\equiv 2\sqrt{C_{\max}/C_{\min}}\kappa^{2}$. We now use the Bernstein-type inequality for centered sub-exponential random variables [Ver12] to get

[TABLE]

Choosing ${\varepsilon}=a\sqrt{(\log p)/n}$, and assuming $n\geq[a/(e\kappa^{\prime})]^{2}\log p$, we arrive at

[TABLE]

The result follows by union bounding over all possible pairs $i,j\in[p]$. ∎

A.10 Proof of Lemma 9.4

Define the event

[TABLE]

In other words, ${\mathcal{H}}_{n}(a)$ is the event that $\Sigma^{-1}u_{i}$ is a feasible solution of (12), for $1\leq i\leq k$. By Lemma A.3, ${\mathbb{P}}({\mathcal{H}}_{n}(a))\geq 1-2p^{-c_{2}}$. On this event, letting $g_{i}$ be the solution of the optimization problem (12), we have

[TABLE]

where $V_{j}=u_{i}^{\sf T}\Sigma^{-1}x_{j}$ are i.i.d. random variables with ${\mathbb{E}}(V_{j}^{2})=u_{i}^{\sf T}\Sigma^{-1}u_{i}$ and sub-gaussian norm

[TABLE]

Letting $S_{j}=V_{j}^{2}-u_{i}^{\sf T}\Sigma^{-1}u_{i}$, we have that $S_{j}$ is zero mean and sub-exponential with $\|S_{j}\|_{\psi_{1}}\leq 2\|V_{j}^{2}\|_{\psi_{1}}\leq 4\|V_{j}\|_{\psi_{2}}^{2}\leq 4\kappa^{2}C_{\rm min}^{-1}\equiv\kappa^{\prime}$. Hence, by applying Bernstein inequality for centered sub-exponential random variables [Ver12] (similar to the proof of Lemma A.3), we have, for ${\varepsilon}\leq e\kappa^{\prime}$,

[TABLE]

We can make $c_{2}\geq 2$ by a suitable choice of $a$. The result then follows by letting ${\varepsilon}=e\kappa^{\prime}\sqrt{6c_{2}(\log p)/n}$.

A.11 Proof of Lemma 5.3

By definition of sub-gaussian norm, given by (7), we have ${\mathbb{E}}(|X_{ij}|^{q})\leq C^{q}q^{q/2}$, for all $q\geq 1$. To prove $(i)$, note that $\bar{{\mathbb{E}}}(y_{i}^{2})\leq{\mathbb{E}}(y_{i}^{4})^{1/2}\leq\sqrt{C^{\prime}}$, and ${\mathbb{E}}_{n}[X_{ij}^{2}y_{i}^{2}]\leq({\mathbb{E}}_{n}[X_{ij}^{4}])^{1/2}({\mathbb{E}}_{n}[y_{i}^{4}])^{1/2}\leq 4\sqrt{C^{\prime}}C^{2}$. To prove $(ii)$, note that by Holder’s inequality we have that for any fixed $j\in[p]$, ${\mathbb{E}}_{n}[|X_{ij}^{3}w_{i}^{3}|]\leq({\mathbb{E}}_{n}[X_{ij}^{12}])^{1/4}({\mathbb{E}}_{n}[w_{i}^{4}])^{3/4}\leq 12^{3/2}C^{3}C^{\prime\prime 3/4}=O(1)$ and also by our assumption $\log p=o(n^{1/3})$. Finally to show $(iii)$, we note that by simple union bounds and tail properties of sub-gaussian variables, $\max_{ij}X_{ij}^{2}=O_{p}(\log p)$. Further, $s=o(n/\log^{2}(p))$ and hence the first part of $(iii)$ holds. To prove the second part of $(iii)$, we note that $\max_{j\in[p]}{\mathbb{E}}_{n}[X_{ij}^{4}w_{i}^{4}]\leq{\mathbb{E}}_{n}[w_{i}^{4}]\max_{i\in[n],j\in[p]}X_{ij}^{4}=O_{P}(\log^{2}p)$. Now by an application of maximal inequality (See [BCCH12, Lemma S.4]), we obtain

[TABLE]

Likewise, we have

[TABLE]

Combining these two equations we get the second part of $(iii)$.