2
$\begingroup$

I'm a little bit puzzled by the default standardization of the lasso/elastic net/ridge regression algorithms implemented in the (great!) glmnet package.

In most other applications, people would standardize a data matrix $X$ by centering all columns $j$ around the mean ($\sum_{i=1}^n x_{ij}=0$) and scaling columns to unit variance, meaning:

$\frac{1}{n-1}\sum_{i=1}^nx_{ij}^2=1$. (a)

In the original glmnet and lasso papers, they state instead a normalization with factor $\frac{1}{n}$, that is, the biased estimator of the variance:

$\frac{1}{n}\sum_{i=1}^nx_{ij}^2=1$. (b)

The original ridge regression paper, on the other hand, states that $X^TX$ should have $1$s on the diagonal like a correlation matrix, leading to

$\sum_{i=1}^nx_{ij}^2=1$, (c)

without any multiplying factor. Now comes the most puzzling part, checking the Fortran code of the glmnet package in R (see lines 116 to 140), I cannot shake the feeling that they are actually standardizing such that:

$\frac{n}{n-1}\sum_{i=1}^nx_{ij}^2=1$, (d)

leading to the entries of the diagonal of $X^TX$ being just a little bit below $1$.

Now, I'm not fit in Fortran and I might read something wrong, but I've got two open questions now:

  1. Why do they use the biased estimate of the variance/standard deviation for standardization in glmnet?
  2. Or do they actually use the standardization with $\frac{n}{n-1}$?

The standardization in 2.) would actually make sense, since the diagonal of $X^TX$ would not depend on the number of observations anymore.

$\endgroup$
  • $\begingroup$ Why does it matter? stats.stackexchange.com/questions/117369/… suggests it does not. $\endgroup$ – whuber Aug 19 '19 at 21:14
  • $\begingroup$ Is diferent scaling of the matrix $X$ equivalent to different scalings of the loss function/lambda? I don't see it to be obvious if this is indeed the case... $\endgroup$ – Edgar Aug 19 '19 at 21:29
  • $\begingroup$ Just multiply the entire objective function by any constant you want: it won't change any of the solutions. Thus, the only thing that might matter is the ratio between the normalization coefficients of the design matrix and the multiplier of the penalty. But the link I gave explains why even that doesn't matter (except for comparing solutions based on different conventions). $\endgroup$ – whuber Aug 19 '19 at 21:32
  • 1
    $\begingroup$ I get the first part of your comment, but I don't see the immediate connection to the normalized design matrix. Anyhow, I'm now quite convinced that scaling X by a factor c will lead to some scaling of lambda but won't change the solution, I will try to work this out myself $\endgroup$ – Edgar Aug 19 '19 at 21:58
  • $\begingroup$ Could you tell me why they standardize the variables with the (1/N) variance estimator formula in glmnet package? Is it to mimic the MLE of the variance? If I use scale() function in R to standardize the variables by the 1/(N-1) formula, input the standadized data and use option standardize=FALSE, the de-standardized coefficients will be slight different from what the default method (input data in original scale with standardize=TRUE) would produce. I'm wondering if it is theoretically okay to use scale to standardize the data first and use them as the inputs. Could anyone help? $\endgroup$ – Blain Waan Jan 23 at 13:14
2
$\begingroup$

@whuber's comments led me to the solution: Yes, it doesn't matter. Scaling of the design matrix by a constant does not change the underlying optimization problem.

More precisely:

Suppose we are given some data matrix $X_0$ and a scaled design matrix $X=c\cdot X_0$, where $c>0$ is some real constant. E.g. $X_0$ could be such that (c) holds: $X_0^TX_0$ has $1$s on the diagonal, and with $c=\sqrt{1/n}$ we have $X$ such that (b) holds.

The glmnet optimization problem for $X_0$ would be to find the $\beta_0$ that minimizes $||y-X_0\beta_0||_2^2+\lambda_0||\beta_0||_p$, where $p$ is either $1$ (lasso) or $2$ (ridge). Elastic net is a combination of lasso and ridge, so I'm pretty sure the next step holds true anyway.

We have now the following: $$ ||y-X_0\beta_0||_2^2+\lambda_0||\beta_0||_p\\ = ||y-cX_0\cdot \frac{\beta_0}{c}||_2^2+\lambda_0c||\frac{\beta_0}{c}||_p\\ = ||y-X\beta||_2^2+\lambda||\beta||_p $$ where $\lambda=\lambda_0\cdot c$ and $\beta=\beta_0/c$. That means the optimization problem is the same for both $X_0$ and $X$, but the lambdas and betas will live on a different scale. The solution for either $X_0$ or $X$ is equivalent to the solution for the other one.

$\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.