¿Por qué la regresión de cresta glmnet me da una respuesta diferente al cálculo manual?

Estoy usando glmnet para calcular las estimaciones de regresión de crestas. Obtuve algunos resultados que me hicieron sospechar que glmnet realmente está haciendo lo que creo que hace. Para verificar esto, escribí un script R simple donde comparo el resultado de la regresión de cresta realizada por resolver y el de glmnet, la diferencia es significativa:

n    <- 1000
p.   <-  100
X.   <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y    <- X%*%beta+rnorm(n,0,0.5)

beta1 <- solve(t(X)%*%X+5*diag(p),t(X)%*%Y)
beta2 <- glmnet(X,Y, alpha=0, lambda=10, intercept=FALSE, standardize=FALSE, 
                family="gaussian")$beta@x
beta1-beta2

La norma de la diferencia es generalmente alrededor de 20, lo que no puede deberse a algoritmos numéricamente diferentes, debo estar haciendo algo mal. ¿Cuáles son las configuraciones que tengo que configurar glmnetpara obtener el mismo resultado que con la cresta?

La diferencia que está observando se debe a la división adicional por el número de observaciones, N, que GLMNET utiliza en su función objetivo y la estandarización implícita de Y por su desviación estándar de muestra como se muestra a continuación.

$$\frac{1}{2N}\left\|\frac{y}{s_y}-X\beta\right\|^2_{2}+\lambda\|\beta\|^2_{2}/2$$

donde usamos en lugar de para , 1 / ( n - 1 ) s y s y = ∑ i ( y i - ˉ y ) $1/n$$1/(n-1)$$s_y$

$$s_y=\frac{\sum_i(y_i-\bar{y})^2}{n}$$

Al diferenciar con respecto a beta, establecer la ecuación a cero,

$$X^TX\beta-\frac{X^Ty}{s_y}+N\lambda\beta =0$$

Y resolviendo para beta, obtenemos la estimación,

$$\tilde{\beta}_{GLMNET}= (X^TX+N\lambda I_p)^{-1}\frac{X^Ty}{s_y}$$

Para recuperar las estimaciones (y sus penalizaciones correspondientes) en la métrica original de Y, GLMNET multiplica tanto las estimaciones como las lambdas por y devuelve estos resultados al usuario,$s_y$

lambda u n s t d . = s y

$$\hat{\beta}_{GLMNET}=s_y\tilde{\beta}_{GLMNET}= (X^TX+N\lambda I_p)^{-1}X^Ty$$ $$\lambda_{unstd.}=s_y\lambda$$

Compare esta solución con la derivación estándar de regresión de cresta.

$$\hat{\beta}= (X^TX+\lambda I_p)^{-1}X^Ty$$

Tenga en cuenta que se escala por un factor adicional de N. Además, cuando usamos la función o , la penalización se escalará implícitamente por . Es decir, cuando usamos estas funciones para obtener las estimaciones de coeficientes para algunos , estamos efectivamente obteniendo estimaciones para .1 / s y λ ∗ λ = λ ∗ / s $\lambda$predict()coef()$1/s_y$$\lambda^*$$\lambda=\lambda^*/s_y$

Basándose en estas observaciones, la pena utilizado en GLMNET necesita ser reducido por un factor de .$s_y/N$

set.seed(123)

n    <- 1000
p   <-  100
X   <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y    <- X%*%beta+rnorm(n,0,0.5)

sd_y <- sqrt(var(Y)*(n-1)/n)[1,1]

beta1 <- solve(t(X)%*%X+10*diag(p),t(X)%*%(Y))[,1]

fit_glmnet <- glmnet(X,Y, alpha=0, standardize = F, intercept = FALSE, thresh = 1e-20)
beta2 <- as.vector(coef(fit_glmnet, s = sd_y*10/n, exact = TRUE))[-1]
cbind(beta1[1:10], beta2[1:10])

           [,1]        [,2]
[1,]  0.23793862  0.23793862
[2,]  1.81859695  1.81859695
[3,] -0.06000195 -0.06000195
[4,] -0.04958695 -0.04958695
[5,]  0.41870613  0.41870613
[6,]  1.30244151  1.30244151
[7,]  0.06566168  0.06566168
[8,]  0.44634038  0.44634038
[9,]  0.86477108  0.86477108
[10,] -2.47535340 -2.47535340

Los resultados generalizan a la inclusión de una intercepción y variables X estandarizadas. Modificamos una matriz X estandarizada para incluir una columna de unos y la matriz diagonal para tener una entrada cero adicional en la posición [1,1] (es decir, no penalizar la intersección). A continuación, puede desestandarizar las estimaciones por sus respectivas desviaciones estándar de muestra (nuevamente asegúrese de estar usando 1 / n al calcular la desviación estándar).

$$\hat\beta_{j}=\frac{\tilde{\beta_j}}{s_{x_j}}$$

$$\hat\beta_{0}=\tilde{\beta_0}-\bar{x}^T\hat{\beta}$$

mean_x <- colMeans(X)
sd_x <- sqrt(apply(X,2,var)*(n-1)/n)
X_scaled <- matrix(NA, nrow = n, ncol = p)
for(i in 1:p){
    X_scaled[,i] <- (X[,i] - mean_x[i])/sd_x[i] 
}
X_scaled_ones <- cbind(rep(1,n), X_scaled)

beta3 <- solve(t(X_scaled_ones)%*%X_scaled_ones+1000*diag(x = c(0, rep(1,p))),t(X_scaled_ones)%*%(Y))[,1]
beta3 <- c(beta3[1] - crossprod(mean_x,beta3[-1]/sd_x), beta3[-1]/sd_x)

fit_glmnet2 <- glmnet(X,Y, alpha=0, thresh = 1e-20)
beta4 <- as.vector(coef(fit_glmnet2, s = sd_y*1000/n, exact = TRUE))

cbind(beta3[1:10], beta4[1:10])
             [,1]        [,2]
 [1,]  0.24534485  0.24534485
 [2,]  0.17661130  0.17661130
 [3,]  0.86993230  0.86993230
 [4,] -0.12449217 -0.12449217
 [5,] -0.06410361 -0.06410361
 [6,]  0.17568987  0.17568987
 [7,]  0.59773230  0.59773230
 [8,]  0.06594704  0.06594704
 [9,]  0.22860655  0.22860655
[10,]  0.33254206  0.33254206

Código agregado para mostrar X estandarizado sin intercepción:

set.seed(123)

n <- 1000
p <-  100
X <- matrix(rnorm(n*p,0,1),n,p)
beta <- rnorm(p,0,1)
Y <- X%*%beta+rnorm(n,0,0.5)

sd_y <- sqrt(var(Y)*(n-1)/n)[1,1]

mean_x <- colMeans(X)
sd_x <- sqrt(apply(X,2,var)*(n-1)/n)

X_scaled <- matrix(NA, nrow = n, ncol = p)
for(i in 1:p){
    X_scaled[,i] <- (X[,i] - mean_x[i])/sd_x[i] 
}

beta1 <- solve(t(X_scaled)%*%X_scaled+10*diag(p),t(X_scaled)%*%(Y))[,1]

fit_glmnet <- glmnet(X_scaled,Y, alpha=0, standardize = F, intercept = 
FALSE, thresh = 1e-20)
beta2 <- as.vector(coef(fit_glmnet, s = sd_y*10/n, exact = TRUE))[-1]
cbind(beta1[1:10], beta2[1:10])

             [,1]        [,2]
 [1,]  0.23560948  0.23560948
 [2,]  1.83469846  1.83469846
 [3,] -0.05827086 -0.05827086
 [4,] -0.04927314 -0.04927314
 [5,]  0.41871870  0.41871870
 [6,]  1.28969361  1.28969361
 [7,]  0.06552927  0.06552927
 [8,]  0.44576008  0.44576008
 [9,]  0.90156795  0.90156795
[10,] -2.43163420 -2.43163420

De acuerdo con https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html , cuando la familia está gaussian, glmnet()debe minimizar

$$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda\sum_{j=1}^p(\alpha|\beta_j| +(1-\alpha)\beta_j^2/2). \tag{1}$$

Cuando se usa glmnet(x, y, alpha=1)para ajustar el lazo con las columnas en estandarizadas, la solución para la penalización reportada es la solución para minimizar Sin embargo, al menos en , cuando se usa para ajustar la regresión de cresta, la solución para una penalización reportada es la solución para minimizar donde es la desviación estándar de . Aquí, la penalización debería haberse informado como .$x$$\lambda$

$$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda \sum_{j=1}^p |\beta_j|.$$ glmnet_2.0-13glmnet(x, y, alpha=0)$\lambda$ $$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\lambda \frac{1}{2s_y} \sum_{j=1}^p \beta_j^2.$$ $s_y$$y$$\lambda/s_y$

Lo que podría suceder es que la función primero estandariza a y luego minimiza que efectivamente es minimizar o equivalente, para minimizar $y$$y_0$

$$\frac{1}{2n} \sum_{i=1}^n (y_{0i}-x_i^T\gamma)^2 +\eta \sum_{j=1}^p(\alpha|\gamma_j| +(1-\alpha)\gamma_j^2/2), \tag{2}$$ $$\frac{1}{2n s_y^2} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\eta \frac{\alpha}{s_y} \sum_{j=1}^p |\beta_j| +\eta \frac{1-\alpha}{2s_y^2} \sum_{j=1}^p \beta_j^2,$$ $$\frac{1}{2n} \sum_{i=1}^n (y_i-\beta_0-x_i^T\beta)^2 +\eta s_y \alpha \sum_{j=1}^p |\beta_j| +\eta (1-\alpha) \sum_{j=1}^p \beta_j^2/2.$$

Para el lazo ( ), volver a escalar para informar la penalización como tiene sentido. Entonces, para todos , debe ser reportado como la penalización para mantener la continuidad de los resultados en . Esta es probablemente la causa del problema anterior. Esto se debe en parte al uso de (2) para resolver (1). Solo cuando o existe cierta equivalencia entre los problemas (1) y (2) (es decir, una correspondencia entre en (1) y en (2)). Para cualquier otro$\alpha=1$$\eta$$\eta s_y$$\alpha$$\eta s_y$$\alpha$$\alpha=0$$\alpha=1$$\lambda$$\eta$$\alpha\in(0,1)$, los problemas (1) y (2) son dos problemas de optimización diferentes, y no hay correspondencia uno a uno entre en (1) y en (2).$\lambda$$\eta$