Merge remote-tracking branch 'upstream'

Author: jonathan-taylor

selectiveInference-currentCRAN/DESCRIPTION

@@ -1,12 +1,12 @@
 Package: selectiveInference
 Type: Package
 Title: Tools for Post-Selection Inference
-Version: 1.2.3
+Version: 1.2.4
 Date: 2017-09-18
 Author: Ryan Tibshirani, Rob Tibshirani, Jonathan Taylor,
     Joshua Loftus, Stephen Reid
 Maintainer: Rob Tibshirani <[email protected]>
-Depends: glmnet, intervals, survival
+Depends: glmnet, intervals, survival, R (>= 3.4.0)
 Suggests: Rmpfr
 Description: New tools for post-selection inference, for use
     with forward stepwise regression, least angle regression, the

selectiveInference/R/RcppExports.R

@@ -5,6 +5,10 @@ solve_QP <- function(Sigma, bound, maxiter, theta, linear_func, gradient, ever_a
     .Call('_selectiveInference_solve_QP', PACKAGE = 'selectiveInference', Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active)
 }
 
+solve_QP_wide <- function(X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active) {
+    .Call('_selectiveInference_solve_QP_wide', PACKAGE = 'selectiveInference', X, bound, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, max_active)
+}
+
 update1_ <- function(Q2, w, m, k) {
     .Call('_selectiveInference_update1_', PACKAGE = 'selectiveInference', Q2, w, m, k)
 }

selectiveInference/R/funs.fixed.R

@@ -82,10 +82,11 @@ fixedLassoInf <- function(x, y, beta,
 
     tol.coef = tol.beta * sqrt(n^2 / colSums(x^2))
     # print(tol.coef)
-    vars = which(abs(beta) > tol.coef)
+      vars = which(abs(beta) > tol.coef)
+   #    vars = abs(beta) > tol.coef
     # print(beta)
     # print(vars)
-    if(length(vars)==0){
+    if(sum(vars)==0){
       cat("Empty model",fill=T)
       return()
     }
@@ -96,8 +97,10 @@ fixedLassoInf <- function(x, y, beta,
                     "'thresh' parameter, for a more accurate convergence."))
 
     # Get lasso polyhedral region, of form Gy >= u
-    if (type == 'full' & p > n) out = fixedLassoPoly(x,y,lambda,beta,vars,inactive=TRUE)
-    else out = fixedLassoPoly(x,y,lambda,beta,vars)
+logical.vars=rep(FALSE,p)
+logical.vars[vars]=TRUE
+    if (type == 'full') out = fixedLassoPoly(x,y,lambda,beta,logical.vars,inactive=TRUE)
+    else out = fixedLassoPoly(x,y,lambda,beta,logical.vars)
     A = out$A
     b = out$b
 
@@ -127,7 +130,7 @@ fixedLassoInf <- function(x, y, beta,
     # add additional targets for inference if provided
     if (!is.null(add.targets)) vars = sort(unique(c(vars,add.targets,recursive=T)))
 
-    k = length(vars)
+      k = length(vars)
     pv = vlo = vup = numeric(k)
     vmat = matrix(0,k,n)
     ci = tailarea = matrix(0,k,2)
@@ -154,20 +157,24 @@ fixedLassoInf <- function(x, y, beta,
 
       # Reorder so that active set S is first
       Xordered = Xint[,c(S,notS,recursive=T)]
+      hsigmaS = 1/n*(t(XS)%*%XS) # hsigma[S,S]
+      hsigmaSinv = solve(hsigmaS) # pinv(hsigmaS)
 
-      hsigma <- 1/n*(t(Xordered)%*%Xordered)
-      hsigmaS <- 1/n*(t(XS)%*%XS) # hsigma[S,S]
-      hsigmaSinv <- solve(hsigmaS) # pinv(hsigmaS)
+      FS = rbind(diag(length(S)),matrix(0,pp-length(S),length(S)))
+      GS = cbind(diag(length(S)),matrix(0,length(S),pp-length(S)))
 
-      # Approximate inverse covariance matrix for when (n < p) from lasso_Inference.R
+      is_wide = n < (2 * p) # somewhat arbitrary decision -- it is really for when we don't want to form with pxp matrices
 
-      htheta = debiasingMatrix(hsigma, n, 1:length(S), verbose=FALSE, max_try=linesearch.try, warn_kkt=TRUE)
+      # Approximate inverse covariance matrix for when (n < p) from lasso_Inference.R
+      if (!is_wide) {
+           hsigma = 1/n*(t(Xordered)%*%Xordered)
+           htheta = debiasingMatrix(hsigma, is_wide, n, 1:length(S), verbose=FALSE, max_try=linesearch.try, warn_kkt=TRUE)
+           ithetasigma = (GS-(htheta%*%hsigma))
+      } else {
+           htheta = debiasingMatrix(Xordered, is_wide, n, 1:length(S), verbose=FALSE, max_try=linesearch.try, warn_kkt=TRUE)
+           ithetasigma = (GS-((htheta%*%t(Xordered)) %*% Xordered)/n)
+      }
 
-      FS = rbind(diag(length(S)),matrix(0,pp-length(S),length(S)))
-      GS = cbind(diag(length(S)),matrix(0,length(S),pp-length(S)))
-      ithetasigma = (GS-(htheta%*%hsigma))
-      # ithetasigma = (diag(pp) - (htheta%*%hsigma))
-      
       M <- (((htheta%*%t(Xordered))+ithetasigma%*%FS%*%hsigmaSinv%*%t(XS))/n)
       # vector which is offset for testing debiased beta's
       null_value <- (((ithetasigma%*%FS%*%hsigmaSinv)%*%sign(hbetaS))*lambda/n)
@@ -264,10 +271,11 @@ fixedLassoPoly =
 ## Approximates inverse covariance matrix theta
 ## using coordinate descent 
 
-debiasingMatrix = function(Sigma, 
+debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n depending on is_wide
+                           is_wide,
                            nsample, 
                            rows, 
-		           verbose=FALSE, 
+ 		           verbose=FALSE, 
 		           mu=NULL,             # starting value of mu
    			   linesearch=TRUE,     # do a linesearch?
    		           scaling_factor=1.5,  # multiplicative factor for linesearch
@@ -284,7 +292,7 @@ debiasingMatrix = function(Sigma,
      max_active = max(50, 0.3 * nsample)
   } 
 
-  p = nrow(Sigma);
+  p = ncol(Xinfo);
   M = matrix(0, length(rows), p);
 
   if (is.null(mu)) {
@@ -295,19 +303,19 @@ debiasingMatrix = function(Sigma,
   xp = round(p/10);
   idx = 1;
   for (row in rows) {
-
     if ((idx %% xp)==0){
       xperc = xperc+10;
       if (verbose) {
         print(paste(xperc,"% done",sep="")); }
     }
 
-    output = debiasingRow(Sigma,
+    output = debiasingRow(Xinfo,               # could be X or t(X) %*% X / n depending on is_wide
+                          is_wide,
                           row,
                           mu,
                           linesearch=linesearch,
                           scaling_factor=scaling_factor,
-			  max_active=max_active,
+                          max_active=max_active,
 			  max_try=max_try,
 			  warn_kkt=FALSE,
 			  max_iter=max_iter,
@@ -329,31 +337,32 @@ debiasingMatrix = function(Sigma,
   return(M)
 }
 
-# Find one row of the debiasing matrix
+# Find one row of the debiasing matrix -- assuming X^TX/n is not too large -- i.e. X is tall
 
-debiasingRow = function (Sigma, 
+debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n depending on is_wide
+                         is_wide, 
                          row, 
                          mu, 
-		         linesearch=TRUE,     # do a linesearch?
+	                 linesearch=TRUE,     # do a linesearch?
 		         scaling_factor=1.2,  # multiplicative factor for linesearch
-			 max_active=NULL,     # how big can active set get?
+		         max_active=NULL,     # how big can active set get?
 			 max_try=10,          # how many steps in linesearch?
 			 warn_kkt=FALSE,      # warn if KKT does not seem to be satisfied?
 			 max_iter=100,        # how many iterations for each optimization problem
                          kkt_tol=1.e-4,       # tolerance for the KKT conditions
 			 objective_tol=1.e-8  # tolerance for relative decrease in objective
                          ) {
 
-  p = nrow(Sigma)
+  p = ncol(Xinfo)
 
   if (is.null(max_active)) {
-      max_active = nrow(Sigma)
+      max_active = min(nrow(Xinfo), ncol(Xinfo))
   }
 
   # Initialize variables 
 
   soln = rep(0, p)
-
+  Xsoln = rep(0, n)
   ever_active = rep(0, p)
   ever_active[1] = row      # 1-based
   ever_active = as.integer(ever_active)
@@ -371,17 +380,33 @@ debiasingRow = function (Sigma,
 
   while (counter_idx < max_try) {
 
-      result = solve_QP(Sigma, 
-                        mu, 
-                        max_iter, 
-                        soln, 
-                        linear_func, 
-                        gradient, 
-                        ever_active, 
-                        nactive, 
-                        kkt_tol, 
-                        objective_tol, 
-                        max_active) 
+      if (!is_wide) {
+          result = solve_QP(Xinfo, # this is non-neg-def matrix
+                            mu, 
+                            max_iter, 
+                            soln, 
+                            linear_func, 
+                            gradient, 
+                            ever_active, 
+                            nactive, 
+                            kkt_tol, 
+                            objective_tol, 
+                            max_active) 
+      } else {
+          result = solve_QP_wide(Xinfo, # this is a design matrix
+                                 mu, 
+                                 max_iter, 
+                                 soln, 
+                                 linear_func, 
+                                 gradient, 
+                                 Xsoln,
+                                 ever_active, 
+                                 nactive, 
+                                 kkt_tol, 
+                                 objective_tol, 
+                                 max_active) 
+
+      }
 
       iter = result$iter
 
@@ -439,6 +464,7 @@ debiasingRow = function (Sigma,
 
 }
 
+
 ##############################
 
 print.fixedLassoInf <- function(x, tailarea=TRUE, ...) {
@@ -478,4 +504,3 @@ print.fixedLassoInf <- function(x, tailarea=TRUE, ...) {
  # lambda = 2*mean(apply(t(x)%*%eps,2,max))
  # return(lambda)
 #}
-

selectiveInference/R/funs.inf.R

@@ -209,7 +209,7 @@ TG.limits = function(Z, A, b, eta, Sigma=NULL) {
     target_estimate = sum(as.numeric(eta) * as.numeric(Z))
 
     if (max(A %*% as.numeric(Z) - b) > 0) {
-        warning('Contsraint not satisfied. A %*% Z should be elementwise less than or equal to b')
+        warning('Constraint not satisfied. A %*% Z should be elementwise less than or equal to b')
     }
 
     if (is.null(Sigma)) {

selectiveInference/man/debiasingMatrix.Rd

@@ -11,7 +11,8 @@ Newton step from some consistent estimator (such as the LASSO)
 to find a debiased solution.
 }
 \usage{
-debiasingMatrix(Sigma, 
+debiasingMatrix(Xinfo, 
+                is_wide,			
                 nsample, 
                 rows, 
 		verbose=FALSE, 
@@ -26,8 +27,14 @@ debiasingMatrix(Sigma,
 		objective_tol=1.e-8)
 }
 \arguments{
-\item{Sigma}{
-A symmetric non-negative definite matrix, often a cross-covariance matrix.
+\item{Xinfo}{
+Either a non-negative definite matrix S=t(X) %*% X / n or X itself. If 
+is_wide is TRUE, then Xinfo should be X, otherwise it should be S.
+}
+\item{is_wide}{
+Are we solving for rows of the debiasing matrix assuming it is 
+a wide matrix so that Xinfo=X and the non-negative definite
+matrix of interest is t(X) %*% X / nrow(X).
 }
 \item{nsample}{
 Number of samples used in forming the cross-covariance matrix.
@@ -101,8 +108,9 @@ set.seed(10)
 n = 50
 p = 100
 X = matrix(rnorm(n * p), n, p)
-S = t(X) \%*\% X / n
-M = debiasingMatrix(S, n, c(1,3,5))
-
+S = t(X) %*% X / n
+M = debiasingMatrix(S, FALSE, n, c(1,3,5))
+M2 = debiasingMatrix(X, TRUE, n, c(1,3,5))
+max(M - M2)
 }
 

selectiveInference/src/Makevars

@@ -2,7 +2,7 @@ PKG_CFLAGS= -I.
 PKG_CPPFLAGS= -I.
 PKG_LIBS=-L. 
 
-$(SHLIB): Rcpp Rcpp-matrixcomps.o Rcpp-debias.o RcppExports.o quadratic_program.o
+$(SHLIB): Rcpp Rcpp-matrixcomps.o Rcpp-debias.o RcppExports.o quadratic_program.o quadratic_program_wide.o
 
 clean:
 	rm -f *o