R/autotune_lasso.R

In autotune: Faster and more Efficient Lasso (than 'glmnet' and 'scalreg') with Data-Driven Tuning

#' Data-adaptive automatic and fast tuning of LM with Lasso regularization
#' @name autotune_lasso
#' @param x Matrix of predictors, with one column for each predictor.
#' Dimension will be \code{nobs} \eqn{\times} \code{nvars}; so each row is a new observation
#' and \code{nvars} > 1.
#' @param y Vector of responses.
#' @param alpha Default 0.01, significance level of sequential F-tests 
#' used for estimation of support set.
#' @param standardize Logical flag for standardization of all variables in x, prior to
#' fitting the model sequence. The coefficients are always returned on the
#' original scale. Default value is \code{TRUE}.
#' @param standardize_response Logical flag for demeaning the reponse variable y.
#' Default value is \code{TRUE}.
#' @param intercept Should intercept(s) be fitted (default=\code{TRUE}) or set to zero
#' (\code{FALSE}).
#' @param active Should active set selection be used (\code{TRUE}) or avoided
#' (default = \code{FALSE}). It is under experimentation phase, use it 
#' with care.
#' @param trace_it Logical input, default \code{FALSE}; if \code{TRUE}
#' prints out the iteration number while running, useful for big datasets that 
#' take a long time to fit.
#' @param tolerance Numeric input for an additional stopping criteria on the
#' coordinate descent when sigma is updating. Stops that coordinate
#' descent when the relative change between successive iterates of
#' coefficients' estimates is less than the \code{tolerance}. Default value 
#' is \eqn{10^{-4}}.
#' @param beta_tolerance Numeric input for the stopping criteria on the
#' coordinate descent when sigma is not updating. Stops that coordinate
#' descent when the relative change between successive iterates of
#' coefficients' estimates is less than the \code{beta_tolerance}.
#' Default value is \eqn{10^{-3}}.
#' @param iter_max Maximum number of iterations of coordinate descent
#' allowed when sigma is updating. Default value is 30.
#' @param beta_iter_max Maximum number of iterations of coordinate 
#' descent allowed when sigma is not updating. Default value is 40.
#' @param active_iter_max If \code{active = TRUE}, maximum number of 
#' times the active set is updated as per the violations of KKT 
#' conditions. Default value is 5.
#' @param PR_norm_l2 Logical flag to whether use the l2 norm of partial 
#' residuals for ordering them instead of the default l1 norm. 
#' @param ... Additional arguments passed to the internal fitting routine.
#' 
#' 
#' @return A list with various intermediate and final outputs produced by autotune Lasso in its regularization path.
#' 
#' \item{beta}{ Final estimates of regression coefficients. }
#' \item{a0}{ Intercept of the fit. }
#' \item{lambda}{Final thresholding value \eqn{\lambda =
#' \lambda_0\hat\sigma^2} used in the coordinate descent after noise
#' variance estimate \eqn{\hat\sigma^2} has converged.}
#' \item{sigma_sq}{ Final estimate of noise variance \eqn{\sigma^2}}
#' \item{nobs}{ Number of observations.}
#' \item{nvars}{ Number of variables.}
#' \item{CD.path.details}{ A list of additional details about the coordinate descent path taken by
#' autotune Lasso:
#'   \describe{
#'     \item{\code{sorted_predictors}}{Decreasing ordering of predictors in terms of
#'     their contribution to predicting the response values.}
#'     \item{\code{sigma_sq_seq}}{A \code{no_of_iterations}-length sequence of
#'     noise variance estimates \eqn{\sigma^2}.}
#'     \item{\code{beta_matrix}}{A (\code{length(sigma_sq_seq) + 1})
#'     \eqn{\times} \code{nvars} matrix of estimated coefficients.}
#'     \item{\code{no_of_iter_before_lambda_conv}}{Number of coordinate descent
#'     iterations performed before the noise variance estimate
#'     \eqn{\hat{\sigma}^2} converged.}
#'     \item{\code{no_of_iter_after_lambda_conv}}{After the noise variance estimate
#'     \eqn{\hat{\sigma}^2} has converged, the number of coordinate descent
#'     iterations required for coefficients \eqn{\hat\beta} to converge.}
#'     \item{\code{no_of_iterations}}{Total number of coordinate descent iterations
#'     implemented by autotune Lasso.}
#'     \item{\code{lambda0}}{Value of \eqn{\lambda_0} used in autotune Lasso.
#'     Refer to the original paper for details.}
#'     \item{\code{support_set}}{Final set of predictors included in the support set
#'     for noise variance estimation by autotune Lasso.}
#'     \item{\code{count_sig_beta}}{A \code{no_of_iterations}-length vector containing
#'     the support set sizes across the coordinate descent iterations while the
#'     noise variance estimate is being updated.}
#'     \item{\code{null_support}}{Boolean output indicating whether the final support
#'     set estimate is a null set. If so, \code{autotune_lasso} uses the support set
#'     estimate from the previous iteration to obtain the final noise variance
#'     estimate \eqn{\sigma^2}.}
#'     \item{\code{active_iterations}}{If \code{active = TRUE}, the number of times
#'     the active set is updated according to violations of the KKT conditions.}
#'     \item{\code{active_set_sizes}}{An \code{active_iterations}-length vector
#'     containing the active-set sizes used for coordinate descent across the
#'     active-set iterations.}
#'   }
#' }
#' 
#' 
#' 
#' 
#' @description
#' Fits a linear model via alternative optimization of penalized gaussian maximum likelihood 
#' which is a biconvex function of regression coefficients \eqn{\beta} and noise variance \eqn{\sigma^2}. 
#' The regularization path of \code{autotune_lasso} quickly picks out a good lambda for Lasso and then
#' returns the corresponding linear fit along with various attributes related to the fit.
#' 
#' 
#' @export
#'
#' @examples
#' library(autotune)
#' set.seed(10)
#' n <- 80
#' p <- 400
#' s <- 5
#' snr <- 4
#' betatrue <- c(rep(1,s), rep(0, p - s))
#' x <- matrix(rnorm(n * p), ncol = p)
#' error.sd <- sqrt((betatrue %*% betatrue)/snr)
#' err <- rnorm(n, sd = error.sd)
#' y <- x %*% betatrue + err
#' ans <- autotune_lasso(x, y, trace_it = TRUE)
#' b <- betatrue
#' # The Predictors which are actually significant:
#' which(b != 0)
#' # The Predictors which had nonzero estmated coefficients:
#' which(ans$beta != 0)
#' # Top 10 predictors X_i's in the ranking of X_i's given by autotune:
#' ans$CD.path.details$sorted_predictors[1:10]
#' # Cardinality of autotune's estimated support set in each CD iteration before lambda converged:
#' ans$CD.path.details$count_sig_beta
#' # Sigma estimates in each CD iteration:
#' ans$CD.path.details$sigma_sq_seq
#' # Empirical noise variance:
#' var(err)
#' 
#' 
autotune_lasso <- function(x,
                           y,
                           alpha = 0.01,
                           standardize = TRUE,
                           standardize_response = TRUE,
                           intercept = TRUE,
                           active = FALSE,
                           trace_it = FALSE,
                           tolerance = 1e-4,
                           beta_tolerance = 1e-3,
                           iter_max = 30,
                           beta_iter_max = 40,
                           active_iter_max = 5,
                           PR_norm_l2 = FALSE,
                           ...) {
  cl <- match.call()
  
  fit <- autotune_lasso_cpp(
    xin = x,
    yin = y,
    alpha = alpha,
    standardize = standardize,
    standardize_response = standardize_response,
    intercept = intercept,
    active = active,
    trace_it = trace_it,
    tolerance = tolerance,
    beta_tolerance = beta_tolerance,
    iter_max = iter_max,
    beta_iter_max = beta_iter_max,
    active_iter_max = active_iter_max,
    PR_norm_l2 = PR_norm_l2,
    ...
  )
  fit$call <- cl
  fit$x <- x
  fit$y <- y
  fit
}