我正在使用 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
,你会看到高斯响应的惩罚目标函数是:
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),因此其他系数的估计不会受到显着影响。另一个答案只是展示了这一点。跟进
As for the using
crossprod
andsolve(A, b)
- interesting! Do you by chance have any reference to simulation comparison for that?
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)
为 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/