Replicating results for glmnet linear regression using a generic optimizer

Question

As the title states, I'm trying to replicate the results from glmnet linear using the LBFGS optimizer from library lbfgs. This optimizer allows us to add an L1 regularizer term without having to worry about differentiability, as long as our objective function (without the L1 regularizer term) is convex.

The elastic net linear regression problem in the glmnet paper is given by $$\min_{\beta \in \mathbb{R}^p} \frac{1}{2n}\Vert \beta_0 + X\beta - y \Vert_2^2 + \alpha \lambda \Vert \beta\Vert_1 + \frac{1}{2}(1-\alpha)\lambda\Vert\beta\Vert^2_2$$ where $X \in \mathbb{R}^{n \times p}$ is the design matrix, $y \in \mathbb{R}^p$ is the vector of observations, $\alpha \in [0,1]$ is the elastic net parameter, and $\lambda > 0$ is the regularization parameter. The operator $\Vert x \Vert_p $ denotes the usual Lp norm.

The below code defines the function, and then includes a test to compare the results. As you can see, the results are acceptable when alpha = 1, but are way off for values of alpha < 1. The error gets worse as we go from alpha = 1 to alpha = 0, as the following plot shows (the "comparison metric" is the mean Euclidean distance between the parameter estimates of glmnet and lbfgs for a given regularization path).

Okay, so here's the code. I've added comments wherever possible. My question is: Why are my results different from those of glmnet for values of alpha < 1? It clearly has to do with the L2 regularization term, but as far as I can tell, I've implemented this term exactly as per the paper. Any help would be much appreciated!

library(lbfgs)
linreg_lbfgs <- function(X, y, alpha = 1, scale = TRUE, lambda) {
  p <- ncol(X) + 1; n <- nrow(X); nlambda <- length(lambda)

  # Scale design matrix
  if (scale) {
    means <- colMeans(X)
    sds <- apply(X, 2, sd)
    sX <- (X - tcrossprod(rep(1,n), means) ) / tcrossprod(rep(1,n), sds)
  } else {
    means <- rep(0,p-1)
    sds <- rep(1,p-1)
    sX <- X
  }
  X_ <- cbind(1, sX)

  # loss function for ridge regression (Sum of squared errors plus l2 penalty)
  SSE <- function(Beta, X, y, lambda0, alpha) {
    1/2 * (sum((X%*%Beta - y)^2) / length(y)) +
      1/2 * (1 - alpha) * lambda0 * sum(Beta[2:length(Beta)]^2) 
                    # l2 regularization (note intercept is excluded)
  }

  # loss function gradient
  SSE_gr <- function(Beta, X, y, lambda0, alpha) {
    colSums(tcrossprod(X%*%Beta - y, rep(1,ncol(X))) *X) / length(y) + # SSE grad
  (1-alpha) * lambda0 * c(0, Beta[2:length(Beta)]) # l2 reg grad
  }

  # matrix of parameters
  Betamat_scaled <- matrix(nrow=p, ncol = nlambda)

  # initial value for Beta
  Beta_init <- c(mean(y), rep(0,p-1)) 

  # parameter estimate for max lambda
  Betamat_scaled[,1] <- lbfgs(call_eval = SSE, call_grad = SSE_gr, vars = Beta_init, 
                              X = X_, y = y, lambda0 = lambda[2], alpha = alpha,
                              orthantwise_c = alpha*lambda[2], orthantwise_start = 1, 
                              invisible = TRUE)$par

  # parameter estimates for rest of lambdas (using warm starts)
  if (nlambda > 1) {
    for (j in 2:nlambda) {
      Betamat_scaled[,j] <- lbfgs(call_eval = SSE, call_grad = SSE_gr, vars = Betamat_scaled[,j-1], 
                                  X = X_, y = y, lambda0 = lambda[j], alpha = alpha,
                                  orthantwise_c = alpha*lambda[j], orthantwise_start = 1, 
                                  invisible = TRUE)$par
    }
  }

  # rescale Betas if required
  if (scale) {
    Betamat <- rbind(Betamat_scaled[1,] -
colSums(Betamat_scaled[-1,]*tcrossprod(means, rep(1,nlambda)) / tcrossprod(sds, rep(1,nlambda)) ), Betamat_scaled[-1,] / tcrossprod(sds, rep(1,nlambda)) )
  } else {
    Betamat <- Betamat_scaled
  }
  colnames(Betamat) <- lambda
  return (Betamat)
}

# CODE FOR TESTING
# simulate some linear regression data
n <- 100
p <- 5
X <- matrix(rnorm(n*p),n,p)
true_Beta <- sample(seq(0,9),p+1,replace = TRUE)
y <- drop(cbind(1,X) %*% true_Beta)

library(glmnet)

# function to compare glmnet vs lbfgs for a given alpha
glmnet_compare <- function(X, y, alpha) {
  m_glmnet <- glmnet(X, y, nlambda = 5, lambda.min.ratio = 1e-4, alpha = alpha)
  Beta1 <- coef(m_glmnet)
  Beta2 <- linreg_lbfgs(X, y, alpha = alpha, scale = TRUE, lambda = m_glmnet$lambda)
  # mean Euclidean distance between glmnet and lbfgs results
  mean(apply (Beta1 - Beta2, 2, function(x) sqrt(sum(x^2))) ) 
}

# compare results
alpha_seq <- seq(0,1,0.2)
plot(alpha_seq, sapply(alpha_seq, function(alpha) glmnet_compare(X,y,alpha)), type = "l", ylab = "Comparison metric")

---

@hxd1011 I tried your code, here are some tests (I made some minor tweaks to match the structure of glmnet - note we do not regularize the intercept term, and the loss functions must be scaled). This is for alpha = 0, but you can try any alpha - the results do not match.

rm(list=ls())
set.seed(0)
# simulate some linear regression data
n <- 1e3
p <- 20
x <- matrix(rnorm(n*p),n,p)
true_Beta <- sample(seq(0,9),p+1,replace = TRUE)
y <- drop(cbind(1,x) %*% true_Beta)

library(glmnet)
alpha = 0

m_glmnet = glmnet(x, y, alpha = alpha, nlambda = 5)

# linear regression loss and gradient
lr_loss<-function(w,lambda1,lambda2){
  e=cbind(1,x) %*% w -y
  v= 1/(2*n) * (t(e) %*% e) + lambda1 * sum(abs(w[2:(p+1)])) + lambda2/2 * crossprod(w[2:(p+1)])
  return(as.numeric(v))
}

lr_loss_gr<-function(w,lambda1,lambda2){
  e=cbind(1,x) %*% w -y
  v= 1/n * (t(cbind(1,x)) %*% e) + c(0, lambda1*sign(w[2:(p+1)]) + lambda2*w[2:(p+1)])
  return(as.numeric(v))
}

outmat <- do.call(cbind, lapply(m_glmnet$lambda, function(lambda) 
  optim(rnorm(p+1),lr_loss,lr_loss_gr,lambda1=alpha*lambda,lambda2=(1-alpha)*lambda,method="L-BFGS")$par
))

glmnet_coef <- coef(m_glmnet)
apply(outmat - glmnet_coef, 2, function(x) sqrt(sum(x^2)))

Answer 1

tl;dr version:

The objective implicitly contains a scaling factor $\hat{s} = sd(y)$, where $sd(y)$ is the sample standard deviation.

Longer version

If you read the fine print of the glmnet documentation, you will see:

Note that the objective function for ‘"gaussian"’ is

               1/2  RSS/nobs + lambda*penalty,                  

and for the other models it is

               -loglik/nobs + lambda*penalty.                   

Note also that for ‘"gaussian"’, ‘glmnet’ standardizes y to have unit variance before computing its lambda sequence (and then unstandardizes the resulting coefficients); if you wish to reproduce/compare results with other software, best to supply a standardized y.

Now this means that the objective is actually $$ \frac{1}{2n} \lVert y/\hat{s}-X\beta \rVert_2^2 + \lambda\alpha \lVert \beta \rVert_1 + \lambda(1-\alpha)\lVert \beta \rVert_2^2, $$ and that glmnet reports $\tilde{\beta} = \hat{s} \beta$.

Now, when you were using a pure lasso ($\alpha=1$), then unstandardization of glmnet's $\tilde{\beta}$ means that the answers are equivalent. On the other hand, with a pure ridge, then you need to scale the penalty by a factor $1/\hat{s}$ in order for the glmnet path to agree, because an extra factor of $\hat{s}$ pops out from the square in the $\ell_2$ penalty. For intermediate $\alpha$, there is not an easy way to scale the penalty of coefficients to reproduce glmnets output.

Once I scale the $y$ to have unit variance, I find

which still doesn't match exactly. This seems to be due to two things:

1. The lambda sequence may be too short for the warm-start cyclic coordinate descent algorithm to be fully convergent.
2. There's no error term in your data (the $R^2$ of the regression is 1).
3. Note also there's a bug in the code as provided in which it takes lambda[2] for the initial fit, but that should be lambda[1].

Once I rectify items 1-3, I get the following result (though YMMV depending on the random seed):