---
title: "autotune Lasso"
output: rmarkdown::html_vignette
# output: rmarkdown::pdf_document
vignette: >
  %\VignetteIndexEntry{Autotune Lasso}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

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
```{r install-github, eval = FALSE}
# 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](https://cran.r-project.org/bin/windows/Rtools/) on Windows, or [these tools](https://mac.r-project.org/tools/) on Mac).


### Load Package
After installation, the package can be loaded in the standard way:
```{r setup}
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.

```{r}
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$.

```{r}
ptm <- proc.time()
fit.autotune <- autotune_lasso(x, y, alpha = 0.01)
proc.time() - ptm
```

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
```{r}
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

```{r}
head(coef(fit.autotune), n = 21)
```



## 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.
```{r}
library(glmnet)
```


### Runtime of CV Lasso
```{r}
ptm2 <- proc.time()
fit.glmnet <- cv.glmnet(x, y)
proc.time()-ptm2
```
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.
```{r, fig.width=7.5, fig.height=7.5}
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)
```

```{r, fig.width=7.5, fig.height=7.5}
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.

```{r, fig.width=7.5, fig.height=7.5}
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)
```

```{r}
library(AUC)
```

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


```{r}
auc(roc(abs(b.autotune), as.factor(beta != 0)))
auc(roc(abs(b.glmnet), as.factor(beta != 0)))
```

## 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
```{r}
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)
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
pred_err_scallas
pred_err_cv_min
pred_err_cv_1se

scaled_lasso_observed$nonzero_coefs
sum(ans_autotune$beta != 0)
sum(coef(cv_fit, s = "lambda.min") != 0)
sum(coef(cv_fit, s = "lambda.1se") != 0)
detach(sp500)
```


### Comparing via bootstrap
```{r}
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)))
}
```


```{r, fig.width=7.5, fig.height=7.5}
ans <- real_life_prediction_comparison(
  x_train,
  y_train,
  x_test,
  y_test,
  scaled_lasso_result = scaled_lasso_observed,
  plot_regu = TRUE
)
names(ans)
```

```{r}
library(dplyr)
```



```{r, fig.width=7.5, fig.height=7.5}
#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)
  )
```



```{r, fig.width=7.5, fig.height=7.5}
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()
```