# autotune Lasso

The R package `autotune` implements the Lasso with data-driven tuning for linear models.

## Installation and Loading

### Installation

The development version of the `autotune` package can be installed from GitHub using

    # install.packages("devtools")
    # Ensure that you have the Rcpp package installed with version >=1.0.13
    # devtools::install_github("Tathagata-S/autotune")

When installing from GitHub, in order to build the package from source, you need to have the appropriate R development tools installed (Rtools on Windows, or these tools on Mac).

### Load Package

After installation, the package can be loaded in the standard way:

    library(autotune)

## High Dimensional Regression

`autotune` performs lasso via the `autotune_lasso()` function.

We illustrate autotune lasso on simulated data using a linear model with \(s=10\), \(n=300\), \(p=500\). Reader can also specify different configurations of data generating process in the following code chunk.

    set.seed(10)
    n = 300
    p = 500
    s = 10
    beta = c(rep(1, s), rep(0, p - s))
    x = matrix(rnorm(n * p), ncol = p)
    # Maunal sigma allocation
    # y = x %*% beta + rnorm(n, sd = 1)
    # Dynamic sigma allocation with snr specified
    snr = 2
    y = x %*% beta + rnorm(n, sd = sqrt(var(x%*%beta)/snr))

### Runtime of autotune Lasso

Given data (x, y), run the autotune lasso as follows with default \(\alpha = 0.01\).

    ptm <- proc.time()
    fit.autotune <- autotune_lasso(x, y, alpha = 0.01)
    proc.time() - ptm
    #>    user  system elapsed
    #>   0.013   0.003   0.016

The regression coefficients \(\beta\)s, intercept, final lambda and the sequence of estimated sigmas can be extracted from the fitted autotune lasso object as follows

    b.autotune <- fit.autotune$beta
    intercept.autotune <- fit.autotune$a0
    lambda.autotune <- fit.autotune$lambda
    sigma.seq.autotune <- fit.autotune$CD.path.details$sigma_sq_seq
    sigma.estimate.autotune <- fit.autotune$sigma_sq

Now, `autotune` also supports the S3 function `coef` which returns the estimated regression coefficients \(\beta\)s along with the intercept

    head(coef(fit.autotune), n = 21)
    #> 21 x 1 sparse Matrix of class "dgCMatrix"
    #>
    #> (Intercept) -0.2506709
    #> V1           0.5329995
    #> V2           0.8984378
    #> V3           0.6647993
    #> V4           0.7491217
    #> V5           0.7049080
    #> V6           0.8335613
    #> V7           0.8657125
    #> V8           0.6522736
    #> V9           0.6986307
    #> V10          0.7410557
    #> V11          .
    #> V12          .
    #> V13          .
    #> V14          .
    #> V15          .
    #> V16         -0.1021190
    #> V17          .
    #> V18          .
    #> V19          .
    #> V20          .

## Comparison with Lasso tuned via Cross Validation (CV Lasso)

We contrast our solution with Cross-Validation-tuned Lasso using `cv.glmnet()` in the `glmnet` package.

    library(glmnet)
    #> Loading required package: Matrix
    #> Loaded glmnet 4.1-8

### Runtime of CV Lasso

    ptm2 <- proc.time()
    fit.glmnet <- cv.glmnet(x, y)
    proc.time()-ptm2
    #>    user  system elapsed
    #>   0.191   0.007   0.199

So, autotune Lasso shows faster runtimes as compared to CV Lasso.

## Comparison of quality of tuning between autotune and CV Lasso

Now, we will visualize the quality of lambdas selected by autotune and CV with respect to the 10-fold CV MSE and true Relative MSE.

    plot(fit.glmnet, cex.lab = 1.5, ylab = "10-fold CV MSE")

    # plotting log of lambdas selected by CV
    abline(v = log(c(fit.glmnet$lambda.min, fit.glmnet$lambda.1se)), lty = "dashed",col = c(rgb(0,1,0), rgb(0,0.75,0)), lwd = 3)

    # plotting log of lambda selected by autotune
    abline(v = log(lambda.autotune), col = "blue", lty = "dashed", lwd = 3)
    legend("topleft", inset = 0.05,             legend = c("Autotune LASSO",
                 "CV_LASSO(min)",
                 "CV_LASSO(1se)"),
      col = c("blue", rgb(0,1,0), rgb(0,0.75,0)),
      lty = "dashed", lwd = 2, cex = 1)

    mse_glmnet <- apply(fit.glmnet$glmnet.fit$beta, 2, function(x) mean((x- beta)^2)/mean(beta^2))
    ymin = 0.95 * min(mse_glmnet)
    ymax = max(1, max(mse_glmnet))
    oldpar <- par(mgp = c(3.5, 1, 0), mar = c(5, 4, 2.5, 0.5) + 1.4)
    plot(log(fit.glmnet$lambda), mse_glmnet,
         type = 'b', col = "red", cex.lab = 2, cex=2, cex.axis = 2,
         ylim = c(ymin, ymax),
         ylab = "Relative MSE w.r.t. true coefficients", xlab = expression(paste("log(",lambda,")")))

    # plotting log of lambdas selected by CV
    abline(v = log(c(fit.glmnet$lambda.min, fit.glmnet$lambda.1se)), lty = "dashed",col = c(rgb(0,1,0), rgb(0,0.75,0)), lwd = 3)

    # plotting log of lambdas selected by autotune
    abline(v = log(lambda.autotune), col = "blue", lty = "dashed", lwd = 3)
    legend(
      "topleft", inset = 0.05,
      legend = c("Autotune LASSO",
                 "CV_LASSO(min)",
                 "CV_LASSO(1se)"),
      col = c("blue", rgb(0,1,0), rgb(0,0.75,0)),
      lty = "dashed", lwd = 2, cex = 1)

    par(oldpar)

Across the lambda grid, we plot RMSE of solution path taken by different tuners.

    temp <- fit.autotune$CD.path.details$lambda0
    seq.lambdas.autotune <- temp * c(var(y), sigma.seq.autotune)

    final.rmse.autotune <- (mean(b.autotune - beta)^2) / (mean(beta)^2)
    intermediate.rmse.autotune <- approx(x = fit.glmnet$lambda, y = mse_glmnet, xout = seq.lambdas.autotune)
    intermediate.rmse.autotune$y[1] <- 1
    rmse.path.autotune <- c(intermediate.rmse.autotune$y, final.rmse.autotune)

    oldpar <- par(mgp = c(3.5, 1, 0), mar = c(5, 4, 2.5, 0.5) + 1.4)
    plot(log(fit.glmnet$lambda), mse_glmnet,
         type = 'b', col = "red", cex.lab = 2, cex=2, cex.axis = 2,
         ylim = c(min(ymin, 0.95 * final.rmse.autotune), ymax),
         ylab = "Relative MSE w.r.t. the true coefficients", xlab = expression(paste("log(",lambda,")")))

    lines(c(log(seq.lambdas.autotune), log(fit.autotune$lambda)),
          rmse.path.autotune,
          col = "blue",
          lwd = 2)
    points(log(seq.lambdas.autotune),
           intermediate.rmse.autotune$y,
           col = "blue",
           pch = 1,
           lwd = 2,
           cex = 2)
    points(log(fit.autotune$lambda),
          final.rmse.autotune,
          col = "blue",
          pch = 11,
          lwd = 3,
          cex = 3
        )

    par(oldpar)

    library(AUC)
    #> AUC 0.3.2
    #> Type AUCNews() to see the change log and ?AUC to get an overview.

    b.glmnet = coef(fit.glmnet, s = "lambda.min")[-1]

    auc(roc(abs(b.autotune), as.factor(beta != 0)))
    #> [1] 1
    auc(roc(abs(b.glmnet), as.factor(beta != 0)))
    #> [1] 1

## Comparison on real data provided in Scaled Lasso’s R package `scalreg`

The `scalreg` package was removed from CRAN in July 2026.

### Comparing prediction errors of autotune and benchmarks: CV and Scaled Lasso

    data("sp500")
    attach(sp500)

    scaled_lasso_benchmark <- read.csv(
      system.file(
        "extdata",
        "sp500-scalreg-1.0.1.csv",
        package = "autotune"
      )
    )
    scaled_lasso_observed <- scaled_lasso_benchmark[
      which(!scaled_lasso_benchmark$bootstrap),
      ,
      drop = FALSE
    ]

    n <- 190
    rang <- 252 - n

    X_train = sp500.percent[1:n, 3:(dim(sp500.percent)[2])]
    Y_train = sp500.percent[1:n, 1]

    x_train = scale(X_train)
    y_train = Y_train - mean(Y_train)

    X_test = sp500.percent[(1:rang) + n, 3: (dim(sp500.percent)[2])]
    Y_test = sp500.percent[(1:rang) + n, 1]
    x_test <- scale(X_test)
    y_test <- Y_test - mean(Y_test)
    # Default value of beta_iter_max is 40
    ans_autotune <- autotune_lasso(x_train, y_train, beta_iter_max = 40, trace_it = TRUE)
    #> Iteration: 1Iteration: 2Lambda converged, Iteration: 3Lambda converged, Iteration: 4Lambda converged, Iteration: 5Lambda converged, Iteration: 6Lambda converged, Iteration: 7Lambda converged, Iteration: 8Lambda converged, Iteration: 9Lambda converged, Iteration: 10Lambda converged, Iteration: 11Lambda converged, Iteration: 12Lambda converged, Iteration: 13Lambda converged, Iteration: 14Lambda converged, Iteration: 15Lambda converged, Iteration: 16Lambda converged, Iteration: 17Lambda converged, Iteration: 18Lambda converged, Iteration: 19Lambda converged, Iteration: 20Lambda converged, Iteration: 21Lambda converged, Iteration: 22Lambda converged, Iteration: 23Lambda converged, Iteration: 24Lambda converged, Iteration: 25Lambda converged, Iteration: 26Lambda converged, Iteration: 27Lambda converged, Iteration: 28Lambda converged, Iteration: 29Lambda converged, Iteration: 30Lambda converged, Iteration: 31Lambda converged, Iteration: 32Lambda converged, Iteration: 33Lambda converged, Iteration: 34Lambda converged, Iteration: 35Lambda converged, Iteration: 36Lambda converged, Iteration: 37Lambda converged, Iteration: 38Lambda converged, Iteration: 39Lambda converged, Iteration: 40Lambda converged, Iteration: 41Lambda converged, Iteration: 42
    #>
    #> No of predictors significant for sigma estimation: 2
    pred_err_autotune <- mean( (y_test - x_test %*% ans_autotune$beta)^2 )

    pred_err_scallas <- scaled_lasso_observed$test_mse

    cv_fit <- cv.glmnet(x_train, y_train, alpha = 1, intercept = F)
    pred_err_cv_min <- mean( (y_test - predict(cv_fit, newx = x_test, s = "lambda.min"))^2 )
    pred_err_cv_1se <- mean( (y_test - predict(cv_fit, newx = x_test, s = "lambda.1se"))^2 )

    pred_err_autotune
    #> [1] 0.5001709
    pred_err_scallas
    #> [1] 0.712887
    pred_err_cv_min
    #> [1] 0.7034008
    pred_err_cv_1se
    #> [1] 0.6785715

    scaled_lasso_observed$nonzero_coefs
    #> [1] 64
    sum(ans_autotune$beta != 0)
    #> [1] 47
    sum(coef(cv_fit, s = "lambda.min") != 0)
    #> [1] 60
    sum(coef(cv_fit, s = "lambda.1se") != 0)
    #> [1] 55
    detach(sp500)

### Comparing via bootstrap

    real_life_prediction_comparison <- function(x, y, x_test, y_test, scaled_lasso_result, alpha = 0.01, beta_iter_max = 40, plot_regu = TRUE, trace_it = TRUE){

    ans_autotune <- autotune_lasso(x, y, alpha = alpha, beta_iter_max = beta_iter_max, trace_it = trace_it)
    # ans_glmnet <- glmnet(x, y, alpha = 1, intercept = FALSE)
    cv_fit <- cv.glmnet(x, y, alpha = 1, intercept = F)
    vary <- var(y_test)
    mse_glmnet <- apply(cv_fit$glmnet.fit$beta, 2, function(b) mean((y_test - x_test%*%b)^2)/vary)

    # log_glmnet_lambdas <- log(fit.glmnet$lambda)
    # log_autotune_lambdas <- log(intermediate_lambdas_used)

    if(plot_regu) {
      intermediate_lambdas_used <- c( c(var(y)/2, ans_autotune$CD.path.details$sigma_sq_seq) * ans_autotune$CD.path.details$lambda0)

      autotune_lambdas <- c(intermediate_lambdas_used, ans_autotune$lambda)

      final.rmse.autotune <- mean((y_test - x_test%*%ans_autotune$beta)^2) / vary
      intermediate.rmse.autotune <- approx(x = log(cv_fit$lambda), y = mse_glmnet, xout = log(intermediate_lambdas_used), rule = 2)
      intermediate.rmse.autotune$y[1] <- 1
      rmse.path.autotune <- c(intermediate.rmse.autotune$y, final.rmse.autotune)

      ysquaredmean <- mean(y^2)
      cv_fit$cvm <- cv_fit$cvm/ysquaredmean
      cv_fit$cvsd <- cv_fit$cvsd/ysquaredmean
      cv_fit$cvup <- cv_fit$cvm + cv_fit$cvsd
      cv_fit$cvlo <- cv_fit$cvm - cv_fit$cvsd

      cv_fit_lowerlim <- min(cv_fit$cvm - cv_fit$cvsd, mse_glmnet)
      cv_fit_upperlim <- max(cv_fit$cvm + cv_fit$cvsd, mse_glmnet)

      ymin = 0.95 * min(cv_fit_lowerlim , min(rmse.path.autotune))
      ymax = max(cv_fit_upperlim, rmse.path.autotune)
      ymax_trunc = ymin + 0.8 * (ymax - ymin)  # Truncate vertical lines here

      ydiff = ymax - ymin
      oldpar <- par(mgp = c(3.5, 1, 0), mar = c(6, 5, 4, 2) + 0.1)
      on.exit(par(oldpar), add = TRUE)

      par(mar = c(5, 7, 4, 2))
      plot(cv_fit,
           ylim = c(ymin, max(cv_fit_upperlim, rmse.path.autotune)),
           cex.lab = 2,
           cex.axis = 2,
           ylab = expression(
        "||" * (Y[test] - X[test] * hat(beta)) * "||"[2]^2
        ~ "/" ~
        "||" * Y[test]  * "||"[2]^2
      ))
      # expression(paste(sqrt("RMSE")))

      points(log(cv_fit$lambda), mse_glmnet,
             type = 'b', col = "orange", lwd = 3)

      lines(rep(log(ans_autotune$lambda), 2),
            c(ymin, ymin + 0.92 * ydiff), col = "blue", lwd = 0.5, lty = "dashed")

      lines(log(autotune_lambdas), rmse.path.autotune,
            col = "blue", lwd = 1.5)

      points(log(autotune_lambdas), rmse.path.autotune,
             col = "blue", pch = 1, lwd = 1.5, cex = 1)

      points(log(ans_autotune$lambda), final.rmse.autotune,
             col = "blue", pch = 11, lwd = 2, cex = 2)

      vlines <- c(
        log(cv_fit$lambda.min),
        log(cv_fit$lambda.1se)
      )
      vline_labels <- c("CV(min)", "CV(1se)")
      vline_colors <- c(rgb(0,1,0), rgb(0,0.75,0))
      vline_ltys <- c("dashed", "dashed")
      vline_lwds <- c(3, 3)

      for (i in seq_along(vlines)) {
        lines(rep(vlines[i], 2),
              c(ymin, ymin + (0.55 + i * 0.1) * ydiff),
              col = vline_colors[i], lty = vline_ltys[i],
              lwd = vline_lwds[i])

        text(x = vlines[i], y = ymin + (0.57 + i * 0.1) * ydiff,
             labels = vline_labels[i], srt = -30, adj = 0,
             xpd = TRUE, cex = 1.4, col = vline_colors[i])
      }

      text(x = log(ans_autotune$lambda),
           y =  ymin + 0.95 * ydiff,
           labels = "Autotune", srt = -30, adj = 0,
           xpd = TRUE, cex = 1.4, col = "blue")

      if (!is.null(cv_fit$cvsd)) {
        arrows(x0 = log(cv_fit$lambda),
               y0 = cv_fit$cvm - cv_fit$cvsd,
               x1 = log(cv_fit$lambda),
               y1 = cv_fit$cvm + cv_fit$cvsd,
               angle = 90, code = 3, length = 0.05, col = "red", lwd = 1.5)
      }

      legend("topright",
             inset = c(0.05, 0.02),
             legend = c("CV Error", "Test Error"),
             col = c("red", "orange"),
             pt.cex = c(1.5, 2), pch = c(19, 19),
             lwd = c(1, 3),
             text.col = c("red", "orange"),
             horiz = FALSE,
             bty = "n")

    }

    # mse_autotune <- mean( (y_test - x_test %*% ans_autotune$beta)^2 )
    # mse_cv_min <- mean( (y_test - predict(cv_fit, newx = x_test, s = "lambda.min"))^2 )
    # mse_cv_1se <- mean( (y_test - predict(cv_fit, newx = x_test, s = "lambda.1se"))^2 )

    oos_scallas <- scaled_lasso_result$out_of_sample_r2
    oos_autotune <- 1 -  mean( (y_test - x_test %*% ans_autotune$beta)^2 ) / vary
    oos_cv_min <- 1 - mean( (y_test - predict(cv_fit, newx = x_test, s = "lambda.min"))^2 ) / vary
    oos_cv_1se <- 1 - mean( (y_test - predict(cv_fit, newx = x_test, s = "lambda.1se"))^2 ) / vary

    return(data.frame(out_of_sample_autotune = oos_autotune,
    out_of_sample_scaled_lasso = oos_scallas,
    out_of_sample_cv_min = oos_cv_min,
    out_of_sample_cv_1se = oos_cv_1se,
    nonzero_coefs_autotune = sum(ans_autotune$beta != 0),
    nonzero_coefs_scaled_lasso = scaled_lasso_result$nonzero_coefs,
    nonzero_coefs_cvmin = sum(coef(cv_fit, s = "lambda.min") != 0),
    nonzero_coefs_cv1se = sum(coef(cv_fit, s = "lambda.1se") != 0)))
    }

    ans <- real_life_prediction_comparison(
      x_train,
      y_train,
      x_test,
      y_test,
      scaled_lasso_result = scaled_lasso_observed,
      plot_regu = TRUE
    )
    #> Iteration: 1Iteration: 2Lambda converged, Iteration: 3Lambda converged, Iteration: 4Lambda converged, Iteration: 5Lambda converged, Iteration: 6Lambda converged, Iteration: 7Lambda converged, Iteration: 8Lambda converged, Iteration: 9Lambda converged, Iteration: 10Lambda converged, Iteration: 11Lambda converged, Iteration: 12Lambda converged, Iteration: 13Lambda converged, Iteration: 14Lambda converged, Iteration: 15Lambda converged, Iteration: 16Lambda converged, Iteration: 17Lambda converged, Iteration: 18Lambda converged, Iteration: 19Lambda converged, Iteration: 20Lambda converged, Iteration: 21Lambda converged, Iteration: 22Lambda converged, Iteration: 23Lambda converged, Iteration: 24Lambda converged, Iteration: 25Lambda converged, Iteration: 26Lambda converged, Iteration: 27Lambda converged, Iteration: 28Lambda converged, Iteration: 29Lambda converged, Iteration: 30Lambda converged, Iteration: 31Lambda converged, Iteration: 32Lambda converged, Iteration: 33Lambda converged, Iteration: 34Lambda converged, Iteration: 35Lambda converged, Iteration: 36Lambda converged, Iteration: 37Lambda converged, Iteration: 38Lambda converged, Iteration: 39Lambda converged, Iteration: 40Lambda converged, Iteration: 41Lambda converged, Iteration: 42
    #>
    #> No of predictors significant for sigma estimation: 2

    names(ans)
    #> [1] "out_of_sample_autotune"     "out_of_sample_scaled_lasso"
    #> [3] "out_of_sample_cv_min"       "out_of_sample_cv_1se"
    #> [5] "nonzero_coefs_autotune"     "nonzero_coefs_scaled_lasso"
    #> [7] "nonzero_coefs_cvmin"        "nonzero_coefs_cv1se"

    library(dplyr)
    #>
    #> Attaching package: 'dplyr'
    #> The following objects are masked from 'package:stats':
    #>
    #>     filter, lag
    #> The following objects are masked from 'package:base':
    #>
    #>     intersect, setdiff, setequal, union

    #Number of bootstraps
    B <- 20

    result_list <- list()

    result_list[[1]] <- data.frame(
      n = n,
      bootstrap = FALSE,
      real_life_prediction_comparison(
        x_train,
        y_train,
        x_test,
        y_test,
        scaled_lasso_result = scaled_lasso_observed,
        plot_regu = FALSE,
        trace_it = FALSE
      )[1:4]
    )

    for( j in (1:B)+1) {
      set.seed(j)
      boot.sample <- sample(n, size = n, replace = TRUE)
      boot.x <- x_train[boot.sample, ]
      boot.y <- y_train[boot.sample]
      scaled_lasso_result <- scaled_lasso_benchmark[
        which(scaled_lasso_benchmark$seed == j),
        ,
        drop = FALSE
      ]
      result_list[[j]] <- c(
        n = n,
        bootstrap = TRUE,
        real_life_prediction_comparison(
          boot.x,
          boot.y,
          x_test,
          y_test,
          scaled_lasso_result = scaled_lasso_result,
          plot_regu = FALSE,
          trace_it = FALSE
        )[1:4]
      )
    }

    mat <- do.call(rbind, result_list)
    df <- as.data.frame(mat, stringsAsFactors = FALSE)

    colnames(df) <- c("n", "bootstrap", "Autotune",
                      "Scaled", "CVmin",
                      "CV1se")
    df <- df %>%
      mutate(
        n = as.integer(n),
        bootstrap = as.logical(bootstrap),
        across(!c(n, bootstrap), as.numeric)
      )

    alg_cols <- colnames(df)[3:6]
    boot_only <- df %>% filter(bootstrap == TRUE)
    boxlist <- lapply(alg_cols, function(cn) boot_only[[cn]])
    names(boxlist) <- alg_cols

    # Uncomment the following line and line 448 for saving the boxplot in your working directory

    # pdf(file.path(paste0("OOS_boxplot_for_n=",df[1,1],"_B=", B,".pdf")), width = 8, height = 9)

    oldpar <- par(mar = c(4, 6, 4, 2))
    boxplot(boxlist, main = paste0("Bootstrapped Out-of-Sample R^2 for n = ", df[1,1]), ylab = expression(paste("Out-of-Sample ",R^2)), ylim = c(max(0.98 * min(df[alg_cols]), 0), 1),
            cex.lab = 2,
      cex.axis = 1.8,
      cex.main = 1.5)

    nonboot <- df %>% filter(bootstrap == FALSE)
    if(nrow(nonboot) > 0){
      points(1:4, as.numeric(nonboot[1, alg_cols]), pch = 18, col = "red", cex = 2)
    }

    legend(
      "topright",
      legend = expression(paste("Observed O.O.S. ",R^2)),
      col = "red",
      pch = 18,
      pt.cex = 2,
      cex = 2,
      bty = "n"
    )

    par(oldpar)
    # dev.off()