L0-pseudonorm penalized least squares regression (aka best subset regression) solves $\widehat{\beta}(\lambda)$ as $$\min_\beta \frac{1}{2}||y-X\beta||_2^2 +\lambda||\beta||_0.$$ where $||\beta||_0$ is the number of nonzero coefficients. I was wondering what would be (1) the minimum value of $\lambda$ that would result in no variables being selected and (2) the maximum value of $\lambda$ that would result in the maximum number of variables being selected, either for the case where the coefficients are unconstrained or where they are nonnegativity constrained (ie required to be all zero or positive as in nnls)? For LASSO regression, where we work with the L1-norm penalty $\lambda||\beta||_1$ I understand that (1) is given by $\lambda_1 = \max_j |X_j^Ty|$, but what would be it's value in case of L0-penalized regression (as implemented in the L0Learn package)?

Example in R:

# Simulate some data
data <- GenSynthetic(n=500,p=100,k=10,seed=1)
X = data$X
y = data$y
# make L0 penalized fit:
system.time(L0fit <- L0Learn.fit(x=X, y=y, penalty="L0", algorithm="CDPSI", nLambda=1000, intercept=F, maxSuppSize = 100)) 

Maximum lambda that would result in no variables being selected = 0.0618124 :

unlist(L0fit$lambda)[unlist(L0fit$suppSize)==0][1] # = 0.0618124

Lambda that would result in the maximum number of variables (100 here, ie all variables) being selected = 6.5916e-09 :

unlist(L0fit$lambda)[unlist(L0fit$suppSize)==max(unlist(L0fit$suppSize))][1] # = 6.5916e-09
max(unlist(L0fit$suppSize)) # size of largest model = 100

So I'm looking for a way to calculate those two lambda values - here 0.0618124 and 6.5916e-09 - a priori. For the 0.0618124 I tried with the recipe in the answer below but couldn't quite reproduce this value - instead of 0.0618124 I am getting 677 in my example:

max(diag(1/crossprod(X, X)) * (crossprod(X, y)^2)) # 677.1252

This paper ("Efficient Regularized Regression with L0 Penalty for Variable Selection and Network Construction", 2016, by Liu & Li, page 6) mentions a maximum $\lambda^\star = \max_{i = 1,\dots,p}~ (X^\top_i y)^2 / (4X^\top_i X_i)$ but again that seems to give a different value... Any thoughts?

EDIT: So it seems that L0Learn first centers & L2 norm normalizes both the columns of the design matrix & the outcome variable y. Hence, the maximum lambda that would cause all variables to be penalized, based on the logic in the answers below, in L0Learn is given by

Xcentnorm = apply(X, 2, function (col) (col-mean(col))/norm((col-mean(col)),"2"))
ycentnorm = (y-mean(y))/(norm(y-mean(y),"2"))
max((crossprod(Xcentnorm, ycentnorm)^2)/2) # = 0.06262011

The factor diag(1/crossprod(Xcentnorm, Xcentnorm)) drops out because of the L2 norm normalization (i.e. it would be a vector of 1s).


Assuming the columns have unit L2 norm, the $\lambda^{*}$ which sets all coefficients to zero is given by $\frac{1}{2} \max_{j} (X_j^T y)^2$ (the reasoning in the answer above is correct; but the final answer misses the factor of $\frac{1}{2}$).

L0Learn centers and then normalizes the columns before fitting the model. The $\lambda$'s are reported after centering and normalization. So to reproduce L0Learn's $\lambda^{*}$ you can try centering and then normalizing the columns.

| cite | improve this answer | |
  • $\begingroup$ So in my example would max((crossprod(apply(X, 2, function (col) (col-mean(col))/norm((col-mean(col)),"2")), (y-mean(y))/(norm(y-mean(y),"2")))^2)/2) = 0.06262011 correspond to the max lambda of L0Learn? I only get a result close to 0.06 if I also center & L2 norm normalize the y values - is that what L0Learn does? $\endgroup$ – Tom Wenseleers Jul 12 '19 at 15:07
  • $\begingroup$ Right it also centers/normalizes y. The small difference you're observing is probably because L0Learn uses a lambda slightly smaller than lambda*, so that only 1 variable is non-zero. I will add a note to the documentation to clarify the normalization process. $\endgroup$ – Hussein Hazimeh Jul 12 '19 at 17:07
  • $\begingroup$ Ha that's great - thanks for clarifying this! $\endgroup$ – Tom Wenseleers Jul 13 '19 at 16:30

To answer the first part of your question: No variable is selected, if it is optimal to not change any element of $\beta$ away from zero, i.e. $y^\top y \leq (y-x_i \hat{\beta_i})^\top (y-x_i \hat{\beta_i}) + 2 \lambda ~ \forall i$. For simplicity, I always only consider one of the parameters at a time (which should be fine due to linearity). Now, if we only have the variable $x_i$, the optimal coefficent is given by $\hat{\beta_i} = (X^\top_i X_i)^{-1} X^\top_i y$. Plugging this into the inequality above, results in $\lambda \geq (X^\top_i X_i)^{-1} (X^\top_i y)^2$. If this inequality holds for all $i$, all coefficents should be equal to zero. Therefore: $\lambda^\star = \max_{i = 1,\dots,p}~ (X^\top_i X_i)^{-1} (X^\top_i y)^2$ if I'm not mistaken.

| cite | improve this answer | |
  • $\begingroup$ Thanks for that - makes sense - but when I tried your last formula I couldn't quite reproduce the result that L0Learn gives me in the numerical example that I added above. Any thoughts? $\endgroup$ – Tom Wenseleers Jul 5 '19 at 12:11

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.