使用`glmnet`进行岭回归得到的系数与我通过“教科书定义”计算得到的系数不同？

Question

我正在使用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")

- Marta Karas

---

这实际上是一份关于统计辅导的请求，因此更适合CrossValidated.com。 （我认为答案是岭回归是一种受惩罚的方法，但您可能会从CV群体中得到更权威的答案。） - IRTFM

看起来这实际上是一个编程问题。如果我理解正确，楼主问的是为什么使用glmnet函数对于给定的lambda（惩罚项）值返回的系数与直接使用相同的lambda值解决回归系数得到的系数不一样。 - eipi10

有趣的是，OP使用100*ridge.fit.lambda进行“手动”计算，其结果与使用solve(t(X) %*% X + 100*ridge.fit.lambda * diag(p.tmp)) %*% t(X) %*% Y和glmnet中的ridge.fit.lambda得到的系数（几乎）完全相同。 - eipi10

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使用的是惩罚均方误差。

请注意，我没有使用你原来的代码，而是使用了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; 而tcrossprod(X)的主要用途是计算X %*% t(X)，在这种情况下会调用DSYRK而不是DGEMM。您可以阅读为什么内置的lm函数在R中如此缓慢？的第一部分了解原因和基准测试。

请注意，如果X不是方阵，则crossprod(X)和tcrossprod(X)的速度并不相同，因为它们涉及不同数量的浮点运算，您可以阅读对于对称密集矩阵乘法，是否有比“tcrossprod”更快的R函数？的附注。

关于 solvel(A, b) 和 solve(A) %*% b，请阅读 如何在不取矩阵逆的情况下高效计算 diag(X %% solve(A) %% t(X)) 的第一部分。

- Zheyuan Li

---

李哲远，非常感谢您的回答！关于使用 crossprod 和 solve(A, b) - 很有趣！您是否有任何模拟比较的参考资料？它仅在时间计算方面高效还是会产生更好的底层计算解决方案，从而在计算一些“棘手”的矩阵（非常大或具有高振幅值）时获得更高的精度？我能看到的一个缺点是它使代码不如链式 %*% 操作清晰易懂。 - Marta Karas

李哲远，感谢您的跟进！我一定会在周末查看这些链接。祝一切顺利！ - Marta Karas

Answer 2

在Zheyuan有趣的帖子基础上，我进行了更多的实验，发现我们也可以通过截距得到相同的结果，如下：

# 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

- Sandipan Dey

---

嗨sandipan，感谢您的这篇文章！我甚至惊讶于残差差异是显著的（而不是1e-10或更小阶数的超边际差异）。在这篇文章中，我实际上只关心在同一图上调整两个向量的简单性（我需要使它们的长度一致）。 - Marta Karas

@MartaKaras 但是如果你考虑到RSS的数量级，差异就非常小了，它是（110063.4-110059.9）/110063.4 = 3.179985e-05。 - Sandipan Dey