1 Introduction↩︎

Suppose \(A, B \in \mathbb{R}^{n \times p}\) are two matrices differing only in the \(j\)-th column. Let \(Y \in \mathbb{R}^n\) and let \(\rho \colon \mathbb{R} \to [0, +\infty)\) be a convex function, and define \[\begin{align} \hat{\alpha} &:= \operatorname*{argmin}_{\beta \in \mathbb{R}^p} \left\{ \frac{1}{2n} \|Y - A\beta\|_2^2 + \sum_{k=1}^p \rho(\beta_k) \right\}; \tag{1} \\ \hat{\beta} &:= \operatorname*{argmin}_{\beta \in \mathbb{R}^p} \left\{ \frac{1}{2n} \|Y - B\beta\|_2^2 + \sum_{k=1}^p \rho(\beta_k) \right\}. \tag{2} \end{align}\] Motivated in part by accelerating resampling-based variable selection methods [1], we consider how to efficiently solve the second optimization given the solution to the first. In general, there is no explicit formula for updating \(\hat{\alpha}\) to \(\hat{\beta}\), and approximations based on Taylor expansions are not accurate in high-dimensional settings with correlated features. It is therefore natural to consider alternative representations of the solutions that preserve stability under local updates. The main message of this paper is that debiased estimators satisfy an orthogonality principle that enables stable updates under fairly general conditions. In particular, we show that correlated sub-Gaussian designs are sufficient for the approximation to hold. To be precise, we establish approximation results for a generalized definition of the debiased estimator, which reduces to the standard definition in the literature under more restrictive assumptions.

The standard definition of debiased Lasso in the high-dimensional regime can be found in [2]. We focus on the Lasso case \(\rho(t)=\lambda|t|\) for \(\lambda \in [0,\infty)\), where the objective is convex but not strongly convex and lacks smoothness, leading to additional technical difficulties; the case of smooth and strongly convex \(\rho\) is discussed later and is technically simpler. Suppose that \(A\) has i.i.d.rows following the normal distribution \(\mathcal{N}(0,\Sigma)\), and \(Y=A\alpha +w\), where \(w\) is an independent Gaussian noise vector, and \(\alpha\in\mathbb{R}^p\). The “number of nonzero coefficients” is defined by \[\begin{align} k:=\|\chi^{\alpha}\|_0, \label{e95k} \end{align}\tag{3}\] where \[\begin{align} \psi^{\alpha}&:=\frac{1}{n\lambda}A^{\top}R, \quad R:=Y-A\hat{\alpha}, \tag{4} \\ \chi^{\alpha}_j&:=1\{\psi^{\alpha}_j=1\} - 1\{\psi^{\alpha}_j=-1\}. \tag{5} \end{align}\] Note that \(\psi^{\alpha}\) is the subgradient of the \(\ell_1\) norm, so \(k\ge \|\hat{\alpha}\|_0\), although equality is achieved in most cases. We use the definition 3 instead of \(\|\hat{\alpha}\|_0\) since \(\hat{\alpha}\) may not be unique, due to the lack of strong convexity of the optimization, whereas \(\psi^{\alpha}\), and hence \(\chi^{\alpha}\), is always uniquely defined. Similarly \(\chi^{\alpha}_j\) can be understood as the ‘essential sign’ of \(\hat{\alpha}_j\). Then, the debiased Lasso defined in [2] is \[\begin{align} \hat{\alpha}^u:=\hat{\alpha}+\frac{1}{n-k}\Sigma^{-1}A^{\top}(Y-A\hat{\alpha}). \label{eq95betau} \end{align}\tag{6}\] Under suitable conditions, it has been shown that \(\hat{\alpha}^u\approx \alpha+ \tau \Sigma^{-1/2}z\), where \(z\sim \mathcal{N}(0,I)\) and \(\tau\in(0,\infty)\) is a constant determined by a set of fixed point equations [2]. Rigorously establishing such Gaussian limit properties for general non-Gaussian \(A\) in the proportional growth regime (\(p\), \(n\), and the sparsity level \(s\) have fixed ratios) is open; see discussions Section 2.

In this work, we propose a generalization of 6 that extends its applicability under more relaxed assumptions. This generalization is similar to the orthogonality principle in the semiparametric and double machine learning literature [3], but we will discuss the differences in Section 2. Let \(A_{:\setminus j}\in\mathbb{R}^{n\times p}\) denote the matrix obtained by excluding the \(j\)-th column of \(A\). Define the “residualized column” \[\begin{align} \check{A}_{:j}&:=A_{:j}-\mu_{:j}. \label{eq68951} \end{align}\tag{7}\] where \(\mu_{:j}\in\mathbb{R}^n\) is an arbitrary vector. To control the error of the approximate formula, we will need to choose \(\mu_{:j}\) so that the residualized column is “orthogonal” to other columns in a sense that will be made precise in 19 (specific examples are provided in the itemized below). We then define a generalization of the debiased estimator in 6 :

\[\begin{align} \hat{\alpha}_j^U :=\hat{\alpha}_j+ \left( \frac{1}{n}\check{A}_{:j}^{\top} (I-P_{\mathcal{A}})A_{:j} \right)^{-1} \frac{\check{A}_{:j}^{\top}R}{n}, \label{eq95alpha95U} \end{align}\tag{8}\] where \[\begin{align} \mathcal{A}:=\{l\neq j\colon \chi^{\alpha}_l\neq 0\}, \end{align}\] and \(P_{\mathcal{A}}\) denotes the projection onto the columns of \(A\) corresponding to \(\mathcal{A}\) in the case of Lasso, and \(P_{\mathcal{A}}:=A_{:\setminus j}[A_{:\setminus j}^{\top}A_{:\setminus j}+\mathop{\mathrm{diag}}(\rho''(\hat{\alpha}_{\setminus j}))]^{-1}A_{:\setminus j}^{\top}\) in the case of general \(\rho\). The inverse in 8 is well-defined in the event of interest. Again, we adopt the convention in 3 when \(\hat{\alpha}\) is not unique. Note that the definition of \(\hat{\alpha}^U_j\) uses only \((Y,A)\) and has no reference to \(B\), hence we can use \(\hat{\alpha}^U_j\) to build \(t(j,A,Y)\). Specific examples of \(\mu_{:j}\) and the corresponding \(\hat{\alpha}^U\) include:

• If \(n>p\), we can take \(\mu_{:j}\) to be the projection of \(A_{:j}\) on the column space of \(A_{:\setminus j}\). Then \(\check{A}_{:j}^{\top}P_{\mathcal{A}}=0\), and we have \(\hat{\alpha}_j^U= \hat{\alpha}_j+(\check{A}_{:j}^{\top}A_{:j})^{-1}\check{A}_{:j}^{\top}(Y-A\hat{\alpha}) =\|\check{A}_{:j}^{\top}\|_2^{-2}\check{A}_{:j}^{\top}Y=[(A^{\top}A)^{-1}A^{\top}]_{j:}Y\) which is the ordinary least squares estimator. In the specific case of orthogonal designs, we simply take \(\mu_{:j}=0\) (see Example 1).

• If \(n<p\), the above choice of \(\mu_{:j}\) leads to a degenerate definition of \(\hat{\alpha}_j^U\). However, due to the sparsity of \(\hat{\alpha}\), other choices of \(\mu_{:j}\) may still ensure 19 . For the random design case where the rows of \(A\) are independent, we can choose \(\mu_{:j}=\mathbb{E}[A_{:j}|A_{:\setminus j}]\). Specifically for Gaussian \(A\), we have \(\check{A}_{:j}:=A_{:j}-A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj}\) in terms of the covariance matrix \(\Sigma\).

Contributions

We will show that \(\hat{\alpha}^U_j\approx \hat{\alpha}^u_j\) under suitable assumptions, but the definition of \(\hat{\alpha}^U_j\) is applicable in more general settings. Under increasingly stronger assumptions, we establish the following approximation results:

• Explicit nonasymptotic error bounds for any given design matrix: Define \(\hat{\beta}^U_j\) analogously to \(\hat{\alpha}^U_j\), i.e. \[\begin{align} \hat{\beta}_j^U :=\hat{\beta}_j+ \left( \frac{1}{n}\check{B}_{:j}^{\top} (I-P_{\mathcal{B}})B_{:j} \right)^{-1} \frac{\check{B}_{:j}^{\top}S}{n}, \label{eq95beta95U} \end{align}\tag{9}\] where \(\check{B}_{:j}:=B_{:j}-\mu_{:j}\), \(S:=Y-B\hat{\beta}\), and \(P_{\mathcal{B}}\) denotes the projection onto the columns of \(A\) (equivalently, columns of \(B\) since the two matrices differ only in the \(j\)-th column) corresponding to \(\mathcal{B}:=\{l\neq j\colon \chi^{\beta}_l\neq 0\}\). For any given \((A,B,Y)\), we can show \[\begin{align} \check{B}_{:j}^{\top} (I-P_{\mathcal{B}}) B_{:j} \hat{\beta}^U_j &\approx \check{B}^{\top}_{:j}R + \check{B}^{\top}_{:j}(I-P_{\mathcal{A}})A_{:j} \hat{\alpha}_j \label{e95approx} \end{align}\tag{10}\] with error controlled in Theorem 1. The linear combination on the right side of 9 orthogonalizes the estimator so that the approximation error in 10 becomes the inner product of two incoherent vectors, canceling the leading error term. In Section 2, we discuss the connection between this orthogonalization and the one in the double machine learning literature. If \(B_{:j}\) and \(A_{:j}\) are independent conditioned on \(A_{:\setminus j}\), we can further show that the right side of 10 is approximately \(\check{B}^{\top}_{:j}R\), although the right side of 10 is already computable without using \(\hat{\beta}\). Thus, for any given \(Y\), if we define \[\begin{align} t(j,A,Y) &:=\frac{1}{n}\check{A}_{:j}^{\top} (I-P_{\mathcal{A}}) A_{:j}\hat{\alpha}_j^U; \tag{11} \\ t(j,B,Y) &:=\frac{1}{n}\check{B}_{:j}^{\top} (I-P_{\mathcal{B}}) B_{:j}\hat{\beta}_j^U, \tag{12} \end{align}\] then we can compute both \(t(j,A,Y)\) and approximately \(t(j,B,Y)\) using only \(\hat{\alpha}\) and without \(\hat{\beta}\).

• Convergence for sub-Gaussian designs: Consider the case where \(Y=A\alpha+w\) and design matrices have i.i.d.rows with covariance of bounded max and min eigenvalues and with bounded sub-Gaussian variance proxy. Also assume that \(B_{:j}\) and \(A_{:j}\) are i.i.d.given \(A_{:\setminus j}\). In this setting, we show that the approximation error in 10 vanishes asymptotically for almost all \(j\) (see Definition 1 and Theorem 4). Our proof only uses certain concentration and anti-concentration properties to give order-wise control of quantities, rather than more precise calculation of limits, which often requires stronger assumptions (see Section 2).

• Reduction to \(\hat{\beta}_j^u\): Further assuming that \(\frac{1}{n}\check{B}_{:j}^{\top} (I-P_{\mathcal{A}})B_{:j}\) is bounded away from 0 (which is true when \(\mathbb{E}[|\check{B}_{1j}|^2]\) is bounded away from 0; see 27 ), we have from 10 that \[\begin{align} \hat{\beta}^U_j &\approx \frac{\check{B}^{\top}_{:j}R + \check{B}^{\top}_{:j}(I-P_{\mathcal{A}})A_{:j} \hat{\alpha}_j}{\check{B}_{:j}^{\top} (I-P_{\mathcal{A}}) B_{:j}}. \label{e95approx95na} \end{align}\tag{13}\] If we further have \[\begin{align} \frac{1}{n}\check{B}_{:j}^{\top} (I-P_{\mathcal{A}})B_{:j} \approx (1-\frac{k}{n})\Sigma_{j|\setminus j}, \label{e95app95chi} \end{align}\tag{14}\] then \(\hat{\beta}_j^U\approx \hat{\beta}_j^u\). For example, 14 is true in the case of Gaussian feature vector, as a consequence of concentration of the chi-square distribution; see Theorem 7 and Remark 2). In more general settings, 14 may no longer hold, since \(\mathbb{E}[B_{:j}|A_{:\setminus j}]=A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj}\) may fail (see Remark 10). This is our main motivation for introducing \(\hat{\beta}_j^U\).

We apply the approximation formulas to accelerate resampling-based variable selection procedures. Consider the proportional regime where \(n\) grows linearly with \(p\), and suppose that solving the Lasso requires computational cost \(L\) (for example, \(L = O(p^3)\) via least angle regression [4]). The knockoff filter requires solving a single regression and hence has complexity \(\Theta(L)\), but its statistical power can be reduced due to the increase in the number of variables from \(p\) to \(2p\). We consider two alternative methods that avoid regressing on \(2p\) variables but require solving multiple regression problems, and show that their computational cost can be substantially reduced using the proposed approximation formulas:

• We consider “local knockoff filter”, a variant of the knockoff method that resamples only a single feature in each Lasso run. A concurrent work, [5], proposed this strategy under the name “one-at-a-time knockoffs” and showed experimentally that it is “substantially more powerful.” While [5] focused on least squares and ridge regression, exploiting the explicit solution and rank-one matrix update formulas, our work extends the method to Lasso regression while maintaining computational complexity \(\Theta(L + p^3)\).

• Conditional randomization test (CRT) iteratively updates each feature vector (column of the design matrix) with \(K\) conditionally independent sample and computes the corresponding test statistics (e.g., \(K=\Theta(p)\) in [1]). When CRT is applied to the variable selection problem, its asymptotic power is often observed to be worse than that of CRT [1], [6]. However, it has a complexity of \(\Theta(LpK)\). To reduce the complexity, distilled CRT has been proposed, and we will see in Section 2.4 that \(\rm d_0CRT\) from [7] is asymptotically equivalent to CRT using the debiased coefficients. Implementing \(\rm d_0CRT\) for FDR control requires complexity at least \(\Omega(pL)\), for solving \(p\) leave-one-out regression problems. In contrast, using 10 the complexity is reduced to \(O(L+p^2K)\).

Since previously the debiased Lasso \(\hat{\alpha}^u\) often appears in the literature on asymptotic normality, and asymptotic normality results can be used to directly estimate the \(p\)-value of \(\hat{\alpha}^u\), one might ask what is the benefit of resampling \(B_{:j}\) and using \(\hat{\beta}^U_j\) to estimate the \(p\)-values. The answer is that asymptotic normality requires more stringent conditions than the validity of the update formula. One simple example is the limiting case where the Lasso is reduced to a least square problem (\(n>p\) and \(\lambda\to 0\)). In this case, 10 is in fact exact equality regardless of the distributions and the dimensions, whereas normality is approximate and requires distributional assumptions for asymptotic convergence. For the general \(\lambda\) case, Example 1 ahead shows that asymptotic normality is false when the column of the design matrix contains a non-vanishing component, yet our main result (Theorem 1) is still capable of controlling the error of the update formula. For the case of i.i.d.rows of the design matrix, as mentioned before, our proof of the approximation in 10 only uses certain concentration and anti-concentration properties to control the order of the errors, rather than more precise characterization of limits such as Gaussian convergence. Indeed, our Theorem 4 shows asymptotic approximation assuming that the covariance matrix \(\Sigma\) of \(B_{1:}\) has bounded conditional numbers, and that \(B_{1:}\) is a sub-Gaussian vector. In contrast, a Gaussian limit result for \(\hat{\alpha}^u\) with correlated sub-Gaussian designs has been noted as an open problem (see Section 2).

Organization

Further related work is discussed in Section 2. In Section 3, we present main results on general nonasymptotic error bounds and asymptotic analysis for the sub-Gaussian case. Section 4 proves the nonasymptotic error bound. In Section 5, applies the approximation formulas in accelerating resampling-based variable selection algorithms. Section 6 concludes with an outlook for future directions. In the appendices, we present omitted proofs, experimental validations of the approximation errors in the debiased and the standard Lasso, and FDR control for synthetic and semi-real (riboflavin and HIV) data.

2 Related work↩︎

2.1 Definitions of the debiased estimator↩︎

Debiasing the Lasso for inference was suggested by [8], [9], [10], and [2]. The replica analysis heuristic calculation in [2] was perhaps the first to show that \(\alpha^u\) in 6 satisfies asymptotic normality in the proportional growth regime, with i.i.d.\(\mathcal{N}(0,\Sigma)\) rows in the design matrix. More specifically, in a suitable sense there is the approximation \[\begin{align} \hat{\alpha}^u-\alpha \approx \tau\Sigma^{-1/2}z \label{e95gapprox} \end{align}\tag{15}\] for some \(z\in\mathcal{N}(0,I)\), where \(\alpha\) is the ground truth, \(Y=A\alpha+w\), \(w\sim \mathcal{N}(0,nI)\), and \(\tau\) is the solution to a fixed point equation. The replica calculation, as well as other existing proof techniques for 15 in the proportional regime (see Section 2.3), are limited to the random design settings, and require more restricted assumptions than our proof of 10 .

In the landmark paper [3], a generalized debiased estimator similar in spirit to 10 was proposed, which was dubbed “double machine learning” since it involves regressing the outcome on the controls and regressing the treatment on the controls (\(\mu_{:j}\) in 7 ). To explain the idea, suppose that \(p\) is fixed, \(n\to\infty\), and the log likelihood is given by \(\ell(y,b,\beta)=-\frac{1}{2}(y-b^{\top}\beta)^2\) where \(b\in\mathbb{R}^p\). The approach of [3] (see also [11]) then defines a score \(\psi(y,b,\beta,\mu):=\partial_{\beta_j}\ell(y,b,\beta)-\mu\partial_{\beta_{\setminus j}}\ell(y,b,\beta)\) satisfying a Neyman orthogonality condition, ensuring first-order insensitivity with respect to the parameters \((\mu,\beta_{\setminus j})\). Data points indexed by the complement \(I^c\) of a set \(I\subseteq \{1,\dots,n\}\) are used to learn \(\hat{\beta}_{\setminus j}(I^c)\) and \(\mu(I^c)= \sum_{i\in I^c}B_{ij}^{\top}B_{i\setminus j}(\sum_{i\in I^c}B_{i\setminus j}^{\top}B_{i\setminus j})^{-1}\), the latter following from the orthogonality condition. Then \(\psi(y,b,\hat{\beta}_{\setminus j}(I^c),\beta_j,\mu(I^c)) =(b_j-\mu(I^c)b_{\setminus j})(y-b_{\setminus j}^{\top}\hat{\beta}_{\setminus j}(I^c)-b_j\beta_j)\), and the debiased estimator is defined by solving \(\sum_{i\in I}\psi(Y_i,B_{i:},\hat{\beta}_{\setminus j}(I^c),\beta_j,\mu(I^c))=0\) for \(\beta_j\), which is \[\begin{align} \frac{\sum_{i\in I}(B_{ij}-\mu(I^c)B_{i\setminus j}^{\top})(Y_i-B_{i\setminus j}\hat{\beta}_{\setminus j}(I^c))}{\sum_{i\in I}(B_{ij}-\mu(I^c)B_{i\setminus j}^{\top})B_{ij}}. \label{e17} \end{align}\tag{16}\] If the sample-splitting estimator \(\hat{\beta}_{\setminus j}(I^c)\) in 16 is replaced by the leave-one-out estimator, 16 becomes similar to the right side of 10 (note that \(\hat{\alpha}_j=0\) in 10 when \(A_{:j}=0\)). However, 10 also has a factor \(I-P_{\mathcal{B}}\) arising from the \(\ell_1\) penalization, which is unique for the high-dimensional (proportional) regime and was not discussed in [3]. Our result 10 establishes the approximate equivalence of orthogonalized estimators derived from two seemingly different principles:

• Robustness is the principle behind constructing debiased estimators in double machine learning. In 16 , if the nuisance parameter estimator \(\hat{\beta}_{\setminus j}(I^c)\) has an error, the error in \(Y_i-B_{i\setminus j}\hat{\beta}_{\setminus j}(I^c)\) will be in \(\mathop{\mathrm{col}}(B_{i\setminus j})_{i\in I}\), which is incoherent with \((B_{ij}-\mu(I^c)B_{i\setminus j}^{\top})_{i\in I}\), ensuring that 16 is robust.

• Stability in the literature usually refers to the regularity of a function under the perturbation in one entry. In 10 , the observations in one feature are changed, yet the estimator can be updated in a well-controlled manner, manifesting a stability property.

2.2 Leave-one-out analysis↩︎

Leave-one-out analysis is a fruitful approach for establishing limiting distributions or algorithmic properties of regression [[12]][13][14], and is closely related to techniques of the present paper. In [12], it is shown using the leave-one-out technique that the M-estimator converges asymptotically to a normal distribution (see also [15] and [16]). The problem considered there is different from the distribution of the Lasso considered in the present paper: the M-estimation problem concerns the \(n>p\) regime, and there is no need for debiasing; the asymptotic normality follows immediately from the rotation invariance of the distribution. We remark that a duality between M-estimation estimation and penalized least squares was mentioned in [17]. However, the duality only applies when the design matrix of the lasso has orthonormal rows, which does not cover the setting of the present paper.

A leave-one-out analysis for the Lasso was carried out in [18]. In addition to bounded singular values of \(\Sigma\), their analysis requires bounded \(\ell_1\) norms of the rows of the inverses of the submatrices of \(\Sigma\) (see [18]). The latter condition can be more restrictive than ours in Definition 1: for example a \(k\times k\) random matrix with independent entries of scale \(1/\sqrt{k}\) has spectral norm of order \(O(1)\), yet the \(\ell_1\) norm of each of its row has order \(\Omega(\sqrt{k})\) which is unbounded. Furthermore, [18] requires a sublinear sparsity level \(s_0=o(n/(\ln p)^2)\). In that regime, there is no need for the degrees of freedom adjustment factor \(\frac{1}{n-k}\) in 6 , and in fact in the approximation formula 10 it suffices to replace \(R\) with the noise \(w\) (see [18]). An extension of the analysis was done in [19], where the role of degrees-of-freedom adjustment was highlighted for sparsity level \(s_0\gg n^{2/3}\), but still \(s_0/p\to 0\) is required. In contrast, the present paper considers the regime of proportional sparsity level.

2.3 Asymptotic normality in the high-dimensional regimes↩︎

It appears that the first asymptotic normality result for debiased Lasso estimates in the proportional regime for correlated designs was derived in [20] (see the discussions therein). The technique of [20] (see also [21]) was based on the Second Order Stein theorems bounding the non-Gaussianity of a random variable of the form \(z^{\top}f(z)-{\rm div}f(z)\), where \(z\sim\mathcal{N}(0,I)\). To apply it to the debiased Lasso problem, consider given \(A_{:\setminus j}\) (the submatrix of \(B\) formed by excluding the \(j\)-th column), ground truth \(\alpha\) and noise \(w\). Let \(A_{:j}:=A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj}+z\sqrt{(\Sigma^{-1})_{jj}^{-1}}\) and \(Y=A\alpha+w\), which are both viewed as functions of \(z\). Let \(f(z):=A\hat{\alpha}-Y\). Then it can be verified that \(z^{\top}f(z)-{\rm div}f(z)\) is the debiased Lasso estimate up to a linear transform. The method of [20] made essential uses of the Gaussian random design assumption, e.g. Gaussian integration by parts.

Gaussian comparison is another powerful approach for deriving the asymptotic distribution of the Lasso. Building on an earlier idea of [22] that constructs a simpler but comparable Gaussian process, [23] proved asymptotic normality of 6 (in the Wasserstein distance in \(\mathbb{R}^p\)) for i.i.d.\(\mathcal{N}(0,I)\) rows, and [24] extended the result to i.i.d.\(\mathcal{N}(0,\Sigma)\) rows. By nature, the Gaussian comparison argument strongly relies on the Gaussianity of the design matrix.

Characterizing the asymptotic distribution of the Lasso for dependent non-Gaussian designs is an open challenge (see comments in [25] and [24]). Proof of universality based on the Lindeberg-type argument typically assumes independent entries [26] [27].

Approximate message passing (AMP) is not only an algorithm but also a method of characterizing asymptotic distributions. The most common approach for analyzing the state evolution of AMP is through a conditioning technique, which shows that vector approximate message passing works for design matrices with a general spectrum but satisfying right-rotational invariance [[28]][29][[30]][31]. In particular, rotation invariance implies that the feature distribution is permutation invariant, which does not subsume our setting. Other representative approaches for AMP analysis [[32]][33] assume independent matrix entries.

2.4 Resampling-based variable selection methods↩︎

Traditionally, the most well-known variable selection method with guaranteed false discovery rate (FDR) control is the Benjamini-Hochberg procedure [34], which typically assumes that the \(p\)-values are independent or positively correlated. The knockoff filter [[35]][1] is a recent popular approach that controls the FDR without such restrictive dependency assumptions. Intuitively, the knockoff filter creates knockoff features which have the same distribution as the true features, but are conditionally independent of the response, so that the knockoff statistics can be used as a benchmark/control for understanding the \(p\)-values. Remarkably, the knockoff filter extends such an intuition by offering provable finite sample FDR control via an elegant martingale analysis [35].

The knockoff filter performs regression on \(2p\) features and is often observed to have lower statistical power than methods that resample one feature at a time, such as the conditional randomization test (CRT), and Gaussian mirror [36]. CRT is primarily a hypothesis testing procedure. While it can be used to compute \(p\)-values and achieve FDR control when combined with BHq, it is computationally intensive in this setting. The holdout randomization test (HRT; [37]) improves the efficiency of CRT via sample splitting, at the cost of potential power loss. Distilled CRT [7] is a well-known approach for mitigating the computational cost of CRT. Assuming that the \(j\)-th feature is to be tested, we can express their test statistic (see [7]) as \[\begin{align} \frac{|(Y-d_y)^{\top}(B_{:j}-d_x)|}{\|B_{:j}-d_x\|_2^2} \end{align}\] where we follow the notation of [7] to denote by \(d_x\) and \(d_y\) the distilled information. According to [7] we can take \(d_x=\mathbb{E}[B_{:j}|B_{:\setminus j}]\) and \(d_y=A\hat{\alpha}\), where \(A_{:j}=0\) and \(A_{:\setminus j}=B_{:\setminus j}\) (that is, \(A\hat{\alpha}\) is the leave-one-out estimator). Then \(\hat{\alpha}=0\) and \((Y-d_y)^{\top}(B_{:j}-d_x)\) is identical to the right side of 10 . Thus our results can be interpreted as rigorously proving the asymptotic equivalence between the debiased estimator and the distilled CRT statistics in the nontrivial proportional regime, which was not covered in [7] or [3]. Furthermore, we can take \(A\) in 10 to be the original design matrix (rather than the leave-one-out matrix). This obviates the need of solving the leave-one-out Lasso for each \(j\), reducing the computational complexity of the original \(\rm d_0CRT\) in [7] by a factor of \(p\).

3 Main results on approximate formulas↩︎

Notations

We use the standard Landau notations such as \(O(n)\), \(\Omega(n)\), and \(\omega(1)\). The notation \(\tilde{O}(f(n))=O(f(n)(\ln n)^C)\) indicates an upper bound up to a factor of a polynomial of \(\ln n\). To emphasize the dependence on a set of parameters \(\mathcal{P}\) in the implicit constants, we may write \(O_{\mathcal{P}}()\) or \(\tilde{O}_{\mathcal{P}}()\). The \(L_2\) norm and the operator norm are denoted by \(\|\cdot\|_2\) and \(\|\cdot\|_{\rm op}\), whereas \(\|\cdot \|_0\) denotes the number of nonzero coefficients. \(\lambda_{\max}\) and \(\lambda_{\min}\) denote the largest and smallest singular values of a matrix. \(\mathop{\mathrm{col}}()\) denotes the column space of a matrix. The projection matrix onto a set \(\mathcal{A}\) of columns in the design matrix is denoted \(P_{\mathcal{A}}\). The max and min of two numbers \(a,b\) are denoted \(a\vee b\) and \(a\wedge b\). For \(a,b\in\mathbb{R}^p\), write the empirical distribution \(\hat{P}_{ab}:=\frac{1}{p}\sum_{j=1}^p\delta_{(a_j,b_j)}\). We use \(A_{:\setminus j}\) to denote the submatrix of \(A\) formed by all except the \(j\)-th column. The standard normal distribution in \(\mathbb{R}^d\) is written as \(\mathcal{N}(0,I_n)\). \(\mathop{\mathrm{diag}}(\cdot)\) can denote either a diagonal matrix with diagonals specified by a vector, or the vector formed by the diagonal values of a matrix.

3.1 Nonasymptotic error bounds↩︎

Suppose that \(A,B\in\mathbb{R}^{n\times p}\) are matrices differing only in the \(j\)-th column. Recall that we define \(\hat{\alpha}\) and \(\hat{\beta}\) as the solutions to 1 2 (any choice of the minimizer when the minimizer is not unique), and \(\hat{\alpha}^U_j\) and \(\hat{\beta}^U_j\) in 8 and 9 . Let \(\mu_{:j}\in\mathbb{R}^n\) be an arbitrary vector, and define \[\begin{align} \check{B}_{:j}:=B_{:j}-\mu_{:j}. \end{align}\] Our main nonasymptotic result is the following:

The proof the theorem is given in Section 4. Theorem 1 suggests the formula 13 for fast calculation of the debiased estimator when the design matrix is updated by one-column, since the right side only depends on the result of solving 1 . Note that Theorem 1 applies to any given \(B\) as long as conditions 17 19 are satisfied. For deterministic designs, \(\mu_{:j}\) can be taken as any vector that ensures 19 .

A basic example of Theorem 1 is simply \(A=B\), in which case \(\varepsilon=0\) and \(\Gamma=0\). Then 20 simply recovers our definition of the debiased estimator. As a more advanced example, consider random designs where \(A\) and \(B\) are conditionally independent given \(A_{:\setminus j}\). Then, we can take \(D\) as a constant independent of \(n\), \(\Gamma\) as a slowly growing function (e.g.polylog of \(n\)), and \(\varepsilon\) vanishing in \(n\). Moreover, let \[\begin{align} \mu_{:j}:=\mathbb{E}[B_{:j}|B_{:\setminus j}]=\mathbb{E}[B_{:j}|A_{:\setminus j}], \end{align}\] so that 19 is satisfied with high probability when \(\varepsilon\) is small, because \(\mathcal{A}\) is a function of \((A,Y)\) whereas \(\check{B}_{:j}\) is a zero mean vector conditioned on \((A,Y)\). More precise results will be discussed in Section 3.2.

The update formula is closely related to the leave-one-out analysis and the asymptotic normality of debiased Lasso, e.g. [18], [20] and [19]. To see this, observe that by slightly changing the proof of Theorem 1 to allow different observation vector \(Y\) in the two Lasso problems, we obtain:

The proof can be found in Section 4.3. Observe that if \(A\) and \(B\) are independent conditioned on \(A_{:\setminus j}\), 21 suggests that \[\begin{align} \hat{\alpha}_j^U -\alpha_j \approx \frac{\check{A}_{:j}^{\top}S}{\check{A}_{:j}^{\top}(I-P_{\mathcal{B}})A_{:j}} \approx \frac{\check{A}_{:j}^{\top}S}{(n-k)\Sigma_{j|\setminus j}} \label{e95approxa} \end{align}\tag{22}\] which approximately follows \(\mathcal{N}(0,(n-k)^{-2}\Sigma_{j|\setminus j}^{-1}\|S\|_2^2)\) under appropriate regularity conditions. This recovers the asymptotic normality of the debiased Lasso previously proved by other methods; see Section 1.

In contrast to its debiased version, the Lasso estimator has update formula in a more restricted setting:

The proof the theorem is given in Section 4. Theorem 3 suggests the approximate formula \[\begin{align} \hat{\beta}_j \approx \left(\frac{1}{n}B_{:j}^{\top}(I-P_{\mathcal{A}})B_{:j} \right)^{-1} \mathsf{S}_{\lambda} \left( \frac{1}{n}B_{:j}^{\top}R + \frac{1}{n}B_{:j}^{\top}(I-P_{\mathcal{A}})A_{:j} \hat{\alpha}_j \right). \end{align}\] The approximation is good if \(\tau, \varepsilon\to0\). Note that \(\tau\to 0\) implies that the features are approximately independent. In the proof of Theorem 1, error terms of the form \(\frac{1}{n}\check{B}_{:j}^{\top}(P_{\mathcal{A}\cup\Delta}-P_{\mathcal{A}\setminus\Delta})A_{:j}\) asymptotically vanishes since \(\check{B}_{:j}\) and \((P_{\mathcal{A}\cup\Delta}-P_{\mathcal{A}\setminus\Delta})A_{:j}\) are uncorrelated. On the other hand, in the proof of Theorem 3, error terms of the form \(\frac{1}{n}B_{:j}^{\top}(P_{\mathcal{A}\cup\Delta}-P_{\mathcal{A}\setminus\Delta})A_{:j}\) arise, and it is not necessarily vanishing unless the features are independent. In the experiments, we also observe much better approximation of the debiased estimator under correlated features.

3.2 Asymptotic error bounds↩︎

Using Theorem 1, we can show that asymptotically and for most \(j\), the approximation error is negligible, under the following condition:

For each \(j\), we generate \(B^{(j)}_{:j}\) by setting \(B^{(j)}_{:\setminus j}=A_{:\setminus j}\), and independently sampling \(B^{(j)}_{:j}\) according to \(P_{A_{:j}|A_{:\setminus j}}\). For each \(j\), recall \(P_{\mathcal{A}}\) and \(P_{\mathcal{B}}\) defined around 8 and 9 , and set \[\begin{align} P_j:=P_{\mathcal{A}}; \quad P^{(j)}:=P_{\mathcal{B}}. \label{e95pj} \end{align}\tag{24}\] We then have a general asymptotic approximation result:

The proof is given in Section 9. In Theorem 4 we bound the approximation error in computing \(t(j,B^{(j)},Y) =\check{B}^{(j)\top}_{:j} (I-P^{(j)}) B^{(j)}_{:j}\hat{\beta}_j^{(j)U}\) from the update formula. A natural question is, what about the approximation error for \(\hat{\beta}_j^{(j)U}\) itself? That is, by Theorem 4, do we have \[\begin{align} \hat{\beta}_j^{(j)U} \approx \frac{\frac{1}{n}\check{B}_{:j}^{\top}R + \frac{1}{n}\check{B}_{:j}^{\top}(I-P_j)A_{:j} \hat{\alpha}_j}{\frac{1}{n}\check{B}_{:j}^{\top}(I-P_j)B_{:j}}? \label{e95approxj} \end{align}\tag{26}\] As we will see in Theorem 5, for 26 to hold it suffices make an additional assumption of \[\begin{align} \mathbb{E}[|\check{B}^{(j)}_{1j}|^2] \ge \kappa_2^{-1}, \quad \forall j=1,\dots, p, \label{e95lb} \end{align}\tag{27}\] for some \(\kappa_2>0\) independent of \(n\). 27 is true, for example, if \(\mu_{:j}=A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj}\) (i.e., the linear predictor is optimal), under the assumption of \(\mathcal{P}\). In the meantime, this is essentially also necessary: if \(\mathbb{E}[|\check{B}^{(j)}_{:j}|^2]=0\), the right side of 26 is undefined. We have the following answer to 26 :

The proof is given in Section 10. Theorem 4 and Theorem 5 can be used for variable selection under the false discovery rate (FDR) control, because they bound the approximation error for all but a small fraction of \(j\), and a vanishing fraction of coordinates does not change the asymptotic FDR and power of the selection algorithm. To formalize the notion of “approximation in most coordinates”, recall the notion of Levy-Prokhorov metric which metricizes weak convergence of probability measures (see for example [39]):

Corollary 6 follows directly from Theorem 4, Theorem 5 and the definition of the Levy-Prokhorov metric, so the proof is omitted. An analysis of the asymptotic impact of the approximation error on variable selection algorithms is given in Section 5. Finally, we have the following result that the generalized debiased Lasso coefficients reduces to the traditional debiased Lasso coefficients in the case of Gaussian designs, whose proof is deferred to Section 11.

3.3 Beyond \(\ell_1\) regularizer↩︎

In this section, we extend the previous results to the case where the regularization term \(\lambda\|\beta\|_1\) is replaced by \(\sum\rho(\beta_j)\), where \(\rho\) is a general function satisfying: \[\begin{align} \rho(0)&=\psi(0)=0; \\ \psi'&:=\rho''\in [\zeta_{2-},\zeta_{2+}]; \\ |\psi''|&\le \zeta_3, \end{align}\] where \(\zeta_{2-},\zeta_{2+},\zeta_3>0\) are constants. The proofs for such \(\rho\) is simpler than the Lasso case, due to strong convexity. Define \[\begin{align} \bar{G}_{\mathcal{A}} &:=\frac{1}{n}A_{:\setminus j}^{\top}A_{:\setminus j} +\mathop{\mathrm{diag}}(\psi'(\hat{\alpha}_{\setminus j})), \\ P_{\mathcal{A}}&:=\frac{1}{n}A_{:\setminus j} \bar{G}_{\mathcal{A}}^{-1} A_{:\setminus j}^{\top}. \end{align}\] The following theorems establish the nonasymptotic and asymptotic error bounds analogous to the \(\ell_1\) case, whose proofs can be found in the appendix.

4 Proof of the approximation formula↩︎

The goal of this section is to prove Theorem 1, Theorem 2, and Theorem 3.

4.1 Intuitions↩︎

Recall the optimization problems given in 1 2 . Intuitively an update formula is possible because the Taylor expansion is asymptotically correct; the challenge though lies in the non-differentiability of the objective function and in showing that error is indeed negligible in the high-dimensional setting. Before the proof let us first give a heuristic derivation of the approximate update formula. Let \(R:=Y-A\hat{\alpha}\) and \(S:=Y-B\hat{\beta}\) denote the residuals. From the normal equations we have \[\begin{align} -\frac{1}{n}A_{:\setminus j}^{\top}R +\lambda\psi^{\alpha}_{\setminus j}&=0_{p-1}, \tag{33} \\ -\frac{1}{n}A_{:\setminus j}^{\top}S +\lambda\psi^{\beta}_{\setminus j}&=0_{p-1}, \tag{34} \\ -\frac{1}{n}B_{:j}^{\top}S +\lambda\psi^{\beta}_j&=0, \tag{35} \end{align}\] where the subdifferential \(\psi^{\alpha}\in[-1,1]^p\) is intuitively the derivative of (the non-differentiable function) \(\|\cdot\|_1\) at \(\hat{\alpha}\). Because of the non-differentiability, \(\psi^{\alpha}\) is not a function \(\hat{\alpha}\). But as a heuristic argument, we pretend that \(\psi^{\alpha}=\psi(\hat{\alpha})\), where \(\psi\colon \mathbb{R}^p\to\mathbb{R}^p\) applies element-wise the function of the derivative of the absolute value function. We also pretend that \(\psi(\cdot)\) is a smooth function so that we can Taylor expand \(\psi(\cdot)\) around \(\hat{\alpha}\). Now using \[\begin{align} S-R=A_{:\setminus j}(\hat{\alpha}_{\setminus j}-\hat{\beta}_{\setminus j})+A_{:j}\hat{\alpha}_j-B_{:j}\hat{\beta}_j, \label{eqn3} \end{align}\tag{36}\] and ignoring higher order terms in the Taylor expansion, we obtain \[\begin{align} \bar{G}(\hat{\alpha}_{\setminus j}-\hat{\beta}_{\setminus j}) +u\hat{\alpha}_j -v\hat{\beta}_j&\approx 0_{p-1}, \tag{37} \\ -\frac{1}{n}B_{:j}^{\top}R -v^{\top}(\hat{\alpha}_{\setminus j}-\hat{\beta}_{\setminus j}) -\frac{1}{n}B_{:j}^{\top}A_{:j}\hat{\alpha}_j +\frac{1}{n}B_{:j}^{\top}B_{:j}\hat{\beta}_j+\lambda \psi(\hat{\beta}_j)&\approx 0, \tag{38} \end{align}\] where we defined \[\begin{align} \bar{G}&:=\frac{1}{n}A_{:\setminus j}^{\top}A_{:\setminus j} +\lambda\mathop{\mathrm{diag}}(\psi'(\hat{\alpha}_{\setminus j})); \label{e95barg} \\ u&:=\frac{1}{n}A_{:\setminus j}^{\top}A_{:j}; \quad v:=\frac{1}{n}A_{:\setminus j}^{\top}B_{:j}. \end{align}\tag{39}\] Cancelling \(\hat{\alpha}_{\setminus}-\hat{\beta}_{\setminus j}\), we can solve for \(\hat{\beta}_j\) to obtain \[\begin{align} (\frac{1}{n}\|B_{:j}\|^2-v^{\top}\bar{G}^{-1}v)\hat{\beta}_j +\lambda\psi(\hat{\beta}_j) \approx(\frac{1}{n}B_{:j}^{\top}A_{:j} -v^{\top}\bar{G}^{-1}u)\hat{\alpha}_j +\frac{1}{n}B_{:j}^{\top}R. \label{eq59} \end{align}\tag{40}\] Using the matrix inversion formula, we see that \(\bar{G}^{-1}\) is nonzero only in the principal submatrix corresponding to the nonzeros of \(\hat{\alpha}_{\setminus j}\). In the high dimension setting, supposing that the entries of \(A\) are i.i.d.with unit variance, we have \((1-\frac{k}{n})\hat{\beta}_j +\lambda\psi(\hat{\beta}_j) \approx \frac{1}{n}B_{:j}^{\top}R\), therefore \((1-\frac{k}{n})\hat{\beta}_j^{u} \approx \frac{1}{n}B_{:j}^{\top}R\), which recovers the known formula of the debiased estimator 6 for i.i.d.features.

4.2 Proof of Theorem 1↩︎

To deal with the non-differentiability of \(\psi\) in 39 , define \[\begin{align} \Lambda=\mathop{\mathrm{diag}}\left(\left\{\frac{\psi^{\alpha}_l-\psi^{\beta}_l}{\hat{\alpha}_l-\hat{\beta}_l}\right\}_{l\neq j}\right) \label{eq95d} \end{align}\tag{41}\] where the \(\frac{0}{0}\) case in 41 is resolved by setting \(\frac{\psi^{\alpha}_l-\psi^{\beta}_l}{\hat{\alpha}_l-\hat{\beta}_l}=0\) if \(\psi^{\alpha}_l=\psi^{\beta}_l\notin (-1,1)\), and \(\frac{\psi^{\alpha}_l-\psi^{\beta}_l}{\hat{\alpha}_l-\hat{\beta}_l}=+\infty\) if \(\psi^{\alpha}_l=\psi^{\beta}_l\in (-1,1)\). Then set \[\begin{align} G:=\frac{1}{n}A_{:\setminus j}^{\top}A_{:\setminus j} +\lambda \Lambda. \label{e95g} \end{align}\tag{42}\] Note that by replacing \(\bar{G}\) with \(G\) in 40 , exact equality is achieved, and we have \[\begin{align} (\frac{1}{n}\|B_{:j}\|^2-v^{\top}G^{-1}v)\hat{\beta}_j +\lambda\psi^{\beta}_j = \frac{1}{n}B_{:j}^{\top}R + (\frac{1}{n}B_{:j}^{\top}A_{:j} -v^{\top}G^{-1}u)\hat{\alpha}_j. \label{eq9569} \end{align}\tag{43}\] Note that \[\begin{align} \frac{1}{n}B_{:j}^{\top}R&=\frac{1}{n}\check{B}_{:j}^{\top}R +\frac{1}{n}\mu_{:j}^{\top}S +\frac{1}{n}\mu_{:j}^{\top}(R-S), \label{eq66} \end{align}\tag{44}\] and next we will simplify the last term in 44 using 36 and Proposition 16: \[\begin{align} \frac{1}{n}\mu_{:j}^{\top}(R-S) &= -\frac{1}{n}\mu_{:j}^{\top} [A_{:\setminus j}(\hat{\alpha}_{\setminus j}-\hat{\beta}_{\setminus j}) +A_{:j}\hat{\alpha}_j -B_{:j}\hat{\beta}_j] \\ &=-\frac{1}{n}\mu_{:j}^{\top}A_{:\setminus j}G^{-1}v\hat{\beta}_j + \frac{1}{n}\mu_{:j}^{\top}A_{:\setminus j}G^{-1}u\hat{\alpha}_j - \frac{1}{n}\mu_{:j}^{\top}A_{:j}\hat{\alpha}_j + \frac{1}{n}\mu_{:j}^{\top}B_{:j}\hat{\beta}_j. \end{align}\] Note that even if \(G\) is not invertible, we can show that \(A_{:\setminus j}G^{-1}A_{:\setminus j}^{\top}\) is well-defined via limit. Collecting terms and using 7 and \(\lambda\psi^{\beta}_j -\frac{1}{n}\mu_{j:}^{\top}S=\frac{1}{n}\check{B}_{:j}^{\top}S\), we see 43 becomes: \[\begin{align} &\quad(\frac{1}{n}\|B_{:j}\|^2 -\frac{1}{n}\mu_{:j}^{\top}B_{:j} -\frac{1}{n}\check{B}_{:j}^{\top}A_{:\setminus j}G^{-1}v)\hat{\beta}_j +\frac{1}{n}\check{B}_{:j}^{\top}S \nonumber \\ &=\frac{1}{n}\check{B}_{:j}^{\top}R + (\frac{1}{n}\check{B}_{:j}^{\top}A_{:j} - \frac{1}{n}\check{B}_{:j}^{\top}A_{:\setminus j}G^{-1}u)\hat{\alpha}_j. \label{eq69} \end{align}\tag{45}\] Next we estimate the coefficients on the two sides of 45 . Define \[\begin{align} \bar{G}_{\mathcal{A}}&:=\frac{1}{n}A_{:\setminus j}^{\top}A_{:\setminus j} +\lambda\bar{\Lambda}, \label{e95ga} \end{align}\tag{46}\] where \(\bar{\Lambda}\) is a diagonal matrix, \(\bar{\Lambda}_{ll}:=0\) if \(l\in\mathcal{A}\) and \(\bar{\Lambda}_{ll}:=+\infty\) otherwise. Define \(\bar{G}_{\mathcal{B}}\) analogously but with \(\mathcal{A}\) above replaced by \(\mathcal{B}:=\{l\neq j\colon \chi^{\beta}_l\neq 0\}\). Then the third term in the coefficient of \(\hat{\beta}_j\) is \[\begin{align} \frac{1}{n}\check{B}_{:j}^{\top}A_{:\setminus j}G^{-1}v &=\frac{1}{n^2}\check{B}_{:j}^{\top}A_{:\setminus j}\bar{G}_{\mathcal{B}}^{-1}A_{:\setminus j}^{\top}B_{:j} + \frac{1}{n^2}\check{B}_{:j}^{\top}A_{:\setminus j}(G^{-1}-\bar{G}_{\mathcal{B}}^{-1})A_{:\setminus j}^{\top}B_{:j} \label{eq70} \\ &=\frac{1}{n}\check{B}_{:j}^{\top}P_{\mathcal{B}}B_{:j} + \frac{1}{n^2}\check{B}_{:j}^{\top}A_{:\setminus j}(G^{-1}-\bar{G}_{\mathcal{B}}^{-1})A_{:\setminus j}^{\top}B_{:j} \end{align}\tag{47}\] where \(P_{\mathcal{B}}:=\frac{1}{n}A_{:\setminus j}\bar{G}_{\mathcal{B}}^{-1}A_{:\setminus j}^{\top}\) denotes the projection onto the span of the columns corresponding to the indices \(\mathcal{B}\). Therefore the coefficient of \(\hat{\beta}_j\) on the left side of 45 differs from \(\frac{1}{n}\check{B}_{:j}^{\top}(I-P_{\mathcal{B}})B_{:j}\) by at most \[\begin{align} \frac{1}{n^2}|\check{B}_{:j}^{\top}A_{:\setminus j}(G^{-1}-\bar{G}_{\mathcal{B}}^{-1})A_{:\setminus j}^{\top}B_{:j}| &\le \frac{1}{n^2}\left\| \check{B}_{:j}^{\top}A_{:\setminus j}(G^{-1}-\bar{G}_{\mathcal{B}}^{-1})A_{:\setminus j}^{\top} \right\|_2 \left\|B_{:j}\right\|_2 \\ &\le \frac{1}{n} \left\| \check{B}_{:j}^{\top}(P_{\mathcal{A}\cup\mathcal{J}}-P_{\mathcal{A}\setminus\mathcal{J}}) \right\|_2\left\|B_{:j}\right\|_2 \label{e52} \\ &\le \frac{1}{n}\cdot \Gamma \sqrt{|\mathcal{J}|}\cdot D\sqrt{n}, \end{align}\tag{48}\] where 48 follows from Lemma 15. Similarly, for one term in the coefficient for \(\hat{\alpha}_j\) in 45 , \[\begin{align} \left|\frac{1}{n}\check{B}_{:j}^{\top}A_{:\setminus j}G^{-1}u -\frac{1}{n}\check{B}_{:j}^{\top}P_{\mathcal{A}}A_{:j} \right| &= \frac{1}{n^2}\left|\check{B}_{:j}^{\top}A_{:\setminus j}(G^{-1}-\bar{G}_{\mathcal{A}}^{-1})A_{:\setminus j}^{\top}A_{:j} \right| \\ &\le \frac{1}{n^2}\left\|\check{B}_{:j}^{\top}A_{:\setminus j}(G^{-1}-\bar{G}_{\mathcal{A}}^{-1})A_{:\setminus j}^{\top}\right\|_2 \left\|A_{:j} \right\|_2 \\ &\le \Gamma D\sqrt{\varepsilon}. \label{eq57} \end{align}\tag{49}\] Then 45 yields \[\begin{align} &\quad\left|\frac{1}{n}\check{B}_{:j}^{\top}(I-P_{\mathcal{B}})B_{:j}\hat{\beta}_j^U-\frac{1}{n}\check{B}_{:j}^{\top}R - \frac{1}{n}\check{B}_{:j}^{\top}(I-P_{\mathcal{A}})A_{:j} \hat{\alpha}_j \right| \le \Gamma D\sqrt{\varepsilon} (|\hat{\alpha}_j|+|\hat{\beta}_j|). \end{align}\] Finally the proof is completed by \[\begin{align} \frac{1}{n}|\check{B}_{:j}^{\top}(P_{\mathcal{A}}-P_{\mathcal{B}})B_{:j}|&\le \frac{D}{\sqrt{n}}\|\check{B}_{:j}^{\top}(P_{\mathcal{A}}-P_{\mathcal{B}})\|_2 \\ &\le D\Gamma \sqrt{\varepsilon}. \end{align}\]

4.3 Proof of Theorem 2↩︎

Similar to Theorem 1, but make changes to 36 , 37 , 38 , 40 , and Section 8.1 by replacing \(\hat{\beta}_j\) with \(\hat{\beta}_j-\alpha_j\), and switch the notations \(A\leftrightarrow B\), \(\hat{\alpha}\leftrightarrow \hat{\beta}\).

4.4 Proof of Theorem 3↩︎

Recall that we showed 43 . Now \[\begin{align} v^{\top}G^{-1}v-\frac{1}{n}B_{:j}^{\top}P_{\mathcal{A}}B_{:j} &= \frac{1}{n}B_{:j}^{\top}EB_{:j} \\ &=\frac{1}{n}\mu_{:j}^{\top}E\mu_{:j} +\frac{2}{n}\mu_{:j}^{\top}E\check{B}_{:j} +\frac{1}{n}\check{B}_{:j}^{\top}E\check{B}_{:j}, \end{align}\] where \(E:=\frac{1}{n}A_{:\setminus j}GA_{:\setminus j}^{\top}-P_{\mathcal{A}}\), and recall that \(\mu_{:j}=B_{:j}-\check{B}_{:j}\). As shown in Lemma 15, \[\begin{align} -(P_{\mathcal{A}\cup \mathcal{J}}-P_{\mathcal{A}\setminus \mathcal{J}}) \preceq E,E^2 \preceq P_{\mathcal{A}\cup \mathcal{J}}-P_{\mathcal{A}\setminus \mathcal{J}} \end{align}\] where we note that \(P_{\mathcal{A}\cup \mathcal{J}}-P_{\mathcal{A}\setminus \mathcal{J}}\) is a projection matrix. Therefore, \[\begin{align} \left|v^{\top}G^{-1}v-\frac{1}{n}B_{:j}^{\top}P_{\mathcal{A}}B_{:j}\right| &\le \frac{1}{n}\|\mu_{:j}\|_2^2 + \frac{2}{n}\|\mu_{:j}\|_2\|E\check{B}_{:j}\|_2 + \frac{1}{n}\check{B}_{:j}^{\top}(P_{\mathcal{A}\cup \mathcal{J}}-P_{\mathcal{A}\setminus \mathcal{J}})\check{B}_{:j}, \tag{50} \\ &\le \tau^2+2\tau\Gamma\sqrt{\varepsilon}+\varepsilon\Gamma^2. \tag{51} \end{align}\] Similarly, we have \[\begin{align} \left|v^{\top}G^{-1}u-\frac{1}{n}B_{:j}^{\top}P_{\mathcal{A}}A_{:j}\right| &\le \frac{1}{n}\|\mu_{:j}\|_2^2 + \frac{1}{n}\|\mu_{:j}\|_2 (\|E\check{B}_{:j}\|_2 +\|E\check{A}_{:j}\|_2) +\frac{1}{n}\|E\check{B}_{:j}\|_2\|\check{A}_{:j}\|_2 \\ &\le \tau^2+\tau(\Gamma\sqrt{\varepsilon}+\check{D})+\Gamma\sqrt{\varepsilon}\check{D}. \end{align}\] Then the claim follows from 43 .

5 Application in false discovery rate control↩︎

In this section, we provide a theoretical analysis showing that, in certain settings, the power of the knockoff filter can be strictly improved by the local knockoff filter and the conditional randomization test. Although these methods have been studied empirically in the literature and are typically at least \(p\) times slower than knockoffs, we show that this \(p\)-fold computational overhead can be eliminated by using the approximation formula to accelerate the computation.

5.1 Review of the knockoff filter and its limitation↩︎

Consider the variable selection problem with observations \(Y\in\mathbb{R}^n\) and \(A\in\mathbb{R}^{n\times p}\) The false discovery rate (FDR) is defined as \[\begin{align} {\rm FDR}:=\mathbb{E}\left[ FDP\right], \end{align}\] where \(FDP:=\frac{|H_0\wedge \hat{H}_1|}{|\hat{H}_1|}\), \(H_0\subseteq \{1,2,\dots,p\}\) is the set of true null variables, and \(\hat{H}_1\subseteq \{1,2,\dots,p\}\) is the set of selected variables. A good variable selection algorithm is expected to control the FDR below a given budget, while ensuring a large statistical power: \[\begin{align} {\rm POWER}:=\mathbb{E}\left[ \frac{|H_1\wedge \hat{H}_1|}{|H_1|} \right]. \end{align}\]

The present paper focuses on the Model-X knockoff framework [1], which aligns with our emphasis on resampling-based procedures. In this framework, a knockoff matrix \(\tilde{A}\) is constructed so that \((A, \tilde{A})\) is exchangeable, while \(\tilde{A}\) is conditionally independent of \(Y\). Then we regress \(Y\) on the matrix \([A,\tilde{A}]\), so that the test statistics for \(\tilde{A}\) can be used to estimate the number of false discoveries. For a full description, see [1].

The exchangeability condition may create challenges in the construction of the knockoff distribution, and the increase in the number of features (from \(p\) to \(2p\)) can often induce a power loss [[40]][6]. In this section, we provide a simple example which always suffers from this suboptimality, regardless of the choice of the knockoff distribution (Theorem 13). For simplicity, consider the \(2p<n\) regime, and we simply use the least squares estimator rather than the Lasso. We recall the following result from [41] which gives a necessary and sufficient condition for asymptotic consistency of the knockoff filter in this regime:

The setting of [41] is to use the debiased Lasso coefficients as test statistics, but when the Lasso regularization \(\lambda=0\) we recover the case of least squares statistics. Moreover, [41] assumes the existence of the standard distributional limit defined in [2], so the empirical distribution of \((\underline{\Theta}_{jj}^{(n)})_{j=1}^{2p}\) has a weak limit, and by “\(\lim_{n\to \infty}{\rm POWER}^{(n)}\) to converge to 1” in Proposition 10 we mean \(\lim_{n\to \infty}{\rm POWER}^{(n)}\) is viewed as a function of this limiting distribution.

The intuition for Proposition 10 is as follows: Let \(\underline{\alpha}\in\mathbb{R}^{2p}\) be the true coefficients and zero paddings for the knockoff coefficients. There exists \(\tau>0\), bounded above and below, such that \[\begin{align} \hat{\underline{\alpha}}\approx \underline{\alpha}+\tau\underline{\Theta}^{1/2}z \end{align}\] where \(z\sim \mathcal{N}(0,I_{2p})\). Therefore if variables are selected based on a threshold test for the coefficients of \(\hat{\underline{\alpha}}\), then asymptotic consistency is true only if most diagonal entries of \(\underline{\Theta}\) vanish (equivalently, the empirical distribution of the diagonal entries must converge to zero).

For Gaussian knockoff filters, it is known [1] that selecting a knockoff distribution satisfying exchangeability is equivalent to choosing \(s\) such that the joint covariance matrix of the true and knockoff variables, \[\begin{align} \underline{\Sigma}:= \begin{pmatrix} \Sigma & \Sigma-\mathop{\mathrm{diag}}(s) \\ \Sigma-\mathop{\mathrm{diag}}(s) & \Sigma \end{pmatrix} \label{e95sigmau} \end{align}\tag{52}\] is positive semidefinite. Thus we see that the Schur complement \(2S-S\Sigma^{-1}S\) plays a key role. We now establish an auxiliary result that will help building a suboptimality example.

Proof can be found in Section 15. Note that \(\Sigma^{-1}\) in the lemma is not invertible, but we do not need its inverse; alternatively, we may perturb it to make it invertible, and then pass the final result to a limit.

The oracle threshold algorithm is not realistic since the threshold is not data-driven. Nevertheless, it serves as a natural benchmark, and has been considered in, e.g.[42], under the name prototype method. We then have the following result, whose proof is omitted since it is analogous to Proposition 10.

By Schur’s complement theorem, it is easy to see that the diagonals of \(\underline{\Theta}\) dominate the corresponding diagonals of \(\Theta\), so the knockoff algorithm is asymptotically consistent only if the oracle threshold algorithm is. Using Lemma 11, we can construct an example in the following theorem exhibiting a strict separation: the knockoff filter is not asymptotically consistent, no matter how \(s\) is chosen in 52 , even though the oracle threshold algorithm remains asymptotically consistent.

Proof can be found in Appendix 15.

5.2 Local knockoff filter↩︎

As mentioned in Section 1, a variant of the knockoff method that resamples only a single feature in each Lasso run, referred to as the local knockoff filter (or “one-at-a-time knockoffs” in the concurrent work [5]), often achieves higher statistical power. In fact, its asymptotic power is approximately that of the oracle threshold algorithm in Definition 3. However, this improvement comes at the cost of increased computational complexity, namely \(\Theta(pL)\), where \(L\) denotes the cost of solving a single regression problem. In this section, we show how the approximation formula in Section 3.1 can be used to reduce this complexity by a factor of \(p\).

In the local knockoff method (See Algorithm 1), the exchangeability of \((A,\tilde{A})\) in the knockoff filter is weakened to the conditional exchangeability of columns \((A_{:j},\tilde{A}_{:j})|A_{:\setminus j}\), for each \(j=1,\dots,p\). Then for each \(j\), one regression \(Y\) on \((A_{:j},A_{:,\setminus j})\) and on \((\tilde{A}_{:j},A_{:,\setminus j})\), thus generating a pair of test statistics which have the same distribution if \(j\in H_0\). To see why the algorithm controls the FDR, note that \(|\{j\colon |\hat{\alpha}_j^U|>t\}|\) is the number of selected variables, whereas \(|\{j\colon |\hat{\gamma}_j|>t\}|\ge |\{j\in H_0\colon |\hat{\gamma}_j|>t\}|\) approximately controls the number of false discoveries.

The two debiased estimators \(\hat{\alpha}^u\) and \(\hat{\alpha}^U\) in 6 and 8 have similar asymptotic behaviors under sufficient regularity conditions. We shall use \(\hat{\alpha}^U\) since it arises more naturally in the derivations of the update formula and also appears to induce smaller error in numerical experiments. To analyze the computational complexity, note that

• Although the definition of \(\hat{\alpha}^U\) may appear to require \(O(p^4)\) computation time since each \(P_j\) requires \(O(p^3)\) computation time directly from the formula of the projection matrix, we can actually compute each \(P_j\) in \(O(p^2)\) time by the rank one update formula and hence \(\hat{\alpha}^U\) in \(O(p^3)\) time (see Section 15.4). Thus, the computation time for \(\hat{\alpha}^U\) is the same order as obtaining \(\hat{\gamma}\) using the approximate update formula.

• In the preprocessing step, the conditional means \([\mu_{:1},\dots,\mu_{:p}]\) can be calculated in \(O(p^3)\) time, if it can be approximated using the linear estimator \(-A[\Theta-\mathop{\mathrm{diag}}\mathop{\mathrm{diag}}\Theta]\cdot (\mathop{\mathrm{diag}}\mathop{\mathrm{diag}}\Theta)^{-2}\), where \(\Theta\) denotes the precision matrix. Otherwise, we follow the practice in [7] by regressing each \(A_{:j}\) on \(A_{:\setminus j}\).

In summary, the computational complexity of Algorithm 1 is \(O(p^3 + L)\), which becomes \(O(p^3)\) when \(L = O(p^3)\) (e.g., when solving the Lasso via least angle regression [4]). This matches the computational complexity of the knockoff filter.

We now show that the local knockoff filter guarantees FDR control under certain asymptotic assumptions. In order for the proof to proceed smoothly, we introduce \[\begin{align} T(\epsilon):=\inf\left\{t\colon \frac{|\{j\colon |\hat{\gamma}_j|>t-\epsilon\}|}{|\{j\colon |\hat{\alpha}_j^U|>t\}|} \le q \right\}. \label{e95te} \end{align}\tag{53}\] for any given \(\epsilon\in(0,1)\). Note that \(T(0)=T\) in Algorithm 1, and the numerator in 53 is an overestimate of the number of false discoveries when \(\epsilon>0\). In practical applications though, it may not be necessary to use \(\epsilon>0\), if the empirical measure of \((\hat{\gamma}_j)\) is not concentrated around one point (consider for example, the setting where the Gaussian limit property [2] is true, so that the empirical measure is a Gaussian-smoothed density).

The following consequence of Algorithm 1 is rather direct:

Proof can be found in Section 15.

5.3 Fast conditional randomization test↩︎

As discussed in Section 1, the complexity of exactly solving CRT is \(O(KpL)\), where \(L\) denotes the complexity of solving lasso. In this section, we adopt the approximate update formula 13 to compute \(\hat{\gamma}^b_j\), which we call approx-CRT-db; see Algorithm 2. Following the same argument as in Theorem 14, approx-CRT-db still guarantees asymptotic FDR control. Note that computing \(Y-A\hat{\alpha}\) has complexity \(O(L)\), and then computing the product \(\check{B}_{:j}^{\top}(Y-A\hat{\alpha})\) takes time \(O(n)=O(p)\) in the proportion regime. This implies that running CRT approximately takes only \(O(L+p^2K)\) time. As discussed in Section 2.4, approx-CRT-db is asymptotically equivalent to a version of distilled CRT in [7], but exactly implementing the latter takes \(O(Lp+p^2K)\) times.

6 Conclusion and future work↩︎

For a given design matrix and observation, we derived a formula for updating a debiased Lasso coefficient under perturbation in one column, and upper bounded the nonasymptotic error (Theorem 1). Then by general concentration and anti-concentration machinery, vanishing approximation error in all but a vanishing fraction of coordinates was established under mild assumptions (Theorem 4). The heart of the argument lies in establishing a stability result for the signs of the Lasso coefficients, which implies that a heuristic Taylor expansion argument can be used for updating the debiasing Lasso coefficients. Interestingly, updating the Lasso coefficients (without debiasing) requires much stronger assumptions (Theorem 3). As a consequence, several FDR control algorithms based on feature resampling can be implemented much faster with asymptotically zero impact, if debiased Lasso coefficients are used as the test statistics.

We expect that the asymptotic error control in Theorem 4 can be extended to even more general matrix classes. In the variable selection literature, a common assumption for the design matrix is the factor model [44]. We expect that Theorem 4 can be extended to the setting of \(A=QU\), where \(Q\in\mathbb{R}^{n\times r}\) is a deterministic matrix with bounded singular values, and \(U\) is a random matrix from the class of Definition 1. Furthermore, the proof technique for Theorem 4 should be extendable beyond \(\ell_1\) penalties.

More broadly, we expect that some of our stability type analysis can be adapted to other related problems. For example, algorithmic stability [45] is defined as the stability of a function of the training data when one data point is removed, which implies desirable generalization properties and predictive inference guarantees [46] [47]. Our approximation error bounds for debiased Lasso based on the number of sign changes may be adapted to establish algorithmic stability or differential privacy guarantees. Furthermore, while the original definition of the debiased coefficients may be non-differentiable in \((A,Y)\), its approximation formula is piecewise differentiable. This suggests the possibility of establishing concentration or Gaussian limit results for the empirical distribution of the debiased Lasso through the leave-one-out type analysis leveraging the Poincare inequality or Stein’s method, under conditions more general than (or at least not covered by) existing approaches based on vector approximate message passing [[48]][49] or Lindeberg’s universality argument [26].

Finally, the recent work of [32] also showcased the power and generality of the leave-one-out approach, by applying it in deriving nonasymptotic error bounds for the AMP state evolution, which differs from the previous proofs using the conditioning technique that require rotational invariance. While [32] focused on the case of independent matrix entries, it may be that an extension to the vector sub-Gaussian case is possible by combining techniques in our paper.

7 Errors in the projection matrices↩︎

Recall \(G\), \(\bar{G}_{\mathcal{A}}\) and \(\bar{G}_{\mathcal{B}}\) defined in Section 4. Their inverses may not be well-defined as matrices since the diagonal values may be \(+\infty\). However, the inverses can be defined as linear operators on the column space of \(A_{:\setminus j}\). That is, we observe that the following map \[\begin{align} a \mapsto A_{:\setminus j}(A_{:\setminus j}^{\top} A_{:\setminus j} +\mathop{\mathrm{diag}}(a))^{-1}A_{:\setminus j}^{\top} \label{e95gmap} \end{align}\tag{54}\] is well-defined from \([0,+\infty]^{p-1}\) to the set of positive semidefinite matrices, since the minimum eigenvalue of \(A_{:\setminus j}^{\top} A_{:\setminus j} +\mathop{\mathrm{diag}}(a)\) as an operator on the column space of \(A_{:\setminus j}\) is positive and hence admits an inverse. Details of this continuity proof can be found in Lemma 29. Since our final bounds will only depend on \(G^{-1}\) through \(A_{:\setminus j}G^{-1}A_{:\setminus j}^{\top}\) (and similarly for \(\bar{G}_{\mathcal{A}}\) and \(\bar{G}_{\mathcal{B}}\)), we can write \(G^{-1}\) with the understanding that it is well-defined restricted to the column space of \(A_{:\setminus j}\). We now show that \(\frac{1}{n}A_{:\setminus j}G^{-1}A_{:\setminus j}^{\top}\) is close to the projection matrices \(P_{\mathcal{A}}:=\frac{1}{n}A_{:\setminus j}\bar{G}_{\mathcal{A}}^{-1}A_{:\setminus j}^{\top}\) and \(P_{\mathcal{B}}:=\frac{1}{n}A_{:\setminus j}\bar{G}_{\mathcal{B}}^{-1}A_{:\setminus j}^{\top}\).

8 Sign stability↩︎

8.1 Control of change of residual↩︎

The residual change \(R-S\) is easier to control than \(\hat{\alpha}-\hat{\beta}\). From the normal equations in Section 4, we have \[\begin{align} \frac{1}{n}A_{:\setminus j}^{\top}(S-R)+\lambda \Lambda(\hat{\alpha}_{\setminus j}-\hat{\beta}_{\setminus j})=0. \label{e110} \end{align}\tag{58}\] If \(0\cdot\infty\) is encountered in 58 , we can add perturbations as in [prop95err], and pass to a limit in the end, so the end result still holds. Solving \(\hat{\alpha}_{\setminus j}-\hat{\beta}_{\setminus j}\) and plugging it in the normal equations, we have \[\begin{align} S-R=-\frac{1}{\lambda n}A_{:\setminus j}\Lambda^{-1}A_{:\setminus j}^{\top}(S-R) +A_{:j}\hat{\alpha}_j-B_{:j}\hat{\beta}_j. \label{eq85} \end{align}\tag{59}\] Therefore \[\begin{align} S-R=\left[I+\frac{1}{\lambda n}A_{:\setminus j}\Lambda^{-1}A_{:\setminus j}^{\top}\right]^{-1}(A_{:j}\hat{\alpha}_j-B_{:j}\hat{\beta}_j). \label{eq86} \end{align}\tag{60}\] Since the eigenvalues of \(\left[I+\frac{1}{\lambda n}A_{:\setminus j}\Lambda^{-1}A_{:\setminus j}^{\top}\right]^{-1}\) are upper bounded by 1, we have \[\begin{align} \|S-R\|_2\le \left\|A_{:j}\hat{\alpha}_j-B_{:j}\hat{\beta}_j\right\|_2. \label{eq82} \end{align}\tag{61}\]

8.2 Control of \(\psi^{\alpha}\) near the edge↩︎

In this section we prove sufficient conditions under which \(\psi_m(Y)\in (-1,-1+\epsilon]\cup[1-\epsilon,1)\) with small probability whenever \(\epsilon>0\) is small. In this section, we assume that \(A\) is an arbitrary fixed matrix, \(Y\) is a Gaussian random vector with a certain mean, and define \(\psi()\) and \(\chi()\) as the subgradient and the signs of the Lasso solution as functions of \(Y\), using the definitions in Section 1.

8.3 Control of change of signs↩︎

9 The case of sub-Gaussian designs↩︎

In this section we prove Theorem 4.

9.1 Anticoncentration↩︎

In this section we prove auxiliary results that will be used in justifying the conditions in Corollary 19. The following Lemma is a consequence of [38], which lower bounds the minimum singular values of a random matrix with i.i.d.rows, under the assumption that the projection of a row in each direction has bounded (both from the above and below) \(L_1\) and \(L_2\) norms.

9.2 Expected sign changes↩︎

Recall the sign vector \(\chi^{\alpha}\) defined in Section 1, and define \(\chi^{\beta^{(j)}}\) analogously. We show that when the Lasso regularization parameter \(\lambda\) is sufficiently large (but independent of \(n\)), the expected number of sign changes is \(o(n)\).

9.3 Conditioning sub-Gaussians↩︎

We establish auxiliary results for controlling \(\check{B}_{:j}^{(j)}\) conditioned on \((A,Y)\).

9.4 Concentration bounds↩︎

We need the following auxiliary result to control the projection of \(\check{B}_j^{(j)}\) onto small subspaces.

9.5 Proof of Theorem 4↩︎

Define \[\begin{align} \mathcal{I}_3(j,l) &= 1\{\chi^{\alpha}_l\neq \chi^{\beta^{(j)}}_l\}; \\ \mathcal{I}_4(j) &:=1\left\{\frac{1}{n}\sum_l \mathcal{I}_3(j,l) >\varepsilon \right\}. \end{align}\] It follows from Corollary 22 that \[\begin{align} \frac{1}{p}\sum_{j=1,l}^p\mathbb{E}[\mathcal{I}_3(j,l)] = \tilde{O}_{\mathcal{P}}(n^{2/3}). \end{align}\] Then (writing \(\bar{\mathcal{I}}:=1-\mathcal{I}\)), \[\begin{align} \mathbb{P} \left[ \frac{1}{p}\sum_j \mathcal{I}_4(j) >\delta_1 \right] &\le \frac{1}{p\delta_1} \sum_j \mathbb{P}\left[\frac{1}{n}\sum_l \mathcal{I}_3(j,l) >\varepsilon \right] \\ &\le \frac{1}{\delta_1\varepsilon pn} \sum_{j,l} \mathbb{E}[\mathcal{I}_3(j,l)] \\ &\le \tilde{O}_{\mathcal{P}}(p^{-1/3}\delta_1^{-1}\varepsilon^{-1}) \label{eqq118} \end{align}\tag{88}\] where \(\varepsilon\) and \(\delta_1\) will be optimized later (see 93 ). Let \(c\) be as in Corollary 22, and let \[\begin{align} \mathcal{I}_8 := 1\{ \|A\|_{\rm op}> \sqrt{n}/c \textrm{ or \eqref{e129} fails for some j}\}. \end{align}\] Applying Theorem 1 with \(D=\frac{1}{c}\) and \(\Gamma=10\sqrt{c\ln n}\ln \frac{1}{\varepsilon}\), we see that under the event \(\mathcal{I}_4(j)=\mathcal{I}_8=0\), \[\begin{align} |t(j,B^{(j)},Y)-t_j| \le 10\sqrt{\frac{\ln n}{c}} \ln \frac{1}{\varepsilon} \cdot \sqrt{\varepsilon}(|\hat{\alpha}_j|+|\hat{\beta}^{(j)}_j|). \label{e136} \end{align}\tag{89}\]

To finish the proof, we will control the right side of 89 . Let \[\begin{align} \mathcal{I}_6(j) :=1\{\max\{ |\hat{\alpha}_j|, |\hat{\beta}^{(j)}_j| \}>\delta_3^{-1/2}\}, \end{align}\] and next we will show that \(\mathcal{I}_6(j)=0\) for most \(j\). Using 83 , 84 and the Markov inequality, we have that for any \(\delta_3>0\), \[\begin{align} \max\left\{\frac{1}{p}\sum_{j=1}^p \mathbb{P}[|\hat{\alpha}_j|>\delta_3^{-1/2}] ,\, \frac{1}{p}\sum_{j=1}^p \mathbb{P}[|\hat{\beta}^{(j)}_j|>\delta_3^{-1/2}] \right\} =\tilde{O}_{\mathcal{P}}(\delta_3) \end{align}\] and \[\begin{align} \mathbb{P}[\frac{1}{p}\sum_j\mathcal{I}_6(j)>\delta_1] =\tilde{O}_{\mathcal{P}}(\delta_3/\delta_1). \label{e138} \end{align}\tag{90}\] Returning to 89 , we see that under the event \(\mathcal{I}_4(j) =\mathcal{I}_8 =\mathcal{I}_6(j) =0\) we have \[\begin{align} |t(j,B^{(j)},Y)-t_j| =\tilde{O}_{\mathcal{P}} \left(\sqrt{\frac{\varepsilon}{\delta_3}}\ln \frac{1}{\varepsilon} \right). \label{e143} \end{align}\tag{91}\] From 88 , 90 and Lemma 26, we have \[\begin{align} \mathbb{P}[\frac{1}{p}\sum_j\mathcal{I}_4(j) \vee \mathcal{I}_6(j)>2\delta_1 \textrm{ or }\mathcal{I}_8] =\tilde{O}_{\mathcal{P}}(p^{-1/3}\delta_1^{-1}\varepsilon^{-1}) + \tilde{O}_{\mathcal{P}}(\delta_3/\delta_1) +n^{-8}. \label{e144} \end{align}\tag{92}\] The values of \(\delta_1,\delta_3\), and \(\varepsilon\) can be optimized; here we specify them in a way such that the Levy-Prokhorov error in Corollary 6 will be minimized. Thus we will set 91 , \(2\delta_1\), and the right side of 92 to be the same order (up to log factors), which gives \[\begin{align} \sqrt{\frac{\varepsilon}{\delta_3}} =\delta_1 =p^{-1/3}\delta_1^{-1}\varepsilon^{-1} =\frac{\delta_3}{\delta_1}, \label{e95opt95d1} \end{align}\tag{93}\] yielding \(\delta_1=p^{-1/18}\), \(\varepsilon=p^{-2/9}\), and \(\delta_3=p^{-1/9}\). which yields the claimed convergence rate.

It remains to show the claim about further simplification when \(B^{(j)}_{:j}\) and \(A_{:j}\) are independent of \(A_{:\setminus j}\). Recall that for each \(j\), \((I-P_j)A_{:j}\) is a function of \((A,Y)\), and \(\check{B}^{(j)}_{:j}\) is \(\tilde{O}_{\mathcal{P}}(1)\)-sub-Gaussian conditioned on a set \(\mathcal{G}_j\) of \(A\) with probability \(1-O(n^{-10})\) (by Lemma 24). Therefore conditioned on \(G_j\) and \(\mathcal{I}_8\), \(\frac{1}{n}\check{B}^{(j)\top}_{:j}(I-P_j)A_{:j}\) is sub-Gaussian with variance proxy \(\tilde{O}_{\mathcal{P}}(\|(I-P_j)A_{:j}\|_2^2/n^2)=\tilde{O}_{\mathcal{P}}(\frac{1}{n})\). Hence conditioned on \(\mathcal{G}_j\) and \(\mathcal{I}_8\), \(|\frac{1}{n}\check{B}^{(j)\top}_{:j}(I-P_j)A_{:j}| \le \tilde{O}_{\mathcal{P}}(\frac{1}{\sqrt{n}})\) with probability at least \(1-n^{-20}\). Therefore \[\begin{align} |\frac{1}{n}\check{B}^{(j)\top}_{:j}(I-P_j)A_{:j}| \le \tilde{O}_{\mathcal{P}}(\frac{1}{\sqrt{n}}), \quad \forall j\in\{1,\dots,p\} \end{align}\] with probability at least \(1-pn^{-20}-\sum_{j=1}^p\mathbb{P}[\mathcal{G}_j^c]-\mathbb{E}[\mathcal{I}_8]\ge 1-\tilde{O}_{\mathcal{P}}(n^{-1/18})\). Then \[\begin{align} |\frac{1}{n}\check{B}^{(j)\top}_{:j}(I-P_j)A_{:j} \hat{\alpha}_j|=\tilde{O}_{\mathcal{P}}((n\delta_3)^{-1/2}), \quad \forall j\colon \mathcal{I}_6(j)=0 \end{align}\] with probability at least \(1-\tilde{O}_{\mathcal{P}}(n^{-1/18})\). Thus the \(\tilde{O}_{\mathcal{P}}(n^{-1/18})\) rate bound also holds when the term \(\check{B}^{(j)\top}_{:j}(I-P_j)A_{:j}\hat{\alpha}_j\) is dropped from the definition of \(t_j\) in Theorem 4.

10 The Case of non-vanishing \(\mathbb{E}[|\check{B}^{(j)}_{1j}|^2]\)↩︎

In this section we prove Theorem 5, which is essentially based on showing that the denominator \(\frac{1}{n}\check{B}_{:j}^{\top}(I-P_j)B_{:j}\) is bounded below.

10.1 Auxiliary lemmas↩︎

10.2 Proof of Theorem 5↩︎

Recall that in 91 and with the choice of parameters in 93 , we established that \[\begin{align} |\frac{1}{n} \check{B}_{:j}^{(j)\top}(I-P^{(j)}) B^{(j)}_{:j}\hat{\beta}^{(j)U}_j-\frac{1}{n}\check{B}_{:j}^{(j)\top}R - \frac{1}{n}\check{B}^{(j)\top}_{:j}(I-P_j)A_{:j} \hat{\alpha}_j| =\tilde{O}_{\mathcal{P}} \left(n^{-1/18} \right) \label{e155} \end{align}\tag{98}\] for all \(j\) except for a set of size at most \(2n^{-1/18}p\), under the event \[\begin{align} \frac{1}{p}\sum_j\mathcal{I}_4(j) \vee \mathcal{I}_6(j)\le 2\delta_1 \textrm{ and }\mathcal{I}_8=0 \label{e95event} \end{align}\tag{99}\] which has probability \(1-\tilde{O}_{\mathcal{P}}(n^{-1/18})\). Now by Lemma 24, there exists a set \(\mathcal{E}\) of \(A_{:\setminus j}\), \(\mathbb{P}[\mathcal{E}]\le n^{-10}\), such that conditioned on any \((A,Y)\) satisfying \(A_{:\setminus j}\in \mathcal{E}^c\), \(\frac{1}{n}\check{B}_{:j}^{(j)\top}R\) is sub-Gaussian with variance proxy \(O_{\mathcal{P}}(\frac{1}{n^2}\|R\|_2^2\ln n)\le O_{\mathcal{P}}(\frac{1}{n^2}\|Y\|_2^2\ln n)\). Thus conditioned on any \(A_{:\setminus j}\in \mathcal{E}^c\), \(\frac{1}{n}\check{B}_{:j}^{(j)\top}R =O_{\mathcal{P}}(\ln n)\) with probability \(1-e^{-n}\). Thus unconditionally, \[\begin{align} \frac{1}{n}\check{B}_{:j}^{(j)\top}R =O_{\mathcal{P}}(\ln n) \label{e152} \end{align}\tag{100}\] with probability \(1-O(n^{-10})\). Under the event of 99 and for the same set of \(j\), by a similar argument, we have \(\|\frac{1}{n}(I-P_j)A_{:j}\hat{\alpha}_j\|_2 \le \frac{1}{n}\|A_{:j}\|_2\delta_3^{-1/2}=O_{\mathcal{P}}(\frac{1}{\sqrt{n}}\delta_3^{-1/2})\) and hence \[\begin{align} \frac{1}{n}\check{B}^{(j)\top}_{:j}(I-P_j)A_{:j} \hat{\alpha}_j =\tilde{O}_{\mathcal{P}}(\frac{1}{\sqrt{n}}\delta_3^{-1/2}) \label{e153} \end{align}\tag{101}\] with probability \(1-O(n^{-10})\). Moreover, from 99 and Lemma 26 we see that with probability \(1-\tilde{O}_{\mathcal{P}}(n^{-1/18})\), \[\begin{align} \frac{1}{n} \check{B}_{:j}^{(j)\top}(P^{(j)}-P_j) B^{(j)}_{:j}=o_{\mathcal{P}}(1) \end{align}\] for all \(j\) except for a set of size \(\tilde{O}_{\mathcal{P}}(pn^{-1/18})\). Then using 96 we have that with probability \(1-\tilde{O}_{\mathcal{P}}(n^{-1/18})\), \[\begin{align} \frac{1}{n} \check{B}_{:j}^{(j)\top}(I-P^{(j)}) B^{(j)}_{:j} = \frac{1}{n} \check{B}_{:j}^{(j)\top}(I-P_j) B^{(j)}_{:j} -o_{\mathcal{P}}(1) =\Omega_{\mathcal{P}}(1) \end{align}\] for all \(j\) except for a set of size \(\tilde{O}_{\mathcal{P}}(pn^{-1/18})\). Thus under 99 , 100 and 101 , from 98 we see that \[\begin{align} \hat{\beta}^{(j)U}_j =O_{\mathcal{P}}(\ln n) \end{align}\] and hence \[\begin{align} &\quad|\frac{1}{n} \check{B}_{:j}^{(j)\top}(I-P_j) B^{(j)}_{:j}\hat{\beta}^{(j)U}_j-\frac{1}{n}\check{B}_{:j}^{(j)\top}R - \frac{1}{n}\check{B}^{(j)\top}_{:j}(I-P_j)A_{:j} \hat{\alpha}_j| \nonumber\\ &\le \tilde{O}_{\mathcal{P}} \left(n^{-1/18} \right) + \left|\frac{1}{n} \check{B}_{:j}^{(j)\top}(P^{(j)}-P_j) B^{(j)}_{:j}\hat{\beta}^{(j)U}_j\right| \\ &\le \tilde{O}_{\mathcal{P}} \left(n^{-1/18} \right). \end{align}\] The claim then follows by dividing both sides by \(\frac{1}{n} \check{B}_{:j}^{(j)\top}(I-P_j) B^{(j)}_{:j}\), and noting that 99 , 100 and 101 hold simultaneously with probability at least \(1-\tilde{O}_{\mathcal{P}} \left(n^{-1/18}\right)\).

11 Gaussian designs: Proof of Theorem 7↩︎

In this section we prove that our definition of the generalized debiased estimators \(\hat{\beta}^{U(j)}\), \(\hat{\alpha}^U\) asymptotically matches the traditional debiased estimators \(\hat{\beta}^{u(j)}\), \(\hat{\alpha}^u\) in the case of Gaussian designs. Recall eqrefeqq118 with parameters chosen in 93 states that \[\begin{align} \mathbb{P} \left[ \frac{1}{p}\sum_j \mathcal{I}_4(j) >p^{-1/18} \right] \le \tilde{O}_{\mathcal{P}}(p^{-1/18}). \end{align}\] Thus with probability at least \(1-\tilde{O}_{\mathcal{P}}(p^{-1/18})\), there is a set of \(j\) of cardinality at least \(p(1-p^{-1/18})\), for which \[\begin{align} \frac{1}{n}\sum_l 1\{\chi_l^{\alpha} \neq \chi_l^{\beta^{(j)}} \} &\le p^{-2/9}, \label{e157} \end{align}\tag{102}\] and hence by Lemma 26, \[\begin{align} \max_{\Delta}\|(P_{\mathcal{A}}-P_{\mathcal{B}})\check{B}^{(j)}_{:j}\|_2 &\le \tilde{O}_{\mathcal{P}} (p^{\frac{7}{18}}). \label{e158} \end{align}\tag{103}\] Now in the Gaussian case, \(\check{B}_{:j}\) is independent of the rest of the columns, and \(\mu_{:j}=A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj}\). We have \[\begin{align} \frac{1}{n}\check{B}_{:j}^{\top}(I-P_{\mathcal{B}})B_{:j} &=\frac{1}{n}\check{B}_{:j}^{\top}(I-P_{\mathcal{B}})(\check{B}_{:j}+A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj}). \end{align}\] Then we see that with probability \(1-\tilde{O}_{\mathcal{P}}(p^{-1/18})\), \[\begin{align} \frac{1}{n}\check{B}_{:j}^{\top}(I-P_{\mathcal{A}})(\check{B}_{:j}+A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj}) &=\frac{1}{n}\check{B}_{:j}^{\top}(I-P_{\mathcal{A}})\check{B}_{:j}+O_{\mathcal{P}}(\frac{1}{\sqrt{n}}) \tag{104} \\ &=\Sigma_{j|\setminus j}(1-\frac{1}{n}\|\hat{\alpha}\|_0) +O_{\mathcal{P}}(n^{-2/9}) \tag{105} \\ &=\Sigma_{j|\setminus j}(1-\frac{1}{n}\|\hat{\beta}^{(j)}\|_0) +O_{\mathcal{P}}(n^{-2/9}) \tag{106} \end{align}\] where 104 follows since with probability \(1-e^{-n}\), \((I-P_{\mathcal{A}})A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj}\) has norm \(O_{\mathcal{P}}(\sqrt{n})\), and \(\check{B}_{:j}\) is independent of it with norm \(O_{\mathcal{P}}(\sqrt{n})\). 105 follows from concentration of the Chi-square distribution and the independence of \(P_{\mathcal{A}}\) and \(\check{B}_{:j}\). 106 follows from 102 . Moreover, with probability \(1-\tilde{O}_{\mathcal{P}}(p^{-1/18})\), we also have the error bound \[\begin{align} &\quad\frac{1}{n}|\check{B}_{:j}^{\top}(P_{\mathcal{A}}-P_{\mathcal{B}})(\check{B}_{:j}+A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj})| \nonumber\\ &\le \frac{1}{n}|\check{B}_{:j}^{\top}(P_{\mathcal{A}}-P_{\mathcal{B}})\check{B}_{:j}| + \frac{1}{n}\|\check{B}_{:j}^{\top}(P_{\mathcal{A}}-P_{\mathcal{B}})\|_2 \|A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj}\|_2 \\ &= \tilde{O}_{\mathcal{P}}(n^{-2/9}) + \frac{1}{n}\cdot\tilde{O}_{\mathcal{P}}(n^{7/18}) \cdot O_{\mathcal{P}}(\sqrt{n}) \\ &=\tilde{O}_{\mathcal{P}} (n^{-1/9}) \label{e140461} \end{align}\tag{107}\] where we used 103 . This establishes \(\frac{1}{n}\check{B}_{:j}^{\top}(I-P_{\mathcal{B}})B_{:j} = \Sigma_{j|\setminus j}(1-\|\hat{\beta}^{(j)}\|_0/n) +O_{\mathcal{P}}(n^{-1/9})\). Moreover, from Lemma 28 and 107 , there exists \(\lambda_{\mathcal{P}}\) such that for all \(\lambda \ge\lambda_{\mathcal{P}}\), with probability at least \(1-\tilde{O}_{\mathcal{P}}(p^{-1/18})\), set of \(j\) of cardinality at least \(p(1-p^{-1/18})\), the denominator \(\frac{1}{n}\check{B}^{\top}_{:j}(I-P^{(j)})B_{:j}\) is bounded away from 0, which implies \(|\hat{\beta}^{U(j)}_j -\hat{\beta}^{u(j)}_j| =O_{\mathcal{P}}(n^{-1/9})\).

We now turn to \(\hat{\alpha}^U\), whose proof is similar to the above. With probability at least \(1-\tilde{O}_{\mathcal{P}}(p^{-1/18})\), there is a set of \(j\) of cardinality at least \(p(1-p^{-1/18})\), for which 102 holds, and \[\begin{align} \frac{1}{n}\check{A}_{:j}^{\top}(I-P_{\mathcal{B}})(\check{A}_{:j}+A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj}) &= \Sigma_{j|\setminus j}(1-\frac{1}{n}\|\hat{\alpha}\|_0) +O_{\mathcal{P}}(n^{-2/9}) \label{e141461} \end{align}\tag{108}\] by steps similar to 106 . Now Let \(\mathcal{I}_1(j)\) be the probability that 86 is not satisfied, which, by symmetry, is also the probability that \[\begin{align} \max_{\Delta}\|(P_{\mathcal{B}\cup\Delta}-P_{\mathcal{B}\setminus\Delta})\check{A}^{(j)}_{:j}\|_2 &\le \tilde{O}_{\mathcal{P}}(p^{\frac{7}{18}}) \label{e141} \end{align}\tag{109}\] does not hold, where the max is over the same collection of \(\Delta\) as in 103 . By Lemma 26, \(\sum_{j=1}^p\mathcal{I}_1(j)\le n^{-8}p\). By Markov’s inequality, with probability at least \(1-n^{-4}\), 109 is violated only by a set of \(j\) of size at most \(n^{-4}p\). Then by 102 , 109 , and the union bound, we have \[\begin{align} \max_{\Delta}\|(P_{\mathcal{B}}-P_{\mathcal{A}})\check{A}^{(j)}_{:j}\|_2 &\le \tilde{O}_{\mathcal{P}}(p^{\frac{7}{18}}) \label{e142} \end{align}\tag{110}\] for all \(j\) in a good set of size at least \(p(1-p^{-1/18}-n^{-4})\), with probability at least \(1-\tilde{O}_{\mathcal{P}}(p^{-1/18})\), and hence \[\begin{align} &\quad\frac{1}{n}|\check{A}_{:j}^{\top}(P_{\mathcal{A}}-P_{\mathcal{B}})A_{:j}| \nonumber\\ &\le \frac{1}{n}|\check{A}_{:j}^{\top}(P_{\mathcal{A}}-P_{\mathcal{B}})\check{A}_{:j}| + \frac{1}{n}\|\check{A}_{:j}^{\top}(P_{\mathcal{A}}-P_{\mathcal{B}})\|_2 \|A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj}\|_2 \\ &= \tilde{O}_{\mathcal{P}}(n^{-2/9}) + \frac{1}{n}\cdot\tilde{O}_{\mathcal{P}}(n^{7/18}) \cdot O_{\mathcal{P}}(\sqrt{n}) \\ &=\tilde{O}_{\mathcal{P}} (n^{-1/9}) \label{e146} \end{align}\tag{111}\] where we used 110 . Then \(\frac{1}{n}\check{A}_{:j}^{\top}(I-P_{\mathcal{A}})A_{:j} = \Sigma_{j|\setminus j}(1-\|\hat{\alpha}\|_0/n) +O_{\mathcal{P}}(n^{-1/9})\) follows by 108 and 111 (the cardinality estimate \(p(1-3p^{-1/18})\) arises from taking the union bound for sets of \(j\)). Finally, from the proof of Lemma 28, there exists \(\lambda_{\mathcal{P}}\) such that for all \(\lambda \ge\lambda_{\mathcal{P}}\), with high probability \(1-\|\hat{\alpha}\|_0/n\) is bounded below and away from 0, and hence by 111 , the denominator \(\Sigma_{j|\setminus j}(1-\frac{1}{n}\|\hat{\alpha}\|_0)\) in the definition of \(\hat{\alpha}^u_j\) is bounded away from 0, which implies \(|\hat{\alpha}^U_j -\hat{\alpha}^u_j| =O_{\mathcal{P}}(n^{-1/9})\).

12 Proof of matrix continuity↩︎

13 Convergence for general regularizers↩︎

13.1 Auxiliary lemmas↩︎

13.2 Proof of Theorem 8↩︎

Recall from 61 that \[\begin{align} \lambda\|\psi(\hat{\alpha}_{\setminus j})- \psi(\hat{\beta}_{\setminus j}) \|_2 &\le \frac{1}{n}\|A_{:\setminus j}^{\top}(R-S)\|_2 \\ &\le \frac{1}{n}\|A\|_{\rm op} \|A_{:j}\hat{\alpha}_j-B_{:j}\hat{\beta}_j\|_2 \\ &\le \frac{1}{n}(\|A\|_{\rm op}^2\vee \|B\|_{\rm op}^2) (|\hat{\alpha}_j|+|\hat{\beta}_j|). \end{align}\] Under \(\inf_{t\in\mathbb{R}}\psi'(t)\ge \zeta_{2-}\), the left side above also upper bounds \(\lambda \zeta_{2-}\|\hat{\alpha}_{\setminus j}-\hat{\beta}_{\setminus j}\|_2\). Thus \[\begin{align} \|\frac{1}{n}A_{:\setminus j} (\bar{G}_{\mathcal{A}}^{-1}-\bar{G}_{\mathcal{B}}^{-1}) A_{:\setminus j}^{\top} \check{B}_{:j}\|_2 &\le \frac{c_2^{1.5}\zeta_3}{c_1^2\sqrt{n}} \|\bar{G}_{\mathcal{A}}^{-1} A_{:\setminus j}^{\top} \check{B}_{:j} \|_{\infty} \cdot \frac{1}{n\lambda \zeta_{2-}}(\|A\|_{\rm op}^2\vee \|B\|_{\rm op}^2) (|\hat{\alpha}_j|+|\hat{\beta}_j|) \\ &\le \frac{c_2^{1.5}\zeta_3}{c_1^2\lambda \zeta_{2-}}\Gamma D^2 (|\hat{\alpha}_j|+|\hat{\beta}_j|). \end{align}\] The rest of the proof is the same as Section 4, using Lemma 32 in place of Lemma 15.

13.3 Proof of Theorem 9↩︎

There exists \(c>0\) (depending only on \(\mathcal{P}\) and \(C\)) such that the following holds: there exists a set \(\mathcal{E}\) of \(A\) such that \[\begin{align} \mathbb{P}[A\in \mathcal{E}] \le n^{-C} \label{e9545} \end{align}\tag{117}\] and that for each \(A\in \mathcal{E}^c\), \[\begin{align} \|A\|_{\rm op} &\vee \max_{j=1}^p\|B^{(j)}\|_{\rm op}\le \frac{1}{c}\sqrt{n}; \tag{118} \\ \check{B}_j^{(j)} &:=B_j^{(j)}-\mu_{:j} \textrm{ is }\frac{\ln n}{c}\textrm{-sub-Gaussian.} \tag{119} \end{align}\] Here 118 is a standard results following from references around 80 , and 119 follows from Lemma 24 and the union bound. Then conditioned on any \(A\in\mathcal{E}^c\), we have \(\max_{1\le j \le p}\|\bar{G}_{\mathcal{A}}^{-1} A_{:\setminus j}^{\top}\|_{\rm op} \le \frac{\sqrt{n}}{c}\zeta_{2-}^{-1}\), and each coordinate of \(\bar{G}_{\mathcal{A}}^{-1} A_{:\setminus j}^{\top} \check{B}_{:j}^{(j)}\) is \(\frac{n\ln n}{c^3}\zeta_{2-}^{-2}\)-sub-Gaussian, and we have (by union bound) \[\begin{align} \max_{1\le j\le p}\|\bar{G}_{\mathcal{A}}^{-1} A_{:\setminus j}^{\top} \check{B}_{:j}^{(j)} \|_{\infty} \le O_{\mathcal{P},\zeta,C} (\sqrt{n}\ln n ), \end{align}\] with probability at least \(1-n^{-C}\), where the constants in \(O_{\mathcal{P},\zeta,C}\) may depend on parameters in \(\mathcal{P}\), \(\zeta:=(\zeta_{2-},\zeta_{2+},\zeta_3)\), and \(C\). Therefore, by taking \(D=\frac{1}{c}\) and \(\Gamma=\gamma_{\mathcal{P},\zeta,C}\ln n\) with sufficiently large \(\gamma_{\mathcal{P},\zeta,C}>0\), we have \(\mathcal{G}\) with probability at least \(1-O_{\mathcal{P}}(n^{-C})\), without conditioning on \(A\), where we used the union bound and 117 .

Next, we show that \[\begin{align} \sum_{j=1}^p\mathbb{E}[|\hat{\alpha}_j|^2+|\hat{\beta}_j|^2]=O_{\mathcal{P},\zeta}(n). \label{e48} \end{align}\tag{120}\] We first observe that by the optimality of \(\hat{\alpha}\), \[\begin{align} \frac{1}{2n}\|Y\|_2^2\ge \sum_{j=1}^n\rho(\hat{\alpha}_j) \ge\frac{\zeta_{2-}}{2}\|\hat{\alpha}\|_2^2, \end{align}\] implying \(\sum_{j=1}^p\mathbb{E}[|\hat{\alpha}_j|^2]=O_{\mathcal{P},\zeta}(n)\). To control \(\mathbb{E}[\hat{\beta}_j^2]\), we use arguments similar to the proof of Corollary 22. Define \(A(t)\in\mathbb{R}^{n\times p}\) such that \(A_{:\setminus j}(t)=A_{:\setminus j}\) and \(A_{:j}(t)=(1-t)A_{:j}+tB^{(j)}_{:j}\). Thus \(A(0)=A\) and \(A(1)=B^{(j)}\). Let \(Y(t):=A(t)\alpha+w\), and let \(\hat{\alpha}(t)\) be the solution to the penalized regression problem for the data \((A,Y(t))\). Thus \(\hat{\alpha}(0)=\hat{\alpha}\), and \(\hat{\alpha}(1)\) has the same distribution as \(\hat{\beta}^{(j)}\) (since \(B^{(j)}\) and \(A\) are identically distributed). The normal equation reads as \[\begin{align} -\frac{1}{n}A^{\top}(Y(t)-A\hat{\alpha}(t))+\psi(\hat{\alpha}(t))=0, \end{align}\] from which we obtain the derivative \[\begin{align} -\frac{1}{n}A^{\top} (B_{:j}^{(j)}-A_{:j})\alpha_j +\left(\frac{1}{n}A^{\top}A +\mathop{\mathrm{diag}}(\psi'(\hat{\alpha}(t)))\right) \hat{\alpha}'(t)=0. \end{align}\] Since \(\left\|\left(\frac{1}{n}A^{\top}A +\mathop{\mathrm{diag}}(\psi'(\hat{\alpha}(t)))\right)^{-1}\right\|_{\rm op}\le \zeta_{2-}^{-1}\), we have \(\|\hat{\alpha}'(t)\|_2\le |\alpha_j|O_{\mathcal{P},\zeta}(1)\). Thus \(|\hat{\alpha}_j(0)- \hat{\alpha}_j(1)| \le \|\hat{\alpha}(0)- \hat{\alpha}(1)\|_2\le |\alpha_j|O_{\mathcal{P},\zeta}(1)\) and \(|\hat{\alpha}_j(1)|^2 \le 2|\hat{\alpha}_j(0)|^2 +2|\hat{\alpha}_j(1)-\hat{\alpha}_j(0)|^2 \le 2\hat{\alpha}_j^2+|\alpha_j|O_{\mathcal{P},\zeta}(1)\). Then \(\sum_{j=1}^p\mathbb{E}[|\hat{\beta}_j|^2]= \sum_{j=1}^p\mathbb{E}[|\hat{\alpha}_j(1)|^2] =O_{\mathcal{P},\zeta}(n)\), verifying 120 .

14 Experiments↩︎

Results are generated by the following codes, which can be found in https://github.com/jingboliu1/local_change.git

• Figures 3–8: localk1_fig1and2.R.

• Table 1: err_eval.R.

• Table 2: err_eval1.R.

• Tables 3 and 5: FDR_demo_eq_table1.R.

• Table 4: FDR_demo_eq_table2.R.

• Table 6: FDR_demo_riboflavin_table3.R.

• Table 7: FDR_HIV.R.

14.1 Table of approximation errors in the update formula↩︎

Consider the setting where the rows of the design matrix \(A\in \mathbb{R}^{n\times p}\) are i.i.d.\(\mathcal{N}(0,\Sigma)\), where \[\begin{align} \Sigma_{ij}=\rho^{|i-j|}. \end{align}\] Let \(\alpha\in\mathbb{R}^p\) be such that the first \(s\) entries are generated i.i.d. according to \(\mathcal{N}(0,1)\), and the remaining entries are 0. Take \[\begin{align} Y=A\alpha+w \end{align}\] where \(w\sim \mathcal{N}(0,0.01nI_{n\times n})\) is an independent Gaussian noise. We take \(\require{upgreek} n=800\upalpha\), \(\require{upgreek} p=1000\upalpha\), \(\require{upgreek} s=400\upalpha\), where \(\require{upgreek} \upalpha\) controls the scale of the problem. We uniformly sample \(\require{upgreek} \lfloor \frac{p}{10\upalpha}\rfloor\) coordinates \(j\) in \(\{1,\dots,p\}\), resample the corresponding column in \(A\) with the conditional distribution \(P_{A_{:j}|A_{:\setminus j}}\), and compute the Lasso estimator \(\hat{\beta}^{(j)}_j\) and the debiased estimator \(\hat{\beta}^{(j)U}_j\) (both the exact values as well as the approximate values via the approximation formulas in Section 3.1 in the main paper). We compute the normalized approximation error (mean squared error normalized by the square of the \(\ell_2\) norm) of the approximation formulas across \(\require{upgreek} \upalpha\) and \(\rho\). The results are shown in Table 1. These results confirm that the approximation error for the debiased estimator is consistently smaller than that of the standard Lasso, with the improvement becoming more pronounced as the correlation \(\rho\) increases, in agreement with the theoretical analysis.

In Table 2, we choose \(\rho=0.5\), and use the same method for generating \(A\), but generate \(B\) by \(B_{ij}:=\mathbb{E}[A_{ij}|A_{i\setminus j}] +((\Sigma^{-1})_{jj})^{-1/2}\Xi_{ij}\), where \(\Xi_{ij}\) are i.i.d. random variables satisfying \(\Xi_{ij}=\pm (\xi)^{-1/2}\) with probability \(\xi/2\), and \(\Xi_{ij}=0\) with probability \(1-\xi\). In other words, \(\check{B}_{ij}\) is now a 3-point distribution with the same mean and variance as in Table 1. This demonstrates the validity of the approximation formula beyond Gaussian designs, consistent with the theoretical analysis.

14.2 Plots of approximation errors in the update formula↩︎

We consider the setting in Section 14.1, fixing \(n=1000\), \(p=1200\), \(s=500\), In Fig. 3 - Fig. 5, we plot the values of \(\hat{\beta}^{(j)U}_j\) and \(\hat{\beta}^{(j)U}_j-\tilde{\gamma}_j\) for \(\rho=0,0.5,0.95\), where \(\hat{\beta}^{(j)U}_j\) is the debiased Lasso coefficient computed exactly according to the definition, and \(\tilde{\gamma}_j\) denotes the value computed using the approximation formula. Fig. 6 - Fig. 8 compare \(\hat{\beta}^{(j)}_j\) and its approximation error for \(\rho=0,0.5,0.95\). In these plots, we only uniformly select 1/12 of all the coordinates, to avoid cluttering of the picture. It can be seen that the approximation error for the debiased estimator is better than the plain Lasso for large \(\rho\) (in turns of size of the error relative to the magnitude of the debiased coefficients), which is consistent with the theoretical analysis.

14.3 FDR control experiments with synthetic data↩︎

Consider an FDR control problem similar to the setting in Theorem 13 in the main paper, where \(p=300\), \(n=200\), and the target FDR \(q= 0.1\). We take \(\Sigma^{-1}=a_p(E+\epsilon I)\), where \(E\) denotes the matrix whose entries are all 1, and \(a_p,\epsilon>0\) are parameters to be specified later. Then we have \(\Sigma=\frac{1}{a_p\epsilon} \left( I-\frac{\epsilon^{-1}}{1+\epsilon^{-1}p}E \right)\). We then generate \(\alpha\) with a random set of \(s\) coordinates equal to \(A_{\rm val}/\sqrt{n}\) (\(A_{\rm val}>0\) being a parameter to be specified), and the rest coordinates equal to 0. The observation is \(Y=A\alpha+w\), where \(w\sim \mathcal{N}(0,\sigma I)\).

We compare the performance of 6 variable selection methods in Table 3 and Table 4. For the knockoff method, we will use the eq-knockoff construction of [1], which is natural in this setting since \(\Sigma\) has equal values for the off-diagonals and for the diagonals. The max eigenvalue of \(\Sigma^{-1}\) is \(a_p(p+\epsilon)\), so the condition \(2S^{-1}\succeq \Sigma^{-1}\) in eq-knockoff becomes \(S_{11}=\dots=S_{pp}\le \frac{2}{\lambda_{\rm max}(\Sigma^{-1})}=\frac{2}{a_p(p+\epsilon)}\). To design the knockoff filter, one tries to minimize \(|\Sigma_{11}-S_{11}|=\left|\frac{1+\epsilon^{-1}(p-1)}{a_p(\epsilon+p)}-S_{11}\right|\) and so for \(0<\epsilon<p-1\) the optimal \(S_{11}=\frac{2}{a_p(p+\epsilon)}\). For \(\epsilon>p-1\), the optimal \(S_{11}=\frac{1+\epsilon^{-1}(p-1)}{a_p(\epsilon+p)}\).

In Table 3, we use \(\epsilon=50,a_p=1/p\), number of nonzero coefficients \(s=20\), and noise standard deviation \(\sigma = 1/\sqrt{p}\), and vary the signal strength \(A_{\rm val}\). We see that FDR is controlled at approximately \(q=0.1\) in most cases. As \(s=20\) is relatively small in this setting, there are a few instances of FDR overflow for Knockoff-db, due to fluctuations. Meanwhile, the power achieved by the local knockoff filter and CRT are better than the knockoff filter, with or without debiasing.

In Table 4, we change the sparsity level to \(s=60\). Again the FDR is controlled at approximately \(q=0.1\) in most cases, while local knockoff filter and CRT achieve higher power than the knockoff filter. Note that in contrast to Table 3, the power for the debiased versions are noticeably better than without debiasing in Table 4, which is expected since \(s\) increased. For correlated designs, the debiased coefficients tend to outperform the standard coefficients when \(s\) is not too small.

In Table 5, we increase the dimensions to \(p=1000\), \(n=600\), and sparsity to \(s=200\). Now FDR is well-controlled at approximately \(q=0.1\) in all cases. The local knockoff filter and CRT (especially the debiased versions) achieve high power, whereas the knockoff filters fail. This confirms the theoretical analysis in Theorem 13 in the main paper: when the precision matrix is close to \(E\) (the matrix consisting of 1’s), the knockoff filter fails in the high-dimensional limit, regardless of the choice of the knockoff mechanism, whereas methods based on more relaxed local exchangeability conditions (such as local knockoff and CRT) remains powerful.

14.4 FDR control with Riboflavin data↩︎

We use the riboflavin dataset, available in the supplemental materials of [57], which was widely used in FDR control experiments [[58]][57][59] . It contains a \(p_0=4088\) by \(n=71\) matrix of the logarithm of the expression levels, and a response vector \(Y\) of the logarithm of the riboflavin production rate. The original measurement matrix contains measurements of many similar (highly correlated) genes, so we use the findCorrelation function in R to remove the highly correlated columns with cutoff \(\rho=0.5\). We normalize the means and variance of the columns, and use graphical Lasso function glasso to estimate the covariance matrix \(\Sigma\) of the features. Then we use the best linear estimator \(A_{:\setminus j}\Sigma_{\setminus j}^{-1}\Sigma_{\setminus jj}\) for the \(\mu_{:j}\) in the definition of the debiased estimator. The FDR and power cannot be precisely evaluated since we do not know the ground truth. To tackle this issue, we first use cross-validated Lasso to obtain \(\alpha\) for the observed \(Y\), and then generate new \(Y=A\alpha+w\), where the noise level is estimated using the norm of the residual in the previous Lasso regression, so that we can calculate the FDR and power using the new \(Y\) and \(\alpha\). Previously, a similar approach for testing FDR control methods on real datasets was adopted in the literature; see for example [60]. The results of the FDR and power values are shown in Table 6. We see an increase in the power by using local knockoffs or CRT, while roughly controlling the FDR in most cases.

14.5 FDR control with HIV data↩︎

We use the same procedure as in Section 14.4, but apply it to the HIV Drug Resistance Database, available at https://hivdb.stanford.edu. We use the code from [5] to extract the design matrix and response vector. For this dataset, we have \(n = 2026\), \(p = 163\), and sparsity level \(s = 79\). For the knockoff function, we use the implementation from the official package. The results are shown in Table 7. While all methods kept FDR below approximately \(q=0.1\), the debiased versions of local knockoff and CRT achieve higher powers.

15 Proofs and implementation details for variable se lection↩︎

15.1 Proof of Lemma 11↩︎

Set \(a=\frac{1}{\sqrt{2}}S^{1/2}e\) where \(e=(1,\dots,1)^{\top}\). From \(\|s\|_1=s^{\top}S^{-1}s\ge\frac{1}{2}s^{\top}Es=\frac{1}{2}\|s\|_1^2\) we obtain \(\|a\|^2=\frac{1}{2}\|s\|_1\le 1\). Moreover, \[\begin{align} (I-\frac{1}{2}S^{1/2}ES^{1/2})^{-1} &=(I-aa^{\top})^{-1} \\ &=I+\frac{1}{1-\|a\|_2^2}aa^{\top}. \end{align}\] Therefore, the diagonal values of \((2S-S\Sigma^{-1}S)^{-1}\) are \[\begin{align} \frac{1}{2s_j}(1+\frac{a_j^2}{1-\|a\|_2^2}) &=\frac{1}{2s_j}(1+\frac{s_j/2}{1-\|s\|_1/2}) \\ &=\frac{1}{2s_j}+\frac{1}{4-2\|s\|_1} \\ &>\frac{1}{2s_j}. \end{align}\] However, by the Markov inequality we have \[\begin{align} \frac{1}{p}|\{j\colon \frac{1}{2s_j}>\frac{p}{10}\}| &=1-\frac{1}{p}|\{j\colon s_j\ge\frac{5}{p}\}| \\ &\ge 1-\frac{1}{5}\|s\|_1 \\ &\ge \frac{3}{5}. \end{align}\]

15.2 Proof of Theorem 13↩︎

From the definition of \(\Sigma^{-1}\) it is easy to see that the empirical distribution of its diagonals converges to 0. From the Schur complement theorem we know that the diagonals of \(\underline{\Sigma}^{-1}\) is two copies of the diagonals of \((2S-S\Sigma^{-1}S)^{-1}\), which we denote by \(d_1,\dots,d_p\). From Lemma 11 we see that \(\frac{1}{p}|\{j\colon d_j>\frac{pa_p}{10}\}|\ge \frac{3}{5}\), hence \(\frac{1}{2p} |\{j\colon (\underline{\Sigma}^{-1})_{jj} >\frac{1}{10}\}|\ge \frac{3}{5}\).

15.3 Proof of Theorem 14↩︎

Under Assumption 1, for large enough \(n\) we have the following with high probability: \[\begin{align} {\rm FDR} &=\mathbb{E}\left[\frac{|\{j\in H_0\colon |\hat{\alpha}^U_j| >T(\epsilon)\}|}{|\{j\colon |\hat{\alpha}^U_j|>T(\epsilon)\}|} \right] \\ &\le \frac{|\{j\in H_0\colon |\hat{\alpha}^U_j|>T(\epsilon)-0.1\epsilon\}|}{|\{j\colon |\hat{\alpha}^U_j|>T(\epsilon)\}|} \\ &=\frac{\hat{P}_{\hat{\alpha}^U\alpha}[(T(\epsilon)-0.1\epsilon,\infty)\times \{0\}]}{\hat{P}_{\hat{\alpha}^U}[(T(\epsilon),\infty)]} \\ &\le \frac{P_{\hat{\alpha}^U\alpha}[(T(\epsilon)-0.2\epsilon,\infty)\times \{0\}]}{\hat{P}_{\hat{\alpha}^U}[(T(\epsilon),\infty)]} \\ &=\frac{P_{\hat{\gamma}\alpha}[(T(\epsilon)-0.2\epsilon,\infty)\times \{0\}]}{\hat{P}_{\hat{\alpha}^U}[(T(\epsilon),\infty)]} \\ &=\frac{P_{\tilde{\gamma}\alpha}[(T(\epsilon)-0.3\epsilon,\infty)\times \{0\}]}{\hat{P}_{\hat{\alpha}^U}[(T(\epsilon),\infty)]} \\ &=\frac{\hat{P}_{\tilde{\gamma}\alpha}[(T(\epsilon)-\epsilon,\infty)\times \{0\}]}{\hat{P}_{\hat{\alpha}^U}[(T(\epsilon),\infty)]} \\ &\le q. \end{align}\] Since the FDP is bounded by 1, convergence in probability implies convergence of expectations.

15.4 A note on computation of \(\hat{\alpha}^U\)↩︎

We discuss how to efficiently compute \(\hat{\alpha}^U\); see Algorithm 9. Recall the definition of \(\hat{\alpha}^U\): \[\begin{align} \hat{\alpha}^U_j &=\hat{\alpha}_j +\left(\frac{1}{n}\check{A}_{:j}^{\top} (I-P_j)A_{:j} \right)^{-1} \frac{\check{A}^{\top}_{:j}R}{n} \end{align}\] where \(P_j\) denotes the projection onto the column space of \(A_{:\mathcal{S}_j}\), \(\mathcal{S}_j:=\{l\neq j\colon \hat{\alpha}_l\neq 0\}\), and \(\Theta:=\Sigma^{-1}\) which can be computed in time \(O(n^3)\). If we use the formula \(P_j=A_{:\mathcal{S}_j}(A_{:\mathcal{S}_j}^{\top}A_{:\mathcal{S}_j})^{-1}A_{:\mathcal{S}_j}^{\top}\), then computing \(P_j\) takes time \(O(n^3)\), so that computing \(\hat{\alpha}^U\) takes time \(O(pn^3)=O(n^4)\). However, it is possible to efficiently compute the \(P_j\)’s using a “rank-one-update” formula so that the time of computing \(\hat{\alpha}^U\) is \(O(n^3)\): First set \[\begin{align} \Omega:=A_{:\mathcal{S}}^{\top}A_{:\mathcal{S}} \end{align}\] and \(\Pi=\Omega^{-1}\), which can be computed in time \(O(n^3)\). Then for each \(j\), note that by the Schur complement theorem, \[\begin{align} (A_{:\mathcal{S}_j}^{\top}A_{:\mathcal{S}_j})^{-1} ={\Omega_{\mathcal{S}_j}}^{-1} =\Pi_{\mathcal{S}_j} -\Pi_{\mathcal{S}_jj}\Pi_j^{-1}\Pi_{j\mathcal{S}_j}. \end{align}\] We can compute \(\check{A}_{:j}=A_{:j}+A_{:\setminus j}\Theta_{\setminus jj}\Theta_j^{-1}\) in time \(O(n^2)\). Since \(\Pi_{\mathcal{S}_jj}\) is a vector, we can compute \(\check{A}_{:j}^{\top}\Pi_{\mathcal{S}_jj}\Pi_j^{-1}\Pi_{j\mathcal{S}_j}A_{:j}\) in time \(O(n^2)\). We can also compute \[\begin{align} A_{:\mathcal{S}_j}\Pi_{\mathcal{S}_j}A_{:\mathcal{S}_j}^{\top}=A_{:\mathcal{S}}\Pi A_{:\mathcal{S}}^{\top}-(A_{:\mathcal{S}}\Pi_{\mathcal{S}j})A_{:j}^{\top} -A_{:j}(\Pi_{j\mathcal{S}_j}A_{:\mathcal{S}_j}^{\top}) \end{align}\] in time \(O(n^2)\). Then \(P_j=A_{:\mathcal{S}_j}(A_{:\mathcal{S}_j}^{\top}A_{:\mathcal{S}_j})^{-1}A_{:\mathcal{S}_j}^{\top}\) can be computed in time \(O(n^2)\) for each \(j\). Hence \(\hat{\alpha}^U\) is computed in \(O(n^3)\) time.

Acknowledgment↩︎

This research was supported in part by NSF Grant DMS-2515510.

References↩︎

---

1. Department of Statistics, University of Illinois Urbana-Champaign. Also affiliated with Department of Electrical and Computer Engineering, the Grainger College of Engineering. jingbol@illinois.edu
   Mathematics Subject Classification: Primary 62J07, 62E17; Secondary 62F12, 62F05.↩︎