WIP: function to fit randomized lasso

Author: jonathan-taylor

selectiveInference/R/funs.randomized.R

@@ -0,0 +1,64 @@
+# Functions to fit and "infer" about parameters in the
+# randomized LASSO
+#
+# 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
+                                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)
+    			
+    noise_type = match.arg(noise_type)
+
+    if (noise_type == 'gaussian') {
+       D = Norm(mean=0, sd=noise_scale)
+    }
+    else if (noise_type == 'laplace') {
+       D = DExp(rate = 1 / noise_scale) # D is a Laplace distribution with rate = 1.
+    }
+    perturb_ = distr::r(D)(p)
+			
+    lam = as.numeric(lam)
+    if (length(lam) == 1) {
+       lam = rep(lam, p)
+    }
+    if (length(lam) != p) {
+       stop("Lagrange parameter should be single float or of length ncol(X)")
+    }    
+
+    soln = rep(0, p)
+    Xsoln = rep(0, n)
+    linear_func = (- t(X) %*% y - perturb_)
+    gradient = 1. * linear_func
+    ever_active = rep(0, p)
+    nactive = as.integer(0)
+
+    result = solve_QP_wide(X,                  # design matrix
+    	                   lam,                # vector of Lagrange multipliers
+		           ridge_term / n,     # ridge_term 
+                           max_iter, 
+                           soln, 
+                           linear_func, 
+                           gradient, 
+                           Xsoln,
+                           ever_active, 
+                           nactive, 
+                           kkt_tol, 
+                           objective_tol, 
+                           p,
+		           objective_stop,     # objective_stop
+			   kkt_stop,           # kkt_stop
+			   param_stop)         # param_stop
+    return(result)
+}

tests/test_randomized.R

@@ -0,0 +1,14 @@
+library(selectiveInference)
+
+test = function() {
+
+    n = 100; p = 50
+    X = matrix(rnorm(n * p), n, p)
+    y = rnorm(n)
+    lam = 20 / sqrt(n)
+    noise_scale = 0.01 * sqrt(n)
+    ridge_term = .1 / sqrt(n)
+    fit_randomized_lasso(X, y, lam, noise_scale, ridge_term)
+}
+
+print(test())