working! bugs fixed. new sampler

Author: Jelena Markovic

selectiveInference/R/funs.randomized.R

@@ -33,11 +33,15 @@ randomizedLasso = function(X,
     if (is.null(noise_scale)) {
         noise_scale = 0.5 * sd(y) * sqrt(mean_diag)
     }
-
+    
+    print(paste("ridge term", ridge_term))
+    print(paste("noise scale", noise_scale))
+    
     noise_type = match.arg(noise_type)
 
     if (noise_scale > 0) {
         if (noise_type == 'gaussian') {
+            set.seed(1)
             perturb_ = rnorm(p) * noise_scale
         }
         else if (noise_type == 'laplace') {
@@ -65,8 +69,8 @@ randomizedLasso = function(X,
     nactive = as.integer(0)
 
     result = solve_QP_wide(X,                  # design matrix
-    	                   lam / n,            # vector of Lagrange multipliers
-		           ridge_term / n,     # ridge_term 
+    	                     lam / n,            # vector of Lagrange multipliers
+		                       ridge_term / n,     # ridge_term 
                            max_iter, 
                            soln, 
                            linear_func, 
@@ -76,12 +80,12 @@ randomizedLasso = function(X,
                            nactive, 
                            kkt_tol, 
                            objective_tol, 
-			   parameter_tol,
+			                     parameter_tol,
                            p,
-		           objective_stop,     # objective_stop
-			   kkt_stop,           # kkt_stop
-			   parameter_stop)         # param_stop
-
+		                       objective_stop,     # objective_stop
+			                     kkt_stop,           # kkt_stop
+			                     parameter_stop)         # param_stop
+    
     sign_soln = sign(result$soln)
 
     unpenalized = lam == 0
@@ -96,7 +100,11 @@ randomizedLasso = function(X,
 
     observed_scalings = abs(result$soln)[active]
     observed_unpen = result$soln[unpenalized]
-    observed_subgrad = result$gradient[inactive]
+    observed_subgrad = -n*result$gradient[inactive]
+    
+    if (length(which(abs(observed_subgrad)>lam[1]))){
+      print("subgradient eq not satisfied")
+    }
 
     observed_opt_state = c(observed_unpen, observed_scalings, observed_subgrad)
 
@@ -111,14 +119,15 @@ randomizedLasso = function(X,
     coef_term = L_E
 
     signs_ = c(rep(1, sum(unpenalized)), sign_soln[active])
+    
+    coef_term[active,] = coef_term[active,] + ridge_term * diag(rep(1, sum(active)))  # ridge term
+  
     if (length(signs_) == 1) {
-        coef_term = coef_term * signs_
+      coef_term = coef_term * signs_
     } else {
-        coef_term = coef_term %*% diag(signs_)  # scaligns are non-negative
+      coef_term = coef_term %*% diag(signs_)  # scaligns are non-negative
     }
-
-    coef_term[active,] = coef_term[active,] + ridge_term * diag(rep(1, sum(active)))  # ridge term
-
+    
     subgrad_term = matrix(0, p, sum(inactive)) # for subgrad
     for (i in 1:sum(inactive)) {
         subgrad_term[inactive_set[i], i] = 1
@@ -155,7 +164,8 @@ randomizedLasso = function(X,
     inactive_lam = lam[inactive_set]
     inactive_start = sum(unpenalized) + sum(active)
     active_start = sum(unpenalized)
-
+    
+    
     # XXX only for Gaussian so far
 
     log_optimization_density = function(opt_state) {
@@ -185,9 +195,11 @@ randomizedLasso = function(X,
                 optimization_transform = opt_transform,
                 internal_transform = internal_transform,
                 log_optimization_density = log_optimization_density,
-		observed_opt_state = observed_opt_state,
+		            observed_opt_state = observed_opt_state,
                 observed_raw = observed_raw,
-		noise_scale = noise_scale
+		            noise_scale = noise_scale,
+		            soln = result$soln,
+		            perturb = perturb_
                 ))
 
 }
@@ -314,19 +326,22 @@ conditional_density = function(noise_scale, lasso_soln) {
   lasso_soln$log_optimization_density = log_condl_optimization_density
   lasso_soln$observed_opt_state = observed_opt_state[1:nactive]
   lasso_soln$optimization_transform = opt_transform
-  return(lasso_soln)
+  reduced_opt_transform =list(linear_term = reduced_B, offset_term = reduced_beta_offset)
+  return(list(lasso_soln=lasso_soln, 
+              reduced_opt_transform = reduced_opt_transform))
 }
 
 randomizedLassoInf = function(X, 
                               y, 
                               lam, 
+                              sampler="A",
                               sigma=NULL, 
                               noise_scale=NULL, 
                               ridge_term=NULL, 
                               condition_subgrad=TRUE, 
                               level=0.9,
-			      nsample=10000,
-			      burnin=2000,
+			                        nsample=10000,
+			                        burnin=2000,
                               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
@@ -353,22 +368,35 @@ randomizedLassoInf = function(X,
                                parameter_stop=parameter_stop)
 
   active_set = lasso_soln$active_set
-  if (length(active_set)==0){
+  nactive = length(active_set)
+  print(paste("nactive", nactive))
+  if (nactive==0){
     return (list(active_set=active_set, pvalues=c(), ci=c()))
   }
   inactive_set = lasso_soln$inactive_set
-  nactive = length(active_set)
+  
 
   noise_scale = lasso_soln$noise_scale # set to default value in randomizedLasso
 
  if (condition_subgrad==TRUE){
-   lasso_soln=conditional_density(noise_scale, lasso_soln)
+   condl_lasso=conditional_density(noise_scale, lasso_soln)
+   lasso_soln = condl_lasso$lasso_soln
+   reduced_opt_transform = condl_lasso$reduced_opt_transform
  } 
 
   ndim = length(lasso_soln$observed_opt_state)
-
-  S = sample_opt_variables(lasso_soln, jump_scale=rep(1/sqrt(n), ndim), nsample=nsample)
-  opt_samples = as.matrix(S$samples[(burnin+1):nsample,,drop=FALSE])
+  
+  if (sampler =="R"){
+    S = sample_opt_variables(lasso_soln, jump_scale=rep(1/sqrt(n), ndim), nsample=nsample)
+    opt_samples = as.matrix(S$samples[(burnin+1):nsample,,drop=FALSE])
+  } else if (sampler == "A"){
+    opt_samples = gaussian_sampler(noise_scale, 
+                                 lasso_soln$observed_opt_state, 
+                                 reduced_opt_transform$linear_term,
+                                 reduced_opt_transform$offset_term,
+                                 nsamples=nsample)
+    opt_sample = opt_samples[(burnin+1):nsample,]
+  }
 
   X_E = X[, active_set]
   X_minusE = X[, inactive_set]
@@ -458,3 +486,6 @@ randomizedLassoInf = function(X,
   }
   return(list(active_set=active_set, pvalues=pvalues, ci=ci))
 }
+
+
+

selectiveInference/R/sampler.R

@@ -0,0 +1,96 @@
+
+log_concave_sampler = function(negative_log_density, 
+                               grad_negative_log_density, 
+                               constraints,
+                               observed,
+                               nsamples){
+  #print(constraints)
+  constraints = as.matrix(constraints)
+  dim = nrow(constraints)
+  
+  get_poisson_process = function(state){
+    pos = as.matrix(state$pos)
+    velocity = as.matrix(state$velocity)
+    neg_velocity = velocity<0
+    pos_velocity = velocity>0
+    tau_min = 0
+    tau_max = 10
+    if (sum(neg_velocity)>0){
+      R = (-constraints[neg_velocity,1]+pos[neg_velocity])/(-velocity[neg_velocity])
+      tau_max = min(tau_max, min(R))
+      L = (-constraints[neg_velocity,2]+pos[neg_velocity])/(-velocity[neg_velocity])
+      tau_min = max(tau_min, max(L))
+    }
+    if (sum(pos_velocity)>0){
+      R = (constraints[pos_velocity,2]-pos[pos_velocity])/velocity[pos_velocity]
+      tau_max = min(tau_max, min(R))
+      L = (constraints[pos_velocity,1]-pos[pos_velocity])/velocity[pos_velocity]
+      tau_min = max(tau_min, max(L))
+    }
+    
+    f=function(t){as.numeric(t(velocity) %*% grad_negative_log_density(pos+velocity*t))}
+    tau_star = tau_max
+    if (f(tau_min)*f(tau_max)<0){
+      tau_star = uniroot(f, c(tau_min, tau_max))$root
+    } else{
+      if (negative_log_density(pos+velocity*tau_min)<negative_log_density(pos+velocity*tau_max)){
+        tau_star = tau_min
+      }
+    }
+    
+    tau_min = max(tau_min, tau_star)
+    
+    RHS = negative_log_density(pos+velocity*tau_star)+rexp(1)
+    g = function(t){negative_log_density(pos+velocity*t)-RHS}
+    if (g(tau_min)*g(tau_max)<0){
+      tau = uniroot(g, c(tau_min, tau_max))$root
+    } else{
+      tau = tau_max
+    }
+    return (tau)
+  }
+  
+  update_velocity = function(){
+    Z=rnorm(dim)
+    return(Z/sqrt(t(Z)%*%Z))
+  }
+  
+  compute_next = function(state){
+    bounce_time = get_poisson_process(state)/2
+    #print(paste("bounce time", bounce_time))
+    next_pos = state$pos+state$velocity*bounce_time
+    next_velocity=update_velocity()
+    return(list(pos=next_pos, velocity=next_velocity))
+  }
+  
+  state = list(pos=observed, velocity = update_velocity())
+  samples = matrix(0, nrow = nsamples, ncol = dim)
+  for (i in 1:nsamples){
+    #print(paste("pos", toString(state$pos)))
+    #print(paste("velocity", toString(state$velocity)))
+    samples[i,]=state$pos
+    state = compute_next(state)
+  }
+  return (samples)
+}
+
+gaussian_sampler = function(noise_scale, observed, linear_term, offset_term, nsamples){
+  
+  negative_log_density = function(x) {
+    recon = linear_term %*% x+offset_term
+    return(as.numeric(t(recon)%*%recon/(2*noise_scale^2)))
+  }
+  grad_negative_log_density=function(x){
+    recon = linear_term %*% x+offset_term
+    return(t(linear_term)%*% recon/(noise_scale^2))
+  }
+  dim = length(observed)
+  constraints = matrix(0,dim,2)
+  constraints[,2] = Inf
+  
+  return(log_concave_sampler(negative_log_density, 
+                             grad_negative_log_density,
+                             constraints,
+                             observed,
+                             nsamples))
+}

tests/randomized/test_instances.R

@@ -24,10 +24,46 @@ gaussian_instance = function(n, p, s, sigma=1, rho=0, signal=6, X=NA,
   return(result)
 }
 
+test_randomized_lasso = function(n=100,p=200,s=0){
+  set.seed(1)
+  data = gaussian_instance(n=n,p=p,s=s, rho=0.3, sigma=3)
+  X=data$X
+  y=data$y
+  lam = 2.
+  noise_scale = 0.5
+  ridge_term = 1./sqrt(n)
+  result = selectiveInference:::randomizedLasso(X,y,lam, noise_scale, ridge_term)
+  print(result$soln)
+  print(length(which(result$soln!=0)))
+  print(result$observed_opt_state) # compared with python code
+}
+
+test_KKT=function(){
+  set.seed(1)
+  n=200
+  p=100
+  data = gaussian_instance(n=n,p=p,s=0, rho=0.3, sigma=3)
+  X=data$X
+  y=data$y
+  lam = 2.
+  noise_scale = 0.5
+  ridge_term = 1./sqrt(n)
+  result = selectiveInference:::randomizedLasso(X,y,lam, noise_scale, ridge_term)
+  print("check KKT")
+  opt_linear = result$optimization_transform$linear_term
+  opt_offset = result$optimization_transform$offset_term
+  observed_opt_state=result$observed_opt_state
+  #print(dim(opt_linear))
+  #print(opt_offset)
+  #print(result$perturb)
+  print(opt_linear %*% observed_opt_state+opt_offset+result$observed_raw-result$perturb) ## should be zero
+}
+  
+
 
 collect_results = function(n,p,s, nsim=100, level=0.9, condition_subgrad=TRUE, lam=1.2){
 
-  rho=0.3
+  rho=0.
   sigma=1
   sample_pvalues = c()
   sample_coverage = c()
@@ -36,7 +72,14 @@ collect_results = function(n,p,s, nsim=100, level=0.9, condition_subgrad=TRUE, l
     X=data$X
     y=data$y
     beta=data$beta
-    result = selectiveInference:::randomizedLassoInf(X, y, lam, level=level, burnin=2000, nsample=4000, condition_subgrad=condition_subgrad)
+    result = selectiveInference:::randomizedLassoInf(X, y, 
+                                                     lam=lam, 
+                                                     sigma=sigma,
+                                                     level=level, 
+                                                     sampler = "A",
+                                                     burnin=1000, 
+                                                     nsample=5000, 
+                                                     condition_subgrad=condition_subgrad)
     true_beta = beta[result$active_set]
     coverage = rep(0, nrow(result$ci))
     if (length(result$active_set)>0){
@@ -61,7 +104,7 @@ collect_results = function(n,p,s, nsim=100, level=0.9, condition_subgrad=TRUE, l
 }
 
 set.seed(1)
-collect_results(n=200, p=100, s=0, lam=2)
-
-
+collect_results(n=100, p=2000, s=0, lam=3)
+#test_randomized_lasso()
+#test_KKT()
 

tests/randomized/test_sampler.R

@@ -0,0 +1,16 @@
+
+test_log_concave_sampler = function(){
+  samples = log_concave_sampler(negative_log_density= function(x){x^2/2}, 
+                                grad_negative_log_density=function(x){x},
+                                constraints = t(as.matrix(c(2,3))),
+                                observed = 2, nsamples=10000)
+  mean(samples)
+  hist(samples)
+}
+
+
+test_gaussian_sampler =function(){
+  samples = gaussian_sampler(1, 1, 1, 0,10000)
+  mean(samples)
+  hist(samples)
+}