# Replicating results for glmnet linear regression using a generic optimizer

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)
}

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)))  • I'm not sure your question is on topic (I think it might be, as it's about the underlying optimization technique), and I can't really check your code now, but lbfgs raises a point about the orthantwise_c argument regarding glmnet equivalence. – Firebug Sep 25 '16 at 22:08 • The problem is not really with lbfgs and orthantwise_c, as when alpha = 1, the solution is near exactly the same with glmnet. It has to do with the L2 regularization side of things i.e. when alpha < 1. I think making some kind of modification to the definition of SSE and SSE_gr should fix it, but I'm not sure what the modification should be - as far as I know, those functions are defined exactly as described in the glmnet paper. – user3294195 Sep 25 '16 at 22:25 • This may be more of a stackoverflow, programming question. – Matthew Gunn Sep 27 '16 at 22:03 • I thought it has more to do with optimization & regularization rather than the code itself, which is why I posted it here. – user3294195 Sep 27 '16 at 22:08 • For a pure optimization question, scicomp.stackexchange.com is also an option. I am not sure if language specific questions are on topic there though? (e.g. "do this in R") – GeoMatt22 Sep 28 '16 at 3:57 ## 1 Answer ## 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):