Merge pull request #31 from jonathan-taylor/wide_matrix_QP

Wide matrix qp

Author: tibshirani

selectiveInference/R/RcppExports.R

@@ -1,15 +1,19 @@
-# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# This file was generated by Rcpp::compileAttributes
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 solve_QP <- function(Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active) {
-    .Call('_selectiveInference_solve_QP', PACKAGE = 'selectiveInference', Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, max_active)
+    .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)
+    .Call('selectiveInference_update1_', PACKAGE = 'selectiveInference', Q2, w, m, k)
 }
 
 downdate1_ <- function(Q1, R, j0, m, n) {
-    .Call('_selectiveInference_downdate1_', PACKAGE = 'selectiveInference', Q1, R, j0, m, n)
+    .Call('selectiveInference_downdate1_', PACKAGE = 'selectiveInference', Q1, R, j0, m, n)
 }
 

selectiveInference/R/funs.fixed.R

@@ -154,20 +154,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 +268,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 +289,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 +300,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 +334,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 +377,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 +461,7 @@ debiasingRow = function (Sigma,
 
 }
 
+
 ##############################
 
 print.fixedLassoInf <- function(x, tailarea=TRUE, ...) {

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

selectiveInference/src/Rcpp-debias.cpp

@@ -1,5 +1,5 @@
 #include <Rcpp.h>      // need to include the main Rcpp header file 
-#include <debias.h>    // where find_one_row_void is defined
+#include <debias.h>    // where solve_QP, solve_QP_wide are defined
 
 // Below, the gradient should be equal to Sigma * theta + linear_func!!
 // No check is done on this.
@@ -68,3 +68,98 @@ Rcpp::List solve_QP(Rcpp::NumericMatrix Sigma,
 			    Rcpp::Named("max_active_check") = max_active_check));
 
 }
+
+
+// [[Rcpp::export]]
+Rcpp::List solve_QP_wide(Rcpp::NumericMatrix X,
+			 double bound,
+			 int maxiter,
+			 Rcpp::NumericVector theta,
+			 Rcpp::NumericVector linear_func,
+			 Rcpp::NumericVector gradient,
+			 Rcpp::NumericVector X_theta,
+			 Rcpp::IntegerVector ever_active,
+			 Rcpp::IntegerVector nactive,
+			 double kkt_tol,
+			 double objective_tol,
+			 int max_active
+			 ) {
+
+  int ncase = X.nrow(); // number of cases
+  int nfeature = X.ncol(); // number of features
+
+  // Active set
+
+  int icase, ifeature;
+
+  // A vector to keep track of gradient updates
+
+  Rcpp::IntegerVector need_update(nfeature);
+
+  // Extract the diagonal
+  Rcpp::NumericVector nndef_diag(nfeature);
+  double *nndef_diag_p = nndef_diag.begin();
+
+  for (ifeature=0; ifeature<nfeature; ifeature++) {
+    nndef_diag_p[ifeature] = 0;
+    for (icase=0; icase<ncase; icase++) {
+      nndef_diag_p[ifeature] += X(icase, ifeature) * X(icase, ifeature);
+    }
+    nndef_diag_p[ifeature] = nndef_diag_p[ifeature] / ncase;
+  }
+  
+  // Now call our C function
+
+  int iter = solve_wide((double *) X.begin(),
+			(double *) X_theta.begin(),
+			(double *) linear_func.begin(),
+			(double *) nndef_diag.begin(),
+			(double *) gradient.begin(),
+			(int *) need_update.begin(),
+			(int *) ever_active.begin(),
+			(int *) nactive.begin(),
+			ncase,
+			nfeature,
+			bound,
+			(double *) theta.begin(),
+			maxiter,
+			kkt_tol,
+			objective_tol,
+			max_active);
+  
+  // Check whether feasible
+
+  int kkt_check = check_KKT_wide((double *) theta.begin(),
+				 (double *) gradient.begin(),
+				 (double *) X_theta.begin(),
+				 (double *) X.begin(),
+				 (double *) linear_func.begin(),
+				 (int *) need_update.begin(),
+				 nfeature,
+				 ncase,
+				 bound,
+				 kkt_tol);
+
+  int max_active_check = (*(nactive.begin()) >= max_active);
+
+  // Make sure gradient is updated -- essentially a matrix multiply
+
+  update_gradient_wide((double *) gradient.begin(),
+		       (double *) X_theta.begin(),
+		       (double *) X.begin(),
+		       (double *) linear_func.begin(),
+		       (int *) need_update.begin(),
+		       nfeature,
+		       ncase);
+
+  return(Rcpp::List::create(Rcpp::Named("soln") = theta,
+			    Rcpp::Named("gradient") = gradient,
+			    Rcpp::Named("X_theta") = X_theta,
+			    Rcpp::Named("linear_func") = linear_func,
+			    Rcpp::Named("iter") = iter,
+			    Rcpp::Named("kkt_check") = kkt_check,
+			    Rcpp::Named("ever_active") = ever_active,
+			    Rcpp::Named("nactive") = nactive,
+			    Rcpp::Named("max_active_check") = max_active_check));
+
+}