我在使用 glmnet R 包进行岭回归时发现，通过 glmnet::glmnet 函数得到的系数与我使用相同 lambda 值通过定义计算的系数不同。能有人解释一下这是为什么吗？
数据（响应变量 Y 和设计矩阵 X）均已标准化。
library(MASS)library(glmnet)# 数据维度p.tmp <- 100n.tmp <- 100# 数据对象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)] <- 10Y.true <- X %*% betaY <- scale(Y.true + matrix(rnorm(n.tmp))) # Y.true + 高斯噪声# 运行 glmnet ridge.fit.cv <- cv.glmnet(X, Y, alpha = 0)ridge.fit.lambda <- ridge.fit.cv$lambda.1se# 提取 lambda.1se 的系数值（不含截距）ridge.coef <- (coef(ridge.fit.cv, s = ridge.fit.lambda))[2:(p.tmp+1)]# 通过定义获取系数ridge.coef.DEF <- solve(t(X) %*% X + ridge.fit.lambda * diag(p.tmp)) %*% t(X) %*% Y# 绘制估计值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) %*% Ydrop(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) %*% Ydrop(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)# 数据维度p.tmp <- 100n.tmp <- 100# 数据对象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)] <- 10Y.true <- X %*% betaY <- scale(Y.true + matrix(rnorm(n.tmp)))# 运行 glmnet ridge.fit.cv <- cv.glmnet(X, Y, alpha = 0, intercept = FALSE)ridge.fit.lambda <- ridge.fit.cv$lambda.1se# 提取 lambda.1se 的系数值（不含截距）ridge.coef <- (coef(ridge.fit.cv, s = ridge.fit.lambda))[-1]# 通过定义获取系数ridge.coef.DEF <- drop(solve(crossprod(X) + diag(n.tmp * ridge.fit.lambda, p.tmp), crossprod(X, Y)))# 绘制估计值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。即使使用 intercept = TRUE（glmnet 的默认设置），你也可以检查到截距的估计值为 ~e-17（数值上为0），因此其他系数的估计值不会受到显著影响。另一个答案只是展示了这一点。
后续
关于使用
crossprod和solve(A, b)– 很有趣！你是否碰巧有任何关于此的模拟比较参考？
t(X) %*% Y 将首先对 X 进行转置 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 函数用于对称密集矩阵乘法？
关于 solve(A, b) 和 solve(A) %*% b，请阅读 第一部分 的 如何高效计算 diag(X %*% solve(A) %*% t(X)) 而不需要取矩阵的逆？
