Debiasing Precision Estimates For Edge Inference

Why Correct a Regularized Precision Estimate?

In a Gaussian graphical model, an off-diagonal entry of the precision matrix is zero exactly when the corresponding variables are conditionally independent given the others. This makes precision entries natural targets for edge inference. But in a high-dimensional problem, the unpenalized estimate may be unstable or unavailable, so regularization is used. Regularization stabilizes estimation and encourages sparsity, but it also shrinks estimates. A small regularized entry is therefore not, by itself, evidence that the corresponding population entry is small.

Debiasing uses the score equation of the Gaussian likelihood to remove the leading effect of that shrinkage. It does not make every regularized estimate unbiased, nor does it guarantee valid inference without assumptions. The key is to separate the corrected estimate into a tractable random leading term and a remainder that is negligible at the scale of sampling error.

Model, Estimator, and Exact Decomposition

Let \(X_1,\ldots,X_n\) be independent observations from \(N_p(0,\Sigma)\), where \(\Sigma\) is positive definite, and write \(\Omega=\Sigma^{-1}\). Assume the mean is known to be zero and use the sample second-moment matrix \(S=n^{-1}\sum_{k=1}^n X_kX_k^T\). Let \(\hat\Omega\) be a regularized precision estimate, such as a graphical lasso estimate, and define \(E=S-\Sigma\) and \(\Delta=\hat\Omega-\Omega\). The correction below uses the Gaussian likelihood score through \(S\); it does not require \(\hat\Omega\) to be an unpenalized likelihood estimate.

The following identity shows exactly what the correction accomplishes. It also identifies the condition needed for the leading random term to govern the sampling distribution.

To verify the identity, expand the product using \(S=\Sigma+E\) and \(\hat\Omega=\Omega+\Delta\). Since \(\Omega\Sigma=\Sigma\Omega=I\), \[ (\Omega+\Delta)\Sigma(\Omega+\Delta) =\Omega+2\Delta+\Delta\Sigma\Delta. \] The part involving \(E\) is \[ (\Omega+\Delta)E(\Omega+\Delta) =\Omega E\Omega+\Delta E\Omega+\Omega E\Delta+\Delta E\Delta. \] Substituting both expansions into \(2(\Omega+\Delta)-\hat\Omega S\hat\Omega-\Omega\) gives the identity. Thus \(-\Omega E\Omega\) is the leading random term, and the other four terms are the remainder.

The Leading Term and Its Variance

Write \(X_k\) for observation \(k\), and put \(Y_k=\Omega X_k\). Because \(X_k\sim N_p(0,\Sigma)\), the transformed vector satisfies \(Y_k\sim N_p(0,\Omega)\). The \((i,j)\) entry of the leading term has the exact representation \[ (-\Omega E\Omega)_{ij} =-\frac{1}{n}\sum_{k=1}^n \bigl(Y_{ki}Y_{kj}-\Omega_{ij}\bigr). \] The summands have mean zero. For jointly Gaussian, mean-zero variables, Isserlis’ formula gives \(E[Y_i^2Y_j^2]=\Omega_{ii}\Omega_{jj}+2\Omega_{ij}^2\). Therefore \[ \operatorname{Var}(Y_iY_j) =E[Y_i^2Y_j^2]-\{E[Y_iY_j]\}^2 =\Omega_{ii}\Omega_{jj}+\Omega_{ij}^2. \] Set \(v_{ij}=\Omega_{ii}\Omega_{jj}+\Omega_{ij}^2\). This is positive because \(\Omega\) is positive definite, so \(\Omega_{ii}\) and \(\Omega_{jj}\) are positive.

The condition says that the \((i,j)\) entry of the exact remainder is negligible compared with the \(n^{-1/2}\) sampling scale. It is an assumption to be justified by the estimation rates and structural conditions of the problem; it does not follow merely from using the correction. Here is the complete limit argument. The variables \(Y_{ki}Y_{kj}-\Omega_{ij}\) are independent and identically distributed across \(k\), have mean zero, and have finite variance \(v_{ij}\). The classical central limit theorem therefore implies \[ \sqrt n\,(-\Omega E\Omega)_{ij}\xrightarrow{d}N(0,v_{ij}). \] The stated remainder condition and Slutsky’s theorem give the claimed limit for \(\tilde\Omega\). \(\square\)

Standardized Approximation When the Dimension or Indices Vary

When \(p\) or the target indices vary with \(n\), the variance factor \(v_{ij}\) may itself vary. In that setting it is not enough to say that it is bounded above and away from zero and then claim convergence to one fixed normal distribution: a sequence of normal distributions with alternating variances need not converge to a single limit. A natural statement is instead a standardized normal approximation.

To prove the leading-term claim, define \(U_{k,n}=Y_{k,n,i_n}/\sqrt{(\Omega_n)_{i_ni_n}}\) and \(V_{k,n}=Y_{k,n,j_n}/\sqrt{(\Omega_n)_{j_nj_n}}\), where \(Y_{k,n}=\Omega_nX_{k,n}\). These are standard normal variables with correlation \(\rho_n=(\Omega_n)_{i_nj_n}/\sqrt{(\Omega_n)_{i_ni_n}(\Omega_n)_{j_nj_n}}\), so \(|\rho_n|\leq 1\). Their product has mean \(\rho_n\) and variance \(1+\rho_n^2\). The summands in the standardized leading term are independent across \(k\), have mean zero, and each has variance one after division by \(\sqrt{1+\rho_n^2}\).

For completeness, their fourth moments are bounded uniformly even if \(\Omega_n\) has no uniform eigenvalue bound. If \(Z_n=(U_{k,n}V_{k,n}-\rho_n)/\sqrt{1+\rho_n^2}\), then \((a-b)^4\leq 8(a^4+b^4)\), the Cauchy–Schwarz inequality gives \(E[U_{k,n}^4V_{k,n}^4]\leq\sqrt{E[U_{k,n}^8]E[V_{k,n}^8]}=105\), and \(|\rho_n|\leq1\). Since \(1+\rho_n^2\geq1\), it follows that \(E[Z_n^4]\leq8(105+1)=848\). The Lyapunov ratio for the row of \(n\) standardized summands is consequently at most \(848/n\), which tends to zero. Lyapunov’s central limit theorem proves the standardized leading-term convergence. The stated remainder condition and Slutsky’s theorem prove the final claim. \(\square\)

From the Approximation to Edgewise Inference

For a fixed-dimensional problem, or under suitable high-dimensional consistency conditions, estimate the variance factor by \[ \hat v_{ij}=\tilde\Omega_{ii}\tilde\Omega_{jj}+\tilde\Omega_{ij}^{\,2}. \] For this plug-in standard error to be valid, one needs \(\hat v_{ij}/v_{ij}\xrightarrow{p}1\) in the fixed-target setting. In a varying-dimension setting, require the analogous ratio \(\hat v_n/v_n\xrightarrow{p}1\). These are additional conditions, not consequences of the normal approximation alone. With them and the corresponding negligible-remainder condition, the statistic \[ Z_{ij}=\frac{\sqrt n\,\tilde\Omega_{ij}}{\sqrt{\hat v_{ij}}} \] is asymptotically standard normal under the null hypothesis \(H_0:\Omega_{ij}=0\). The factor \(\sqrt n\) appears because the variance of the estimate is approximately \(v_{ij}/n\).

For example, a two-sided single-entry test rejects at approximate level \(\alpha\) when \(|Z_{ij}|>z_{1-\alpha/2}\), where \(z_{1-\alpha/2}\) is the corresponding standard normal quantile. Equivalently, an approximate interval for \(\Omega_{ij}\) is \[ \tilde\Omega_{ij}\ \pm\ z_{1-\alpha/2}\sqrt{\hat v_{ij}/n}. \] These statements concern one target entry at a time. They do not, by themselves, provide simultaneous error control for testing many edges.

Common Pitfall: Treating the Penalized Estimate as Normally Centered

A tempting argument is that a regularized estimate is close to \(\Omega\), so its centered entries must have the same normal limit as an unregularized estimate. Closeness alone does not justify this: the shrinkage error can be larger than the \(n^{-1/2}\) sampling error. The exact decomposition shows why debiasing is useful: it isolates the leading term \(-\Omega E\Omega\), while making the remaining terms explicit. In a particular regime, those remainder terms still have to be shown negligible at the relevant scale. If that step fails, the normal approximation and the resulting edgewise test are not justified.