added a warning about feasibility -- but it didn't print in gist?

Author: jonathan-taylor

selectiveInference/R/funs.fixed.R

@@ -334,13 +334,15 @@ InverseLinfty <- function(sigma, n, e, resol=1.5, mu=NULL, maxiter=50, threshold
   return(M)
 }
 
-InverseLinftyOneRowC <- function ( sigma, i, mu, maxiter=50, threshold=1e-2 ) {
+InverseLinftyOneRowC <- function (Sigma, i, mu, maxiter=50, threshold=1e-2 ) {
 
-         p = nrow(sigma)		
+         p = nrow(Sigma)		
+         basis_vector = rep(0, p)
+	 basis_vector[i] = 1.
          theta = rep(0, p)     
 
 	 val = .C("find_one_row",
-          	 Sigma=as.double(sigma),
+          	 Sigma=as.double(Sigma),
 		 nrow=as.integer(p),
    		 bound=as.double(mu),
 		 theta=as.double(theta),
@@ -350,6 +352,13 @@ InverseLinftyOneRowC <- function ( sigma, i, mu, maxiter=50, threshold=1e-2 ) {
 		 dup=FALSE,
 		 package="selectiveInference")
 
+	# Check feasibility
+
+	if (max(abs(Sigma %*% val$theta - basis_vector)) > 1.01 * mu) {
+           print("Solution for row of M does not seem to be feasible")
+	   warning("Solution for row of M does not seem to be feasible")
+	}
+
 	return(val$theta)
 }
 

selectiveInference/src/debiasing_matrix.c

@@ -3,8 +3,12 @@
 
 // Find an approximate row of \hat{Sigma}^{-1}
 
-// Problem (4) of ????
+// Solves a dual version of problem (4) of https://arxiv.org/pdf/1306.3171.pdf
 
+// Dual problem: \text{min}_{\theta} 1/2 \theta^T \Sigma \theta - e_i^T\theta + \mu \|\theta\|_1
+
+// This is the "negative" of the problem as in https://gist.github.com/jonathan-taylor/07774d209173f8bc4e42aa37712339bf
+// Therefore we don't have to negate the answer to get theta.
 // Update one coordinate 
 
 double update_one_coord(double *Sigma,     /* A covariance matrix: X^TX/n */
@@ -36,19 +40,22 @@ double update_one_coord(double *Sigma,     /* A covariance matrix: X^TX/n */
   }
 
   if (row == coord) {
-    linear_term += 1;
+    linear_term -= 1;
   }
 
   // Now soft-threshold the coord entry of theta 
 
   // Objective is t \mapsto q/2 * t^2 + l * t + bound |t|
   // with q=quadratic_term and l=linear_term
 
+  // With a negative linear term, solution should be
+  // positive
+
   if (linear_term < -bound) {
-    value = - (-linear_term - bound) / quadratic_term;
+    value = (-linear_term - bound) / quadratic_term;
   }
   else if (linear_term > bound) {
-    value = (linear_term - bound) / quadratic_term;
+    value = -(linear_term - bound) / quadratic_term;
   }
 
   theta_ptr = ((double *) theta + coord);