BF: computing inverse correctly after centering

Author: jonathan-taylor

selectiveInference/R/funs.ROSI.R

@@ -276,7 +276,14 @@ approximate_BN = function(X, active_set){
   nactive=length(active_set)
 
   svdX = svd(X)
-  inv = solve(svdX$u %*% diag(svdX$d^2) %*% t(svdX$u))
+  
+  rank = sum(svdX$d > max(svdX$d) * 1.e-9)
+  inv_d = 1. / svdX$d
+  if (rank < length(svdX$d)) { 
+     inv_d[(rank+1):length(svdX$d)] = 0.
+  }
+  inv = svdX$u %*% diag(inv_d^2) %*% t(svdX$u)
+
   D = rep(0, nactive)
   for (i in 1:nactive){
     var = active_set[i]

selectiveInference/man/ROSI.Rd

@@ -111,9 +111,13 @@ p is comparable to n. Ignored when n<p and set to TRUE.
 
 \references{
 
+
 Keli Liu, Jelena Markovic, Robert Tibshirani. More powerful post-selection 
 inference, with application to the Lasso. arXiv:1801.09037
 
+Tom Boot, Didier Nibbering. Inference in high-dimensional linear regression models.
+arXiv:1703.03282
+
 }
 \author{Jelena Markovic, Jonathan Taylor}
 
@@ -122,6 +126,7 @@ inference, with application to the Lasso. arXiv:1801.09037
 library(selectiveInference)
 library(glmnet)
 set.seed(43)
+
 n = 200
 p = 100
 s = 2
@@ -153,5 +158,58 @@ out = ROSI(x,
            dispersion=sigma^2)
 out
 
+# an alternate approximate inverse from Boot and Nibbering
+
+out = ROSI(x,
+           y,
+           beta,
+           lambda,
+           dispersion=sigma^2,
+           debiasing_method="BN")
+out
+
+# wide matrices work too
+
+n = 100
+p = 200
+s = 2
+sigma = 1
+
+x = matrix(rnorm(n*p),n,p)
+x = scale(x,TRUE,TRUE)
+
+beta = c(rep(10, s), rep(0,p-s)) / sqrt(n)
+y = x \%*\% beta + sigma*rnorm(n)
+
+# first run glmnet
+gfit = glmnet(x,y,standardize=FALSE)
+
+# extract coef for a given lambda; note the 1/n factor!
+# (and we don't save the intercept term)
+lambda = 4 * sqrt(n)
+lambda_glmnet = 4 / sqrt(n)
+beta = selectiveInference:::solve_problem_glmnet(x, 
+                                                 y, 
+                                                 lambda_glmnet, 
+                                                 penalty_factor=rep(1,p),
+                                                 family="gaussian")
+# compute fixed lambda p-values and selection intervals
+out = ROSI(x,
+           y,
+           beta,
+           lambda,
+           dispersion=sigma^2)
+out
+
+# an alternate approximate inverse
+
+out = ROSI(x,
+           y,
+           beta,
+           lambda,
+           dispersion=sigma^2,
+           debiasing_method="BN")
+out
+
 }