So I decided to run nested cross-validation using the specialized mlr package in R to see what's actually coming from the modelling approach.
#Code (it takes a few minutes to run on an ordinary notebook)
library(mlr)
daf = read.csv("https://pastebin.com/raw/p1cCCYBR", sep = " ", header = FALSE)
tsk = list(
tsk1110 = makeRegrTask(id = "tsk1110", data = daf, target = colnames(daf)[1]),
tsk500 = makeRegrTask(id = "tsk500", data = daf[, c(1,sample(ncol(daf)-1, 500)+1)], target = colnames(daf)[1]),
tsk100 = makeRegrTask(id = "tsk100", data = daf[, c(1,sample(ncol(daf)-1, 100)+1)], target = colnames(daf)[1]),
tsk50 = makeRegrTask(id = "tsk50", data = daf[, c(1,sample(ncol(daf)-1, 50)+1)], target = colnames(daf)[1]),
tsk10 = makeRegrTask(id = "tsk10", data = daf[, c(1,sample(ncol(daf)-1, 10)+1)], target = colnames(daf)[1])
)
lrn = makeLearner("regr.glmnet", alpha = 0)
rdesc = makeResampleDesc("CV", iters = 10)
msrs = list(mse, rsq)
ctrl = makeTuneControlGrid(resolution = 15L)
pset = makeParamSet( makeNumericParam("lambda", lower = -10, upper = 10, trafo = function(x) exp(x)))
cvl = makeTuneWrapper(learner = lrn, resampling = rdesc, measures = msrs, control = ctrl, par.set = pset)
bm = benchmark(learners = cvl, tasks = tsk, resamplings = rdesc, measures = msrs)
#Results
> bm task.id learner.id mse.test.mean rsq.test.mean 1 tsk10 regr.glmnet.tuned 0.9029739 -0.34385207 2 tsk100 regr.glmnet.tuned 0.9860033 -0.20159950 3 tsk1110 regr.glmnet.tuned 0.6836474 0.08630208 4 tsk50 regr.glmnet.tuned 0.9159190 -0.15177725 5 tsk500 regr.glmnet.tuned 0.8874529 0.04095012
They basically got the same slightly better than random performance, but glmnet is good enough to squeeze some information the more features you feed it.
So, what about the optimal lambdas?
getBMRTuneResults(bm, as.df = TRUE, task.ids = "tsk1110")
# task.id learner.id iter lambda mse.test.mean rsq.test.mean
#1 tsk1110 regr.glmnet.tuned 1 0.239651 0.7234743 0.06932993
#2 tsk1110 regr.glmnet.tuned 2 0.239651 0.6318552 0.10130075
#3 tsk1110 regr.glmnet.tuned 3 0.239651 0.7202225 0.09277869
#4 tsk1110 regr.glmnet.tuned 4 4.172734 0.7169810 -0.08232901
#5 tsk1110 regr.glmnet.tuned 5 0.239651 0.7163615 -0.05643168
#6 tsk1110 regr.glmnet.tuned 6 0.239651 0.6655465 0.13894617
#7 tsk1110 regr.glmnet.tuned 7 0.239651 0.7125394 0.14249063
#8 tsk1110 regr.glmnet.tuned 8 0.239651 0.7312063 -0.56975408
#9 tsk1110 regr.glmnet.tuned 9 0.239651 0.6895064 -0.03308516
#10 tsk1110 regr.glmnet.tuned 10 0.239651 0.6370360 0.02805295
Notice the lambdas are already transformed. Not a single fold picked the minimal lambda $\exp(-10)\approx 4.54 \mathbf E-5$ (this is the same behavior with cv.glmnet in R).
I fiddled a bit more with glmnet and discovered neither there the minimal lambda is picked. Check:
cvfit = cv.glmnet(x = x, y = y, alpha = 0, lambda = exp(seq(-10, 10, length.out = 150)))
plot(cvfit)
#Conclusion
So, basically, $\lambda>0$ really improves the fit.
How is it possible and what does it say about my dataset? Am I missing something obvious or is it indeed counter-intuitive?
We are likely nearer the true distribution of the data setting $\lambda$ to a small value larger than zero. There's nothing counter-intuitive about it though.
EDIT: The lambda selection does say something more about your data. As larger lambdas decrease performance, it means there are preferential, i.e. larger, coefficients in your model, as large lambdas shrink all coefficients towards zero. Though $\lambda\neq0$ means that the effective degrees of freedom in your model is smaller than the apparent degrees of freedom, $p$.
How can there be any qualitative difference between p=100 and p=1000 given that both are larger than n?
$p=1000$ invariably contains at least the same of information or even more than $p=100$.
#Comments
It seems you are getting a tiny minimum for some non-zero lambda (I am looking at your figure), but the curve is still really really flat to the left of it. So my main question remains as to why λ→0 does not noticeably overfit. I don't see an answer here yet. Do you expect this to be a general phenomenon? I.e. for any data with n≪p, lambda=0 will perform [almost] as good as optimal lambda? Or is it something special about these data? If you look above in the comments, you'll see that many people did not even believe me that it's possible.
I think you're conflating validation performance with test performance, and such comparison is not warranted.
We picked the best lambda in the most parsimonious way, and we expect it to be the more performant than any other arbitrary selection. Repeating the benchmark with the addition of a glmnet with lambda set to a ludicrously small value (1E-100), we see that our tuning really improves the result.
bm2 = benchmark(learners = list(cvl, makeLearner("regr.glmnet", alpha = 0, lambda = 1E-100)), tasks = tsk, resamplings = rdesc, measures = msrs)
getBMRAggrPerformances(bm2, as.df = TRUE)
# task.id learner.id mse.test.mean rsq.test.mean
#1 tsk10 regr.glmnet.tuned 1.0302354 -0.36840610
#2 tsk10 regr.glmnet 1.2542921 -0.58709603
#3 tsk100 regr.glmnet.tuned 0.8721003 -0.01025020
#4 tsk100 regr.glmnet 2.6139448 -2.56834215
#5 tsk1110 regr.glmnet.tuned 0.6716137 0.26026566
#6 tsk1110 regr.glmnet 1.4569493 -0.76811536
#7 tsk50 regr.glmnet.tuned 0.7851516 0.02185949
#8 tsk50 regr.glmnet 1.8970910 -1.69729756
#9 tsk500 regr.glmnet.tuned 0.8556580 -0.05821649
#10 tsk500 regr.glmnet 1.6703070 -1.02238746
There's an element of randomness involved, but the tuned model always outperform the asymptotically-OLS ridge, often by quite a marging.
Also, I don't quite understand your last line. Look at the cv.glmnet output for p=100. It will have very different shape. So what affects this shape (asymptote on the left vs. no asymptote) when p=100 or p=1000?
Let's compare the regularization paths for both:
fit1000 = glmnet(x, y, alpha = 0, lambda = exp(seq(-10,10, length.out = 1001)))
fit100 = glmnet(x[, sample(1000, 100)], y, alpha = 0, lambda = exp(seq(-10,10, length.out = 1001)))
plot(fit1000, "lambda")
x11()
plot(fit100, "lambda")
It becomes clear $p=1000$ affords larger coefficients at increasing $\lambda$, even though it has smaller coefficients for asymptotically-OLS ridge, at the left of both plots. So, basically, $p=100$ overfits at the left of the graph, and that probably explains the difference in behavior between them.
It's harder for $p=1000$ to overfit because, even though Ridge shrinks coefficients to zero, they are never reach zero. This mean that the predictive power of the model is shared between many more components, making it easier to predict around the mean instead of being carried away by noise.