R code now produces almost identical answer when assuming wide or not

jonathan-taylor · jonathan-taylor · commit 76f3e7e2d3d2 · 2017-09-28T00:14:02.000-07:00
diff --git a/selectiveInference/R/funs.fixed.R b/selectiveInference/R/funs.fixed.R
@@ -300,7 +300,6 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
   xp = round(p/10);
   idx = 1;
   for (row in rows) {
-
     if ((idx %% xp)==0){
       xperc = xperc+10;
       if (verbose) {
@@ -354,13 +353,13 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
   p = ncol(Xinfo)
 
   if (is.null(max_active)) {
-      max_active = nrow(Xinfo)
+      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)
@@ -379,8 +378,7 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
   while (counter_idx < max_try) {
 
       if (!is_wide) {
-          Sigma = Xinfo
-          result = solve_QP(Sigma, 
+          result = solve_QP(Xinfo, # this is non-neg-def matrix
                             mu, 
                             max_iter, 
                             soln, 
@@ -392,10 +390,7 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
                             objective_tol, 
                             max_active) 
       } else {
-          X = Xinfo
-	  n = nrow(X)
-          Xsoln = rep(0, n)
-          result = solve_QP_wide(X,
+          result = solve_QP_wide(Xinfo, # this is a design matrix
                                  mu, 
                                  max_iter, 
                                  soln, 
@@ -409,6 +404,7 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
                                  max_active) 
 
       }
+
       iter = result$iter
 
       # Logic for whether we should continue the line search
diff --git a/selectiveInference/man/debiasingMatrix.Rd b/selectiveInference/man/debiasingMatrix.Rd
@@ -108,8 +108,9 @@ set.seed(10)
 n = 50
 p = 100
 X = matrix(rnorm(n * p), n, p)
-S = t(X) \%*\% X / n
+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)
 }
  
diff --git a/selectiveInference/src/quadratic_program.c b/selectiveInference/src/quadratic_program.c
@@ -17,7 +17,7 @@ double objective_qp(double *nndef_ptr,       /* A non-negative definite matrix *
 		    int *nactive_ptr,        /* Size of ever active set */
 		    int nrow,                /* how many rows in nndef */
 		    double bound,            /* Lagrange multipler for \ell_1 */
-		    double *theta)           /* current value */
+		    double *theta_ptr)           /* current value */
 {
   int irow, icol;
   double value = 0;
@@ -28,20 +28,20 @@ double objective_qp(double *nndef_ptr,       /* A non-negative definite matrix *
   int active_row, active_col;
   int nactive = *nactive_ptr;
 
-  theta_row_ptr = theta;
-  theta_col_ptr = theta;
+  theta_row_ptr = theta_ptr;
+  theta_col_ptr = theta_ptr;
 
   for (irow=0; irow<nactive; irow++) {
 
     active_row_ptr = ((int *) ever_active_ptr + irow);
     active_row = *active_row_ptr - 1;          // Ever-active is 1-based
-    theta_row_ptr = ((double *) theta + active_row);
+    theta_row_ptr = ((double *) theta_ptr + active_row);
 
     for (icol=0; icol<nactive; icol++) {
       
       active_col_ptr = ((int *) ever_active_ptr + icol);
       active_col = *active_col_ptr - 1;          // Ever-active is 1-based
-      theta_col_ptr = ((double *) theta + active_col);
+      theta_col_ptr = ((double *) theta_ptr + active_col);
 
       nndef_ptr_tmp = ((double *) nndef_ptr + nrow * active_col + active_row); // Matrices are column-major order
 
@@ -90,7 +90,7 @@ int update_ever_active_qp(int coord,
   return(0);
 }
 
-int check_KKT_qp(double *theta,        /* current theta */
+int check_KKT_qp(double *theta_ptr,        /* current theta */
 		 double *gradient_ptr, /* nndef times theta + linear_func */
 		 int nrow,             /* how many rows in nndef */
 		 double bound,         /* Lagrange multipler for \ell_1 */
@@ -99,22 +99,22 @@ int check_KKT_qp(double *theta,        /* current theta */
   // First check inactive
 
   int irow;
-  double *theta_ptr, *gradient_ptr_tmp;
+  double *theta_ptr_tmp, *gradient_ptr_tmp;
   double gradient;
 
   for (irow=0; irow<nrow; irow++) {
-    theta_ptr = ((double *) theta + irow);
+    theta_ptr_tmp = ((double *) theta_ptr + irow);
     gradient_ptr_tmp = ((double *) gradient_ptr + irow);
 
     // Compute this coordinate of the gradient
 
     gradient = *gradient_ptr_tmp;
 
-    if (*theta_ptr != 0) { // these coordinates of gradients should be equal to -bound
-      if ((*theta_ptr > 0) &&  (fabs(gradient + bound) > tol * bound)) {
+    if (*theta_ptr_tmp != 0) { // these coordinates of gradients should be equal to -bound
+      if ((*theta_ptr_tmp > 0) &&  (fabs(gradient + bound) > tol * bound)) {
 	return(0);
       }
-      else if ((*theta_ptr < 0) && (fabs(gradient - bound) > tol * bound)) {
+      else if ((*theta_ptr_tmp < 0) && (fabs(gradient - bound) > tol * bound)) {
 	return(0);
       }
     }
@@ -128,6 +128,57 @@ int check_KKT_qp(double *theta,        /* current theta */
   return(1);
 }
 
+int check_KKT_qp_active(int *ever_active_ptr,           /* Ever active set: 0-based */ 
+		        int *nactive_ptr,               /* Size of ever active set */
+			double *theta_ptr,        /* current theta */
+			double *gradient_ptr, /* nndef times theta + linear_func */
+			int nrow,             /* how many rows in nndef */
+			double bound,         /* Lagrange multipler for \ell_1 */
+			double tol)           /* precision for checking KKT conditions */        
+{
+  // First check inactive
+
+  int iactive;
+  double *theta_ptr_tmp;
+  double gradient;
+  double *gradient_ptr_tmp;
+  int nactive = *nactive_ptr;
+  int active_feature;
+  int *active_feature_ptr;
+
+  for (iactive=0; iactive<nactive; iactive++) {
+
+    active_feature_ptr = ((int *) ever_active_ptr + iactive);
+    active_feature = *active_feature_ptr - 1;          // Ever-active is 1-based
+    theta_ptr_tmp = ((double *) theta_ptr + active_feature);
+
+    gradient_ptr_tmp = ((double *) gradient_ptr + active_feature);
+
+    // Compute this coordinate of the gradient
+
+    gradient = *gradient_ptr_tmp;
+
+    if (*theta_ptr_tmp != 0) { // these coordinates of gradients should be equal to -bound
+
+      if ((*theta_ptr_tmp > 0) &&  (fabs(gradient + bound) > tol * bound)) {
+	return(0);
+      }
+      else if ((*theta_ptr_tmp < 0) && (fabs(gradient - bound) > tol * bound)) {
+	return(0);
+      }
+
+    }
+    else {
+      if (fabs(gradient) > (1. + tol) * bound) {
+	return(0);
+      }
+    }
+  }
+
+  return(1);
+}
+
+
 double update_one_coord_qp(double *nndef_ptr,           /* A non-negative definite matrix */
 			   double *linear_func_ptr,     /* Linear term in objective */
 			   double *nndef_diag_ptr,      /* Diagonal of nndef */
@@ -136,7 +187,7 @@ double update_one_coord_qp(double *nndef_ptr,           /* A non-negative defini
 			   int *nactive_ptr,            /* Size of ever active set */
 			   int nrow,                    /* How many rows in nndef */
 			   double bound,                /* feasibility parameter */
-			   double *theta,               /* current value */
+			   double *theta_ptr,               /* current value */
 			   int coord,                   /* which coordinate to update: 0-based */
 			   int is_active)               /* Is this coord in ever_active */     
 {
@@ -147,7 +198,7 @@ double update_one_coord_qp(double *nndef_ptr,           /* A non-negative defini
   double old_value;
   double *nndef_ptr_tmp;
   double *gradient_ptr_tmp;
-  double *theta_ptr;
+  double *theta_ptr_tmp;
   int icol = 0;
 
   double *quadratic_ptr = ((double *) nndef_diag_ptr + coord);
@@ -156,8 +207,8 @@ double update_one_coord_qp(double *nndef_ptr,           /* A non-negative defini
   gradient_ptr_tmp = ((double *) gradient_ptr + coord);
   linear_term = *gradient_ptr_tmp;
 
-  theta_ptr = ((double *) theta + coord);
-  old_value = *theta_ptr;
+  theta_ptr_tmp = ((double *) theta_ptr + coord);
+  old_value = *theta_ptr_tmp;
 
   // The coord entry of gradient_ptr term has a diagonal term in it:
   // nndef[coord, coord] * theta[coord]
@@ -200,8 +251,8 @@ double update_one_coord_qp(double *nndef_ptr,           /* A non-negative defini
       nndef_ptr_tmp += 1;
     }
 
-    theta_ptr = ((double *) theta + coord);
-    *theta_ptr = value;
+    theta_ptr_tmp = ((double *) theta_ptr + coord);
+    *theta_ptr_tmp = value;
 
   }
 
@@ -230,6 +281,8 @@ int solve_qp(double *nndef_ptr,          /* A non-negative definite matrix */
   int *active_ptr;
 
   int check_objective = 1;
+  int iter_active;
+  int niter_active=5;
 
   double old_value, new_value; 
 
@@ -248,23 +301,38 @@ int solve_qp(double *nndef_ptr,          /* A non-negative definite matrix */
 
   for (iter=0; iter<maxiter; iter++) {
 
-    // Update the active variables first
-
-    active_ptr = (int *) ever_active_ptr;
-
-    for (iactive=0; iactive < *nactive_ptr; iactive++) {
-      update_one_coord_qp(nndef_ptr,
-			  linear_func_ptr,
-			  nndef_diag_ptr,
-			  gradient_ptr,
-			  ever_active_ptr,
-			  nactive_ptr,
-			  nrow,
-			  bound,
-			  theta,
-			  *active_ptr - 1,   // Ever-active is 1-based
-			  1);
-      active_ptr++;
+    // Update the active variables first -- do this niter_active times
+
+    for (iter_active=0; iter_active<niter_active; iter_active++) { 
+
+        active_ptr = (int *) ever_active_ptr;
+        for (iactive=0; iactive < *nactive_ptr; iactive++) {
+
+	  update_one_coord_qp(nndef_ptr,
+			      linear_func_ptr,
+			      nndef_diag_ptr,
+			      gradient_ptr,
+			      ever_active_ptr,
+			      nactive_ptr,
+			      nrow,
+			      bound,
+			      theta,
+			      *active_ptr - 1,   // Ever-active is 1-based
+			      1);
+	  active_ptr++;
+	}
+
+	// Check KKT of active subproblem
+
+	if (check_KKT_qp_active(ever_active_ptr,
+				nactive_ptr,
+				theta,
+				gradient_ptr,
+				nrow,
+				bound,
+				kkt_tol) == 1) {
+	  break;
+	}
     }
 
     // Check KKT
diff --git a/selectiveInference/src/quadratic_program_wide.c b/selectiveInference/src/quadratic_program_wide.c
@@ -1,6 +1,6 @@
 #include <math.h> // for fabs
 
-// Find an approximate row of \hat{Sigma}^{-1}
+// Find an approximate row of \hat{nndef}^{-1}
 
 // Solves a dual version of problem (4) of https://arxiv.org/pdf/1306.3171.pdf
 
@@ -222,7 +222,6 @@ int check_KKT_wide_active(int *ever_active_ptr,           /* Ever active set: 0-
   int iactive;
   double *theta_ptr_tmp;
   double gradient;
-  int ever_active_ptr_tmp;
   int nactive = *nactive_ptr;
   int active_feature;
   int *active_feature_ptr;
@@ -259,16 +258,16 @@ int check_KKT_wide_active(int *ever_active_ptr,           /* Ever active set: 0-
 
 double update_one_coord_wide(double *X_ptr,               /* A design matrix*/
 			     double *linear_func_ptr,     /* Linear term in objective */
-			     double *nndef_diag_ptr,          /* Diagonal entries of Sigma */
-			     double *gradient_ptr,     /* X^TX/ncase times theta + linear_func*/
+			     double *nndef_diag_ptr,      /* Diagonal of nndef */
+			     double *gradient_ptr,        /* X^TX/ncase times theta + linear_func*/
 			     int *ever_active_ptr,        /* Ever active set: 1-based */ 
 			     int *nactive_ptr,            /* Size of ever active set */
 			     double *X_theta_ptr,         /* X\theta -- fitted values */
 			     int *need_update_ptr,        /* Whether a gradient coordinate needs update or not */
-			     int ncase,                    /* How many rows in X */
-			     int nfeature,                    /* How many rows in X */
+			     int ncase,                   /* How many rows in X */
+			     int nfeature,                /* How many rows in X */
 			     double bound,                /* feasibility parameter */
-			     double *theta_ptr,               /* current value */
+			     double *theta_ptr,           /* current value */
 			     int coord,                   /* which coordinate to update: 0-based */
 			     int is_active)               /* Is this coord in ever_active */     
 {
