使用“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")

请注意，当我调用cv.glmnet（或glmnet）时，我已经设置了intercept = FALSE。这比它在实践中的影响具有更多的概念意义。从概念上讲，我们的教科书计算没有截距，所以我们想在使用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；类似地，X %*% t(X) 的 tcrossprod(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?的第一部分

Answer 2

在哲元那篇有趣的帖子之上，又做了一些实验，发现我们也可以用intercept得到同样的结果，如下：

# ridge with intercept glmnet
ridge.fit.cv.int <- cv.glmnet(X, Y, alpha = 0, intercept = TRUE, family="gaussian")
ridge.fit.lambda.int <- ridge.fit.cv.int$lambda.1se
ridge.coef.with.int <- as.vector(as.matrix(coef(ridge.fit.cv.int, s = ridge.fit.lambda.int)))

# ridge with intercept definition, use the same lambda obtained with cv from glmnet
X.with.int <- cbind(1, X)
ridge.coef.DEF.with.int <- drop(solve(crossprod(X.with.int) + ridge.fit.lambda.int * diag(n.tmp, p.tmp+1), crossprod(X.with.int, Y)))

ggplot() + geom_point(aes(ridge.coef.with.int, ridge.coef.DEF.with.int))

# comupte residuals
RSS.fit.cv.int <- sum((Y.true - predict(ridge.fit.cv.int, newx=X))^2) # predict adds inter
RSS.DEF.int <- sum((Y.true - X.with.int %*% ridge.coef.DEF.with.int)^2)

RSS.fit.cv.int
[1] 110059.9
RSS.DEF.int
[1] 110063.4