¿Es inútil la regresión de crestas en dimensiones altas ( $n \ll p$ )? ¿Cómo puede fallar el OLS en el sobreajuste?

Question

Consideremos un viejo problema de regresión con $p$ predictores y tamaño de la muestra $n$ . La sabiduría habitual es que el estimador OLS se ajustará en exceso y generalmente será superado por el estimador de regresión de cresta: $$\hat\beta = (X^\top X + \lambda I)^{-1}X^\top y.$$ Es habitual utilizar la validación cruzada para encontrar un parámetro de regularización óptimo $\lambda$ . Aquí utilizo el CV de 10 veces. Actualización de la aclaración: cuando $n<p$ Por "estimador OLS" entiendo el "estimador OLS de norma mínima" dado por $$\hat\beta_\text{OLS} = (X^\top X)^+X^\top y = X^+ y.$$

Tengo un conjunto de datos con $n=80$ y $p>1000$ . Todos los predictores están estandarizados, y hay bastantes que (por sí solos) pueden hacer un buen trabajo de predicción $y$ . Si yo al azar seleccione una pequeña, digamos $p=50<n$ número de predictores, obtengo una curva de CV razonable: los valores grandes de $\lambda$ producen un R-cuadrado nulo, los valores pequeños de $\lambda$ producen un R-cuadrado negativo (debido al sobreajuste) y hay algún máximo en el medio. En $p=100>n$ la curva es similar. Sin embargo, para $p$ mucho más grande que eso, por ejemplo $p=1000$ no obtengo ningún máximo: la curva se estabiliza, lo que significa que OLS con $\lambda\to 0$ tiene un rendimiento tan bueno como la regresión de cresta con $\lambda$ .

¿Cómo es posible y qué dice sobre mi conjunto de datos? ¿Me estoy perdiendo algo obvio o es realmente contraintuitivo? ¿Cómo puede haber alguna diferencia cualitativa entre $p=100$ y $p=1000$ dado que ambos son mayores que $n$ ?

¿En qué condiciones la solución OLS de norma mínima para $n<p$ no ¿sobreajuste?

---

Actualización: Hubo cierta incredulidad en los comentarios, así que aquí hay un ejemplo reproducible utilizando glmnet . Yo uso Python pero los usuarios de R adaptarán fácilmente el código.

%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()

Answer 1

Una regularización natural ocurre debido a la presencia de muchos componentes pequeños en el PCA teórico de $x$ . Estos pequeños componentes se utilizan implícitamente para ajustar el ruido mediante coeficientes pequeños. Cuando se utiliza OLS de norma mínima, se ajusta el ruido con muchos componentes independientes pequeños y esto tiene un efecto regularizador equivalente a la regularización de Ridge. Esta regularización es a menudo demasiado fuerte, y es posible compensarla utilizando una "antirregularización" conocida como Ridge negativo . En ese caso, verá que el mínimo de la curva MSE aparece para valores negativos de $\lambda$ .

Por PCA teórico, quiero decir:

Dejemos que $x\sim N(0,\Sigma)$ una distribución normal multivariante. Existe una isometría lineal $f$ como $u=f(x)\sim N(0,D)$ donde $D$ es diagonal: los componentes de $u$ son independientes. $D$ se obtiene simplemente diagonalizando $\Sigma$ .

Ahora el modelo $y=\beta.x+\epsilon$ se puede escribir $y=f(\beta).f(x)+\epsilon$ (una isometría lineal preserva el producto punto). Si se escribe $\gamma=f(\beta)$ el modelo se puede escribir $y=\gamma.u+\epsilon$ . Además $\|\beta\|=\|\gamma\|$ por lo que métodos de ajuste como el Ridge o el OLS de norma mínima son perfectamente isomórficos: el estimador de $y=\gamma.u+\epsilon$ es la imagen por $f$ del estimador de $y=\beta.x+\epsilon$ .

El ACP teórico transforma los predictores no independientes en predictores independientes. Sólo está vagamente relacionado con el PCA empírico, en el que se utiliza la matriz de covarianza empírica (que difiere mucho de la teórica con un tamaño de muestra pequeño). El PCA teórico no es computable en la práctica, pero sólo se utiliza aquí para interpretar el modelo en un espacio de predictores ortogonal.

Veamos qué ocurre cuando añadimos muchos predictores independientes de pequeña varianza a un modelo:

Teorema

Regularización de cresta con coeficiente $\lambda$ es equivalente (cuando $p\rightarrow\infty$ ) a:

• añadiendo $p$ falsos predictores independientes (centrados e idénticamente distribuidos) cada uno con varianza $\frac{\lambda}{p}$
• ajuste del modelo enriquecido con el estimador OLS de norma mínima
• manteniendo sólo los parámetros de los verdaderos predictores

(esbozo de) Prueba

Vamos a demostrar que las funciones de coste son asintóticamente iguales. Vamos a dividir el modelo en predictores reales y falsos: $y=\beta x+\beta'x'+\epsilon$ . La función de coste de Ridge (para el verdadero predictores) puede escribirse:

$$\mathrm{cost}_\lambda=\|\beta\|^2+\frac{1}{\lambda}\|y-X\beta\|^2$$

Cuando se utiliza el método OLS de norma mínima, la respuesta se ajusta perfectamente: el término de error es 0. término de error es 0. La función de coste sólo se refiere a la norma de los parámetros. Se puede dividir en los parámetros verdaderos y los falsos falsos:

$$\mathrm{cost}_{\lambda,p}=\|\beta\|^2+\inf\{\|\beta'\|^2 \mid X'\beta'=y-X\beta\}$$

En la expresión correcta, la solución de norma mínima viene dada por:

$$\beta'=X'^+(y-X\beta )$$

Ahora usando SVD para $X'$ :

$$X'=U\Sigma V$$

$$X'^{+}=V^\top\Sigma^{+} U^\top$$

Vemos que la norma de $\beta'$ depende esencialmente del singular valores singulares de $X'^+$ que son los recíprocos de los valores singulares de $X'$ . La versión normalizada de $X'$ es $\sqrt{p/\lambda} X'$ . He mirado la literatura y los valores singulares de las grandes matrices aleatorias son bien conocidos. Para $p$ y $n$ lo suficientemente grande, mínimo $s_\min$ y el máximo $s_\max$ los valores singulares se aproximan por (ver teorema 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)$$

Ya que, para los grandes $p$ , $\sqrt{n/p}$ tiende a 0, podemos decir simplemente que todos los valores singulares se aproximan por $\sqrt p$ . Así:

$$\|\beta'\|\approx\frac{1}{\sqrt\lambda}\|y-X\beta\|$$

Finalmente:

$$\mathrm{cost}_{\lambda,p}\approx\|\beta\|^2+\frac{1}{\lambda}\|y-X\beta\|^2=\mathrm{cost}_\lambda$$

Nota no importa si mantiene los coeficientes de los predictores falsos en su modelo. predictores en su modelo. La varianza introducida por $\beta'x'$ es $\frac{\lambda}{p}\|\beta'\|^2\approx\frac{1}{p}\|y-X\beta\|^2\approx\frac{n}{p}MSE(\beta)$ . De este modo, se incrementa el MSE en un factor $1+n/p$ sólo que tiende hacia 1 de todos modos. De alguna manera no es necesario tratar el predictores falsos de forma diferente a los reales.

Ahora, volviendo a los datos de @amoeba. Después de aplicar el PCA teórico a $x$ (se supone que es normal), $x$ se transforma mediante una isometría lineal en una variable $u$ cuyos componentes son independientes y ordenados por orden de varianza decreciente. El problema $y=\beta x+\epsilon$ es equivalente el problema transformado $y=\gamma u+\epsilon$ .

Ahora imagina que la varianza de los componentes es así:

Considere muchos $p$ de los últimos componentes, llamamos a la suma de su varianza $\lambda$ . Cada uno de ellos tiene una varianza aproximadamente igual a $\lambda/p$ y son independientes. Desempeñan el papel de falsos predictores en el teorema.

Este hecho es más claro en el modelo de @jonny: sólo el primer componente del PCA teórico está correlacionado con $y$ (es proporcional $\overline{x}$ ) y tiene una enorme variabilidad. Todos los demás componentes (proporcionales a $x_i-\overline{x}$ ) tienen una varianza comparativamente muy pequeña (escriba la matriz de covarianza y diagonalícela para ver esto) y juegan el papel de falsos predictores. He calculado que la regularización en este caso corresponde (aproximadamente) a la previa $N(0,\frac{1}{p^2})$ en $\gamma_1$ mientras que el verdadero $\gamma_1^2=\frac{1}{p}$ . Esto definitivamente se sobre-encoge. Esto es visible por el hecho de que el MSE final es mucho mayor que el MSE ideal. El efecto de regularización es demasiado fuerte.

A veces es posible mejorar esta regularización natural mediante Ridge. En primer lugar, a veces es necesario $p$ en el teorema realmente grande (1000, 10000...) para rivalizar seriamente con Ridge y la finitud de $p$ es como una imprecisión. Pero también muestra que Ridge es una regularización adicional sobre una regularización implícita naturalmente existente y que, por tanto, sólo puede tener un efecto muy pequeño. A veces esta regularización natural ya es demasiado fuerte y Ridge puede no ser ni siquiera una mejora. Más que esto, es mejor utilizar una antirregularización: Ridge con coeficiente negativo. Esto muestra el MSE del modelo de @jonny ( $p=1000$ ), utilizando $\lambda\in\mathbb{R}$ :

Answer 2

Gracias a todos por el gran debate en curso. El quid de la cuestión parece ser que El método OLS de norma mínima realiza efectivamente una contracción similar a la regresión de cresta. Esto parece ocurrir siempre que $p\gg n$ . Irónicamente, añadir predictores de ruido puro puede incluso utilizarse como una forma muy extraña de regularización.

---

Parte I. Demostración con datos artificiales y CV analítico

@Jonny (+1) propuso un ejemplo artificial muy sencillo que adaptaré ligeramente aquí. $X$ de $n\times p$ tamaño y $y$ se generan de forma que todas las variables son gaussianas con varianza unitaria, y la correlación entre cada predictor y la respuesta es $\rho$ . Voy a arreglar $\rho=.2$ .

Utilizaré el CV de izquierda a derecha porque existe una expresión analítica para el error al cuadrado: se conoce como PRESS "suma de cuadrados prevista". $$\text{PRESS} = \sum_i \left( \frac{e_i}{1-H_{ii}}\right)^2,$$ donde $e_i$ son residuos $$e = y - \hat y = y - Hy,$$ y $H$ es la matriz del sombrero $$H = X (X^\top X + \lambda I)^{-1} X^\top=U\frac{S^2}{S^2+\lambda} U^\top$$ en términos de SVD $X=USV^\top$ . Esto permite replicar los resultados de @Jonny sin utilizar glmnet y sin realizar realmente la validación cruzada (estoy trazando la relación entre PRESS y la suma de cuadrados de $y$ ):

Este enfoque analítico permite calcular el límite en $\lambda\to 0$ . Simplemente conectando $\lambda=0$ en la fórmula PRESS no funciona: cuando $n<p$ y $\lambda=0$ Los residuos son todos cero y la matriz del sombrero es la matriz de identidad con unos en la diagonal, lo que significa que las fracciones en la ecuación de PRESS son indefinidas. Pero si calculamos el límite en $\lambda \to 0$ entonces corresponderá a la solución OLS de norma mínima con $\lambda=0$ .

El truco es hacer la expansión de Taylor de la matriz del sombrero cuando $\lambda\to 0$ : $$H=U\frac{1}{1+\lambda/S^2} U^\top\approx U(1-\lambda/S^2) U^\top = I - \lambda US^{-2}U^\top = I-\lambda G^{-1}.$$ Aquí introduje la matriz Gram $G=XX^\top = US^2U^\top$ .

Ya casi hemos terminado: $$\text{PRESS} = \sum_i\Big( \frac{\lambda [G^{-1}y]_i}{\lambda G^{-1}_{ii}}\Big)^2 = \sum_i\Big( \frac{ [G^{-1}y]_i}{G^{-1}_{ii}}\Big)^2.$$ Lambda se anuló, así que aquí tenemos el valor límite. Lo he representado con un gran punto negro en la figura anterior (en los paneles donde $p>n$ ), y se combina perfectamente.

Actualización del 21 de febrero. La fórmula anterior es exacta, pero podemos obtener alguna información haciendo más aproximaciones. Parece que $G^{-1}$ tiene valores aproximadamente iguales en la diagonal aunque $S$ tiene valores muy desiguales (probablemente porque $U$ mezcla bastante bien todos los valores propios). Así que para cada $i$ tenemos que $G^{-1}_{ii}\approx \langle S^{-2} \rangle$ donde los paréntesis angulares denotan el promedio. Usando esta aproximación, podemos reescribir: $$\text{PRESS}\approx \Big\lVert \frac{S^{-2}}{\langle S^{-2} \rangle}U^\top y\Big\rVert^2.$$ Esta aproximación se muestra en la figura anterior con círculos rojos abiertos.

Si esto será más grande o más pequeño que $\lVert y \rVert^2 = \lVert U^\top y \rVert^2$ depende de los valores singulares $S$ . En esta simulación $y$ está correlacionada con el primer PC de $X$ así que $U_1^\top y$ es grande y todos los demás términos son pequeños. (En mis datos reales, $y$ también está bien predicho por los PC principales). Ahora, en el $p\gg n$ caso, si las columnas de $X$ son suficientemente aleatorios, entonces todos los valores singulares estarán bastante cerca unos de otros (filas aproximadamente ortogonales). El término "principal $U_1^\top y$ se multiplicará por un factor menor que 1. Los términos hacia el final se multiplicarán por factores mayores que 1 pero no mucho mayores. En general, la norma disminuye. En cambio, en el $p\gtrsim n$ caso, habrá algunos valores singulares muy pequeños. Después de la inversión se convertirán en factores grandes que aumentarán la norma general.

[Este argumento es muy manoseado; espero que se pueda precisar más].

Como comprobación de cordura, si intercambio el orden de los valores singulares por S = diag(flipud(diag(S))); entonces el MSE previsto es superior a $1$ por todas partes en los paneles 2 y 3.

figure('Position', [100 100 1000 300])
ps = [10, 100, 1000];

for pnum = 1:length(ps)
    rng(42)
    n = 80;
    p = ps(pnum);
    rho = .2;
    y = randn(n,1);
    X = repmat(y, [1 p])*rho + randn(n,p)*sqrt(1-rho^2);

    lambdas = exp(-10:.1:20);
    press = zeros(size(lambdas));
    [U,S,V] = svd(X, 'econ');
    % S = diag(flipud(diag(S)));   % sanity check

    for i = 1:length(lambdas)
        H = U * diag(diag(S).^2./(diag(S).^2 + lambdas(i))) * U';
        e = y - H*y;
        press(i) = sum((e ./ (1-diag(H))).^2);
    end

    subplot(1, length(ps), pnum)
    plot(log(lambdas), press/sum(y.^2))
    hold on
    title(['p = ' num2str(p)])
    plot(xlim, [1 1], 'k--')

    if p > n
        Ginv = U * diag(diag(S).^-2) * U';
        press0 = sum((Ginv*y ./ diag(Ginv)).^2);
        plot(log(lambdas(1)), press0/sum(y.^2), 'ko', 'MarkerFaceColor', [0,0,0]);

        press0approx = sum((diag(diag(S).^-2/mean(diag(S).^-2)) * U' * y).^2);
        plot(log(lambdas(1)), press0approx/sum(y.^2), 'ro');
    end
end

---

Parte II. Añadir predictores de ruido puro como forma de regularización

@Jonny, @Benoit, @Paul, @Dikran y otros han presentado buenos argumentos de que el aumento del número de predictores reducirá la solución OLS de norma mínima. De hecho, una vez que $p>n$ cualquier nuevo predictor sólo puede disminuir la norma de la solución de norma mínima. Por lo tanto, añadir predictores empujará la norma hacia abajo, de forma similar a cómo la regresión de cresta está penalizando la norma.

Entonces, ¿se puede utilizar como estrategia de regularización? Comenzamos con $n=80$ y $p=40$ y luego seguir añadiendo $q$ predictores de ruido puro como intento de regularización. Haré el LOOCV y lo compararé con el LOOCV de la cresta (calculado como arriba). Tenga en cuenta que después de obtener $\hat\beta$ en el $p+q$ predictores, lo estoy "truncando" en $p$ porque sólo me interesan los predictores originales.

¡¡¡FUNCIONA!!!

De hecho, no es necesario "truncar" la beta; incluso si utilizo la beta completa y la $p+q$ predictores, puedo obtener un buen rendimiento (línea discontinua en el subgrupo de la derecha). Creo que esto imita los datos reales de la pregunta: sólo unos pocos predictores predicen realmente $y$ La mayoría de ellos son puro ruido, y sirven de regularización. En este régimen, la regularización adicional de las crestas no ayuda en absoluto.

rng(42)
n = 80;
p = 40;
rho = .2;
y = randn(n,1);
X = repmat(y, [1 p])*rho + randn(n,p)*sqrt(1-rho^2);

lambdas = exp(-10:.1:20);
press = zeros(size(lambdas));
[U,S,V] = svd(X, 'econ');

for i = 1:length(lambdas)
    H = U * diag(diag(S).^2./(diag(S).^2 + lambdas(i))) * U';
    e = y - H*y;
    press(i) = sum((e ./ (1-diag(H))).^2);
end

figure('Position', [100 100 1000 300])
subplot(121)
plot(log(lambdas), press/sum(y.^2))
hold on
xlabel('Ridge penalty (log)')
plot(xlim, [1 1], 'k--')
title('Ridge regression (n=80, p=40)')
ylim([0 2])

ps = [0 20 40 60 80 100 200 300 400 500 1000];
error = zeros(n, length(ps));
error_trunc = zeros(n, length(ps));
for fold = 1:n
    indtrain = setdiff(1:n, fold);
    for pi = 1:length(ps)
        XX = [X randn(n,ps(pi))];
        if size(XX,2) < size(XX,1)
            beta = XX(indtrain,:) \ y(indtrain,:);
        else
            beta = pinv(XX(indtrain,:)) * y(indtrain,:);
        end
        error(fold, pi) = y(fold) - XX(fold,:) * beta;
        error_trunc(fold, pi) = y(fold) - XX(fold,1:size(X,2)) * beta(1:size(X,2));
    end
end

subplot(122)
hold on
plot(ps, sum(error.^2)/sum(y.^2), 'k.--')
plot(ps, sum(error_trunc.^2)/sum(y.^2), '.-')
legend({'Entire beta', 'Truncated beta'}, 'AutoUpdate','off')
legend boxoff
xlabel('Number of extra predictors')
title('Extra pure noise predictors')
plot(xlim, [1 1], 'k--')
ylim([0 2])

Answer 3

He aquí una situación artificial en la que esto ocurre. Supongamos que cada variable de predicción es una copia de la variable objetivo con una gran cantidad de ruido gaussiano aplicado. El mejor modelo posible es una media de todas las variables predictoras.

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 se comportan de forma "normal": Algún valor positivo de lambda minimiza el error fuera de la muestra.

Pero aumenta num.vars en el código anterior a 1000, y aquí está la nueva ruta de MSE. (He ampliado a log(Lambda) = -100 para convencerme.

Lo que creo que está sucediendo

Cuando se ajustan muchos parámetros con una regularización baja, los coeficientes se distribuyen aleatoriamente alrededor de su valor real con una varianza alta.

A medida que el número de predictores se hace muy grande, el "error medio" tiende a cero, y es mejor dejar que los coeficientes caigan donde puedan y sumar todo que sesgarlos hacia 0.

Estoy seguro de que esta situación de que la predicción verdadera sea una media de todos los predictores no es la única vez que ocurre, pero no sé cómo empezar a señalar la mayor condición necesaria aquí.

EDITAR:

El comportamiento "plano" para lambda muy baja siempre ocurrirá, ya que la solución está convergiendo a la solución OLS de norma mínima. Del mismo modo, la curva será plana para lambda muy alta, ya que la solución converge a 0. No habrá mínimo si una de esas dos soluciones es óptima.

¿Por qué la solución OLS de norma mínima es tan (comparativamente) buena en este caso? Creo que está relacionado con el siguiente comportamiento que me pareció muy contraintuitivo, pero que al reflexionar tiene mucho sentido.

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)

Con predictores generados aleatoriamente y no relacionados con la respuesta, a medida que p aumenta los coeficientes se hacen más grandes, pero una vez que p es mucho mayor que N se reducen hacia cero. Esto también ocurre en mi ejemplo. Por lo tanto, las soluciones no regularizadas para estos problemas no necesitan reducirse porque ya son muy pequeñas.

Esto sucede por una razón trivial. $y$ puede expresarse exactamente como una combinación lineal de columnas de $X$ . $\hat{\beta}$ es el vector de coeficientes de norma mínima. A medida que se añaden más columnas, la norma de $\hat{\beta}$ debe disminuir o permanecer constante, ya que una posible combinación lineal es mantener los coeficientes anteriores iguales y fijar los nuevos coeficientes en $0$ .

Answer 4

Así que decidí realizar una validación cruzada anidada utilizando el sistema especializado mlr en R para ver lo que realmente viene del enfoque de modelado.

Código (tarda unos minutos en ejecutarse en un portátil normal)

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])
)

rdesc = makeResampleDesc("CV", iters = 10)
msrs = list(mse, rsq)
configureMlr(on.par.without.desc = "quiet")
bm3 = benchmark(learners = list(
    makeLearner("regr.cvglmnet", alpha = 0, lambda = c(0, exp(seq(-10, 10, length.out = 150))),
    makeLearner("regr.glmnet", alpha = 0, lambda = c(0, exp(seq(-10, 10, length.out = 150))), s = 151)
    ), tasks = tsk, resamplings = rdesc, measures = msrs)

Resultados

getBMRAggrPerformances(bm3, as.df = TRUE)
#   task.id    learner.id mse.test.mean rsq.test.mean
#1    tsk10 regr.cvglmnet     1.0308055  -0.224534550
#2    tsk10   regr.glmnet     1.3685799  -0.669473387
#3   tsk100 regr.cvglmnet     0.7996823   0.031731316
#4   tsk100   regr.glmnet     1.3092522  -0.656879104
#5  tsk1110 regr.cvglmnet     0.8236786   0.009315037
#6  tsk1110   regr.glmnet     0.6866745   0.117540454
#7    tsk50 regr.cvglmnet     1.0348319  -0.188568886
#8    tsk50   regr.glmnet     2.5468091  -2.423461744
#9   tsk500 regr.cvglmnet     0.7210185   0.173851634
#10  tsk500   regr.glmnet     0.6171841   0.296530437

Hacen básicamente lo mismo en todas las tareas.

¿Y qué pasa con las lambdas óptimas?

sapply(lapply(getBMRModels(bm3, task.ids = "tsk1110")[[1]][[1]], "[[", 2), "[[", "lambda.min")
# [1] 4.539993e-05 4.539993e-05 2.442908e-01 1.398738e+00 4.539993e-05
# [6] 0.000000e+00 4.539993e-05 3.195187e-01 2.793841e-01 4.539993e-05

Observa que las lambdas ya están transformadas. Algunos pliegues incluso eligieron la lambda mínima $\lambda = 0$ .

He jugado un poco más con glmnet y descubrió que tampoco allí se escoge el lambda mínimo. Compruébalo:

EDITAR:

Tras los comentarios de ameba, quedó claro que la vía de regularización es un paso importante en la glmnet estimación, por lo que el código ahora lo refleja. De este modo, la mayoría de las discrepancias desaparecieron.

cvfit = cv.glmnet(x = x, y = y, alpha = 0, lambda = exp(seq(-10, 10, length.out = 150)))
plot(cvfit)

Conclusión

Así que, básicamente, $\lambda>0$ realmente mejora el ajuste ( edit: ¡pero no por mucho! ).

¿Cómo es posible y qué dice sobre mi conjunto de datos? ¿Me estoy perdiendo algo obvio o es realmente contraintuitivo?

Es probable que estemos más cerca de la verdadera distribución de la configuración de los datos $\lambda$ a un valor pequeño mayor que cero. Sin embargo, no hay nada contraintuitivo en ello.

Editar: Hay que tener en cuenta, sin embargo, que la ruta de regularización de crestas hace uso de las estimaciones previas de los parámetros cuando llamamos a glmnet pero esto va más allá de mis conocimientos. Si fijamos un valor realmente bajo lambda de forma aislada, es probable que se reduzca el rendimiento.

EDIT: La selección de lambdas dice algo más sobre sus datos. Como las lambdas más grandes disminuyen el rendimiento, significa que hay preferentes, es decir mayores, los coeficientes en su modelo, ya que las lambdas grandes reducen todos los coeficientes hacia cero. Aunque $\lambda\neq0$ significa que los grados de libertad efectivos de su modelo son menores que los grados de libertad aparentes, $p$ .

¿Cómo puede haber alguna diferencia cualitativa entre p=100 y p=1000 dado que ambos son mayores que n?

$p=1000$ contiene invariablemente al menos la misma información o incluso más que $p=100$ .

---

Comentarios

Parece que estás obteniendo un mínimo para alguna lambda no nula (estoy mirando tu figura), pero la curva sigue siendo realmente muy plana a la izquierda de la misma. Así que mi pregunta principal sigue siendo por qué 0 no sobreajusta notablemente. No veo una respuesta aquí todavía. ¿Espera que esto sea un fenómeno general? Es decir, ¿para cualquier dato con np, lambda=0 funcionará [casi] tan bien como lambda óptima? ¿O es algo especial de estos datos? Si miras arriba en los comentarios, verás que mucha gente ni siquiera me creía que fuera posible.

Creo que estás confundiendo el rendimiento de la validación con el rendimiento de las pruebas, y esa comparación no está justificada.

Editar: notar, sin embargo, que cuando establecemos lambda ¡a 0 después de ejecutar toda la ruta de regularización el rendimiento no se degrada como tal, por lo tanto la ruta de regularización es clave para entender lo que está pasando!

Además, no entiendo bien tu última frase. Mira la salida de cv.glmnet para p=100. Tendrá una forma muy diferente. Entonces, ¿qué afecta a esta forma (asíntota a la izquierda vs. sin asíntota) cuando p=100 o p=1000?

Comparemos las vías de regularización de ambos:

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")

Queda claro $p=1000$ permite obtener mayores coeficientes al aumentar $\lambda$ aunque tiene coeficientes más pequeños para la cresta asintótica-OLS, a la izquierda de ambos gráficos. Así que, básicamente, $p=100$ sobreajustes a la izquierda del gráfico, y eso probablemente explica la diferencia de comportamiento entre ellos.

Es más difícil para $p=1000$ para sobreajustar porque, aunque Ridge reduce los coeficientes a cero, nunca llegan a cero. Esto significa que el poder de predicción del modelo se reparte entre muchos más componentes, lo que facilita la predicción en torno a la media en lugar de dejarse llevar por el ruido.

Answer 5

¿Cómo es posible que el método OLS (norma mínima) no se ajuste en exceso?

En resumen:

Los parámetros experimentales que se correlacionan con los parámetros (desconocidos) del modelo verdadero tendrán más probabilidades de ser estimados con valores altos en un procedimiento de ajuste OLS de norma mínima. Esto se debe a que se ajustarán al "modelo+ruido", mientras que los demás parámetros sólo se ajustarán al "ruido" (por lo que se ajustarán a una parte mayor del modelo con un valor más bajo del coeficiente y será más probable que tengan un valor alto en el MCO de norma mínima).

Este efecto reducirá la cantidad de sobreajuste en un procedimiento de ajuste OLS de norma mínima. El efecto es más pronunciado si se dispone de más parámetros, ya que entonces es más probable que se incorpore a la estimación una parte mayor del "modelo verdadero".

La parte más larga:
(No estoy seguro de qué poner aquí, ya que la cuestión no me queda del todo clara, o no sé con qué precisión hay que responder a la pregunta)

A continuación se muestra un ejemplo que se puede construir fácilmente y que demuestra el problema. El efecto no es tan extraño y los ejemplos son fáciles de hacer.

• Tomé $p=200$ las funciones sin (porque son perpendiculares) como variables
• creó un modelo aleatorio con $n=50$ medidas.
  • El modelo es construido con sólo $tm=10$ de las variables por lo que 190 de las 200 variables están creando la posibilidad de generar un sobreajuste.
  • los coeficientes del modelo se determinan al azar

En este caso de ejemplo observamos que hay un cierto sobreajuste pero los coeficientes de los parámetros que pertenecen al modelo verdadero tienen un valor más alto. Por lo tanto, el R^2 puede tener algún valor positivo.

La imagen siguiente (y el código para generarla) demuestran que el sobreajuste es limitado. Los puntos que se refieren al modelo de estimación de 200 parámetros. Los puntos rojos se refieren a los parámetros que también están presentes en el "modelo real" y vemos que tienen un valor más alto. Por lo tanto, hay un cierto grado de aproximación al modelo real y de obtener el R^2 por encima de 0.

• Obsérvese que he utilizado un modelo con variables ortogonales (las funciones senoidales). Si los parámetros están correlacionados, entonces pueden aparecer en el modelo con un coeficiente relativamente muy alto y ser más penalizados en la norma mínima OLS.
• Nótese que las "variables ortogonales" no son ortogonales cuando consideramos los datos. El producto interno de $sin(ax) \cdot sin(bx)$ sólo es cero cuando integramos todo el espacio de $x$ y no cuando sólo tenemos unas pocas muestras $x$ . La consecuencia es que incluso con ruido cero se producirá un sobreajuste (y el valor de R^2 parece depender de muchos factores, aparte del ruido. Por supuesto, existe la relación $n$ y $p$ Pero también es importante saber cuántas variables hay en el modelo verdadero y cuántas hay en el modelo de ajuste).

library(MASS)

par(mar=c(5.1, 4.1, 9.1, 4.1), xpd=TRUE)

p <- 200       
l <- 24000
n <- 50
tm <- 10

# generate i sinus vectors as possible parameters
t <- c(1:l)
xm <- sapply(c(0:(p-1)), FUN = function(x) sin(x*t/l*2*pi))

# generate random model by selecting only tm parameters
sel <- sample(1:p, tm)
coef <- rnorm(tm, 2, 0.5)

# generate random data xv and yv with n samples
xv <- sample(t, n)
yv <- xm[xv, sel] %*% coef + rnorm(n, 0, 0.1)

# generate model
M <- ginv(t(xm[xv,]) %*% xm[xv,])

Bsol <- M %*% t(xm[xv,]) %*% yv
ysol <- xm[xv,] %*% Bsol

# plotting comparision of model with true model
plot(1:p, Bsol, ylim=c(min(Bsol,coef),max(Bsol,coef)))
points(sel, Bsol[sel], col=1, bg=2, pch=21)
points(sel,coef,pch=3,col=2)

title("comparing overfitted model (circles) with true model (crosses)",line=5)
legend(0,max(coef,Bsol)+0.55,c("all 100 estimated coefficients","the 10 estimated coefficients corresponding to true model","true coefficient values"),pch=c(21,21,3),pt.bg=c(0,2,0),col=c(1,1,2))

Técnica de beta truncada en relación con la regresión de cresta

He transformado el código python de Amoeba en R y he combinado los dos gráficos juntos. Para cada estimación OLS de norma mínima con variables de ruido añadidas hago coincidir una estimación de regresión ridge con la misma (aproximadamente) $l_2$ -norma para el $\beta$ vectorial.

• Parece que el modelo de ruido truncado hace más o menos lo mismo (sólo computa un poco más lento, y tal vez un poco más a menudo menos bueno).

• Sin embargo, sin el truncamiento el efecto es mucho menos fuerte.

• Esta correspondencia entre la adición de parámetros y la penalización de la cresta no es necesariamente el mecanismo más fuerte detrás de la ausencia de sobreajuste. Esto se puede ver especialmente en la curva de 1000p (en el imagen de la pregunta) llegando casi a 0,3 mientras que las otras curvas con diferentes p, no alcanzan este nivel, sea cual sea el parámetro de de regresión. Los parámetros adicionales, en ese caso práctico, no son lo mismo que un desplazamiento del parámetro de la cresta (y supongo que esto se debe a que los parámetros adicionales crearán un modelo mejor, más completo).

• Los parámetros de ruido reducen la norma por un lado (al igual que la regresión de cresta) pero también introducen ruido adicional. Benoit Sánchez muestra que en el límite, añadiendo muchos parámetros de ruido con una desviación menor, se convertirá finalmente en lo mismo que la regresión de cresta (el número creciente de parámetros de ruido se anulan entre sí). Pero al mismo tiempo, requiere muchos más cálculos (si aumentamos la desviación del ruido, para permitir usar menos parámetros y acelerar el cálculo, la diferencia se hace mayor).

Rho = 0,2

Rho = 0,4

Rho = 0,2 aumentando la varianza de los parámetros de ruido a 2

ejemplo de código

# prepare the data
set.seed(42)
n = 80
p = 40
rho = .2
y = rnorm(n,0,1)
X = matrix(rep(y,p), ncol = p)*rho + rnorm(n*p,0,1)*(1-rho^2)

# range of variables to add
ps = c(0, 5, 10, 15, 20, 40, 45, 50, 55, 60, 70, 80, 100, 125, 150, 175, 200, 300, 400, 500, 1000)
#ps = c(0, 5, 10, 15, 20, 40, 60, 80, 100, 150, 200, 300) #,500,1000)

# variables to store output (the sse)
error   = matrix(0,nrow=n, ncol=length(ps))
error_t = matrix(0,nrow=n, ncol=length(ps))
error_s = matrix(0,nrow=n, ncol=length(ps))

# adding a progression bar
pb <- txtProgressBar(min = 0, max = n, style = 3)

# training set by leaving out measurement 1, repeat n times 
for (fold in 1:n) {
    indtrain = c(1:n)[-fold]

    # ridge regression
    beta_s <- glmnet(X[indtrain,],y[indtrain],alpha=0,lambda = 10^c(seq(-4,2,by=0.01)))$beta
    # calculate l2-norm to compare with adding variables
    l2_bs <- colSums(beta_s^2)

    for (pi in 1:length(ps)) {
        XX = cbind(X, matrix(rnorm(n*ps[pi],0,1), nrow=80))
        XXt = XX[indtrain,]

        if (p+ps[pi] < n) {
            beta = solve(t(XXt) %*% (XXt)) %*% t(XXt) %*% y[indtrain]
        }
        else {
            beta = ginv(t(XXt) %*% (XXt)) %*% t(XXt) %*% y[indtrain]
        }

        # pickout comparable ridge regression with the same l2 norm      
        l2_b <- sum(beta[1:p]^2)
        beta_shrink <- beta_s[,which.min((l2_b-l2_bs)^2)] 

        # compute errors
        error[fold, pi] = y[fold] - XX[fold,1:p] %*% beta[1:p]
        error_t[fold, pi] = y[fold] - XX[fold,] %*% beta[]
        error_s[fold, pi] = y[fold] - XX[fold,1:p] %*% beta_shrink[]
    }
    setTxtProgressBar(pb, fold) # update progression bar
}

# plotting
plot(ps,colSums(error^2)/sum(y^2) , 
     ylim = c(0,2),
     xlab ="Number of extra predictors",
     ylab ="relative sum of squared error")
lines(ps,colSums(error^2)/sum(y^2))
points(ps,colSums(error_t^2)/sum(y^2),col=2)
lines(ps,colSums(error_t^2)/sum(y^2),col=2)
points(ps,colSums(error_s^2)/sum(y^2),col=4)
lines(ps,colSums(error_s^2)/sum(y^2),col=4)

title('Extra pure noise predictors')

legend(200,2,c("complete model with p + extra predictors",
               "truncated model with p + extra predictors",
               "ridge regression with similar l2-norm",
               "idealized model uniform beta with 1/p/rho"),
       pch=c(1,1,1,NA), col=c(2,1,4,1),lt=c(1,1,1,2))

# idealized model (if we put all beta to 1/rho/p we should theoretically have a reasonable good model)
error_op <- rep(0,n)
for (fold in 1:n) {
  beta = rep(1/rho/p,p)
    error_op[fold] = y[fold] - X[fold,] %*% beta
}
id <- sum(error_op^2)/sum(y^2)
lines(range(ps),rep(id,2),lty=2)