Merge pull request #26 from jonathan-taylor/using_Rcpp

Using rcpp to speed up debiased LASSO

Still to decide: strategy to limit size of active set

Author: jonathan-taylor

.travis.yml

@@ -7,8 +7,10 @@ r:
   - devel
 addons:
    apt: 
-      packages: libmpfr-dev
+      packages: libmpfr-dev 
 warnings_are_errors: true
 before_install:
   - tlmgr install index # for texlive and vignette?
+  - R -e 'install.packages("Rcpp", repos="http://cloud.r-project.org")'
+  - make Rcpp
   - cd selectiveInference

Makefile

@@ -0,0 +1,14 @@
+Rcpp: 
+	- rm -f selectiveInference/src/RcppExports.cpp
+	- rm -f selectiveInference/R/RcppExports.R
+	Rscript -e "library(Rcpp); Rcpp::compileAttributes('selectiveInference')"
+
+install: Rcpp
+	R CMD install selectiveInference
+
+build: 
+	R CMD build selectiveInference
+
+check: Rcpp build
+	R CMD build selectiveInference
+	R CMD check selectiveInference_1.2.2.tar.gz # fix this to be a script variable

README.md

@@ -21,3 +21,11 @@ The latest release of the package can be installed through CRAN:
 install.packages("selectiveInference")
 ```
 Code in repo is under development and may be unstable.
+
+## For development
+
+You will have to run
+
+```
+make Rcpp
+```

selectiveInference/DESCRIPTION

@@ -18,3 +18,5 @@ Description: New tools for post-selection inference, for use with forward
     models.
 License: GPL-2
 RoxygenNote: 5.0.1
+LinkingTo: Rcpp
+Imports: Rcpp

selectiveInference/NAMESPACE

@@ -44,5 +44,5 @@ importFrom("stats", dnorm, lsfit, pexp, pnorm, predict,
            qnorm, rnorm, sd, uniroot, dchisq, model.matrix, pchisq)
 importFrom("stats", "coef", "df", "lm", "pf")
 importFrom("stats", "glm", "residuals", "vcov")
-
+importFrom("Rcpp", "sourceCpp")
 

selectiveInference/R/funs.common.R

@@ -152,13 +152,7 @@ updateQR <- function(Q1,Q2,R,col) {
   n = ncol(Q1)
   k = ncol(Q2)
 
-  a = .C("update1",
-    Q2=as.double(Q2),
-    w=as.double(t(Q2)%*%col),
-    m=as.integer(m),
-    k=as.integer(k),
-    dup=FALSE,
-    package="selectiveInference")
+  a = update1_(as.matrix(Q2), t(Q2)%*%col, m, k) # Rcpp call
 
   Q2 = matrix(a$Q2,nrow=m)
   w = c(t(Q1)%*%col,a$w)

selectiveInference/R/funs.fixed.R

@@ -3,9 +3,9 @@
 # min 1/2 || y - \beta_0 - X \beta ||_2^2 + \lambda || \beta ||_1
 
 fixedLassoInf <- function(x, y, beta, lambda, family=c("gaussian","binomial","cox"),intercept=TRUE, add.targets=NULL, status=NULL,
-sigma=NULL, alpha=0.1,
-                     type=c("partial","full"), tol.beta=1e-5, tol.kkt=0.1,
-                     gridrange=c(-100,100), bits=NULL, verbose=FALSE) {
+                          sigma=NULL, alpha=0.1,
+                          type=c("partial","full"), tol.beta=1e-5, tol.kkt=0.1,
+                          gridrange=c(-100,100), bits=NULL, verbose=FALSE, linesearch.try=10) {
 
   family = match.arg(family)
   this.call = match.call()
@@ -26,8 +26,6 @@ sigma=NULL, alpha=0.1,
 
   else{
 
-    
-    
     checkargs.xy(x,y)
     if (missing(beta) || is.null(beta)) stop("Must supply the solution beta")
     if (missing(lambda) || is.null(lambda)) stop("Must supply the tuning parameter value lambda")
@@ -159,8 +157,8 @@ sigma=NULL, alpha=0.1,
       hsigmaSinv <- solve(hsigmaS) # pinv(hsigmaS)
 
       # Approximate inverse covariance matrix for when (n < p) from lasso_Inference.R
-      useC = TRUE
-      htheta <- InverseLinfty(hsigma, n, length(S), verbose=FALSE, useC=useC)
+
+      htheta <- InverseLinfty(hsigma, n, length(S), verbose=FALSE, max.try=linesearch.try)
       # htheta <- InverseLinfty(hsigma, n, verbose=FALSE)
 
       FS = rbind(diag(length(S)),matrix(0,pp-length(S),length(S)))
@@ -270,8 +268,7 @@ fixedLasso.poly=
 ### Functions borrowed and slightly modified from lasso_inference.R
 
 ## Approximates inverse covariance matrix theta
-InverseLinfty <- function(sigma, n, e, resol=1.5, mu=NULL, maxiter=50, threshold=1e-2, verbose = TRUE, useC = FALSE) {
-  # InverseLinfty <- function(sigma, n, resol=1.5, mu=NULL, maxiter=50, threshold=1e-2, verbose = TRUE) {
+InverseLinfty <- function(sigma, n, e, resol=1.2, mu=NULL, maxiter=50, threshold=1e-2, verbose = TRUE, max.try=10) {
     isgiven <- 1;
   if (is.null(mu)){
     isgiven <- 0;
@@ -293,16 +290,15 @@ InverseLinfty <- function(sigma, n, e, resol=1.5, mu=NULL, maxiter=50, threshold
     mu.stop <- 0;
     try.no <- 1;
     incr <- 0;
-    while ((mu.stop != 1)&&(try.no<10)){
+
+    output = NULL
+
+    while ((mu.stop != 1) && (try.no<max.try) ){
       last.beta <- beta
-      if (useC == FALSE) {
-            output <- InverseLinftyOneRow(sigma, i, mu, maxiter=maxiter, threshold=threshold)
-      } else {
-            output <- InverseLinftyOneRowC(sigma, i, mu, maxiter=maxiter)
-      }
-      beta <- output$optsol
+      output <- InverseLinftyOneRow(sigma, i, mu, maxiter=maxiter, soln_result=output) # uses a warm start
+      beta <- output$soln
       iter <- output$iter
-      if (isgiven==1){
+      if (isgiven==1) {
         mu.stop <- 1
       }
       else{
@@ -339,98 +335,54 @@ InverseLinfty <- function(sigma, n, e, resol=1.5, mu=NULL, maxiter=50, threshold
   return(M)
 }
 
-InverseLinftyOneRowC <- function (Sigma, i, mu, maxiter=50) {
-
-         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_diag=as.double(diag(Sigma)),
-		 Sigma_theta=as.double(rep(0, p)),
-		 nrow=as.integer(p),
-   		 bound=as.double(mu),
-		 theta=as.double(theta),
-		 maxiter=as.integer(50),
-		 row=as.integer(i-1),
-		 coord=as.integer(i-1),
-		 dup=FALSE,
-		 package="selectiveInference")
-
-	# Check feasibility
-
-	if (max(abs(Sigma %*% val$theta - basis_vector)) > 1.01 * mu) {
-	   warning("Solution for row of M does not seem to be feasible")
-	}
-
-	return(val$theta)
-}
+InverseLinftyOneRow <- function (Sigma, i, mu, maxiter=50, soln_result=NULL, kkt_tol=1.e-6, objective_tol=1.e-6,
+		                 use_QP=TRUE) {
 
-InverseLinftyOneRow <- function ( sigma, i, mu, maxiter=50, threshold=1e-2 ) {
-  p <- nrow(sigma);
-  rho <- max(abs(sigma[i,-i])) / sigma[i,i];
-  mu0 <- rho/(1+rho);
-  beta <- rep(0,p);
-  
-  #if (mu >= mu0){
-  #  beta[i] <- (1-mu0)/sigma[i,i];
-  #  returnlist <- list("optsol" = beta, "iter" = 0);
-  #  return(returnlist);
-  #}
-  
-  diff.norm2 <- 1;
-  last.norm2 <- 1;
-  iter <- 1;
-  iter.old <- 1;
-  beta[i] <- (1-mu0)/sigma[i,i];
-  beta.old <- beta;
-  sigma.tilde <- sigma;
-  diag(sigma.tilde) <- 0;
-  vs <- -sigma.tilde%*%beta;
-  
-  while ((iter <= maxiter) && (diff.norm2 >= threshold*last.norm2)){
-    
-    for (j in 1:p){
-      oldval <- beta[j];
-      v <- vs[j];
-      if (j==i)
-        v <- v+1;
-      beta[j] <- SoftThreshold(v,mu)/sigma[j,j];
-      if (oldval != beta[j]){
-        vs <- vs + (oldval-beta[j])*sigma.tilde[,j];
-      }
-    }
-    
-    iter <- iter + 1;
-    if (iter==2*iter.old){
-      d <- beta - beta.old;
-      diff.norm2 <- sqrt(sum(d*d));
-      last.norm2 <-sqrt(sum(beta*beta));
-      iter.old <- iter;
-      beta.old <- beta;
-      if (iter>10)
-        vs <- -sigma.tilde%*%beta;
-    }
-  }
-  
-  returnlist <- list("optsol" = beta, "iter" = iter)
-  return(returnlist)
-}
+  # If soln_result is not Null, it is used as a warm start.
+  # It should be a list
+  # with entries "soln", "gradient", "ever_active", "nactive"
+
+  p = nrow(Sigma)
 
-SoftThreshold <- function( x, lambda ) {
-  #
-  # Standard soft thresholding
-  #
-  if (x>lambda){
-    return (x-lambda);}
+  if (is.null(soln_result)) {
+     soln = rep(0, p)
+     ever_active = rep(0, p)
+     ever_active[1] = i      # 1-based
+     ever_active = as.integer(ever_active)
+     nactive = as.integer(1)
+     if (use_QP) {
+          linear_func = rep(0, p)
+	  linear_func[i] = -1
+	  linear_func = as.numeric(linear_func)
+	  gradient = 1. * linear_func 
+     } else {
+          gradient = rep(0, p)
+     }
+  }
   else {
-    if (x< (-lambda)){
-      return (x+lambda);}
-    else {
-      return (0); }
+     soln = soln_result$soln
+     gradient = soln_result$gradient  
+     ever_active = as.integer(soln_result$ever_active)
+     nactive = as.integer(soln_result$nactive)
+     if (use_QP) { 
+         linear_func = soln_result$linear_func
+     }
+  }
+
+  if (use_QP) {
+      result = solve_QP(Sigma, mu, maxiter, soln, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol) 
+  } else {
+      result = find_one_row_debiasingM(Sigma, i, mu, maxiter, soln, gradient, ever_active, nactive, kkt_tol, objective_tol) # C function uses 1-based indexing for active set
   }
+
+  # Check feasibility
+
+  if (!result$kkt_check) {
+     warning("Solution for row of M does not seem to be feasible")
+  } 
+
+  return(result)
+
 }
 
 ##############################

selectiveInference/R/funs.lar.R

@@ -254,14 +254,7 @@ downdateQR <- function(Q1,Q2,R,col) {
   m = nrow(Q1)
   n = ncol(Q1)
 
-  a = .C("downdate1",
-    Q1=as.double(Q1),
-    R=as.double(R),
-    col=as.integer(col-1),
-    m=as.integer(m),
-    n=as.integer(n),
-    dup=FALSE,
-    package="selectiveInference")
+  a = downdate1_(as.matrix(Q1), R, col, m, n) # Rcpp call
 
   Q1 = matrix(a$Q1,nrow=m)
   R = matrix(a$R,nrow=n)

selectiveInference/man/fixedLassoInf.Rd

@@ -10,9 +10,9 @@ fixed value of the tuning parameter lambda
 }
 \usage{
 fixedLassoInf(x, y, beta, lambda, family = c("gaussian", "binomial", 
-    "cox"),intercept=TRUE, add.targets=NULL, status=NULL, sigma=NULL, alpha=0.1,
+         "cox"),intercept=TRUE, add.targets=NULL, status=NULL, sigma=NULL, alpha=0.1,
          type=c("partial","full"), tol.beta=1e-5, tol.kkt=0.1,
-         gridrange=c(-100,100), bits=NULL, verbose=FALSE) 
+         gridrange=c(-100,100), bits=NULL, verbose=FALSE, linesearch.try=10) 
 }
 \arguments{
   \item{x}{
@@ -85,7 +85,11 @@ helpful (though computationally more costly). In particular, extra precision mig
 if the values in the output columns of \code{tailarea} differ noticeably from alpha/2.
 }
 \item{verbose}{
-Print out progress along the way? Default is FALSE}
+Print out progress along the way? Default is FALSE
+}
+\item{linesearch.try}{
+When running type="full" (i.e. debiased LASSO) how many attempts in the line search?
+}
 }
 
 \details{
@@ -145,7 +149,7 @@ p = 10
 sigma = 1
 
 x = matrix(rnorm(n*p),n,p)
-x=scale(x,TRUE,TRUE)
+x = scale(x,TRUE,TRUE)
 
 beta = c(3,2,rep(0,p-2))
 y = x\%*\%beta + sigma*rnorm(n)
@@ -165,10 +169,10 @@ out
 
 ## as above, but use lar function instead to get initial 
 ## lasso fit (should get same results)
- lfit = lar(x,y,normalize=FALSE)
- beta = coef(lfit,s=lambda,mode="lambda")
- out2 = fixedLassoInf(x,y,beta,lambda,sigma=sigma)
- out2
+lfit = lar(x,y,normalize=FALSE)
+beta = coef(lfit,s=lambda,mode="lambda")
+out2 = fixedLassoInf(x,y,beta,lambda,sigma=sigma)
+out2
 
 ## mimic different penalty factors by first scaling x
  set.seed(43)
@@ -249,5 +253,30 @@ status=sample(c(0,1),size=n,replace=TRUE)
      # compute fixed lambda p-values and selection intervals
      out = fixedLassoInf(x,tim,beta_hat,lambda,status=status,family="cox")
      out
+
+# Debiased lasso or "full"
+
+n = 50
+p = 100
+sigma = 1
+
+x = matrix(rnorm(n*p),n,p)
+x = scale(x,TRUE,TRUE)
+
+beta = c(3,2,rep(0,p-2))
+y = x\%*\%beta + sigma*rnorm(n)
+
+# first run glmnet
+gfit = glmnet(x, y, standardize=FALSE, intercept=FALSE)
+
+# extract coef for a given lambda; note the 1/n factor!
+# (and we don't save the intercept term)
+lambda = 2.8
+beta = coef(gfit, s=lambda/n, exact=TRUE)[-1]
+
+# compute fixed lambda p-values and selection intervals
+out = fixedLassoInf(x, y, beta, lambda, sigma=sigma, type='full', intercept=FALSE)
+out
+
 }
 

selectiveInference/src/Makevars

@@ -0,0 +1,11 @@
+PKG_CFLAGS= -I.
+PKG_CPPFLAGS= -I.
+PKG_LIBS=-L. 
+
+$(SHLIB): Rcpp Rcpp-matrixcomps.o Rcpp-debias.o RcppExports.o debias.o
+
+clean:
+	rm -f *o
+
+Rcpp:
+	Rscript -e "library(Rcpp); Rcpp::compileAttributes('..')"