# How is the intercept computed in GLMnet?

I've been implementing the GLMNET version of elastic net for linear regression with another software than R. I compared my results with the R function glmnet in lasso mode on diabetes data.

The variable selection is ok when varying the value of the parameter (lambda) but I obtain slightly different values of coefficients. For this and other reasons I think it comes from the intercept in the update loop, when I compute the current fit, because I don't vary the intercept (which I take as the mean of the target variable) in the whole algorithm : as explained in Trevor Hastie's article ( Regularization Paths for Generalized Linear Models via Coordinate Descent, Page 7, section 2.6):

the intercept is not regularized, [...] for all values of [...] lambda [the L1-constraint parameter]

But despite the article, the R function glmnet does provide different values for the intercept along the regularization path (the lambda different values). Does anyone has a clue about how the values of the Intercept are computed?

I found that the intercept in GLMnet is computed after the new coefficients updates have converged. The intercept is computed with the means of the $y_i$'s and the mean of the $x_{ij}$'s. The formula is siimilar to the previous one I gave but with the $\beta_j$'s after the update loop : $\beta_0=\bar{y}-\sum_{j=1}^{p} \hat{\beta_j} \bar{x_j}$.

In python this gives something like :

        self.intercept_ = ymean - np.dot(Xmean, self.coef_.T)


which I found here on scikit-learn page.

EDIT : the coefficients have to be standardized before :

        self.coef_ = self.coef_ / X_std


$\beta_0=\bar{y}-\sum_{j=1}^{p} \frac{\hat{\beta_j} \bar{x_j}}{\sum_{i=1}^{n} x_{ij}^2}$.

• I should add that this is a standard way of calculating the intercept, assuming that model is linear and the errors have zero mean. Aug 3, 2011 at 12:36
• Indeed, nervertheless the authors explicitly said in their article : $\beta_0=\bar{y}$ for all values of $\alpha$ and $\lambda$, and moreover they don't say at which point of the algorithm it is computed
– yelh
Aug 3, 2011 at 12:44
• Since authors assume that $\bar x=0$ then this is true. Note at the beginning in the page 3, they say that we assume that predictors are centered, but this is not restrictive, since the "... results generalize naturally..." to unstandartized case. The formula you find is this natural generalisation. Aug 3, 2011 at 12:52
• However, even when standardization is applied (and therefore mean centering) on the predictors in the algorithm, they take the unstandardized data (therefore $\bar{x} \ne 0$, in general) to fit the different intercepts displayed by glmnet. They do use the same $\beta_0=constant=\bar{y}$ for updating the coefficients but fit an intercept with the raw data, and they do it a posteriori.
– yelh
Aug 3, 2011 at 13:24

which I take as the mean of the target variable

I think this may be where you're going wrong: unlike the linear model, you can't reparameterise the predictors such that they will always be orthogonal to the intercept, hence the intercept cannot just be calculated as the mean.

• I took the mean of the explained variable because in the article i quoted, the authors of this method write that they do use the mean of the Y_i (the explained variable observations) for all values of alpha and lambda. Now, looking at the output of glmnet function, i guess it's not the case for all lambdas. So that does not tell me how to compute the intercepts for each regularization parameter lambda.
– yelh
Jul 29, 2011 at 13:33
• My guess is that they use the weighted mean (as the weights will change with $\lambda$). Jul 29, 2011 at 14:15
• Could you precise your idea please ? Which expression of weights should I consider ?
– yelh
Jul 29, 2011 at 15:47
• In the paper, each "inner loop" is a penalised, weighted least-squares problem, with the weights being defined by the previous fitted value via equation 17. Jul 30, 2011 at 10:05
• It's ok I found the answer looking at scikit-learn python code (because the glmnet source code is in Fortran and it is not my cup of tea). I will share it later if anyone is interested. Thanks anyway !
– yelh
Aug 2, 2011 at 14:16