使用`glmnet`的岭回归给出的系数与我通过"教科书定义"计算的系数不同？

Question

我使用glmnet R包运行Ridge回归.我注意到我从glmnet::glmnet函数中获得的系数与通过定义计算系数得到的系数不同(使用相同的lambda值).有人可以解释一下为什么吗？

数据(两者:响应Y和设计矩阵X)都是按比例缩放的.

library(MASS)
library(glmnet)

# Data dimensions
p.tmp <- 100
n.tmp <- 100

# Data objects
set.seed(1)
X <- scale(mvrnorm(n.tmp, mu = rep(0, p.tmp), Sigma = diag(p.tmp)))
beta <- rep(0, p.tmp)
beta[sample(1:p.tmp, 10, replace = FALSE)] <- 10
Y.true <- X %*% beta
Y <- scale(Y.true + matrix(rnorm(n.tmp))) # Y.true + Gaussian noise

# Run glmnet 
ridge.fit.cv <- cv.glmnet(X, Y, alpha = 0)
ridge.fit.lambda <- ridge.fit.cv$lambda.1se

# Extract coefficient values for lambda.1se (without intercept)
ridge.coef <- (coef(ridge.fit.cv, s = ridge.fit.lambda))[2:(p.tmp+1)]

# Get coefficients "by definition"
ridge.coef.DEF <- solve(t(X) %*% X + ridge.fit.lambda * diag(p.tmp)) %*% t(X) %*% Y

# Plot estimates
plot(ridge.coef, type = "l", ylim = range(c(ridge.coef, ridge.coef.DEF)),
     main = "black: Ridge `glmnet`\nred: Ridge by definition")
lines(ridge.coef.DEF, col = "red")

Answer 1

如果您阅读?glmnet,您将看到高斯响应的惩罚目标函数是:

1/2 * RSS / nobs + lambda * penalty

如果使用脊罚1/2 * ||beta_j||_2^2,我们有

1/2 * RSS / nobs + 1/2 * lambda * ||beta_j||_2^2

这与...成正比

RSS + lambda * nobs * ||beta_j||_2^2

这与我们通常在有关岭回归的教科书中看到的不同:

RSS + lambda * ||beta_j||_2^2

你写的公式:

##solve(t(X) %*% X + ridge.fit.lambda * diag(p.tmp)) %*% t(X) %*% Y
drop(solve(crossprod(X) + diag(ridge.fit.lambda, p.tmp), crossprod(X, Y)))

用于教科书的结果; 因为glmnet我们应该期待:

##solve(t(X) %*% X + n.tmp * ridge.fit.lambda * diag(p.tmp)) %*% t(X) %*% Y
drop(solve(crossprod(X) + diag(n.tmp * ridge.fit.lambda, p.tmp), crossprod(X, Y)))

因此,教科书使用惩罚的最小二乘法,但glmnet使用惩罚均方误差.

注意我没有使用你的原始代码t(),"%*%"和solve(A) %*% b; 使用crossprod并且solve(A, b)效率更高!请参阅最后的后续部分.

现在让我们进行一个新的比较:

library(MASS)
library(glmnet)

# Data dimensions
p.tmp <- 100
n.tmp <- 100

# Data objects
set.seed(1)
X <- scale(mvrnorm(n.tmp, mu = rep(0, p.tmp), Sigma = diag(p.tmp)))
beta <- rep(0, p.tmp)
beta[sample(1:p.tmp, 10, replace = FALSE)] <- 10
Y.true <- X %*% beta
Y <- scale(Y.true + matrix(rnorm(n.tmp)))

# Run glmnet 
ridge.fit.cv <- cv.glmnet(X, Y, alpha = 0, intercept = FALSE)
ridge.fit.lambda <- ridge.fit.cv$lambda.1se

# Extract coefficient values for lambda.1se (without intercept)
ridge.coef <- (coef(ridge.fit.cv, s = ridge.fit.lambda))[-1]

# Get coefficients "by definition"
ridge.coef.DEF <- drop(solve(crossprod(X) + diag(n.tmp * ridge.fit.lambda, p.tmp), crossprod(X, Y)))

# Plot estimates
plot(ridge.coef, type = "l", ylim = range(c(ridge.coef, ridge.coef.DEF)),
     main = "black: Ridge `glmnet`\nred: Ridge by definition")
lines(ridge.coef.DEF, col = "red")

请注意,我intercept = FALSE在调用cv.glmnet(或glmnet)时设置了.这比在实践中会产生更多的概念意义.从概念上讲,我们的教科书计算没有拦截,因此我们希望在使用时删除拦截glmnet.但实际上,由于你X和它Y是标准化的,拦截的理论估计是0.即使使用intercepte = TRUE(glment默认),你可以检查截距的估计是~e-17(数字0),因此其他系数的估计不会受到显着影响.另一个答案就是表明这一点.

跟进

至于使用crossprod和solve(A, b)- 有趣!您是否偶然参考了模拟比较？

t(X) %*% Y将首先采取转置X1 <- t(X),然后做X1 %*% Y,而crossprod(X, Y)不会进行转置."%*%"是一个DGEMM案例的包装op(A) = A, op(B) = B,crossprod而是一个包装op(A) = A', op(B) = B.同样tcrossprod的op(A) = A, op(B) = B'.

主要用途crossprod(X)是t(X) %*% X; 类似地,tcrossprod(X)for X %*% t(X),在这种情况下DSYRK代替DGEMM被调用.您可以阅读第一部分中为什么内置LM功能是如此缓慢R中？为理由和基准.

请注意,如果X不是一个方阵,crossprod(X)并tcrossprod(X)没有同样快,因为它们涉及不同数量的浮点运算,以便您可以读取侧通知的不是"tcrossprod"对称稠密矩阵乘法的任何快一个R函数？

关于solvel(A, b)和solve(A) %*% b,请阅读第一部分的如何计算诊断(X% %解决(A)%%T(X))有效地而不采取矩阵逆？