Elastic net: How to improve unstable cross validation of lambda I want to tune lambda of an elastic net. It is possible to optimize lambda based on a cross validation with the glmnet package. However, I recognized that the estimated lambda varies a lot, if I run the cross validation several times. See code below, for a reproducible illustration of this situation:
library("glmnet")

set.seed(1234)

# Some example data
N <- 1000
y <- rnorm(N, 5, 10)
x1 <- y + rnorm(N, 2, 10)
x2 <- y + rnorm(N, - 5, 20)
x3 <- y + rnorm(N, 10, 200)
x4 <- rnorm(N, 20, 50)
x5 <- rnorm(N, - 7, 200)
x6 <- rbinom(N, 1, exp(x1) / (exp(x1) + 1))
x7 <- rbinom(N, 1, exp(x2) / (exp(x2) + 1))
x8 <- rbinom(N, 1, exp(x3) / (exp(x3) + 1))
x9 <- rbinom(N, 1, exp(x4) / (exp(x4) + 1))
x10 <- rbinom(N, 1, exp(x5) / (exp(x5) + 1))

data <- data.frame(y, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10)

# Estimate lambda.1se several times
est_lambda.1se <- numeric()

for(i in 1:1000) {
  # Tuning in glmnet
  mod_cv_test <- cv.glmnet(x = as.matrix(data[ , colnames(data) %in% "y" == FALSE]), 
                       y = y, alpha = 0.5, family = "gaussian")
  # Store lambda.1se
  est_lambda.1se[i] <- mod_cv_test$lambda.1se
}

# Distribution of the estimated lambdas
hist(est_lambda.1se, breaks = 50) # Could the median be taken?


In my example, lambda ranges approximately between 1 and 2.3. I am wondering if there is room for improvements. My idea is to take the median of all runs as best lambda. However, I am afraid that this procedure results in any problems which I can not see right now.
Question: How could the tuning of lambda be improved?
 A: I think what you're really asking here is about how sensitive the results are to your choice of $\lambda$, so I suggest you investigate exactly that. When I'm tuning an elastic net, if computation time isn't an issue I like to make a heatmap of the CV loss as a function of $\lambda_1$ and $\lambda_2$ (or $\lambda$ and $\alpha$, depending on your software's parameterization). This exactly tells me the stability of my solution and I can see how much better my optimized point is than a 'random' point. You can also 'zoom in' until the noise of cross validation overpowers the actual increase in model performance (this is like doing a heuristic simplex search). Statements like 'this $\lambda$ is large' aren't really meaningful -- it's all data dependent, so for your data the difference between $\lambda=1$ and $\lambda=2$ might be completely negligible, while for other data sets that's a massive difference.
Here's an example of the kind of heatmap that I'm talking about, taken from this paper (which I suggest you read) by our very own @DikranMarsupial:

A: Good question! 
You can try the train function in the caret package, it can easily tuning the alpha and lambda parameters at the same time.
What you caught at this unstable situation, it telling you maybe the falling in the local minimum.
I have run this example multiple times, this is what i got:
The local minimum of lambda is 0.201, then I use this new lambda seq(0.01,0.4,0.01) which I got the 'global' minimum is 0.07!
Run the code bellow, it's very impressive!

eGrid <- expand.grid(alpha = (1:10) * 0.1, 
                 lambda = seq(0.01,0.4,0.01))

Control <- trainControl(method = "repeatedcv",repeats = 3,verboseIter =TRUE)

netFit <- train(x = as.matrix(data[ , colnames(data) %in% "y" == FALSE]), 
            y = y,
            method = "glmnet",
            tuneGrid = eGrid,
            trControl = Control)
plot(netFit)

