R/LogisticSlope.R

In dkucharc/glmSlope: Sorted L1 Penalized Estimation (SLOPE) for generalized linear models (GLM)

LogisticSlope <- function(X, y, lambda, options=list())
{

  # -------------------------------------------------------------
  # Start times
  # -------------------------------------------------------------
  t0 <- proc.time()[3]
  
  # -------------------------------------------------------------
  # Define function for retrieving option fields with defaults
  # -------------------------------------------------------------
  getDefaultField <- function(options,name,default)
  {  if (!is.null(options[[name]]))
  {  return(options[[name]])  }
    else
    {  return(default)  }
  }
  
  # -------------------------------------------------------------
  # Parse parameters
  # -------------------------------------------------------------
  iterations <- getDefaultField(options,"iterations", 10000)
  verbosity  <- getDefaultField(options,"verbosity" , 1)
  optimIter  <- getDefaultField(options,"optimIter" , 1)
  tolInfeas  <- getDefaultField(options,"tolInfeas" , 1e-6)
  tolRelGap  <- getDefaultField(options,"tolRelGap" , 1e-6)
  wInit      <- getDefaultField(options,"wInit"     , vector())
  
  # -------------------------------------------------------------
  # Ensure that lambda is non-increasing
  # -------------------------------------------------------------
  n <- length(lambda)
  matrix(lambda,c(1,n))
  if ((n > 1) && any(lambda[2:n] > lambda[1:n-1]))
  {  stop("Lambda must be non-increasing.")
  }
  
  # -------------------------------------------------------------
  # Initialize
  # -------------------------------------------------------------
  # Get problem dimension
  n <- ncol(X)
  
  # Get initial lower bound on the Lipschitz constant
  rndKind  <- RNGkind()
  rndState <- tryCatch({.Random.seed},error=function(m) {} )
  set.seed(0,kind="Mersenne-Twister",normal.kind="Inversion")
  w <- matrix(rnorm(n),c(n,1))
  w <- w / sqrt(sum(w^2))
  w <- t(X) %*% (X %*% w)
  L <- sqrt(sum(w^2))/4
  if (!is.null(rndState))
    set.seed(rndState[-1],kind=rndKind[1],normal.kind=rndKind[2])
  
  # Set constants
  STATUS_RUNNING    <- 0
  STATUS_OPTIMAL    <- 1
  STATUS_ITERATIONS <- 2
  STATUS_MSG        <- c('Optimal','Iteration limit reached')
  
  # Initialize parameters and iterates
  if (length(wInit) == 0) wInit <- matrix(0,n,1)
  t       <- 1
  eta     <- 2
  lambda  <- matrix(lambda,nrow=length(lambda))
  y       <- matrix(y,nrow=length(y))
  Y       <- diag(as.vector(y), length(y))
  YX      <- Y %*% X
  w       <- wInit
  w0      <- 0
  v       <- w
  iter    <- 0
  status  <- STATUS_RUNNING
  
  # Deal with Lasso case
  modeLasso <- (length(lambda) == 1)
  if (modeLasso)
    proxFunction <- function(x,lambda) { return(sign(x) * pmax(abs(x) - lambda,0)) }
  else
    proxFunction <- function(v1,v2) { return(prox_sorted_L1(v1,v2)) }
  
  if (verbosity > 0)
  {  printf <- function(...) invisible(cat(sprintf(...)))
  printf('%5s %9s  %9s  %9s\n','Iter','Gap','Infeas.','Rel. gap')
  }
  
  # -------------------------------------------------------------
  # Main loop
  # -------------------------------------------------------------
  while (TRUE)
  {
    # Compute the gradient at f(v)
    r <- 1 / (1 + exp(YX %*% v))
    g <- -as.double(crossprod(YX, r))
    log1mr <- log(1-r)
    f <- -sum(log1mr)

    # Increment iteration count
    iter <- iter + 1
    
    # Check optimality conditions
    if ((iter %% optimIter) == 0)
    {  
      # Compute 'dual', check infeasibility and gap
      if (modeLasso) { 
        infeas    <- max(abs(g)-lambda,0)
        objPrimal <- f + lambda * norm(v,'1')
      } else {  
        gs     <- sort(abs(g), decreasing=TRUE)
        vs     <- sort(abs(v), decreasing=TRUE)
        infeas <- max(max(cumsum(gs-lambda)),0)
        # Compute primal and dual objective
        objPrimal <-  f + as.double(crossprod(lambda,vs))
      }
      objDual   <- as.double(crossprod(r-1, log1mr)) - as.double(crossprod(r, log(r)))
      
      # Format string
      if (verbosity > 0)
        str <- sprintf('   %9.2e  %9.2e  %9.2e',objPrimal - objDual, infeas/lambda[[1]], abs(objPrimal - objDual) / max(1,objPrimal))
      
      # Check primal-dual gap
      if ((abs(objPrimal - objDual)/max(1,objPrimal) < tolRelGap) &&
          (infeas < tolInfeas * lambda[[1]]))
        status <- STATUS_OPTIMAL
    }
    else
    {   str <- ''
    }
    
    if (verbosity > 0)
    {  if ((verbosity == 2) ||
           ((verbosity == 1) && ((iter %% optimIter) == 0)))
    {  printf('%5d %s\n', iter, str)
    }
    }
    
    # Stopping criteria
    if ((status == 0) && (iter >= iterations))
      status <- STATUS_ITERATIONS
    
    if (status != 0)
    {  if (verbosity > 0)
      printf('Exiting with status %d -- %s\n', status, STATUS_MSG[[status]])
      break
    }
    
    # Keep copies of previous values
    wPrev  <- w
    fPrev  <- f
    tPrev  <- t
    
    # Lipschitz search
    while (TRUE)
    {  # Compute prox mapping
      w <- proxFunction(v - (1/L)*g, lambda/L)
      d <- w - v
      
      r  <- 1/(1 + exp(YX %*% w))
      f  <- -sum(log(1-r))
      q  <- fPrev + as.double(crossprod(d,g)) + (L/2)*as.double(crossprod(d))
      
      if (q >= f*(1-1e-12))
        break
      else
        L <- L * eta
    } # Lipschitz search
    
    # Update
    t <- (1 + sqrt(1 + 4*t^2)) / 2
    v <- w + ((tPrev - 1) / t) * (w - wPrev)
  } # While (TRUE)
  
  
  # Information structure
  info <- c()
  info$runtime   <- proc.time()[3] - t0
  info$objPrimal <- objPrimal
  info$objDual   <- objDual
  info$infeas    <- infeas
  info$status    <- status
  
  info$w         <- v
  info$L         <- L
  
  return(info)
} # Function Adlas