0
$\begingroup$

Doing problems from ISLR and I've taken up the task of trying to do linear regression (and by extension lasso and ridge regression) using R's optim function. Recall equation 6.5 of page 215, which states that ridge regression is the minimization of:

$$ \sum_{i = 1}^n \left( y_i - \beta_0 - \sum_{j = 1}^p \beta_j x_{ij} \right)^2 + \lambda \sum_{j = 1}^p \beta_j^2 $$

I created a minimization function based on this. The coefficients it generates do not match up very well with those given by glmnet.

In frustration, I continued searching and I found an existing working implementation:

https://github.com/m-clark/Miscellaneous-R-Code/blob/master/ModelFitting/ridge.R

The part that confuses me is the minimization function:

ridge <- function(w, X, y, lambda = .1) {
  # X: model matrix
  # y: target; lambda: penalty parameter
  # w: the weights/coefficients

  crossprod(y - X%*%w) + lambda*length(y)*crossprod(w)
}

This implementation works as expected, but notice the length(y) term used in the penalty. How did that get there? It is not present in the original equation!

Am I missing something in the original equation? Is this extra term necessary when using optim? If so, why?

============== UPDATE ==============

Here is the simulated_data code:

simulated_data <- function(seed, sample_size, predictors) {
  set.seed(seed)

  X <- scale(matrix(rnorm(sample_size*predictors), ncol = predictors))
  b <- rep(1, predictors)
  y <- scale(X %*% b + rnorm(sample_size))

  return(list(X = X, y = y))
}

Here is the ridge regression code I've written which uses glmnet:

tuple <- simulated_data(1, 100, 3)

# Load the glmnet library
library(glmnet)

# Build the ridge regression model.
# Note that alpha = 0 is for regression while alpha = 1 is for lasso.
ridge.model <- glmnet(tuple$X, tuple$y, 
                      alpha = 0, 
                      lambda = c(10, 1, .1),
                      thresh = 1e-12,
                      intercept = F)
ridge.cf <- coef(ridge.model, s = 0.1)
print(ridge.cf)

This outputs:

## 4 x 1 sparse Matrix of class "dgCMatrix"
##                     1
## (Intercept) .        
## V1          0.4357297
## V2          0.4322341
## V3          0.5070639

And here is the optim code for ridge regression:

minimize.ridge <- function(par, X, y, lambda = 0.1) {
  rss <- crossprod((X %*% par) - y)
  #penalty <- lambda*length(y)*crossprod(par)
  penalty <- lambda * (par %*% par)
  return(rss + penalty)
}

op.result <- optim(rep(0, 3), 
                   fn = minimize.ridge,
                   method = 'BFGS',
                   X = tuple$X,
               y = tuple$y,
                   lambda = 0.1)
op.cf <- op.result$par
print(op.cf)

This outputs

## [1] 0.4781309 0.4778045 0.5589391
$\endgroup$
  • 5
    $\begingroup$ The documentation for glmnet is excellent and it pays to read it all. I suspect the difference might be due to the standardization of the data, but it's impossible to tell because "do not match up very well" doesn't provide enough information. As far as your remaining question goes, if you let the symbol "$\lambda$" refer to lambda*length(y) in the code, then the code implements the formula verbatim. $\endgroup$ – whuber Apr 28 '19 at 16:35
  • $\begingroup$ I've added the full code of what i'm attempting to do. Hopefully that clears up any confusion. I will attempt to read more about glmnet but it still doesn't quite answer my question since this is actually the first step in a series of other problems I'm trying to solve. $\endgroup$ – Farley Knight Apr 28 '19 at 16:49
  • $\begingroup$ It struck me just now that perhaps if the $\lambda$ term was inside the summation, this version of the code would make sense. I just can't see why that would be the case. $\endgroup$ – Farley Knight Apr 28 '19 at 17:29
  • $\begingroup$ You were faster... Yes, if one would misread the formula as $\sum_{i = 1}^n \left( \left( y_i - \beta_0 - \sum_{j = 1}^p \beta_j x_{ij} \right)^2 + \lambda \sum_{j = 1}^p \beta_j^2\right)$ it would be as implemented. But it doesn't matter as whuber remarked. Their lambda is your $\lambda$ divided by $n$. $\endgroup$ – Kornel Apr 28 '19 at 17:33
0
$\begingroup$

I continued to dig and I found the glmnet paper. As it turns out, this version of the minimization function is not for the vanilla ridge regression that is stated in ISLR but for the elastic-net version of ridge regression that is implemented.

full version of elastic net equation

Kudos to Kornel for catching that length(y) is actually $n$ not $p$.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.