Is ridge regression useless in high dimensions ($n \ll p$)? How can OLS fail to overfit? Consider a good old regression problem with $p$ predictors and sample size $n$. The usual wisdom is that OLS estimator will overfit and will generally be outperformed by the ridge regression estimator: $$\hat\beta = (X^\top X + \lambda I)^{-1}X^\top y.$$ It is standard to use cross-validation to find an optimal regularization parameter $\lambda$. Here I use 10-fold CV. Clarification update: when $n<p$, by "OLS estimator" I understand "minimum-norm OLS estimator" given by $$\hat\beta_\text{OLS} = (X^\top X)^+X^\top y = X^+ y.$$
I have a dataset with $n=80$ and $p>1000$. All predictors are standardized, and there are quite a few that (alone) can do a good job in predicting $y$. If I randomly select a small-ish, say $p=50<n$, number of predictors, I get a reasonable CV curve: large values of $\lambda$ yield zero R-squared, small values of $\lambda$ yield negative R-squared (because of overfitting) and there is some maximum in between. For $p=100>n$ the curve looks similar. However, for $p$ much larger than that, e.g. $p=1000$, I do not get any maximum at all: the curve plateaus, meaning that OLS with $\lambda\to 0$ performs as good as ridge regression with optimal $\lambda$.

How is it possible and what does it say about my dataset? Am I missing something obvious or is it indeed counter-intuitive? How can there be any qualitative difference between $p=100$ and $p=1000$ given that both are larger than $n$?
Under what conditions does minimal-norm OLS solution for $n<p$ not overfit?

Update: There was some disbelief in the comments, so here is a reproducible example using glmnet. I use Python but R users will easily adapt the code. 
%matplotlib notebook

import numpy as np
import pylab as plt
import seaborn as sns; sns.set()

import glmnet_python    # from https://web.stanford.edu/~hastie/glmnet_python/
from cvglmnet import cvglmnet; from cvglmnetPlot import cvglmnetPlot

# 80x1112 data table; first column is y, rest is X. All variables are standardized
mydata = np.loadtxt('../q328630.txt')   # file is here https://pastebin.com/raw/p1cCCYBR
y = mydata[:,:1]
X = mydata[:,1:]

# select p here (try 1000 and 100)
p = 1000

# randomly selecting p variables out of 1111
np.random.seed(42)
X = X[:, np.random.permutation(X.shape[1])[:p]]

fit = cvglmnet(x = X.copy(), y = y.copy(), alpha = 0, standardize = False, intr = False, 
               lambdau=np.array([.0001, .001, .01, .1, 1, 10, 100, 1000, 10000, 100000]))
cvglmnetPlot(fit)
plt.gcf().set_size_inches(6,3)
plt.tight_layout()



 A: A natural regularization happens because of the presence of many small components in the theoretical PCA of $x$. These small components are implicitly used to fit the noise using small coefficients. When using minimum norm OLS, you fit the noise with many small independent components and this has a regularizing effect equivalent to Ridge regularization. This regularization is often too strong, and it is possible to compensate it using "anti-regularization" know as negative Ridge. In that case, you will see the minimum of the MSE curve appears for negative values of $\lambda$.
By theoretical PCA, I mean:

Let $x\sim N(0,\Sigma)$ a multivariate normal distribution. There is a
  linear isometry $f$ such as $u=f(x)\sim N(0,D)$ where $D$ is diagonal:
  the components of $u$ are independent. $D$ is simply obtained by diagonalizing $\Sigma$.
Now the model $y=\beta.x+\epsilon$ can be written
  $y=f(\beta).f(x)+\epsilon$ (a linear isometry preserves dot product).
  If you write $\gamma=f(\beta)$, the model can be written
  $y=\gamma.u+\epsilon$. Furthermore $\|\beta\|=\|\gamma\|$ hence
  fitting methods like Ridge or minimum norm OLS are perfectly
  isomorphic: the estimator of $y=\gamma.u+\epsilon$  is the image by $f$
  of the estimator of $y=\beta.x+\epsilon$.

Theoretical PCA transforms non independent predictors into independent predictors. It is only loosely related to empirical PCA where you use the empirical covariance matrix (that differs a lot from the theoretical one with small sample size). Theoretical PCA is not practically computable but is only used here to interpret the model in an orthogonal predictor space.
Let's see what happens when we append many small variance independent predictors to a model:
Theorem
Ridge regularization with coefficient $\lambda$ is equivalent (when $p\rightarrow\infty$) to:


*

*adding $p$ fake independent predictors (centred and identically distributed) each with variance $\frac{\lambda}{p}$ 

*fitting the enriched model with minimum norm OLS estimator

*keeping only the parameters for the true predictors



(sketch of) Proof
We are going to prove that the cost functions are asymptotically
  equal. Let's split the model into real and fake predictors: $y=\beta x+\beta'x'+\epsilon$. The cost function of Ridge (for the true
  predictors) can be written:
$$\mathrm{cost}_\lambda=\|\beta\|^2+\frac{1}{\lambda}\|y-X\beta\|^2$$
When using minimum norm OLS, the response is fitted perfectly: the
  error term is 0. The cost function is only about the norm of the
  parameters. It can be split into the true parameters and the fake
  ones:
$$\mathrm{cost}_{\lambda,p}=\|\beta\|^2+\inf\{\|\beta'\|^2 \mid X'\beta'=y-X\beta\}$$
In the right expression, the minimum norm solution is given by:
$$\beta'=X'^+(y-X\beta )$$
Now using SVD for $X'$:
$$X'=U\Sigma V$$
$$X'^{+}=V^\top\Sigma^{+} U^\top$$
We see that the norm of $\beta'$ essentially depends on the singular
  values of $X'^+$ that are the reciprocals of the singular values of
  $X'$. The normalized version of $X'$ is $\sqrt{p/\lambda} X'$. I've
  looked at literature and singular values of large random matrices are
  well known. For $p$ and $n$ large enough, minimum $s_\min$ and maximum
  $s_\max$ singular values are approximated by (see theorem 1.1):
$$s_\min(\sqrt{p/\lambda}X')\approx \sqrt p\left(1-\sqrt{n/p}\right)$$ 
  $$s_\max(\sqrt{p/\lambda}X')\approx \sqrt p \left(1+\sqrt{n/p}\right)$$
Since, for large $p$, $\sqrt{n/p}$ tends towards 0, we can just say
  that all singular values are approximated by $\sqrt p$. Thus:
$$\|\beta'\|\approx\frac{1}{\sqrt\lambda}\|y-X\beta\|$$
Finally:
$$\mathrm{cost}_{\lambda,p}\approx\|\beta\|^2+\frac{1}{\lambda}\|y-X\beta\|^2=\mathrm{cost}_\lambda$$
Note: it does not matter if you keep the coefficients of the fake
  predictors in your model. The variance introduced by $\beta'x'$ is
  $\frac{\lambda}{p}\|\beta'\|^2\approx\frac{1}{p}\|y-X\beta\|^2\approx\frac{n}{p}MSE(\beta)$.
  Thus you increase your MSE by a factor $1+n/p$ only which tends
  towards 1 anyway. Somehow you don't need to treat the
  fake predictors differently than the real ones.

Now, back to @amoeba's data. After applying theoretical PCA to $x$ (assumed to be normal), $x$ is transformed by a linear isometry into a variable $u$ whose components are independent and sorted in decreasing variance order. The problem $y=\beta x+\epsilon$ is equivalent the transformed problem $y=\gamma u+\epsilon$.
Now imagine the variance of the components look like:

Consider many $p$ of the last components, call the sum of their variance $\lambda$. They each have a variance approximatively equal to $\lambda/p$ and are independent. They play the role of the fake predictors in the theorem.
This fact is clearer in @jonny's model: only the first component of theoretical PCA  is correlated to $y$ (it is proportional $\overline{x}$) and has huge variance. All the other components (proportional to $x_i-\overline{x}$) have comparatively very small variance (write the covariance matrix and diagonalize it to see this) and play the role of fake predictors. I calculated that the regularization here corresponds (approx.) to prior $N(0,\frac{1}{p^2})$ on $\gamma_1$ while the true $\gamma_1^2=\frac{1}{p}$. This definitely over-shrinks. This is visible by the fact that the final MSE is much larger than the ideal MSE. The regularization effect is too strong. 
It is sometimes possible to improve this natural regularization by Ridge. First you sometimes need $p$ in the theorem really big (1000, 10000...) to seriously rival Ridge and the finiteness of $p$ is like an imprecision. But it also shows that Ridge is an additional regularization over a naturally existing implicit regularization and can thus have only a very small effect. Sometimes this natural regularization is already too strong and Ridge may not even be an improvement. More than this, it is better to use anti-regularization: Ridge with negative coefficient. This shows MSE for @jonny's model ($p=1000$), using $\lambda\in\mathbb{R}$:

A: Here is an artificial situation where this occurs. Suppose each predictor variable is a copy of the target variable with a large amount of gaussian noise applied. The best possible model is an average of all predictor variables.
library(glmnet)
set.seed(1846)
noise <- 10
N <- 80
num.vars <- 100
target <- runif(N,-1,1)
training.data <- matrix(nrow = N, ncol = num.vars)
for(i in 1:num.vars){
  training.data[,i] <- target + rnorm(N,0,noise)
}
plot(cv.glmnet(training.data, target, alpha = 0,
               lambda = exp(seq(-10, 10, by = 0.1))))


100 variables behave in a "normal" way: Some positive value of lambda minimizes out of sample error.
But increase num.vars in the above code to 1000, and here is the new MSE path. (I extended to log(Lambda) = -100 to convince myself.

What I think is happening
When fitting a lot of parameters with low regularization, the coefficients are randomly distributed around their true value with high variance.
As the number of predictors becomes very large, the "average error" tends towards zero, and it becomes better to just let the coefficients fall where they may and sum everything up than to bias them toward 0. 
I'm sure this situation of the true prediction being an average of all predictors isn't the only time this occurs, but I don't know how to begin pinpoint the biggest necessary condition here.
EDIT:
The "flat" behavior for very low lambda will always happen, since the solution is converging to the minimum-norm OLS solution. Similarly the curve will be flat for very high lambda as the solution converges to 0. There will be no minimum iff one of those two solution is optimal.
Why is the minimum-norm OLS solution so (comparably) good in this case? I think it is related to the following behavior that I found very counter-intuitive, but on reflection makes a lot of sense.
max.beta.random <- function(num.vars){
  num.vars <- round(num.vars)
  set.seed(1846)
  noise <- 10
  N <- 80
  target <- runif(N,-1,1)
  training.data <- matrix(nrow = N, ncol = num.vars)

  for(i in 1:num.vars){
    training.data[,i] <- rnorm(N,0,noise)
  }
  udv <- svd(training.data)

  U <- udv$u
  S <- diag(udv$d)
  V <- udv$v

  beta.hat <- V %*% solve(S) %*% t(U) %*% target

  max(abs(beta.hat))
}


curve(Vectorize(max.beta.random)(x), from = 10, to = 1000, n = 50,
      xlab = "Number of Predictors", y = "Max Magnitude of Coefficients")

abline(v = 80)


With randomly generated predictors unrelated to the response, as p increases the coefficients become larger, but once p is much bigger than N they shrink toward zero. This also happens in my example. So very loosely, the unregularized solutions for those problems don't need shrinkage because they are already very small!
This happens for a trivial reason. $y$ can be expressed exactly as a linear combination of columns of $X$. $\hat{\beta}$ is the minimum-norm vector of coefficients. As more columns are added the norm of $\hat{\beta}$ must decrease or remain constant, because a possible linear combination is to keep the previous coefficients the same and set the new coefficients to $0$.
A: If you're familiar with linear operators then you may like my answer as most direct path to understanding the phenomenon: why doesn't least norm regression fail outright? The reason is that your problem ($n\ll p$) is the ill posed inverse problem and pseudo-inverse is one of the ways of solving it. Regularization is an improvement though.
This paper is probably the most compact and relevant explanation: Lorenzo Rosasco et al, Learning, Regularization and Ill-Posed Inverse Problems. They set up your regression problem as learning, see Eq.3., where the number of parameters exceeds the number of observations:
$$Ax=g_\delta,$$ where $A$ is a linear operator on Hilbert space and $g_\delta$ - noisy data. 
Obviously, this is an ill-posed inverse problem. So, you can solve it with SVD or Moore-Penrose inverse, which would render the least norm solution indeed. Thus it should not be surprising that your least norm solution is not failing outright.
However, if you follow the paper you can see that the ridge regression would be an improvement upon the above. The improvement is really a better behavior of the estimator, since Moore-Penrose solution is not necessarily bounded.
UPDATE
I realized that I wasn't making it clear that ill-posed problems lead to overfitting. Here's the quote from the paper Gábor A, Banga JR. Robust and efficient parameter estimation in dynamic models of biological systems. BMC Systems Biology. 2015;9:74. doi:10.1186/s12918-015-0219-2:

The ill-conditioning of these problems typically arise from (i) models
  with large number of parameters (over-parametrization), (ii)
  experimental data scarcity and (iii) significant measurement errors
  [19, 40]. As a consequence, we often obtain overfitting of such
  kinetic models, i.e. calibrated models with reasonable fits to the
  available data but poor capability for generalization (low predictive
  value)

So, my argument can be stated as follows:


*

*ill posed problems lead to overfitting

*(n < p) case is an extremely ill-posed inverse problem

*Moore-Penrose psudo-inverse (or other tools like SVD), which you refer to in the question as $X^+$, solves an ill-posed problem

*therefore, it takes care of overfitting at least to some extent, and it shouldn't be surprising that it doesn't completely fail, unlike a regular OLS should 


Again, regularization is a more robust solution still.
