without linesearch we know agree with Adel's code at fixed mu

Author: jonathan-taylor

selectiveInference/DESCRIPTION

@@ -9,7 +9,7 @@ Maintainer: Rob Tibshirani <[email protected]>
 Depends:
     glmnet,
     intervals,
-    survival
+    survival,
 Suggests:
     Rmpfr
 Description: New tools for post-selection inference, for use with forward

selectiveInference/NAMESPACE

@@ -44,4 +44,4 @@ importFrom("stats", dnorm, lsfit, pexp, pnorm, predict,
 importFrom("stats", "coef", "df", "lm", "pf")
 importFrom("stats", "glm", "residuals", "vcov")
 importFrom("Rcpp", "sourceCpp")
-importFrom("distr", "Norm", "DExp")
+

selectiveInference/R/funs.fixed.R

@@ -319,7 +319,7 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
                            nsample, 
                            rows, 
  		           verbose=FALSE, 
-		           mu=NULL,             # starting value of mu
+		           bound=NULL,             # starting value of bound
    			   linesearch=TRUE,     # do a linesearch?
    		           scaling_factor=1.5,  # multiplicative factor for linesearch
 			   max_active=NULL,     # how big can active set get?
@@ -342,8 +342,8 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
   p = ncol(Xinfo);
   M = matrix(0, length(rows), p);
 
-  if (is.null(mu)) {
-      mu = (1/sqrt(nsample)) * qnorm(1-(0.1/(p^2)))
+  if (is.null(bound)) {
+      bound = (1/sqrt(nsample)) * qnorm(1-(0.1/(p^2)))
   }
 
   xperc = 0;
@@ -359,7 +359,7 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
     output = debiasingRow(Xinfo,               # could be X or t(X) %*% X / n depending on is_wide
                           is_wide,
                           row,
-                          mu,
+                          bound,
                           linesearch=linesearch,
                           scaling_factor=scaling_factor,
                           max_active=max_active,
@@ -393,7 +393,7 @@ debiasingMatrix = function(Xinfo,               # could be X or t(X) %*% X / n d
 debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n depending on is_wide
                          is_wide, 
                          row, 
-                         mu, 
+                         bound, 
 	                 linesearch=TRUE,     # do a linesearch?
 		         scaling_factor=1.5,  # multiplicative factor for linesearch
 		         max_active=NULL,     # how big can active set get?
@@ -414,9 +414,11 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
       max_active = min(nrow(Xinfo), ncol(Xinfo))
   }
 
+   
   # Initialize variables 
 
   soln = rep(0, p)
+  soln = as.numeric(soln)
   ever_active = rep(0, p)
   ever_active[1] = row      # 1-based
   ever_active = as.integer(ever_active)
@@ -432,13 +434,16 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
 
   last_output = NULL
 
-  Xsoln = rep(0, n)
+  if (is_wide) {
+     n = nrow(Xinfo)
+     Xsoln = as.numeric(rep(0, n))
+  }
 
   while (counter_idx < max_try) {
 
       if (!is_wide) {
           result = solve_QP(Xinfo, # this is non-neg-def matrix
-                            mu, 
+                            bound, 
                             max_iter, 
                             soln, 
                             linear_func, 
@@ -453,9 +458,9 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
 			    kkt_stop,
 			    parameter_stop)
       } else {
-          result = solve_QP_wide(Xinfo, # this is a design matrix
-                                 rep(mu, p),  # vector of Lagrange multipliers
-				 0,           # ridge_term 
+          result = solve_QP_wide(Xinfo,                      # this is a design matrix
+                                 as.numeric(rep(bound, p)),  # vector of Lagrange multipliers
+				 0,                          # ridge_term 
                                  max_iter, 
                                  soln, 
                                  linear_func, 
@@ -493,13 +498,13 @@ debiasingRow = function (Xinfo,               # could be X or t(X) %*% X / n dep
          if ((iter < (max_iter+1)) && (counter_idx > 1)) { 
            break;      # we've found a feasible point and solved the problem            
          }
-         mu = mu * scaling_factor;
+         bound = bound * scaling_factor;
       } else {         # trying to drop the bound parameter further
          if ((iter == (max_iter + 1)) && (counter_idx > 1)) {
             result = last_output; # problem seems infeasible because we didn't solve it
    	    break;                # so we revert to previously found solution
          }
-         mu = mu / scaling_factor;
+         bound = bound / scaling_factor;
       }
 
       # If the active set has grown to a certain size

selectiveInference/R/funs.randomized.R

@@ -3,19 +3,19 @@
 #
 # min 1/2 || y - \beta_0 - X \beta ||_2^2 + \lambda || \beta ||_1 - \omega^T\beta + \frac{\epsilon}{2} \|\beta\|^2_2
 
-fit_randomized_lasso = function(X, 
-                                y, 
-                                lam, 
-                                noise_scale, 
-                                ridge_term, 
-                                noise_type=c('gaussian', 'laplace'),
-                                max_iter=100,        # how many iterations for each optimization problem
-                                kkt_tol=1.e-4,       # tolerance for the KKT conditions
-                                parameter_tol=1.e-8, # tolerance for relative convergence of parameter
-                                objective_tol=1.e-8, # tolerance for relative decrease in objective
-                                objective_stop=FALSE,
-                                kkt_stop=TRUE,
-                                param_stop=TRUE)
+randomizedLASSO = function(X, 
+                           y, 
+                           lam, 
+                           noise_scale, 
+                           ridge_term, 
+                           noise_type=c('gaussian', 'laplace'),
+                           max_iter=100,        # how many iterations for each optimization problem
+                           kkt_tol=1.e-4,       # tolerance for the KKT conditions
+                           parameter_tol=1.e-8, # tolerance for relative convergence of parameter
+                           objective_tol=1.e-8, # tolerance for relative decrease in objective
+                           objective_stop=FALSE,
+                           kkt_stop=TRUE,
+                           param_stop=TRUE)
 {
 
     n = nrow(X); p = ncol(X)
@@ -24,12 +24,11 @@ fit_randomized_lasso = function(X,
 
     if (noise_scale > 0) {
         if (noise_type == 'gaussian') {
-            D = Norm(mean=0, sd=noise_scale)
+            perturb_ = rnorm(p) * noise_scale
         }
         else if (noise_type == 'laplace') {
-            D = DExp(rate = 1 / noise_scale) # D is a Laplace distribution with rate = 1.
+            perturb_ = rexp(p) * (2 * rbinom(p, 1, 0.5) - 1) * noise_scale
         }
-        perturb_ = distr::r(D)(p)
     } else {
         perturb_ = rep(0, p)
     }

tests/test_debiasing.R

@@ -2,9 +2,9 @@ library(selectiveInference)
 
 
 ## Approximates inverse covariance matrix theta
-InverseLinfty <- function(sigma, n, resol=1.5, mu=NULL, maxiter=50, threshold=1e-10, verbose = TRUE) {
+InverseLinfty <- function(sigma, n, resol=1.5, bound=NULL, maxiter=50, threshold=1e-10, verbose = TRUE) {
   isgiven <- 1;
-  if (is.null(mu)){
+  if (is.null(bound)){
     isgiven <- 0;
   }
 
@@ -19,43 +19,43 @@ InverseLinfty <- function(sigma, n, resol=1.5, mu=NULL, maxiter=50, threshold=1e
         print(paste(xperc,"% done",sep="")); }
     }
     if (isgiven==0){
-      mu <- (1/sqrt(n)) * qnorm(1-(0.1/(p^2)));
+      bound <- (1/sqrt(n)) * qnorm(1-(0.1/(p^2)));
     }
-    mu.stop <- 0;
+    bound.stop <- 0;
     try.no <- 1;
     incr <- 0;
-    while ((mu.stop != 1)&&(try.no<10)){
+    while ((bound.stop != 1)&&(try.no<10)){
       last.beta <- beta
-      output <- InverseLinftyOneRow(sigma, i, mu, maxiter=maxiter, threshold=threshold)
+      output <- InverseLinftyOneRow(sigma, i, bound, maxiter=maxiter, threshold=threshold)
       beta <- output$optsol
       iter <- output$iter
       if (isgiven==1){
-        mu.stop <- 1
+        bound.stop <- 1
       }
       else{
         if (try.no==1){
           if (iter == (maxiter+1)){
             incr <- 1;
-            mu <- mu*resol;
+            bound <- bound*resol;
           } else {
             incr <- 0;
-            mu <- mu/resol;
+            bound <- bound/resol;
           }
         }
         if (try.no > 1){
           if ((incr == 1)&&(iter == (maxiter+1))){
-            mu <- mu*resol;
+            bound <- bound*resol;
           }
           if ((incr == 1)&&(iter < (maxiter+1))){
-            mu.stop <- 1;
+            bound.stop <- 1;
           }
           if ((incr == 0)&&(iter < (maxiter+1))){
-            mu <- mu/resol;
+            bound <- bound/resol;
           }
           if ((incr == 0)&&(iter == (maxiter+1))){
-            mu <- mu*resol;
+            bound <- bound*resol;
             beta <- last.beta;
-            mu.stop <- 1;
+            bound.stop <- 1;
           }
         }
       }
@@ -66,14 +66,14 @@ InverseLinfty <- function(sigma, n, resol=1.5, mu=NULL, maxiter=50, threshold=1e
   return(M)
 }
 
-InverseLinftyOneRow <- function ( sigma, i, mu, maxiter=50, threshold=1e-10) {
+InverseLinftyOneRow <- function ( sigma, i, bound, maxiter=50, threshold=1e-10) {
   p <- nrow(sigma);
   rho <- max(abs(sigma[i,-i])) / sigma[i,i];
-  mu0 <- rho/(1+rho);
+  bound0 <- rho/(1+rho);
   beta <- rep(0,p);
 
-  #if (mu >= mu0){
-  #  beta[i] <- (1-mu0)/sigma[i,i];
+  #if (bound >= bound0){
+  #  beta[i] <- (1-bound0)/sigma[i,i];
   #  returnlist <- list("optsol" = beta, "iter" = 0);
   #  return(returnlist);
   #}
@@ -82,7 +82,7 @@ InverseLinftyOneRow <- function ( sigma, i, mu, maxiter=50, threshold=1e-10) {
   last.norm2 <- 1;
   iter <- 1;
   iter.old <- 1;
-  beta[i] <- (1-mu0)/sigma[i,i];
+  beta[i] <- (1-bound0)/sigma[i,i];
   beta.old <- beta;
   sigma.tilde <- sigma;
   diag(sigma.tilde) <- 0;
@@ -95,7 +95,7 @@ InverseLinftyOneRow <- function ( sigma, i, mu, maxiter=50, threshold=1e-10) {
       v <- vs[j];
       if (j==i)
         v <- v+1;
-      beta[j] <- SoftThreshold(v,mu)/sigma[j,j];
+      beta[j] <- SoftThreshold(v,bound)/sigma[j,j];
       if (oldval != beta[j]){
         vs <- vs + (oldval-beta[j])*sigma.tilde[,j];
       }
@@ -112,7 +112,7 @@ InverseLinftyOneRow <- function ( sigma, i, mu, maxiter=50, threshold=1e-10) {
       #  vs <- -sigma.tilde%*%beta;
     }
 
-    # print(c(iter, maxiter, diff.norm2, threshold * last.norm2, threshold, mu))
+    # print(c(iter, maxiter, diff.norm2, threshold * last.norm2, threshold, bound))
 
   }
 
@@ -142,19 +142,21 @@ n = 100; p = 50
 X = matrix(rnorm(n * p), n, p)
 S = t(X) %*% X / n
 
-mu = 7.791408e-02
+debiasing_bound = 7.791408e-02
 
 tol = 1.e-12
 
-rows = c(1:2)
-A1 = debiasingMatrix(S, FALSE, n, rows, mu=mu, max_iter=1000, kkt_tol=tol, objective_tol=tol, parameter_tol=tol)
-A2 = debiasingMatrix(S / n, FALSE, n, rows, mu=mu, max_iter=1000, kkt_tol=tol, objective_tol=tol, parameter_tol=tol)
+rows = as.integer(c(1:2))
+print('here')
+print(rows)
+A1 = debiasingMatrix(S, FALSE, n, rows, bound=debiasing_bound, max_iter=1000, kkt_tol=tol, objective_tol=tol, parameter_tol=tol, linesearch=FALSE)
 
-B1 = debiasingMatrix(X, TRUE, n, rows, mu=mu, max_iter=1000, kkt_tol=tol, objective_tol=tol, parameter_tol=tol)
-B2 = debiasingMatrix(X / sqrt(n), TRUE, n, rows, mu=mu, max_iter=1000, kkt_tol=tol, objective_tol=tol, parameter_tol=tol)
+A2 = debiasingMatrix(S / n, FALSE, n, rows, bound=debiasing_bound, max_iter=1000, kkt_tol=tol, objective_tol=tol, parameter_tol=tol, linesearch=FALSE)
+B1 = debiasingMatrix(X, TRUE, n, rows, bound=debiasing_bound, max_iter=1000, kkt_tol=tol, objective_tol=tol, parameter_tol=tol, linesearch=FALSE)
+B2 = debiasingMatrix(X / sqrt(n), TRUE, n, rows, bound=debiasing_bound, max_iter=1000, kkt_tol=tol, objective_tol=tol, parameter_tol=tol, linesearch=FALSE)
 
-C1 = InverseLinfty(S, n, mu=mu, maxiter=1000)[rows,]
-C2 = InverseLinfty(S / n, n, mu=mu, maxiter=1000)[rows,]
+C1 = InverseLinfty(S, n, bound=debiasing_bound, maxiter=1000)[rows,]
+C2 = InverseLinfty(S / n, n, bound=debiasing_bound, maxiter=1000)[rows,]
 
 par(mfrow=c(2,3))
 
@@ -172,30 +174,30 @@ print(c('C', sum(C1[1,] == 0)))
 
 ## Are our points feasible
 
-feasibility = function(S, soln, j, mu) {
+feasibility = function(S, soln, j, debiasing_bound) {
      p = nrow(S)
      E = rep(0, p)
      E[j] = 1
      G = S %*% soln - E
-     return(c(max(abs(G)), mu))
+     return(c(max(abs(G)), debiasing_bound))
 }
 
-print(c('feasibility A', feasibility(S, A1[1,], 1, mu)))
-print(c('feasibility B', feasibility(S, B1[1,], 1, mu)))
-print(c('feasibility C', feasibility(S, C1[1,], 1, mu)))
+print(c('feasibility A', feasibility(S, A1[1,], 1, debiasing_bound)))
+print(c('feasibility B', feasibility(S, B1[1,], 1, debiasing_bound)))
+print(c('feasibility C', feasibility(S, C1[1,], 1, debiasing_bound)))
 
-active_KKT = function(S, soln, j, mu) {
+active_KKT = function(S, soln, j, debiasing_bound) {
      p = nrow(S)
      E = rep(0, p)
      E[j] = 1
      G = S %*% soln - E
      print(which(soln != 0))
      print(G[j])
-     return(c(G[soln != 0] * sign(soln)[soln != 0], mu))
+     return(c(G[soln != 0] * sign(soln)[soln != 0], debiasing_bound))
 }
 
-print(c('active_KKT A', active_KKT(S, A1[1,], 1, mu)))
-print(c('active_KKT B', active_KKT(S, B1[1,], 1, mu)))
-print(c('active_KKT C', active_KKT(S, C1[1,], 1, mu)))
+print(c('active_KKT A', active_KKT(S, A1[1,], 1, debiasing_bound)))
+print(c('active_KKT B', active_KKT(S, B1[1,], 1, debiasing_bound)))
+print(c('active_KKT C', active_KKT(S, C1[1,], 1, debiasing_bound)))