Регрессия хребта с "glmnet" дает другие коэффициенты, чем то, что я вычисляю по "определению учебника"?

Если вы прочтете ?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 default), вы можете проверить, что оценка перехвата равна ~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) не одинаково быстры, поскольку они включают различное количество операций с плавающей запятой, для которых вы можете прочитать боковое уведомление любой более быстрой функции R, чем "tcrossprod" для симметричного умножения плотной матрицы?

Относительно solvel(A, b) и solve(A) %*% b, пожалуйста, прочитайте первый раздел из как вычислить diag(X %% solve(A) %% t(X)) эффективно, не принимая матрицу обратной?