Scalable simultaneous inference in high-dimensional linear regression models

TL;DR

This paper introduces a computationally efficient method for simultaneous inference in high-dimensional linear regression models using the Moore-Penrose pseudoinverse, enabling scalable analysis with reliable coverage.

Contribution

It proposes a novel, fast inference method based on the Moore-Penrose pseudoinverse that is free of tuning parameters and effective under certain symmetry assumptions.

Findings

Achieves near-nominal coverage in simulations

Computational complexity is significantly reduced

Regularization can improve efficiency

Abstract

The computational complexity of simultaneous inference methods in high-dimensional linear regression models quickly increases with the number variables. This paper proposes a computationally efficient method based on the Moore-Penrose pseudoinverse. Under a symmetry assumption on the available regressors, the estimators are normally distributed and accompanied by a closed-form expression for the standard errors that is free of tuning parameters. We study the numerical performance in Monte Carlo experiments that mimic the size of modern applications for which existing methods are computationally infeasible. We find close to nominal coverage, even in settings where the imposed symmetry assumption does not hold. Regularization of the pseudoinverse via a ridge adjustment is shown to yield possible efficiency gains.

Tables

Table 1. Table 1 : Computational complexity comparison

Method	Complexity	Method	Complexity
MPI	$n\cdot p^2$	Van de Geer et al., (2014)	$p^4$
RID	$n\cdot p^2$	Javanmard and Montanari, (2014)	$p^4$

Table 2. Table 2 : Simulation designs in the Monte Carlo experiments

1: Independent predictors
${\beta }_j=\{\begin{array}{ccc}{\left(-1\right)}^{u_j}\left(\left|z_j\right|+4\frac{\log \left(n\right)}{\sqrt{n}}\right)& & \\ j=1,\dots ,5& & \\ 0\text{ otherwise.}& & \end{array}$	$\begin{array}{c}{x}_i\sim N(0,\Sigma )\\ \Sigma =I_p\end{array}$
$u_j\sim \text{B}(1,0.4),z_j\sim N(0,1)$
2-4: Equicorrelated predictors
	$\begin{array}{c}{x}_i\sim N(0,\Sigma )\\ \Sigma =\left(1-\rho \right)I_p+\rho \iota {\iota }^{\prime }\end{array}$	$\rho =\{0.3,0.6,0.9\}$
5-7: Autoregressive correlation
	$\begin{array}{c}{x}_i\sim N(0,\Sigma )\\ {\Sigma }_{jk}={\rho }^{\left|j-k\right|}\end{array}$	$\rho =\{0.3,0.6,0.9\}$
8-10: Factor models
	$\begin{array}{c}{x}_{ij}={\sum }_{l=1}^k{f}_{il}{\varphi }_{lj}+{\eta }_{ij}\end{array}$	$k=\{2,10,20\}$
11-13: Group structure
		${\delta }^2=\{0.01,0.05,0.1\}$
14: Extreme correlation

Table 3. Table 3: Results Monte Carlo experiments with simulation design 3, 6, and 9

		$(p,n)=(200,100)$			$(p,n)=(1000,200)$			$(p,n)=(10000,400)$
		MAE	CR	power	MAE	CR	power	MAE	CR	power
method	$S$	Equicorrelated with $\rho =0.6$
MPI	$S$	3.66	0.94	0.19	2.11	0.94	0.43	1.39	0.94	0.79
	$S^c$	3.51	0.95		1.90	0.96		1.23	0.96
RID	$S$	2.86	0.94	0.26	2.00	0.94	0.47	1.38	0.94	0.80
	$S^c$	2.60	0.96		1.76	0.96		1.21	0.96
GBRD	$S$	2.68	0.93	0.28	1.96	0.93	0.48
	$S^c$	2.25	0.96		1.61	0.96
JM	$S$	2.99	0.31	0.39	2.80	0.25	0.51
	$S^c$	0.80	0.95		0.61	0.98
		Factor model with $k=10$
MPI	$S$	3.87	0.97	0.13	2.26	0.96	0.37	1.47	0.96	0.68
	$S^c$	3.76	0.98		2.10	0.97		1.37	0.97
RID	$S$	2.57	0.94	0.27	1.98	0.94	0.44	1.44	0.95	0.70
	$S^c$	2.06	0.98		1.66	0.97		1.28	0.97
GBRD	$S$	2.66	0.92	0.27	2.09	0.88	0.44
	$S^c$	1.99	0.98		1.45	0.97
JM	$S$	3.40	0.15	0.43	3.35	0.10	0.50
	$S^c$	0.36	0.96		0.33	0.98
		Autoregressive with $\rho =0.6$
MPI	$S$	0.64	0.97	0.65	0.34	0.95	0.98	0.21	0.91	1.00
	$S^c$	0.61	0.98		0.29	0.97		0.17	0.97
RID	$S$	0.49	0.95	0.83	0.31	0.95	0.99	0.21	0.91	1.00
	$S^c$	0.43	0.98		0.27	0.97		0.17	0.97
GBRD	$S$	0.47	0.95	0.86	0.31	0.95	0.98
	$S^c$	0.39	0.98		0.28	0.97
JM	$S$	0.48	0.67	0.96	0.33	0.67	1.00
	$S^c$	0.20	0.97		0.15	0.97

Table 4. Table 4 : Results independent and equicorrelated designs

		(p,n)=(200,100)			(p,n)=(1000,200)			(p,n)=(10000,400)
		MAE	CR	power	MAE	CR	power	MAE	CR	power
method	$S$	Independent
MPI	$S$	0.70	0.98	0.63	0.38	0.96	0.97	0.22	0.95	1.00
	$S^c$	0.64	0.99		0.30	0.99		0.18	0.98
RID	$S$	0.62	0.97	0.77	0.37	0.95	0.97	0.22	0.95	1.00
	$S^c$	0.52	0.99		0.29	0.99		0.17	0.98
GBRD	$S$	0.58	0.96	0.85	0.36	0.94	0.98
	$S^c$	0.43	0.99		0.26	0.99
JM	$S$	0.66	0.84	0.89	0.41	0.83	0.99
	$S^c$	0.33	0.98		0.21	0.99
		Equicorrelated with $\rho =0.3$
MPI	$S$	2.31	0.95	0.37	1.38	0.93	0.81	0.90	0.91	0.99
	$S^c$	2.11	0.97		1.15	0.97		0.74	0.96
RID	$S$	1.94	0.93	0.51	1.34	0.92	0.84	0.90	0.91	0.99
	$S^c$	1.65	0.97		1.08	0.97		0.73	0.96
GBRD	$S$	1.82	0.91	0.59	1.32	0.90	0.86
	$S^c$	1.37	0.97		0.98	0.97
JM	$S$	2.00	0.59	0.72	1.71	0.51	0.91
	$S^c$	0.82	0.95		0.57	0.97
		Equicorrelated with $\rho =0.6$
MPI	$S$	3.66	0.94	0.19	2.11	0.94	0.43	1.39	0.94	0.79
	$S^c$	3.51	0.95		1.90	0.96		1.23	0.96
RID	$S$	2.86	0.94	0.26	2.00	0.94	0.47	1.38	0.94	0.80
	$S^c$	2.60	0.96		1.76	0.96		1.21	0.96
GBRD	$S$	2.68	0.93	0.28	1.96	0.93	0.48
	$S^c$	2.25	0.96		1.61	0.96
JM	$S$	2.99	0.31	0.39	2.80	0.25	0.51
	$S^c$	0.80	0.95		0.61	0.98
		Equicorrelated with $\rho =0.9$
MPI	$S$	8.37	0.95	0.08	4.60	0.96	0.12	2.90	0.96	0.23
	$S^c$	8.28	0.96		4.53	0.96		2.92	0.96
RID	$S$	5.70	0.95	0.09	4.07	0.96	0.12	2.83	0.96	0.23
	$S^c$	5.56	0.96		3.99	0.96		2.85	0.96
GBRD	$S$	5.49	0.95	0.09	3.96	0.95	0.12
	$S^c$	5.28	0.96		3.85	0.96
JM	$S$	4.17	0.10	0.12	3.57	0.47	0.13
	$S^c$	0.80	0.96		1.28	0.97

Table 5. Table 5 : Results auto-regressive designs

		(p,n)=(200,100)			(p,n)=(1000,200)			(p,n)=(10000,400)
		MAE	CR	power	MAE	CR	power	MAE	CR	power
method	$S$	Auto-regressive with $\rho =0.3$
MPI	$S$	0.48	0.97	0.81	0.27	0.96	1.00	0.17	0.94	1.00
	$S^c$	0.45	0.98		0.24	0.98		0.15	0.97
RID	$S$	0.40	0.96	0.91	0.26	0.95	1.00	0.17	0.94	1.00
	$S^c$	0.36	0.98		0.23	0.98		0.15	0.97
GBRD	$S$	0.37	0.95	0.95	0.25	0.96	1.00
	$S^c$	0.30	0.98		0.21	0.98
JM	$S$	0.41	0.80	0.96	0.30	0.78	1.00
	$S^c$	0.21	0.97		0.15	0.98
		Auto-regressive with $\rho =0.6$
MPI	$S$	0.64	0.97	0.65	0.34	0.95	0.98	0.21	0.91	1.00
	$S^c$	0.61	0.98		0.29	0.97		0.17	0.97
RID	$S$	0.49	0.95	0.83	0.31	0.95	0.99	0.21	0.91	1.00
	$S^c$	0.43	0.98		0.27	0.97		0.17	0.97
GBRD	$S$	0.47	0.95	0.86	0.31	0.95	0.98
	$S^c$	0.39	0.98		0.28	0.97
JM	$S$	0.48	0.67	0.96	0.33	0.67	1.00
	$S^c$	0.20	0.97		0.15	0.97
		Auto-regressive with $\rho =0.9$
MPI	$S$	1.60	0.95	0.23	0.78	0.90	0.63	0.49	0.74	0.98
	$S^c$	1.60	0.96		0.65	0.96		0.27	0.96
RID	$S$	0.99	0.89	0.49	0.68	0.83	0.79	0.48	0.72	0.99
	$S^c$	0.78	0.96		0.47	0.96		0.25	0.96
GBRD	$S$	1.08	0.91	0.46	0.76	0.91	0.66
	$S^c$	0.89	0.97		0.64	0.97
JM	$S$	0.89	0.38	0.75	0.69	0.57	0.85
	$S^c$	0.23	0.97		0.27	0.96

Table 6. Table 6 : Results factor model designs

		(p,n)=(200,100)			(p,n)=(1000,200)			(p,n)=(10000,400)
		MAE	CR	power	MAE	CR	power	MAE	CR	power
method	$S$	Factor model with $k=2$
MPI	$S$	2.14	0.97	0.40	1.25	0.95	0.79	0.88	0.93	0.95
	$S^c$	2.00	0.98		1.09	0.98		0.67	0.98
RID	$S$	1.70	0.96	0.58	1.19	0.95	0.83	0.87	0.93	0.96
	$S^c$	1.44	0.98		0.99	0.98		0.66	0.98
GBRD	$S$	1.70	0.91	0.63	1.29	0.87	0.84
	$S^c$	1.13	0.98		0.81	0.98
JM	$S$	2.47	0.34	0.68	2.09	0.30	0.85
	$S^c$	0.45	0.96		0.37	0.98
		Factor model with $k=10$
MPI	$S$	3.87	0.97	0.13	2.26	0.96	0.37	1.47	0.96	0.68
	$S^c$	3.76	0.98		2.10	0.97		1.37	0.97
RID	$S$	2.57	0.94	0.27	1.98	0.94	0.44	1.44	0.95	0.70
	$S^c$	2.06	0.98		1.66	0.97		1.28	0.97
GBRD	$S$	2.66	0.92	0.27	2.09	0.88	0.44
	$S^c$	1.99	0.98		1.45	0.97
JM	$S$	3.40	0.15	0.43	3.35	0.10	0.50
	$S^c$	0.36	0.96		0.33	0.98
		Factor model with $k=20$
MPI	$S$	5.08	0.98	0.06	2.97	0.97	0.18	1.99	0.96	0.43
	$S^c$	4.96	0.98		2.89	0.97		1.89	0.97
RID	$S$	2.68	0.91	0.23	2.41	0.92	0.28	1.89	0.93	0.47
	$S^c$	1.85	0.98		1.84	0.97		1.59	0.97
GBRD	$S$	3.00	0.95	0.16	2.52	0.91	0.26
	$S^c$	2.43	0.98		1.88	0.98
JM	$S$	3.13	0.22	0.49	3.61	0.07	0.44
	$S^c$	0.38	0.96		0.20	0.99

Table 7. Table 7 : Results group structure and extreme designs

		(p,n)=(200,100)			(p,n)=(1000,200)			(p,n)=(10000,400)
		MAE	CR	power	MAE	CR	power	MAE	CR	power
method	$S$	Group structure with ${\delta }^2=0.01$
MPI	$S$	4.56	0.85	0.23	3.28	0.76	0.51	2.92	0.70	0.92
	$S^c$	2.85	0.98		1.50	0.98		0.98	0.97
RID	$S$	3.78	0.81	0.27	3.22	0.75	0.54	2.91	0.70	0.92
	$S^c$	2.24	0.98		1.43	0.98		0.97	0.97
GBRD	$S$	6.20	0.90	0.14	5.65	0.88	0.18
	$S^c$	1.86	0.98		1.34	0.98
JM	$S$	3.00	0.79	0.27	3.72	0.79	0.26
	$S^c$	1.44	0.98		1.07	0.98
		Group structure with ${\delta }^2=0.05$
MPI	$S$	4.73	0.90	0.18	3.14	0.73	0.52	2.80	0.62	0.92
	$S^c$	2.76	0.98		1.46	0.98		0.96	0.97
RID	$S$	3.75	0.84	0.27	3.08	0.72	0.55	2.80	0.62	0.92
	$S^c$	2.18	0.98		1.40	0.98		0.95	0.97
GBRD	$S$	5.11	0.94	0.12	4.16	0.90	0.19
	$S^c$	1.81	0.98		1.31	0.98
JM	$S$	2.75	0.77	0.33	3.15	0.79	0.33
	$S^c$	1.41	0.98		1.05	0.98
		Group structure with ${\delta }^2=0.1$
MPI	$S$	4.55	0.92	0.16	2.99	0.74	0.53	2.69	0.59	0.92
	$S^c$	2.68	0.98		1.43	0.98		0.94	0.97
RID	$S$	3.63	0.86	0.26	2.95	0.73	0.57	2.69	0.59	0.93
	$S^c$	2.13	0.98		1.37	0.98		0.93	0.97
GBRD	$S$	4.28	0.94	0.13	3.42	0.91	0.21
	$S^c$	1.77	0.98		1.28	0.98
JM	$S$	2.51	0.77	0.39	2.64	0.75	0.43
	$S^c$	1.37	0.98		1.03	0.98
		Extreme correlation
MPI	$S$	3.05	0.85	0.19	2.66	0.36	0.41	2.62	0.33	0.53
	$S^c$	3.06	0.98		1.63	0.98		1.06	0.98
RID	$S$	2.82	0.42	0.31	2.76	0.33	0.40	2.69	0.31	0.52
	$S^c$	2.09	0.97		1.46	0.98		1.04	0.98
GBRD	$S$	4.03	0.69	0.16	3.75	0.75	0.19
	$S^c$	1.91	0.98		1.36	0.98
JM	$S$	3.62	0.17	0.34	3.48	0.29	0.37
	$S^c$	0.35	0.95		0.48	0.98

Equations

$$y=X\beta +\epsilon ,\epsilon \sim N(0,{\sigma }^2{I}_n),$$

$$\begin{array}{rl}{\hat{\beta }}^b& =\hat{M}y\\ & =\hat{M}X\beta +\hat{M}\epsilon \\ & =\beta +\left(\hat{M}X-I_p\right)\beta +\hat{M}\epsilon ,\end{array}$$

$$\begin{array}{rl}\hat{\beta }& ={\hat{\beta }}^b-\left(\hat{M}X-I_p\right){\hat{\beta }}^{\text{init}}\\ & =\beta +\left(\hat{M}X-I_p\right)\left(\beta -{\hat{\beta }}^{\text{init}}\right)+\hat{M}\epsilon .\end{array}$$

$${\hat{\beta }}^{\text{lasso}}=\arg \underset{b}{\min }\left[\frac{1}{n}(y-Xb{)}^{\top }(y-Xb)+\lambda ||b|{|}_1\right].$$

$$\hat{M}=\hat{D}\tilde{M},{\hat{D}}_{jj}=(e_j^{\top }\tilde{M}X{e}_j{)}^{-1}.$$

$$\hat{M}=\hat{D}{X}^{\top }(X{X}^{\top }{)}^{-1}.$$

$${\hat{D}}_{jj}={\left[e_j^{\top }{X}^{\top }(X{X}^{\top }{)}^{-1}X{e}_j\right]}^{-1}.$$

$$V(\hat{\beta })={\sigma }^2\hat{D}{X}^{\top }(X{X}^{\top }{)}^{-2}X\hat{D}.$$

$$\hat{M}=\hat{D}(X^{\top }X+\gamma I_p{)}^{-1}{X}^{\top },$$

$${\hat{D}}_{jj}={\left(e_j^{\top }(X^{\top }X+\gamma I_p{)}^{-1}{X}^{\top }X{e}_j\right)}^{-1}.$$

$$\begin{array}{rl}{X}^{\top }(X{X}^{\top }{)}^{-1}& =\underset{\gamma \to 0}{\lim }{\left(X^{\top }X+\gamma I_p\right)}^{-1}{X}^{\top }\\ & =\underset{\gamma \to 0}{\lim }{X}^{\top }{\left(X{X}^{\top }+\gamma I_n\right)}^{-1},\end{array}$$

$$\left[{\hat{\beta }}_j-z_{\alpha /2}\sqrt{{\hat{\Omega }}_{jj}/n},{\hat{\beta }}_j+z_{\alpha /2}\sqrt{{\hat{\Omega }}_{jj}/n}\right],$$

$$({\hat{\beta }}^{\text{scaled}},{\hat{\sigma }}^{\text{scaled}})=\arg \underset{b,s}{\min }\left[\frac{1}{2ns}(y-Xb{)}^{\top }(y-Xb)+\frac{s}{2}+{\lambda }_0||b|{|}_1\right],$$

$$\sqrt{n}(\hat{\beta }-\beta )$$

$$Z|X$$

$$\sqrt{n}(\hat{\beta }(\gamma )-\beta )$$

$$Z|X$$

$${\hat{\Omega }}_{jj}^{\text{RID}}(\hat{D})-{\hat{\Omega }}_{jj}^{\text{MPI}}(\hat{D})<0.$$

$$\begin{array}{rl}\text{MAE}& =\frac{1}{s_0}\sum _{j\in S}|{\hat{\beta }}_j-{\beta }_j|,\\ \text{Coverage}& =\frac{1}{s_0}\sum _{j\in S}\text{mathbbm}1[{\beta }_j\in ({\hat{\beta }}_j-1.96\sqrt{n^{-1}{\hat{\Omega }}_{jj}},{\hat{\beta }}_j+1.96\sqrt{n^{-1}{\hat{\Omega }}_{jj}})],\\ \text{Power}& =\frac{1}{s_0}\sum _{j\in S}\text{mathbbm}1[H_{0,j}\text{ is rejected}],\\ \text{FWER}& =\frac{1}{p-s_0}\sum _{j\in S^c}\text{mathbbm}1[H_{0,j}\text{ is rejected}],\end{array}$$

$$y=X\beta +\epsilon ,\epsilon \sim N(0,{\sigma }^2{I}_n),$$

$$X=VS{U}^{\top },$$

$$X=V{S}_n{U}_n^{\top },$$

$$(1-{ϵ}_u)\frac{n}{p}\le \frac{v^{\top }{U}_n{U}_n^{\top }v}{v^{\top }v}\le (1+{ϵ}_u)\frac{n}{p}.$$

$$|e_1^{\top }{U}_n{U}_n^{\top }{e}_2|>(1+{ϵ}_f)\frac{t}{\sqrt{n}}\frac{n}{p}.$$

$$P\left(|e_1^{\top }{U}_n{U}_n^{\top }{e}_2|>(1+{ϵ}_f)\frac{1}{\sqrt{n}}\frac{n}{p}|W|\right)\le 3\exp (-C_{f}n),$$

$$P\left(|e_1^{\top }{U}_n{U}_n^{\top }{e}_2|>(1+{ϵ}_f)\frac{t}{\sqrt{n}}\frac{n}{p}\right)\le 3\exp (-C_{f}n)+2\frac{\exp (-t^2/2)}{t\sqrt{2\pi }}.$$

$$\begin{array}{rl}P({\hat{\Sigma }}_{ii}<c)& \le P\left(\frac{1}{n}{e}_i^{\top }{R}^{\top }R{e}_i-\frac{2}{n}|e_i^{\top }{R}^{\top }F\Lambda e_i|<c\right)\\ & \le \underset{(I)}{\underbrace{P\left(\frac{1}{n}{e}_i^{\top }{R}^{\top }R{e}_i<c+\tilde{c}\right)}}+\underset{(II)}{\underbrace{P\left(\frac{2}{n}|e_i^{\top }{R}^{\top }F{\Lambda }^{\top }{e}_i|>\tilde{c}\right)}}.\end{array}$$

$$\begin{array}{rl}(I)& \le P\left(\frac{1}{n}\sum _{j=1}^n{u}_j\ge c_1\right)=O(-\exp (-n)).\end{array}$$

$$\begin{array}{rl}P\left(\frac{2}{n}|\sum _{j=1}^n{e}_i^{\top }{R}^{\top }{e}_j{e}_j^{\top }F{e}_l{e}_l^{\top }{\Lambda }^{\top }{e}_i|>\tilde{c}|\Lambda \right)& =O\left(\exp \left(-\frac{n}{|e_l^{\top }{\Lambda }^{\top }{e}_i|}\right)\right).\end{array}$$

$$\begin{array}{rl}P\left(\frac{2}{n}|e_i^{\top }{R}^{\top }F{e}_l{e}_l^{\top }{\Lambda }^{\top }{e}_i|>\tilde{c}\right)& =O(\exp (-n)).\end{array}$$

$$\begin{array}{rl}P\left(\frac{2}{n}|e_i^{\top }{R}^{\top }F{\Lambda }^{\top }{e}_i|>\tilde{c}\right)& =O(\exp (-n)).\end{array}$$

Full text

Keywords: high-dimensional regression, confidence intervals, Moore-Penrose pseudoinverse, ridge regression

1 Introduction

Modern data sets require inferential procedures that are valid in high-dimensional settings where the number of variables $p$ exceeds the sample size $n$. Typical examples occur in genome wide association studies, such as Sesia et al., (2018) where $n=O(10^{3})$ and $p=O(10^{5})$, or Rietveld et al., (2013) where $n=O(10^{5})$ and $p=O(10^{6})$.

Debiasing, or desparsification, strategies allow simultaneous inference on full parameters vector in high-dimensional linear regression models.Debiasing relies on an approximate inverse of the singular high-dimensional covariance matrix of the regressors, and the use of the lasso to correct for the bias resulting from this approximation. The remaining bias, generally factorized in a part due to the approximate inverse and one due to the lasso estimation error, is then negligible compared to the variance of the estimator.

One downside of current debiasing strategies is that the computational burden increases rapidly with the number of included variables, as the approximate inverse is constructed by solving a series of $p$ optimization problems. Javanmard and Montanari, (2014) rely on direct numerical optimization for each variable to find an approximate inverse, while Van de Geer et al., (2014); Zhang and Zhang, (2014); Ning and Liu, (2017); Dezeure et al., (2017) solve $p$ lasso problems to estimate the inverse covariance matrix. Even with modern computing power, this puts a strain on computational resources when $p$ is large. As a result, existing simulation studies are typically limited to settings where $(n,p)=O(10^{2})$. This drives a wedge between the scale demanded by modern applications and the scale at which the methods can be tested.

This paper proposes a scalable debiasing method that allows for simultaneous inference in high-dimensional linear regression models. The required approximate inverse is constructed by a scaled Moore-Penrose pseudoinverse, where the scaling ensures that the approximation bias becomes negligible. Replacing the sequence of optimization problems with a single pseudoinverse lowers the computational complexity with several orders of magnitude in the number of variables. Moreover, standard errors are available in closed form and free of tuning parameters.

The confidence intervals of the Moore-Penrose pseudoinverse estimator are shown to be valid under a symmetry assumption on the variables, which for example allows variables to be independent, equicorrelated, or to follow certain factor structures. We show that a ridge adjustment of the Moore-Penrose pseudoinverse estimator can lead to a power gain.

We carry out a simulation study where we consider dimensions $n$ and $p$ for which existing methods have substantial computational costs. We find that coverage is close to nominal for settings where the variables satisfy the symmetry assumption. The proposed ridge adjustment ensures that the method is competitive to existing alternatives in terms of power. We find undercoverage on the set of nonzero coefficients for some, but not all, of the settings that violate the symmetry assumption.

The proposed estimator is related to the corrected ridge estimator of Bühlmann, (2013), who finds that under a fixed design an additional bias correction is required that leads to conservative inference. Our estimator differs in the scaling factor that renders the approximation bias due to the diagonal elements of the approximate inverse exactly equal to zero. Under a random design and a slightly stronger sparsity assumption, we show that the additional bias correction is no longer required.

The outline of the paper is as follows. Section 2 introduces the estimation approach and the proposed estimators. The theoretical properties of the Moore-Penrose pseudoinverse and ridge regression are presented in Section 3. Section 4 illustrates these results through Monte Carlo simulations. Section 5 concludes.

**Notation   ** We use the following notation throughout the paper: For any $n\times 1$ vector $a=(a_{1},\dots,a_{n})^{\top}$, the $l_{q}$-norm is defined as $||a||_{q}:=(\sum^{n}_{i=1}|a_{i}|^{q})^{1/q}$ for $q>0$ and $||a||_{0}$ denotes the number of nonzero elements of $a$. The maximum norm is written as $||a||_{\infty}=\max(|a_{1}|,\dots,|a_{n}|)$. For a $p\times n$ matrix $A$, the $l_{q}$-norm is defined as $||A||_{q}:=\sup_{x,||x||_{q}=1}\left\{||Ax||_{q}\right\}$ and the maximum norm is written as $||A||_{\max}=\max_{i=1,\dots,n,j=1,\dots,p}|A_{ij}|$. When $A$ is $n\times n$, we write $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ for the minimum and maximum eigenvalues of $A$. The $n\times n$ identity matrix is denoted by $I_{n}$. The vector $e_{i}$ has its $i$-th entry equal to 1 and zeros everywhere else. For the regressor matrix $X$, we index the rows with the subscript $i=1,\ldots,n$ and the columns with the subscript $j=1,\ldots,p$. If $U$ is a $p\times p$ orthogonal matrix, we write $U\in\mathcal{O}(p)$. When two random variables $X$ and $Y$ follow the same distribution, this is denoted as $X\overset{(d)}{=}Y$.

2 High-dimensional linear regression

Consider the high-dimensional linear model

[TABLE]

where $y$ is an $n\times 1$ response vector, $X$ an $n\times p$ regressor matrix, $\beta=(\beta_{1},\dots,\beta_{p})^{\top}$ a $p\times 1$ vector of unknown regressor coefficients, and $\varepsilon$ an $n\times 1$ vector of errors which are $N(0,\sigma^{2}I_{n})$ and independent of $X$.

We develop a computationally efficient debiased estimator $\hat{\beta}$ for $\beta$ when $p>n$, with a closed-form expression for the covariance matrix of the estimator. Standard errors and confidence intervals for $\hat{\beta}$ are free of tuning parameters, closed-form, and scalable to a very large number of variables. Below we develop the estimator. Section 2.3 provides an algorithmic overview of the proposed inference method.

2.1 A debiased estimator

We consider estimators for $\beta$ of the form

[TABLE]

where $\hat{M}$ is a $p\times n$ matrix. The second term of (2) represents a bias which depends on the choice of $\hat{M}$. When $p\leq n$ and $\text{E}[\varepsilon|X]=0$, ordinary least squares yields unbiased estimators by choosing $\hat{M}=(X^{\top}X)^{-1}X^{\top}$. When $p>n$, the matrix $X^{\top}X$ is singular, and any choice of $\hat{M}$ will induce bias.

Given an initial estimator $\hat{\beta}^{\text{init}}$, we can reduce the bias in (2) by applying the correction

[TABLE]

For the initial estimator $\hat{\beta}^{\text{init}}$ we use the lasso estimator of Tibshirani, (1996) defined by

[TABLE]

This leads to an interpretation of $\hat{\beta}$ as a ‘desparsified’ version of the lasso estimator (Van de Geer et al.,, 2014), since we can write $\hat{\beta}=\hat{\beta}^{\text{lasso}}+\hat{M}(y-X\hat{\beta}^{\text{lasso}})$.

2.2 The approximate inverse $\hat{M}$

The goal of this paper is to introduce choices of $\hat{M}$ that are computationally efficient and for which the bias of the corrected estimator is small. For the bias to be small, $\hat{M}X$ has to be close to the $p\times p$ identity matrix $I_{p}$. Hence, we refer to $\hat{M}$ as an approximate inverse for $X$.

We ensure that the diagonal terms of $\hat{M}X-I_{p}$ are identically equal to zero by introducing a $p\times p$ diagonal matrix $\hat{D}$, with diagonal elements $\hat{d}_{j}$, and taking

[TABLE]

The effectiveness of $\hat{M}$ as an approximate inverse is then determined by the magnitude of the off-diagonal elements of $\hat{M}X$.

2.2.1 Moore-Penrose pseudoinverse

As a tuning parameter free choice for $\tilde{M}$ in (5) we consider the Moore-Penrose pseudoinverse (MPI) given by $X^{\top}(XX^{\top})^{-1}$. That is,

[TABLE]

The diagonal elements $\hat{D}_{jj}$ of the diagonal scaling matrix $\hat{D}$ equal

[TABLE]

This provides a closed-form expression for the approximate inverse. Since the bias term of the estimator is of lower order compared to the variance, the covariance of $\hat{\beta}$ is available in closed form as well,

[TABLE]

2.2.2 Ridge regularization

We also consider the use of a ridge adjustment (RID), which in finite samples can lead to a power gain. Define

[TABLE]

where $\gamma$ denotes the ridge penalty and the elements of the diagonal scaling matrix $D$ equal

[TABLE]

The ridge-adjustment is the natural regularization procedure as we can define the Moore-Penrose pseudoinverse as

[TABLE]

see e.g. Albert, (1972).

2.2.3 Alternative specifications

Zhang and Zhang, (2014) develop the estimator $\hat{\beta}$ in (3) with $\hat{M}=\frac{1}{n}\bar{M}X^{\top}$, where $\bar{M}$ is an approximate inverse to the empirical covariance matrix $\frac{1}{n}X^{\top}X$. This inverse is found by a series of lasso regressions. Van de Geer et al., (2014) generalize this method to deal with nonlinear models and prove semiparametric optimality. Javanmard and Montanari, (2014) use direct numerical optimization to obtain $\bar{M}$.

The ridge adjusted estimator, but with a different scaling factor, is studied by Bühlmann, (2013) for a fixed design and a less restrictive sparsity assumption then considered in Section 3. In this case an additional bias correction is necessary that results in conservative inference.

2.3 Inference

With a specification for the approximate inverse $\hat{M}$, confidence intervals can be constructed for the $j$-th element of $\beta$ in (1) as

[TABLE]

where $\hat{\beta}_{j}$ is the $j$-th element of $\hat{\beta}$ in (3), $z_{\alpha/2}$ is the $\alpha/2$ critical value for the standard normal distribution, and $\sqrt{\hat{\Omega}_{jj}/n}$ is the standard error of $\hat{\beta}_{j}$.

The covariance matrix of the estimator $\hat{\beta}$ is equal to $\hat{\Omega}=\hat{\sigma}^{2}n\hat{M}\hat{M}^{\top}$. The noise level $\sigma^{2}$ is estimated with the scaled lasso estimator of Sun and Zhang, (2012), defined as

[TABLE]

with $\lambda_{0}=\sqrt{2\log(p)/n}$. This scaled lasso estimator is commonly used for noise estimation in high-dimensional regression settings.

Algorithm 1 outlines the proposed inference method. The algorithm requires a value for the ridge penalty. The numerical results in Section 4 indicate that the ridge adjustment only has an effect when $n$ and $p$ are of the same order of magnitude. In this case, following our theoretical results, setting $\gamma=p\sqrt{\log p/n}$ improves the accuracy compared to the MPI estimator, while keeping coverage at the nominal level.

Algorithm 1 uses the lasso as the initial estimator and the scaled lasso for noise estimation. However, alternative initial estimators, such as proposed by Caner and Kock, (2014), can also be used as long as they satisfy a sufficiently tight accuracy bound on $||\beta-\hat{\beta}^{\text{init}}||_{1}$. Similarly, any consistent estimator for $\sigma$ can be used for noise estimation. Instead of the scaled lasso, Reid et al., (2016) estimate the error variance from the lasso residuals with a degree of freedom correction. Yu and Bien, (2019) take the minimizing value of the lasso regression problem.

2.4 Computational complexity

Table 1 shows that the computational complexity of the proposed methods is several orders of magnitude in the number of variables lower compared to existing methods.

For the Moore-Penrose pseudoinverse, the leading order term in the computation complexity is the matrix multiplication $X(XX^{\top})^{-1}X^{\top}$, which costs $O(n\cdot p^{2})$. The same holds for the ridge regularized version.

The method of Van de Geer et al., (2014) solves a penalized regression problem with $p-1$ explanatory variables for every column of the empirical covariance matrix. A fast solver is provided by the lars algorithm of Efron et al., (2004), which has a complexity of $O(p^{3})$ for each column. Therefore, the complexity for obtaining the approximate inverse covariance matrix equals $O(p^{4})$ in total. Javanmard and Montanari, (2014) state that this is equivalent to the complexity of their method.

All methods use an initial estimator, which is provided by the lasso. With $r=\min\left\{n,p\right\}$, this costs $O(r^{3})$ at a fixed value of the penalization parameter. For the initial estimator we apply $K$-fold cross-validation, and $r=n$, such that the computational complexity of all methods is at least $O(K\cdot n^{3})$. However, for none of the methods this is the leading order term.

3 Theoretical results

3.1 Assumptions

We make the following assumptions on the expansion rates of $n$ and $p$, the regressor matrix $X$ and the errors $\varepsilon$.

Assumption 1**.**

The dimensions $n$ and $p$ are such that $C_{1}\log p<n<C_{2}p$ for some sufficiently large constant $C_{1}>0$ and the constant $C_{2}<1$. The sparsity $s_{0}\equiv||\beta||_{0}$ satisfies $s_{0}=o\left(\frac{\sqrt{n}}{\log p}\right)$.

Assumption 2a**.**

Let $R$ be an $n\times p$ matrix with independent subgaussian rows $r_{i}=(\Sigma^{R})^{1/2}\tilde{r}_{i}$, where $\tilde{r}_{i}$ has independent subgaussian elements. The matrix $\Sigma^{R}=\emph{E}[RR^{\top}]$ satisfies $\max_{i=1,\ldots,p}\Sigma^{R}_{ii}=O(1)$, and $1/\lambda_{\min}(\Sigma^{R})=O(1)$. Assume $X=R$.

Assumption 2b**.**

Let $F$ be an $n\times k$ matrix with independent rows and identically subgaussian distributed columns. Let $\Lambda$ be an $p\times k$ matrix with identically subgaussian distributed columns, and $\Sigma^{\Lambda}=\emph{E}[\Lambda\Lambda^{\top}]$ has finite eigenvalues bounded away from zero by a constant. Assume the factor structure $X=F\Lambda^{\top}+R$, with $k$ fixed, the matrices $R$, $F$ and $\Lambda$ independent, and $R$ as defined in 2a with the additional assumption that $\lambda_{\max}(\Sigma^{R})<M$.

Assumption 3**.**

The rows of $X$ are spherically symmetric: The rows of $R$ satisfy $r_{i}\overset{(d)}{=}Tr_{i}$, and in case of a factor structure the columns of $\Lambda$ also satisfy $\lambda_{i}\overset{(d)}{=}T\lambda_{i}$, for any fixed $T\in\mathcal{O}(p)$.

Assumption 4**.**

The errors $\varepsilon\sim N(0,\sigma^{2}I_{n})$ and $\varepsilon$ is independent of $X$.

1 imposes a condition on the sample size relative to the number of variables. We operate under the assumption that $p>n$, but we assume that $n$ is sufficiently larger than $\log p$. The sparsity constraint restricts the number of non-zero coefficients in $\beta$ by $s_{0}=||\beta||_{0}$. The assumption is standard for inference in high-dimensional regression models (Van de Geer et al.,, 2014; Javanmard and Montanari,, 2014), but slightly stronger than the $s_{0}^{2}=o\left(n/\log p\right)$ required for lasso consistency Van de Geer, (2008). This discrepancy has been studied in more detail in Javanmard and Montanari, (2018). Since our results only depend on the $l_{1}$ norm of the lasso estimation error, the sparsity assumption can be relaxed to allow for approximate sparsity (Chernozhukov et al.,, 2015).

2a is taken from Theorem 2.4 in Javanmard and Montanari, (2014). This assumption is necessary to bound the bias of the initial estimator. It allows for a random $X$ with subgaussian independent rows, but the subgaussian assumption rules out a factor structure for $X$, which is allowed by 2b.

3 is necessary to bound the bias of the approximate inverse. It is for example satisfied when the columns of $X$ have an equicorrelated covariance matrix, but also covers classes of factor models. Fan and Lv, (2008) use this assumption to derive accuracy bounds on the Moore-Penrose pseudoinverse in high-dimensional variable screening. These results depend on the behavior of the right singular vectors under spherical symmetry. In a slightly different context, Shah et al., (2020) show that under spherical symmetry, the right singular vectors can be used to estimate high-dimensional covariance matrices that contain a factor component. We find in Section 4 that our methods continue to work well in a range of settings in which the spherical symmetry assumption is violated, for example when the columns of $X$ have a covariance matrix with a Toeplitz structure.

Finally, 4 ensures that the estimated parameters are normally distributed. This assumption can be relaxed as shown by Van de Geer et al., (2014), in which case valid inference can be achieved for a fixed number of coefficients.

3.2 Theorems

The following theorem presents the main result of this paper.

Theorem 1**.**

Suppose 1 – 4 hold. Let $\hat{\beta}=\hat{M}y-\left(\hat{M}X-I_{p}\right)\hat{\beta}^{\text{lasso}}$, with $\hat{M}$ as defined in (6) and $\hat{D}$ as in (7). Then,

[TABLE]

where $\hat{\Omega}=\hat{\sigma}^{2}n\hat{M}\hat{M}^{\top}$ and $\hat{\Omega}_{jj}=O_{p}(1)$.

The proof is given in Section B.1. Theorem 1 shows that the estimator $\hat{\beta}$ in (3) has a lower order bias compared to the variance. Standard errors decrease at the usual $n^{-1/2}$ rate. Theorem 1 allows for the construction of confidence intervals that are uniformly valid when $p$ is finite. Uniformity is guaranteed since the bound on the lasso estimator given in Lemma B.2 holds uniformly over all sets $S_{0}$ of size $s_{0}=o(\sqrt{n}/\log p)$, see Van de Geer et al., (2014) for a discussion.

Since the resulting covariance matrix of the estimator is available in closed form, efficient multiple testing procedures as in Bühlmann, (2013) can be employed, together with joint tests on estimated coefficients.

Theorem 2 shows that for an appropriate choice of the regularization parameter, the ridge adjusted estimator can also be used for inference.

Theorem 2**.**

Suppose 1 – 4 hold. Let $\hat{\beta}(\gamma)=\hat{M}y-\left(\hat{M}X-I_{p}\right)\hat{\beta}^{\text{lasso}}$, with $\hat{M}$ as in (9) and $\hat{D}$ as in (10). Suppose $\gamma=O(p\sqrt{\log p/n})$. Then,

[TABLE]

where $\hat{\Omega}=\hat{\sigma}^{2}n\hat{M}\hat{M}^{\top}$ and $\hat{\Omega}_{jj}=O_{p}(1)$.

The proof is given in Section B.2.

The reason one would opt for the regularized variants despite the additional tuning parameter $\gamma$ is provided by the following theorem.

Theorem 3**.**

Suppose 1 – 4 hold. Denote the variance of the estimator $\hat{\beta}_{j}$ under a diagonal scaling matrix $\hat{D}$ by $\hat{\Omega}_{jj}^{\text{MPI}}(\hat{D})$ when $\hat{M}$ is as in (6) and $\hat{\Omega}_{jj}^{\text{RID}}(\hat{D})$ when $\hat{M}$ is as in (9). For the choice of $\gamma$ as in Theorem 2, we have

[TABLE]

The proof is given in Section B.3. The proof shows that the difference between variances of the RID and MPI estimators increases in $\gamma$. Since the proof of Theorem 2 in Section B.2 suggests that $\gamma$ can attain a higher value when $p$ is close to $n$, we expect a larger potential efficiency gain of ridge regularization in these settings.

Note that Theorem 3 requires the regularized estimator and the estimator based on the Moore-Penrose pseudoinverse to use the same diagonal scaling matrix. Using $\hat{D}$ from (7) for the Moore-Penrose inverse and $\hat{D}$ from (10) for the ridge regularized inverse, does not yield an ordering in terms of power. However, in all cases we have encountered, the inequality in Theorem 3 is satisfied when using the diagonal matrix specific to the estimator under consideration. This is also evident from the Monte Carlo results in Section 4.

4 Monte Carlo Experiments

This section examines the finite sample behaviour of the proposed estimators in a Monte Carlo experiment. The performance of the estimators is compared with Van de Geer et al., (2014) (GBRD) and Javanmard and Montanari, (2014) (JM), two existing methods for constructing confidence intervals in high-dimensional regression for all coefficients. All estimators use the lasso as initial estimator with a penalty term that minimizes the mean squared error under tenfold cross-validation. The series of lasso regressions in GBRD also use tenfold cross-validation, and we set the tuning parameter $\mu=2\sqrt{n^{-1}\log p}$ in JM, which is equal to the value used in their simulation studies.

The estimators are evaluated on mean absolute bias (MAE), coverage rate, power, and family-wise error rate (FWER). For each Monte Carlo replication, we calculate

[TABLE]

where MAE and coverage are similarly defined on the set of zero coefficients $S^{c}$, and $H_{0,j}\colon\beta_{j}=0$ for $j=1,\ldots,p$. We report the average MAE, coverage and power over the Monte Carlo replications.

4.1 Monte Carlo design

The data generating process takes the form

[TABLE]

where $X$ is an $n\times p$ regressor matrix, and $\beta$ a $p\times 1$ vector of coefficients. Table 2 shows the 14 simulation designs with different specifications for $\beta$ and $X$, as used in Wang and Leng, (2015).

For each simulation design, 1000 data sets are generated with the number of predictors and sample size $(p,n)$ equal to $(200,100)$, $(1000,200)$, and $(10000,400)$. The error variance $\sigma^{2}$ is set to the value that satisfies $R^{2}=\text{var}(X\beta)/(\text{var}(X\beta)+\sigma^{2})=50\%$. Note that due to the computational costs of GBRD and JAM, we only include MPI and RID in the Monte Carlo experiments with $(p,n)=(10000,400)$.

4.2 Results

The results of the Monte Carlo experiments confirm the theory from Section 3. Figure 1 shows the coverage rates of MPI and GBRD for non-zero ($S$) and zero ($S^{c}$) coefficients. The coverage of MPI on $S$ and $S^{c}$ is close to nominal for all independent, equicorrelated, and factor model designs. These simulation designs satisfy the spherical symmetry assumption on the regressor matrix. For the factor model, GBRD attains a slightly lower coverage on the nonzero coefficients, and this effect increases when moving from $p=200$ to $p=1000$.

Some, but not all, of the simulation designs that violate the symmetry assumption report undercoverage for both MPI and GBRD. Figure 1 shows undercoverage on $S$ for autoregressive regressors with correlation parameter $\rho=0.9$, which increases with the dimensionality of the data. For both the group structure and the extreme dependence simulation designs, the undercoverage on $S$ is quite severe.

Table 3 reports the average MAE, coverage and power for experiments with equicorrelated predictors with $\rho=0.6$, autoregressive correlation with $\rho=0.6$, and a factor model with $k=10$. The results for the remaining designs are presented in Appendix D. All results corroborate the undercoverage of JM on $S$ found by Dezeure et al., (2015).

Since the equicorrelated and factor design are spherically symmetric, the coverage rates of MPI and RID are robust to the choice of $p$ and $n$. For the autoregressive model structure, we find a drop in coverage on $S$ when moving to $(p,n)=(10000,400)$.

For all simulation designs with $p=200$ and $n=100$, the noise in the pseudoinverse translates in high MAE for MPI. This can be largely remedied by incorporating the ridge adjustment. In terms of inference, this has the benefit that the coverage rate is not distorted, while power increases. For the two higher-dimensional settings, the findings for MPI and RID are quite similar.

When the number of variables is of the same order as the number of observations, the discussion below Theorem 3 suggests that regularization may improve the efficiency of the confidence intervals. Figure 2 shows histograms of the FWER across experiments for MPI and RID under an equicorrelated design with $\rho=0.6$. In the case with $(p,n)=(200,100)$, the RID rejection rates are indeed more tightly concentrated around the nominal rate of 5% than the MPI rejection rates. For $(p,n)=(1000,200)$ and $(10000,400)$, the number of variables is an order of magnitude larger than the number of observations and the benefit of regularization disappears.

Figure 3 shows the logarithm of the computation times of MPI, GBRD, and JM for each simulation design. MPI demands computation time that is negligible compared to the time needed by GBRD and JM with $p=200$ and $p=1000$. Due to the computational costs of GBRD and JM, they are not included in the experiments with $(p,n)=(10000,400)$. We find that MPI is still computationally feasible in this setting, making them scalable to ultra high-dimensional data.

5 Conclusion

This paper proposes computationally efficient methods for constructing confidence intervals in high-dimensional linear regression models. We employ a debiasing strategy where the Moore-Penrose pseudoinverse is used as an approximate inverse of the singular empirical covariance matrix of the regressors. This strategy is shown to be valid under a symmetry assumption on the regressors. The covariance matrix of the estimates is available in closed form and free of tuning parameters. Confidence intervals can then be constructed using standard procedures. A ridge-based regularization can yield power improvements in finite samples.

We consider large scale Monte Carlo experiments that show that the proposed estimators provide valid confidence intervals with correct coverage rates, even when the symmetry assumption on the regressors does not hold.

Appendix A Preliminary results

Throughout, $C$ and $C_{i}$ denote generic positive constants that can differ between lines. We abbreviate $1-C_{1}\exp(-C_{2}n)$ as $1-O(\exp(-n))$.

A.1 Properties of spherical distributions

Decompose $X$ by a singular value decomposition as

[TABLE]

where $V\in\mathcal{O}(n)$, $S$ the $n\times p$ matrix of singular values, and $U\in\mathcal{O}(p)$. Under 3, $X$ is invariant under right multiplication with an orthogonal matrix, and therefore $U$ is uniformly distributed on $\mathcal{O}(p)$. When $n<p$,

[TABLE]

where $S_{n}$ is an $n\times n$ matrix with the non-zero singular values on its diagonal, and $U_{n}$ is a $p\times n$ matrix that satisfies $U_{n}^{\top}=[I_{n},O_{n,p-n}]U^{\top}$. Since $U$ is uniformly distributed over $\mathcal{O}(p)$, $U_{n}$ is uniformly distributed over the Stiefel manifold $V_{n,p}$, defined as $V_{n,p}=\left\{A\in R^{p\times n}:A^{\top}A=I_{n}\right\}$. See for example Chikuse, (1990).

We will use the following results.

Lemma A.1** (Lemma 5.3.2 of Vershynin, (2018)).**

Let $v$ be a fixed $p\times 1$ vector, and $U_{n}$ a $p\times n$ matrix that is distributed uniformly over $V_{n,p}$. Let $(C_{u},\epsilon_{u})>0$, then with probability at least $1-2\exp(-C_{u}\epsilon_{u}^{2}n)$, we have that

[TABLE]

Lemma A.2** (Fan and Lv, (2008)).**

For any $(C_{f},t)>0$, there exists an $\epsilon_{f}>0$ such that with probability at most $3\exp(-C_{f}n)+2\frac{\exp(-t^{2}/2)}{t\sqrt{2\pi}}$, we have that

[TABLE]

Proof: Lemma 5 of Fan and Lv, (2008) shows that for any $C_{f}>0$, there exists an $\epsilon_{f}>0$, such that

[TABLE]

where $W\sim N(0,1)$ independent of $U_{n}$. For random variables $X$ and $Y$ and a constant $t$, we have $P(X>t)\leq P(X>Y,Y>t)\leq P(X>Y)+P(Y>t)$. Set $t=(1+\epsilon_{f})\frac{1}{\sqrt{n}}\frac{n}{p}$, $X=|e_{1}^{\top}U_{n}U_{n}^{\top}e_{2}|/t$ and $Y=|W|$. Using (A.5) and the standard tail bound bound $P(|W|>t)\leq 2\frac{\exp(-t^{2}/2)}{t\sqrt{2\pi}}$, we have

[TABLE]

Lemma A.3**.**

Let $\hat{\Sigma}_{ii}=e_{i}^{\top}\frac{1}{n}X^{\top}Xe_{i}$. Under 2b there exists constants $0<c<C$, such that $c<\hat{\Sigma}_{ii}<C$ with probability $1-O(\exp(-n))$.

Proof: Notice that

[TABLE]

We first bound (I). By 2a, $Re_{i}$ is a subgaussian vector with independent elements. Under this assumption, $u_{j}=e_{i}^{\top}R^{\top}e_{j}e_{j}^{\top}Re_{i}-\Sigma_{ii}^{R}$ is subexponential with mean zero. From 2a, we also have that $\Sigma_{ii}^{R}=e_{i}^{\top}\Sigma^{R}e_{i}\geq\lambda_{\min}(\Sigma^{R})>0$. Pick $c,\tilde{c}>0$, such that $\lambda_{\min}(\Sigma^{R})-(c+\tilde{c})=c_{1}>0$. Then,

[TABLE]

Continuing with (II), we first condition on $\Lambda$. From 2b, $e_{j}F\Lambda^{\top}e_{i}$ is a subgaussian random variable, and by independence between $R$ and $(F,\Lambda)$, $v_{jl}=e_{i}^{\top}R^{\top}e_{j}e_{j}^{\top}Fe_{l}e_{l}^{\top}\Lambda^{\top}e_{i}$ where $j=1,\ldots,n$ and $l=1,\ldots,k$, is subexponential with mean zero. Then,

[TABLE]

Taking the expectation over $\Lambda$, using Jensen’s inequality and the fact that $\text{E}[|e_{l}^{\top}\Lambda^{\top}e_{i}|]\leq C<\infty$,

[TABLE]

Since $l=1,\ldots,k$ with $k$ finite, we then also have that

[TABLE]

What remains to be shown is that there exists a $C$ such that

[TABLE]

For this it is sufficient to show that

[TABLE]

where $\Lambda$ has rows $\Lambda_{i}$. Starting with the first, notice that

[TABLE]

Conditioning on $\Lambda_{il}$ and using that the elements of $f_{l}$ are independent subgaussian,

[TABLE]

Taking the expectation over $\Lambda_{il}$, using Jensen’s inequality and the fact that $\text{E}[\Lambda_{il}^{2}]\leq C<\infty$, we then find that for fixed $k$, there exists a constant $C_{1}>0$ such that

[TABLE]

By the same arguments as in Section 6.2.2 of Javanmard and Montanari, (2014), we also have that there exists a constant $C_{2}>0$, such that

[TABLE]

This completes the proof.

Lemma A.4**.**

Suppose 1 and 2b hold. Denote by $\hat{\lambda}_{\min\neq 0}$ the smallest nonzero eigenvalue of $p^{-1}X^{\top}X$. Then for some constant $C>0$, $\hat{\lambda}_{\min\neq 0}\geq C$ with probability $1-O(\exp(-n))$.

To determine the minimum nonzero eigenvalue, define $v$ such that $v^{\top}v=1$ and write

[TABLE]

The second to last line holds uniformly over $v$ with probability $1-O(\exp(-n))$ by (A.11) replacing $e_{i}$ by $v$ and noting in applying Jensen’s inequality that $\Lambda$ has subgaussian columns by 2b. This shows that $\lambda_{\min\neq 0}(p^{-1}X^{\top}X)\geq\lambda_{\min\neq 0}(p^{-1}R^{\top}R)-C_{1}(n/p)$. We will now show that the minimum nonzero eigenvalue of $p^{-1}R^{\top}R$ is bounded away from zero by a constant. This implies that for $p$ sufficiently larger than $n$, Lemma A.4 holds.

Since $p>n$, the minimum nonzero eigenvalue of $p^{-1}R^{\top}R$ is equal to

[TABLE]

with $\tilde{r}_{i}=e_{i}^{\top}\tilde{R}$ defined in 2a and the last line follows from 2a. Since $\tilde{R}$ is an $n\times p$ matrix with independent subgaussian elements, by Vershynin, (2010) Theorem 5.39, $\lambda_{\min}(p^{-1}\tilde{R}\tilde{R}^{\top})\geq 1-C_{2}\sqrt{n/p}$ with probability $1-O(\exp(-n))$.

We now have shown that $\lambda_{\min\neq 0}(p^{-1}X^{\top}X)\geq 1-C_{2}\sqrt{n/p}-C_{1}(n/p)$. From 1, we can assume that $p>\max(nC_{2}^{2}\epsilon^{-2},nC_{1}\epsilon^{-1})$ for some positive constant $\epsilon<1$. Then with probability $1-O(\exp(-n))$, $\hat{\lambda}_{\min\neq 0}(p^{-1}X^{\top}X)\geq C>0$.

Appendix B Proofs of the main theory

B.1 Proof of Theorem 1

Rewrite the estimator $\hat{\beta}$ in (3) as

[TABLE]

The bias term of the estimator in (B.1) is bounded by the norm inequality,

[TABLE]

We then have that

[TABLE]

Lemma B.1 provides a probability bound for $\left|\left|\hat{M}X-I_{p}\right|\right|_{\max}$, and Lemma B.2 provides a probability bound for $||\beta-\hat{\beta}^{\text{lasso}}||_{1}$.

Lemma B.1**.**

Suppose 1 and 3 holds. Define $\hat{M}$ as in (6) and $\hat{D}$ as in (7). Then, there exists an $c>0$ such that

[TABLE]

Proof: Section C.1.

Lemma B.2**.**

Suppose 1, 4 and either 2a or 2b hold. Consider the lasso estimator (4) with $\lambda\geq 4C_{1}\sigma\sqrt{\log p/n}$ for some sufficiently large constant $C_{1}$,

[TABLE]

Proof: Section C.2.

Combining 1, Lemma B.1, and Lemma B.2, it follows that the bias can be bounded by

[TABLE]

From 4 and (B.1) it follows that $Z|X\sim N(0,\hat{\Omega})$ with $\hat{\Omega}=\hat{\sigma}^{2}n\hat{M}\hat{M}^{\top}$. The following lemma shows that $\hat{\Omega}_{jj}=O_{p}(1)$, which completes the proof of Theorem 1.

Lemma B.3**.**

Suppose 1 – 4 hold. Define $Z_{j}=\sqrt{n}\hat{D}_{jj}x_{j}^{\top}(XX^{\top})^{-1}\varepsilon$ and $\hat{\Omega}_{jj}=n\hat{\sigma}^{2}e_{j}^{\top}\hat{M}\hat{M}^{\top}e_{j}$ with $\hat{M}$ as before. For $j=1,\ldots,p$,

[TABLE]

Proof: Section C.3.

B.2 Proof of Theorem 2

We show that for a sufficiently small penalty parameter $\lambda$, the results under a Moore-Penrose inverse carry over to a ridge adjusted estimator. We adjust Lemma B.1 and Lemma B.3 to apply to the ridge adjusted estimator. Theorem 2 then follows.

Lemma B.4**.**

Suppose 1 and 3 holds. Define $\hat{M}$ and $\hat{D}$ as in (9) and (10). Choose $\gamma=O\left(p\sqrt{\frac{\log p}{n}}\right)$. Then, there exists an $a>0$ such that

[TABLE]

Proof: Section C.4.

Lemma B.5**.**

Suppose 1 – 4 hold. Define $Z_{j}=\sqrt{n}\hat{D}_{jj}e_{j}^{\top}(X^{\top}X+\gamma I_{p})^{-1}X^{\top}\varepsilon$ and $\hat{\Omega}_{jj}=n\hat{\sigma}^{2}e_{j}^{\top}\hat{M}\hat{M}^{\top}e_{j}$ with $\hat{M}$ as before. For $j=1,\ldots,p$,

[TABLE]

Proof: Section C.5.

B.3 Proof of Theorem 3

The relative efficiency is determined by the difference

[TABLE]

with $\hat{\Omega}_{ii}^{\text{MPI}}$ as in Theorem 1 and $\hat{\Omega}_{ii}^{\text{RID}}$ as in Theorem 2, but both estimators are assumed to either use $\hat{D}$ as in (7) or $\hat{D}$ as in (10). Define the $p\times p$ matrix ${Q}$ as a diagonal matrix with containing the eigenvalues of $X^{T}X$ and by ${Q}_{+}^{-1}$ the Moore-Penrose inverse of ${Q}$. For the ridge regularized inverse we have

[TABLE]

where $A_{\text{RID}}^{2}$ is a diagonal matrix with the diagonal elements satisfying $0<A_{ii}^{2}<1$.

For the Moore-Penrose pseudoinverse,

[TABLE]

where $A_{\text{MPI}}^{2}=I$ is a diagonal matrix with the diagonal elements $A_{ii}^{2}=1$.

Appendix C Proofs of the lemmas

C.1 Proof of Lemma B.1

As the diagonal elements of $\hat{M}X-I_{p}$ are identically zero, we need to establish a probability bound on the off-diagonal elements of $\hat{M}X-I_{p}$. These elements are of the form $e_{i}^{\top}U_{n}U_{n}^{\top}e_{j}/e_{i}^{\top}U_{n}U_{n}^{\top}e_{i}$. Fix $i\neq j$. Then,

[TABLE]

where we use Lemma A.1 and Lemma A.2.

Taking the union bound over all $i\neq j$, we get

[TABLE]

Now choose $t=\frac{1-\epsilon_{u}}{1+\epsilon_{f}}c\sqrt{\log p}$. Suppose $c>3\frac{1+\epsilon_{f}}{1-\epsilon_{u}}$. According to 1, $n>3\max(C_{f}^{-1},C_{u}^{-1}\epsilon_{u}^{-2})\log p$. Then,

[TABLE]

C.2 Proof of Lemma B.2

For the bound in Lemma B.2 we first show that $\max_{i=1,\ldots,p}|e_{i}^{\top}X^{\top}\varepsilon|/\sqrt{n\sigma^{2}}\leq C_{1}\sqrt{\log p}$ with high probability. Denoting by $S_{0}$ the index set of nonzero coefficients, this means that with high probability the lasso estimator satisfies $||\hat{\beta}_{S_{0}^{c}}||_{1}\leq 3||\hat{\beta}_{S_{0}}-\beta_{S_{0}}||_{1}$ by Lemma 6.3 of Bühlmann and Van De Geer, (2011).

We start by decomposing

[TABLE]

First consider (II). Denote by $Z_{i}$ a standard normal random variable and $\hat{\Sigma}_{ii}=e_{i}^{\top}\frac{1}{n}R^{\top}Re_{i}$. By 4, $\varepsilon\sim N(0,\sigma^{2}I_{n})$, and (II) is distributed as

[TABLE]

This term can be bounded in probability as

[TABLE]

where we first use Lemma A.3, and subsequently the tail bound for a standard normal random variable.

Choosing $a=2\sqrt{C}$, the first term on the left hand-side is $O(p^{-1})$, while the second term is $O(\exp(-n))$ when $n$ is sufficiently larger than $\log p$ (1). We conclude that,

[TABLE]

For (I), notice that

[TABLE]

Define $\hat{\Sigma}^{F}_{jj}=e_{j}^{\top}\frac{1}{n}F^{\top}Fe_{j}$ and note again that $|f_{j}^{\top}\varepsilon|/\sqrt{n\sigma^{2}}$ is distributed as $\sqrt{\hat{\Sigma}^{F}_{jj}}|Z_{j}|$. Conditional on $\Lambda$, we can use the same arguments as before to show that

[TABLE]

Taking now the expectation of $\Lambda$ and using Jensen’s inequality, we obtain

[TABLE]

We conclude that for sufficiently large $a$ and $n$, we have that

[TABLE]

From these results, we conclude that there exists an $C_{1}$, such that with high probability

[TABLE]

Choosing then the lasso penalty parameter as $\lambda=2C_{1}\sigma\sqrt{\frac{\log p}{n}}$ for some sufficiently large constant $C_{1}$, the results in Bühlmann and Van De Geer, (2011) show that with high probability $||\hat{\beta}_{S_{0}^{c}}||_{1}\leq 3||\hat{\beta}_{S_{0}}-\beta_{S_{0}}||_{1}$.

The standard $\ell_{1}$ norm accuracy bounds for the lasso follow if the following restricted eigenvalue condition holds (Bickel et al.,, 2009).

Definition 1**.**

The restricted eigenvalue condition holds if

[TABLE]

for all $\beta$ for which $||\beta_{S_{0}^{c}}||_{1}\leq 3||\beta_{S_{0}}||_{1}$ and $\phi_{0}>0$.

We start with the following lower bound.

[TABLE]

where $\tilde{R}$ has rows $\tilde{r}_{i}$ defined in 2a, $\tilde{\beta}=\Sigma^{1/2}\beta$, and $\Lambda_{j}$ denotes the $j$th column of $\Lambda$. (I) satisfies Definition 1 with probability exceeding $1-O(\exp(-n))$ under 2a by the proof in Section 6.2.1 of Javanmard and Montanari, (2014). For (II), since $\Lambda$ and $F$ have identically distributed columns,

[TABLE]

Taking the expectation over $\Lambda_{j}$, using Jensen’s inequality and the fact that $\Lambda$ has subgaussian columns by 2b, we have that, uniformly over $\beta$,

[TABLE]

Since under 2b, $\lambda_{\max}(\Sigma^{R})$ is finite, we conclude that $(II)$ is $O(\log p/n)$ with probability $1-O(p^{-1})$.

C.3 Proof of Lemma B.3

Normality follows from 4. What remains is to show that $\hat{\Omega}_{jj}=O_{p}(1)$.

[TABLE]

By Cauchy-Schwarz,

[TABLE]

It follows that

[TABLE]

By Lemma A.3, with probability at least $1-O(\exp(-n))$, we have $\hat{\Omega}_{jj}\geq C^{-1}\sigma^{2}$. Taking a union bound over $j=1,\ldots,p$, we find that for $n$ sufficiently larger than $\log p$, $\min_{i=1,\ldots,p}\hat{\Omega}_{jj}\geq C_{1}$ for some positive constant $C_{1}$ with probability $1-O(\exp(-n))$.

An upper bound follows similarly. By Lemma A.4 we have that with probability at least $1-O(\exp(-n)$,

[TABLE]

Then, by Lemma A.1, with probability at least $1-O(\exp(-n))$, $\hat{\Omega}_{jj}\leq C\frac{1}{1-\epsilon_{u}}$. By 1 $n$ is sufficiently larger than $\log p$, so applying a union bound over $p$ gives $\max_{j=1,\ldots,p}\hat{\Omega}_{jj}<C$ for some positive constant $C$ with probability $1-O(\exp(-n))$. We conclude that $\hat{\Omega}_{jj}$ is $O_{p}(1)$ for $j=1,\ldots,p$.

C.4 Proof of Lemma B.4

Consider the $(i,j)$-th element of $\hat{M}X$,

[TABLE]

By substituting $X=VSU^{\top}=VS_{n}U_{n}^{\top}$, where $U^{\top}U=UU^{\top}=I_{p}$, $U_{n}^{\top}U_{n}=I_{p}$, ${Q}=S^{\top}S$, and $Q_{n}=S_{n}^{\top}S_{n}$ we have

[TABLE]

where $U_{n}U_{n}^{\top}=X^{\top}(XX^{\top})^{-1}X^{\top}$, $[W]_{ii}=\frac{{q}_{i}}{{q}_{i}+\gamma}$, and $[W]_{ii}-1=-\frac{\gamma}{\gamma+{q}_{i}}$, for $i=1,\ldots,n$.

Now,

[TABLE]

with the last line holding with probability $1-O(\exp(-n))$ by Lemma A.4. Choosing $\gamma=O\left(p\sqrt{\frac{\log(p)}{n}}\right)$, we have that with probability $1-O(\exp(-n))$,

[TABLE]

Then, by Cauchy-Schwarz, $e_{i}^{\top}U_{n}(W-I_{n})U_{n}^{\top}e_{j}/(e_{i}^{\top}U_{n}U_{n}^{\top}e_{i})=O(\sqrt{\log p/n})$ with probability $1-O(\exp(-n))$. In that case, $[MX]_{i,j}=\frac{e_{i}^{\top}U_{n}U_{n}^{\top}e_{j}}{e_{i}^{\top}U_{n}U_{n}^{\top}e_{i}}+O(\sqrt{\log p/n})$ with probability $1-O(\exp(-n))$, and Lemma B.4 then follows from Lemma B.1.

C.5 Proof of Lemma B.5

Normality follows from 4. What remains is to show that $\hat{\Omega}_{jj}=O_{p}(1)$.

[TABLE]

First the lower bound. By Cauchy-Schwarz

[TABLE]

We therefore have

[TABLE]

and the lower bound follows from Lemma A.3.

We now continue with an upper bound. Define $\bar{W}$ as a diagonal matrix with $[\bar{W}]_{ii}=\frac{1}{\gamma+{q}_{i}}$. Note that by Lemma A.4, $[\bar{W}]_{ii}\leq Cp^{-1}$ with probability $1-O(\exp(-n))$. The numerator of (C.21) is then equal to

[TABLE]

with probability $1-O(\exp(-n))$.

From this we have that with probability $1-O(\exp(-n))$,

[TABLE]

where we used that $I_{n}-W=({Q}_{n}+\gamma I_{n})^{-1}{Q}_{n}$ and

[TABLE]

with the last inequality following from 1. Under 3, we can now apply Lemma A.1 and find that with probability $1-O(\exp(-n))$,

[TABLE]

The proof now follows from a union bound.

Appendix D Results Monte Carlo experiments