using active sets

Author: jonathan-taylor

selectiveInference/R/funs.fixed.R

@@ -350,6 +350,8 @@ InverseLinftyOneRowC <- function (Sigma, i, mu, maxiter=50) {
           	 Sigma=as.double(Sigma),
 		 Sigma_diag=as.double(diag(Sigma)),
 		 Sigma_theta=as.double(rep(0, p)),
+                 ever_active=as.integer(i),
+		 nactive_ptr=as.integer(1),
 		 nrow=as.integer(p),
    		 bound=as.double(mu),
 		 theta=as.double(theta),

selectiveInference/src/debiasing_matrix.c

@@ -12,6 +12,8 @@
 // Update one coordinate 
 
 double objective(double *Sigma,       /* A covariance matrix: X^TX/n */
+		 int *ever_active,    /* Ever active set   */ 
+		 int *nactive_ptr,    /* Size of ever active set */
 		 int nrow,            /* how many rows in Sigma */
 		 int row,             /* which row: 0-based */
 		 double bound,        /* Lagrange multipler for \ell_1 */
@@ -21,33 +23,106 @@ double objective(double *Sigma,       /* A covariance matrix: X^TX/n */
   double value = 0;
   double *Sigma_ptr = Sigma;
   double *theta_row_ptr, *theta_col_ptr;
+  int *active_row_ptr, *active_col_ptr;
+  int active_row, active_col;
+  int nactive = *nactive_ptr;
 
   theta_row_ptr = theta;
   theta_col_ptr = theta;
 
+  for (irow=0; irow<nactive; irow++) {
+
+    active_row_ptr = ((int *) ever_active + irow);
+    active_row = *active_row_ptr;
+    theta_row_ptr = ((double *) theta + active_row);
+
+    for (icol=0; icol<nactive; icol++) {
+      
+      active_col_ptr = ((int *) ever_active + icol);
+      active_col = *active_col_ptr;
+      theta_col_ptr = ((double *) theta + active_col);
+      
+      fprintf(stderr, "%d %d \n", active_row, active_col);
+
+      Sigma_ptr = ((double *) Sigma + nrow * active_row + active_col);
+
+      value += 0.5 * (*Sigma_ptr) * (*theta_row_ptr) * (*theta_col_ptr);
+    }
+    value = value + bound * fabs((*theta_row_ptr)); // the \ell_1 term
+  }
+
+  theta_row_ptr = ((double *) theta + row);
+  value -= (*theta_row_ptr); // the elementary basis vector term
+
+  return(value);
+}
+
+int is_active(int coord,
+	      int *nactive_ptr,
+	      int *ever_active) {
+  int iactive;
+  int active_var;
+  int nactive = *nactive_ptr;
+  int *ever_active_ptr = ever_active;
+
+  for (iactive=0; iactive<nactive; iactive++) {
+    active_var = (*ever_active_ptr);
+    if (active_var == coord) {
+      return(1);
+    }
+  }
+  return(0);
+}
+
+int check_KKT(double *theta,       /* current theta */
+	      double *Sigma_theta, /* Sigma times theta */
+	      int nrow,            /* how many rows in Sigma */
+	      int row,             /* which row: 0-based */
+	      double bound)        /* Lagrange multipler for \ell_1 */
+{
+  // First check inactive
+
+  int irow;
+  int fail = 0;
+  double tol = 1.e-4;
+  double *theta_ptr, *Sigma_theta_ptr;
+  double gradient;
+
   for (irow=0; irow<nrow; irow++) {
-    double *theta_col_ptr = theta;
-    if (*theta_row_ptr != 0) {
-      for (icol=0; icol<nrow; icol++) {
-        value += 0.5 * (*Sigma_ptr) * (*theta_row_ptr) * (*theta_col_ptr);
-        Sigma_ptr++;
-        theta_col_ptr++;
+    theta_ptr = ((double *) theta + irow);
+    Sigma_theta_ptr = ((double *) Sigma_theta + irow);
+
+    // Compute this coordinate of the gradient
+
+    gradient = *Sigma_theta_ptr;
+    if (row == irow) {
+      gradient -= 1;
+    }
+
+    if (*theta_ptr != 0) { // these coordinates of gradients should be equal to -bound
+      if ((*theta_ptr > 0) &&  (fabs(gradient + bound) > (1. + tol) * bound)) {
+	fail += 1;
+      }
+      else if ((*theta_ptr < 0) && (fabs(gradient - bound) > (1. + tol) * bound)) {
+	fail += 1;
       }
     }
-    if (irow == row) {
-      value -= (*theta_row_ptr); // the elementary basis vector term
+    else {
+      if (fabs(gradient) > (1. + tol) * bound) {
+	fail += 1;
+      }
     }
-    value = value + bound * fabs((*theta_row_ptr)); // the \ell_1 term
-    theta_row_ptr++;
   }
 
-  return(value);
+  return(fail == 0);
 }
 
 
 double update_one_coord(double *Sigma,           /* A covariance matrix: X^TX/n */
                         double *Sigma_diag,      /* Diagonal entries of Sigma */
                         double *Sigma_theta,     /* Sigma times theta */
+			int *ever_active,        /* Ever active set   */ 
+			int *nactive_ptr,        /* Size of ever active set */
 			int nrow,                /* How many rows in Sigma */
 			double bound,            /* feasibility parameter */
 			double *theta,           /* current value */
@@ -67,6 +142,8 @@ double update_one_coord(double *Sigma,           /* A covariance matrix: X^TX/n
   double *quadratic_ptr = ((double *) Sigma_diag + coord);
   double quadratic_term = *quadratic_ptr;
 
+  int *ever_active_ptr;
+
   Sigma_theta_ptr = ((double *) Sigma_theta + coord);
   linear_term = *Sigma_theta_ptr;
 
@@ -97,7 +174,17 @@ double update_one_coord(double *Sigma,           /* A covariance matrix: X^TX/n
     value = -(linear_term - bound) / quadratic_term;
   }
 
-  if (fabs(old_value - value) > 1.e-6 * (fabs(value) + fabs(old_value))) { // Update the linear term
+  // Add to active set if necessary
+
+  if ((value != 0) && (is_active(coord, ever_active, nactive_ptr) == 0)) {
+    ever_active_ptr = ((int *) ever_active + *nactive_ptr);
+    *ever_active_ptr = coord;
+    *nactive_ptr += 1;
+  }
+
+  // Update the linear term
+
+  if (fabs(old_value - value) > 1.e-6 * (fabs(value) + fabs(old_value))) { 
 
     delta = value - old_value;
     Sigma_ptr = ((double *) Sigma + coord * nrow);
@@ -121,6 +208,8 @@ double update_one_coord(double *Sigma,           /* A covariance matrix: X^TX/n
 void find_one_row(double *Sigma,          /* A covariance matrix: X^TX/n */
                   double *Sigma_diag,     /* Diagonal entry of covariance matrix */
                   double *Sigma_theta,    /* Sigma times theta */
+		  int *ever_active,       /* Ever active set   */ 
+		  int *nactive_ptr,       /* Size of ever active set */
                   int *nrow_ptr,          /* How many rows in Sigma */
 		  double *bound_ptr,      /* feasibility parameter */
                   double *theta,          /* current value */
@@ -135,13 +224,17 @@ void find_one_row(double *Sigma,          /* A covariance matrix: X^TX/n */
   double bound = *bound_ptr;
   int nrow = *nrow_ptr;
 
+  fprintf(stderr, "starting now\n");
+
   double old_value = objective(Sigma,
+			       ever_active,
+			       nactive_ptr,
 			       nrow,
 			       row,
 			       bound,
 			       theta);
   double new_value; 
-  double tol=1.e-6;
+  double tol=1.e-5;
 
   for (iter=0; iter<maxiter; iter++) {
 
@@ -150,31 +243,56 @@ void find_one_row(double *Sigma,          /* A covariance matrix: X^TX/n */
     update_one_coord(Sigma,
 		     Sigma_diag,
 		     Sigma_theta,
+		     ever_active,
+		     nactive_ptr,
 		     nrow,
 		     bound,
 		     theta,
 		     row,
 		     row);
 
+    if (check_KKT(theta,
+		  Sigma_theta,
+		  nrow,
+		  row,
+		  bound) == 1) {
+      fprintf(stderr, "ending in first KKT check\n");
+      break;
+    }
+					  
     for (icoord=0; icoord<nrow; icoord++) {
 
       update_one_coord(Sigma,
 		       Sigma_diag,
 		       Sigma_theta,
+		       ever_active,
+		       nactive_ptr,
 		       nrow,
 		       bound,
 		       theta,
 		       row,
 		       icoord);
     }
 
+    if (check_KKT(theta,
+		  Sigma_theta,
+		  nrow,
+		  row,
+		  bound) == 1) {
+      fprintf(stderr, "ending in second KKT check\n");
+      break;
+    }
+					  
     new_value = objective(Sigma,
+			  ever_active,
+			  nactive_ptr,
 			  nrow,
 			  row,
 			  bound,
 			  theta);
 
     if (((old_value - new_value) < tol * fabs(new_value)) && (iter > 0)) {
+      fprintf(stderr, "ending in objective value check\n");
       break;
     }