WIXLIB

whose properties have now been well understood [see e.g. Bühlmann and van de Geer (2011);Hastie et al. (2015)]. We also try the alternative Lasso formulation without penalizing βj in ournumerical studies and find that it does not improve the finite sample performance. Now given aprojection vector vj (vj,1 , . . . , vj,n ) Rn 1 such that vj Xj n, we define the projection-basedestimator for βj asβ̃j (vj ) : 11 v η̂j βj vj R(vj , β j ),n jn(4)where R(vj , β j ) n 1 vj X j (β j β̂ j ) is the bias term caused by the estimation effect. (4)implies that 1n(β̃j (vj ) βj ) vj nR(vj , β j ).nTo ensure that β̃j (vj ) has asymptotically tractable limiting distribution, we require the bias term nR(vj , β j ) to be dominated by the leading term n 1/2 vj , which converges to a normal limit un der suitable assumptions. In other words, the bias term nR(vj , β j ) controls the non-Gaussianity of β̃j (vj ). A practical challenge here is that the bias nR(vj , β j ) can be hardly estimated directly from the data. It is common in the literature to replace nR(vj , β j ) by a conservative estimatorusing the l1 l bound, i.e., k n(β j β̂ j )k1 kn 1 vj X j k .(5)See Zhang and Zhang (2014), van de Geer et al. (2014), Javanmard and Montanari (2014). Wenote that the variance of n 1/2 vj is equal to σ 2 n 1 kvj k2 . To achieve efficiency, we shall also try to minimize σ 2 n 1 kvj k2 given that the bias nR(vj , β j ) is properly controlled. Because the firstterm in (5) is independent of vj , we can seek a projection direction to minimize a linear combinationof kn 1 vj X j k2 and the variance σ 2 n 1 kvj k2 . However, the l1 l bound on the whole bias termcould be conservative as it does not take into account the specific form of the bias term. We notethat the bias term can be written as 1 X nR(vj , β j ) vj Xk (βk β̂k )nk6 j1 nX(1)k Sj (ν)1vj Xk (βk β̂k ) nXvj Xk (βk β̂k )(6)(2)k Sj (ν) nR(1) (vj , β j ) nR(2) (vj , β j ),(1)(2)where Sj (ν) : S(ν) \ {j} and Sj (ν) : S(ν){ \ {j} denote the index sets (except j) associatedwith the strong and weak signals respectively for S(ν) : {k : βk ν} and both R(1) (vj , β j )and R(2) (vj , β j ) are defined accordingly. Here ν is a threshold that separates the coefficients intotwo-groups namely the group with strong signals and the group with weak or zero signal. For5

pexample, one can set ν c0 log(p)/n for some large enough constant c0 , which is the minimaxrate for support recovery.The formulation (6) using the decomposition associated with signal strengths can be empricallymotivated. Specifically, it generally provides a smaller bias than the one without such decomposition with the simulated data. Figure 4 illustrates one such representative case where we make acomparison of the biases for projection vectors calculated based on two different methods: the onesolves (8) by using the estimated set of strong signals as in Section 3.1 (denoted by “With Decom(1)position”) and the other one solves the same problem but with Aj (denoted by “WithoutDecomposition”). It can be seen that “With Decomposition” shows a smaller bias than “WithoutDecomposition.” Similar results could be observed in various simulation settings.2.2A new projection directionIn this subsection, we propose a novel formulation to find the projection direction. When n, we have the freedom to choose vj to make the term kn 1 vj XS (1) (ν) k arbitrarily(1) Sj (ν) j(1)small. In fact, we can always choose vj such that it is orthogonal to all Xk with k Sj (ν). Thebasic idea here is to find a projection direction vj such that it is “more orthogonal” to the spacespanned by {Xk }k S (1) (ν) as compared to the space spanned by {Xk }k S (2) (ν) . With this intuitionjjin our mind and the goal to balance the squared bias with the variance, we formulate the followingoptimization problem min γ1 max n 1 vj Xk 2 γ2 max n 1 vj Xk 2 σ 2 n 1 kvj k2 ,vj(1)(2)k Sj (ν)k Sj (ν)s.t. vj Xj n,(7)where γ1 , γ2 0 are tuning parameters which control the trade-off between the squared bias andthe variance. The term γ1 maxk S (1) (ν) n 1 vj Xk 2 (γ2 maxk S (2) (ν) n 1 vj Xk 2 ) corresponds to thejj2 (R2 ). By introducing two ancillary variables u , u , (7) can be cast intol1 l bound for R(1)j1 j2(2)the following quadratic programming problemminuj1 ,uj2 ,vj(γ1 u2j1 γ2 u2j2 σ 2 n 1 kvj k2 ),s.t. vj Xj n,(1) uj1 n 1 vj Xk uj1 ,k Sj (ν), uj2 n 1 vj Xk uj2 ,k Sj (ν),(2)which can be solved efficiently using existing quadratic programming solver.(1)(1)The set Sj (ν) is generally unknown and needs to be replaced by a surrogate set Aj with(1)(1) Aj n. In Section 3.1, we describe a method to select Aj6based on a preliminary projection-

(1)based estimators. We show that Ajconverges asymptotically to a nonrandom limit, i.e., (1)(1)P Aj Bj 1,(1)(1)(1)for a nonrandom subset Bj of [p]. We remark that Bj does not need to agree with Sj (ν) forour procedure to be valid. To ensure that the remainder term is negligible, the theoretical analysis in Section 3.2 suggests that γ1 and γ2 should both be of the order O σ 2 n/ log p . Combining theabove discussions, we now state the optimization problem for obtaining the optimal projectiondirection n 2n 2 12u C2u n kvj k ,minC1uj1 ,uj2 ,vjlog p j1log p j2s.t. vj Xj n, uj1 n 1 vj Xk(8) uj1 ,k uj2 n 1 vj Xk uj2 ,k (1)Aj ,(2)Aj , (2)(1) {where Aj : Aj\ {j} and C1 , C2 0 are tuning parameters whose choice will be discussedin Section 5.Remark 1. A related method is the refitted Lasso by Liu and Yu (2013). The idea is to refit themodel selected by the Lasso and conduct inference based on the refitted least squares estimator.Such an estimator fits into the framework of the projection-based estimators. To see this, let Ŝbe the set of active variables selected by the Lasso and note that β̂k 0 for k / Ŝ. For eachj Ŝ, let ŵj be the projection of Xj onto the orthogonal space of XŜ\{j} . Then the refittedleast squares estimator is given by ŵj (Y X j β̂ j )/(ŵj Xj ). It is easy to see that the bias forPthe refitted least squares estimator is proportional to k / Ŝ ŵj Xk βk , which disappears when theselected model contains all significant variables. However, when the model selection consistencyfails, such a procedure is no longer valid due to the nonnegligible bias.33.1MethodologySurrogate setWe describe a procedure to estimate the set of strong signals based on a preliminary projectionbased estimator. It should be noted that the estimator here is different from the original debiasedLasso because it is based on the novel formulation (8). Specifically, for some τ 0, we define ourestimate for the set of strong signals as pA(τ ) : {l : Tl τ log p}7whereTl nβ̃l (v̂l )σ̂n 1/2 v̂l (9)

where σ̂ is an estimator of the noise level σ and β̃l (v̂l ) is a projection-based estimator with v̂l beingthe solution to the following optimization problem n 2 12C0u n kvl k ,minul ,vllog p l(10)s.t. vl Xl n, ul n 1 vl Xk ul ,k 6 l.In practice, both C0 and τ need to be appropriately chosen. The details for the selection arediscussed in Section 8.1. Note that (10) is a special case of (8) when we have no knowledge about(1)the set of strong signals, that is, Al . We define the surrogate sets to be(1)(2)Aj (τ ) : A(τ ){ \ {j}.Aj (τ ) : A(τ ) \ {j},(11)Throughout the paper, we consider the variance estimatorσ̂ 2 1kY Xβ̂k2n(12)which appears to outperform an alternative estimator kY Xβ̂k2 /(n kβ̂k0 ) studied in Reid etal. (2016), see Figure 22 in Supplementary Material for a comparison. Before presenting the mainresult of this subsection, we introduce some assumptions.Assumption 1. There exist a set B [p] {1, 2, . . . , p} and 0 d0 d1 such that nβl pmax d0 log p,σl B{ nβl pmin d1 log p.l BσAssumption 2. The error is a mean-zero sub-gaussian random vector with the sub-gaussiannorm κ .Assumption 3. The preliminary estimator satisfies that nkβ̂ βk1 Op (s0plog(p)).pAssumption 4. The variance estimator σ̂ 2 is consistent in the sense that σ̂/σ 1.Assumption 5. Suppose the design matrix X Rn p has i.i.d. rows with zero population meanand covariance matrix Σ (Σi,j )pi,j 1 . Assume that1. maxj Σj,j ;2. λmin (Σ) Λmin 0;3. The rows of X are sub-gaussian with the sub-gaussian norm κ .8

Assumption 6. n, p and s0 satisfy the rate condition s0 log p/ n o(1).Assumption 1 allows the strengths of strong and weak signals to be the same order and thus ismuch weaker than the “beta-min” condition which requires the weak signals to be of smaller order.Assumptions 3 and 4 are satisfied for the Lasso estimator and the variance estimator σ̂ in (12) undersuitable regularity conditions [Bühlmann and van de Geer (2011)]. Assumptions 2 and 5 requirethe error and design to be sub-gaussian. Similar assumptions have been made in van de Geer etal. (2014). Like Javanmard and Montanari (2014), the validity of our method does not rely on thesparsity of the precision matrix of the design, which is required in the nodewise Lasso regressionfor the original debiased Lasso. In view of Cai and Guo (2017), the rate condition in Assumption6 cannot be relaxed without extra information. Zhu and Bradic (2018a, 2018b) proposed testingprocedures in high-dimensional linear models which impose much weaker restrictions on modelsparsity or the loading vector representing the hypothesis. However, their methods require certainauxiliary sparse models, which are not needed for our procedure. q 12ΣΛDefine Σj\ j Σj,j Σj, j Σ 1Σandκ 21 0j j, j j,jmin j,j κ for 1 j p. The(1)following proposition shows that the surrogate set Aj (τ ) with a properly chosen τ converges toB \ {j}.(1)(2)Proposition 1. Define Aj (τ ) and Aj (τ ) as in (11) and let v̂l be the solution to (10) for l 6 j.Suppose d0 , d1 and τ satisfyqσ2 (τ d0 max Σl,l )2 1l32eκ2 p and d1 /M τ 0 where 2 M min Σl\ l 2C01 l p11min1 l p 8e2 (κ0l )2! 1 max Σl\ l1 l p.Then under Assumptions 1-6, we have pP max Tl τ log p 1, (2)l Bj pP min Tl τ log p 1,(1)l Bj(1)where Bj(2): B \ {j} and Bj (1) {(1)(1): Bj\ {j}. As a consequence, P Aj (τ ) Bj 1.Remark 2. As shown in Proposition 1, the surrogate set in (11) has an asymptotic (nonrandom)limit, which implies that the projection direction obtained in (8) is asymptotically independent ofthe random error . This fact is useful in the proof of Theorem 1 later. To ensure the independencebetween the projection direction and the random error, we can also employ the sample splittingstrategy, i.e., we split the samples into two subsamples, estimate the set of strong signals based on9

the first subsample and construct the projection-based estimator based on another subsample. Aswe use all samples in building the projection-based estimator, our method is more efficient thanthe sample splitting strategy.Remark 3. When d0 0, B coincides with the support of β. Proposition 1 suggests that one canconsistently recover the support of β by thresholding the projection-based estimator.3.2Bias reducing projection (BRP) estimatorIn this subsection, we introduce the bias reducing projection (BRP) estimator and study itsasymptotic behavior. Let ṽj be the solution to (8) based the surrogate sets in (11). Then the BRPestimator β̃j (ṽj ) is defined asβ̃j (ṽj ) 11 ṽj η̂j ṽj (Y X j β̂ j ).nnIn the following, we introduce the two asymptotic results depending on whether the surrogate set isestimated from the same data set used to find the projection direction. We first state the followingtheorem on the asymptotic normality when the surrogate set is estimated via (11).(1)(2)Theorem 1. Denote by ṽj the solution to (8) with Aj (τ ) and Aj (τ ) in (11). Suppose theassumptions in Proposition 1 hold and further assume that for some δ 0,kṽj k2 δ oa.s. (kṽj k).Then we have n β̃j (ṽj ) βjσ̂n 1/2 kṽj k(13)d N (0, 1).(14)Thus an asymptotic 100(1 α)% confidence interval for βj is given by (CI(1 α) b R:n(β̃j (ṽj ) b) z1 α/2σ̂n 1/2 kṽj k),(15)where z1 α/2 is the 1 α/2 quantile of N (0, 1).(13) is a Lyapunov type condition which implies the central limit theorem. This type of assumption regarding the projection direction has also been imposed in Dezeure et al. (2017). It canbe dropped under the Gaussian assumption on the errors. If the surrogate set is chosen based onprior knowledge or estimated from an independent data set (e.g., based on sample splitting), thenAssumptions 1-2 can be relaxed and we have the following result.(1)Corollary 1. Suppose the surrogate set Aj is independent of the data. Under Assumptions 3-6and further assuming that for some δ 0, E[ i 2 δ ] and kṽj k2 δ oa.s. (kṽj k), then (14) stillholds.10

3.3Modified bias reducing projection (MBRP) estimatorWe introduce a modified bias reducing projection (MBRP) estimator which is motivated byProposition 1 and the refitted Lasso idea. This new estimator would lead to a potentially smallerorder of bias compared to that of the original debiased Lasso estimator under suitable assumptionsas shown in Proposition 2. Thus, it is expected to provide better empirical coverage probability.See more details in Section 6. To motivate the MBRP estimator, we note that the bias associatedwith the BRP estimator based on some estimator β̌ for β can be written as 1 X nR(vj , β j ) vj Xk (βk β̌k )nk6 j1 X 1 X vj Xk (βk β̌k ) vj Xk (βk β̌k )nn(1)(2)k Bj(1)k Bj(2)(1)where Bj , Bj are the same as in Proposition 1. When Bj n, we can always require vj to beexactly orthogonal to XB(1) . So, the bias associated with the set of strong signals becomes zero.j(2)Thus it suffices to control the bias term associated with Bj by properly choosing vj and β̌, whichwill be clarified below.To find the projection direction for the MBRP estimator, we consider the optimization problem minuj2 ,vj n 2 12u n kvj k ,C2log p j2s.t. vj Xj n,n 1 vj Xk(16)k 0,(1)Aj , uj2 n 1 vj Xk uj2 ,(2)k Aj .Different from (8), we require the projection direction to be orthogonal to the column space ofXA(1) in (16). Instead of using the Lasso estimator β̂, we shall adopt the refitted least squaresjestimator β̌ as our preliminary estimator, i.e.,β̌A(1) argminjβ̃1kY XA(1) β̃k2 ,2njβ̌A(2) 0.(17)jThe MBRP estimator is then defined asβ̃j (v̄j ) 1 1v̄ (Y X j β̌ j ) βj v̄j R(v̄j , β j )n jn(18)where R(v̄j , β j ) n 1 v̄j X j (β j β̌ j ) and v̄j is the solution to problem (16). The MBRPestimator can be viewed as an intermediate estimator between the refitted Lasso and the BRPestimator based on (8). While (16) is a variant of (8) seeking for a projection direction that isexactly orthogonal to the column space of XA(1) , the modified procedure uses the refitted estimatorj11

for β as the refitted Lasso does as noted in Remark 1. We argue that the bias term nR(v̄j , β j ) which controls non-Gaussianity could have a potnetially smaller order compared to that of the original debiased Lasso estimator in the following.(1)(2)Proposition 2. Denote by v̄j the solution to (16) with Aj (τ ) and Aj (τ ) defined in (11). Let β̌(2)be the refitted least square estimator in (17). Conditional on the event {Aj(2) Bj }, we have p log p nR(v̄j , β j ) Opd0 kβB(2) k0 nj(19)under Assumptions 1 and 5. If we further assume thatpd0 kβB(2) k0 o(s0 ),(20)j the bias nR(v̄j , β j ) is asymptotically negligible with smaller order than that of the original debiased Lasso given by Op (s0 log p/ n).In particular, (20) holds if d0 o(1) and d1 O(1), i.e., the strength of weak signals is ofsmaller order compared to the strong signals. It is more stringent than Assumption 1 where themagnitudes of the set of strong signals and weak signals are allowed to be of the same order.However, it should be mentioned that Proposition 2 is not necessary for the asymptotic normalityin Corollary 2 to be achieved. The following result shows the asymptotic normality of (18) whichcan be proved by using similar arguments as those for Theorem 1.Corollary 2. Under the assumptions in Theorem 1, we have n β̃j (v̄j ) βjσ̂n 1/2 kv̄j kd N (0, 1),where β̃j (v̄j ) is defined in (18) and v̄j is the solution to (16).4Inference on a sparse linear combination of parametersIn some applications, one may be interested in conducting inference on a β for a (sparse) loadingvector a (a1 , . . . , ap ) Rp with kak0 s n. Denote by S S(a) {1 j p : aj 6 0} thesupport set of a. Our method can be generalized to construct estimator and conduct inference fora β a S βS . Recall that β̂ is the preliminary estimator of β. DefineηS Y X S β S XS βS andη̂S Y X S β̂ S XS βS X S (β S β̂ S ).12

We construct an estimator for a β in the form of n 1 va η̂S , where va (va,1 , . . . , va,n ) is aprojection direction such that n 1 va η̂S has tractable asymptotic limit. Notice thatn 1 va η̂S n 1 va XS βS n 1 va n 1 va X S (β S β̂ S ) 1 1 1 a S βS (n va XS aS )βS n va n va X S (β S β̂ S ).Under the equality constraint that n 1 va XS a S 0 and by rearranging the above terms, wehave 1/2 n(n 1 va η̂S a va S βS ) n nR(va , β S ),(21)where R(va , β S ) n 1 va X S (β S β̂ S ). Similar to (6), the bias term can be decomposed into(1)two parts corresponding to different strengths of the signals. Let AS be the surrogate set forthe set of strong signals (excluding the elements in S), which can be obtained in a similar way asdescribed in Section 3.1. Following the derivations in Section 2, we can formulate the followingoptimization problem to find va n 2n 2 12u C2u n kva k ,C1log p a1log p a2minua1 ,ua2 ,vas.t. va XS na S,(22)(1) ua1 n 1 va Xk ua1 ,k AS , ua2 n 1 va Xk ua2 ,k AS ,(2) {(2)(1)where AS : AS S . Denote by (ũa1 , ũa2 , ṽa ) the solution to (22). Our estimato

For statistical inference after model selection, Berk et al. (2013) developed a post-selection inference procedure by reducing the problem to one of simultaneous inference. Lockhart et al. (2014) constructed a statistic from the Lasso solution path and showed that it converges to a standard exponential distribution.