R/getGradientLARS.R

In rRAP: Real-Time Adaptive Penalization for Streaming Lasso Models

getGradientLARS <-
function(XtX, Xty, beta, Eps=0.01){
  # calculate gradient of lasso solution as suggested by LARS formulae
  #
  #
  # INPUT:
  #      - XtX: current estimate of covariance
  #      - Xty: current estimate of t(X)%*%y - only used if all beta entries are 0!
  #      - beta: current estimate of LASSO solution
  
  gradVector = rep(0, length(beta))
  
  s = matrix(sign(beta), ncol=1)
  ActiveID = which(s!=0)
  if (length(ActiveID)>0){
    G_active = (XtX * (s %*% t(s)))[ActiveID, ActiveID]
    
    Ginv_active = try(solve(G_active), silent = TRUE)
    if (class(Ginv_active)=="try-error") {
      Ginv_active = ginv(G_active)
      #cat("gen inv used...\n")
    }
    A_active =  1/sqrt(t(rep(1, length(ActiveID))) %*% Ginv_active %*% rep(1, length(ActiveID)))
    
    gradVector[ActiveID] = as.numeric(A_active) * Ginv_active %*% rep(1, length(ActiveID))
    
    return(gradVector * s)
  } else {
    # here the active set is empty! The direction of the gradient is therefore
    # in the direction of the covariate that is maximally correlated with the response (in absolute terms of course!)
    
    gradVector[which.max(abs(Xty))] = Eps*sign(Xty[which.max(Xty)])
    
    return(gradVector)
  }
}