我正在使用 glmnet R 包运行岭回归。我注意到我从 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")

`glmnet` 的岭回归给出的系数与我通过  &#34;textbook definition&#34;计算的系数不同?-LMLPHP

最佳答案

如果你阅读 ?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 一起使用;使用 crossprodsolve(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")

`glmnet` 的岭回归给出的系数与我通过  &#34;textbook definition&#34;计算的系数不同?-LMLPHP

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

跟进


t(X) %*% Y 将首先进行转置 X1 <- t(X) ,然后进行 X1 %*% Y ,而 crossprod(X, Y) 不会进行转置。 "%*%" DGEMM 的包装器 op(A) = A, op(B) = B ,而 crossprodop(A) = A', op(B) = B 的包装器。 tcrossprodop(A) = A, op(B) = B' 类似。
crossprod(X) 的一个主要用途是用于 t(X) %*% X ;与 tcrossprod(X)X %*% t(X) 类似,在这种情况下,将调用 DSYRK 而不是 DGEMM。您可以阅读 Why the built-in lm function is so slow in R? 的第一部分 以了解原因和基准。

请注意,如果 X 不是方阵,则 crossprod(X)tcrossprod(X) 的速度并不相同,因为它们涉及不同数量的浮点运算,为此您可以阅读 Any faster R function than “tcrossprod” for symmetric dense matrix multiplication? 的旁注

关于 solvel(A, b)solve(A) %*% b ,请阅读 How to compute diag(X %% solve(A) %% t(X)) efficiently without taking matrix inverse? 的第一节

关于 `glmnet` 的岭回归给出的系数与我通过 "textbook definition"计算的系数不同?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/49991887/

10-12 19:57