7
$\begingroup$

I'm struggling to grasp the intuition behind why LASSO can only select at most $n$ features when $n << p$, where $n$ is the number of samples and $p$ is the number of features.

I've read through the other questions asking this problem on SE, e.g., 38299 and 386116, but their answers don't complete the picture for me.

In essence, I'm struggling to find an intuitive answer. I understand the geometric constraints of LASSO in the context of feature selection itself, but reasoning behind a maximum of $n$ features being selected evades me.

The comments on 376318 inspired me to solve a minimal simulation with $x = [1, 2]$ and $y = 5$ using $\lambda = 1$. A grid search of solutions to this problem doesn't suggest that the optimal solution is found with only one coefficient. Instead, the optimal solution includes both coefficients for the two features in $x$: $\beta \approx [0.102, 2.347]$. This breaks the $n$ feature selection logic (although I'm likely missing an important detail in my implementation; Solution PDF and R Code).

enter image description here

Any help or source to an intuitive explanation behind the the maximum of $n$ feature selection aspect of LASSO would be greatly appreciated.

Perhaps a starting point for a specific question would be: is the max $n$ feature selection "baked into" the penalty function of LASSO, or is it a component of the algorithms used to solve it?

$\endgroup$
6
  • $\begingroup$ Welcome to Cross Validated! Are you sure the claim is true, even if the penalty is very small? $\endgroup$
    – Dave
    Commented Nov 20, 2023 at 20:28
  • $\begingroup$ Thanks @Dave. Just to clarify, do you mean even if the $\lambda$ is small? $\endgroup$ Commented Nov 20, 2023 at 20:32
  • $\begingroup$ Different people use different conventions, though when I have seen people use $\lambda$, I have tended to see a small penalty correspond with small $\lambda$ (and large penalty with large $\lambda$), yes. $\endgroup$
    – Dave
    Commented Nov 20, 2023 at 20:34
  • $\begingroup$ Doesn't the Lasso automatically include an intercept? $\endgroup$
    – whuber
    Commented Nov 20, 2023 at 20:37
  • $\begingroup$ @Dave, even as I modify the $\lambda$ parameter to be very low (1E-6) and very high (1E6), the solutions still have both coefficients set as non-zero. $\endgroup$ Commented Nov 20, 2023 at 20:38

2 Answers 2

6
$\begingroup$

I think you can understand why LASSO selects at most $n$-features intuitively if you think about the contours of $rss = |X\beta - y|_2^2$ and $l_1 = |\beta|_1$ on the parameter space $\mathbb R^p$. You can think about finding the minimal $rss + l_1$ solution as looking at the 2 contours of larger and larger values for $rss$ and $l_1$ until they intersect for the first time. Their growth with respect to each other needs a specific ratio, but that's irrelevant for the intuition, in fact in your head you might want to keep the $l_1$-contour constant.

The $l_1$-contour will be a generalized octahedron (https://en.wikipedia.org/wiki/Octahedron) with the corners on the axis. For your example with $p = 2$ that's just a square with sides 45 degrees rotated against the axis. For $rss$ let's first think about $rss = 0$, so the solution space to $X\beta = y$. That's a $p-n$-dimensional affine subspace. In your example that's the line $\beta_2 = -0.5\beta_1 + 2.5$. $rss$ is just distance squared so constant $rrs \implies$ constant distance. Therefor the contours will be parallel(or a elliptic, symmetric arrangement of parallel spaces) to the original subspace.

Let's check it out in your example (R-code, at the bottom, I improved some stuff and choose $\lambda = 10$ for visibility)enter image description here:

The solution is at $\beta_1 = 0, \beta_2 = 1.2$. In red you can see the contours of $rss$ and in orange the contours of $l_1$. The question comes down to how can the red lines be tangent to an orange square?

  1. By going through one of the corners, which is on an axis, so the other parameter will be $0$.

  2. By being parallel to one of the sides and touching the whole side, in which case you have non unique solutions (https://www.stat.cmu.edu/~ryantibs/papers/lassounique.pdf) This is partially discussed in some of the threads you linked. If you change $x_2 = 1$ in your example, then you will have an obvious case.

Generalizing: the $l_1$-contour is made up of $d$-dimensional edges ($d = 0$ is a corner, $d = 1$ is an edge, $d = 2$ is a 2d-surface like the face for $p = 3$, ...) and hitting that $d$-dimensional edge correspond to selecting $d+1$ features. The $rss$-contour covers $p-n$-dimension and the $n$-dimensional edge, that would select $n+1$ features, covers the remaining dimensions, so there is no space to not be parallel.

In $p = 3$ you can check this out with an 8-sided-die for the $l_1$-contour and for $n = 1$ a large piece of paper or for $n = 2$ a pencil, to stand on for the plane/lines that make up the $rss$-contours. A 6-sided-die is also fine for getting the idea.

# Implemented the fix from Eoins answer
# Comment on your Code: you are loading a library, use that for 
# readability instead of andless ::
library(tidyverse)
data <- data.frame(x1 = c(1),
                   x2 = c(2),
                   y = c(5))
range <- 5
b1 <- seq(-range, range, length.out = 51)
b2 <- seq(-range, range, length.out = 51)
model <- expand.grid(b1 = b1, b2 = b2) 

compute_loss <- function(model, data, lambda) {
  y_hat <- model["b1"] * data$x1 + model["b2"] * data$x2
  rss <- mean((data$y - y_hat)^2)
  l1 <- lambda * sum(abs(model))
  
  data.frame(rss = rss,
             l1 = l1,
             lasso = rss + l1,
             b1 = model["b1"],
             b2 = model["b2"],
             y_hat = y_hat)
}
lambda <- 10 # l1 closer to rss values
costs <- apply(model, 1, compute_loss, data, lambda) %>% 
  bind_rows()
costs %>%filter(lasso == min(lasso))
# rss  l1 lasso  b1 b2   y_hat
# 6.76 12 18.76  0 1.2   2.4 
# making about with the breaks argument, to show 
ggplot(costs, aes(x = b1, y = b2, fill = lasso))+
  geom_raster() +
  geom_contour(aes(z = rss), breaks = c(0.01, sqrt(6.76/2), 6.76), 
                     color = "red") +
  geom_abline(intercept = 2.5, slope = -0.5, color = "red", ) +
  geom_contour(aes(z = l1), breaks = c(0.01, 6, 12), 
                   color = "orange") +
  geom_point(data = costs %>%filter(lasso == min(lasso)), size = 3, 
              shape = 1) +
  scale_fill_viridis_c(alpha = 0.7) 
$\endgroup$
8
  • $\begingroup$ I think you can understand why LASSO selects at most n -features intuitively Doesn't the answer by Eoin show that more than $n$ features can selected? $\endgroup$
    – Dave
    Commented Nov 21, 2023 at 13:11
  • 3
    $\begingroup$ Keep in mind that no matter what the sample size, the probability that lasso selects the right features is zero. $\endgroup$ Commented Nov 21, 2023 at 13:15
  • 1
    $\begingroup$ @Dave, no, matter of fact Eoins answer explains why the Code linked in the question couldn't find the correct solution with $b_1 = 0$. $\endgroup$ Commented Nov 21, 2023 at 13:22
  • 2
    $\begingroup$ @LukasLohse the geometry is beautiful because it solves a problem by oversimplification. Contrast that with full Bayesian modeling using more reasonable shrinkage priors, or using data reduction/unsupervised learning to help you live within your sample size budget. $\endgroup$ Commented Nov 21, 2023 at 22:28
  • 3
    $\begingroup$ @FrankHarrell harsh words, but i watched your talk(fharrell.com/talk/stratos19) and read some stuff and I definitely underestimated the ways LASSO could fail, never mind questioned the idea of feature selection as a whole. Thank you for fighting the good fight, i guess:) $\endgroup$ Commented Nov 22, 2023 at 14:36
7
$\begingroup$

This is an excellent question, and approach to building intuition. Unfortunately, the answer is silly: your grid search doesn't include cases where $b_1$ or $b_2 = 0$:

range <- 5
seq(-range, range, length.out = 50)
#  [1] -5.0000000 -4.7959184 -4.5918367 -4.3877551 -4.1836735 -3.9795918 -3.7755102 -3.5714286
#  [9] -3.3673469 -3.1632653 -2.9591837 -2.7551020 -2.5510204 -2.3469388 -2.1428571 -1.9387755
# [17] -1.7346939 -1.5306122 -1.3265306 -1.1224490 -0.9183673 -0.7142857 -0.5102041 -0.3061224
# [25] -0.1020408  0.1020408  0.3061224  0.5102041  0.7142857  0.9183673  1.1224490  1.3265306
# [33]  1.5306122  1.7346939  1.9387755  2.1428571  2.3469388  2.5510204  2.7551020  2.9591837
# [41]  3.1632653  3.3673469  3.5714286  3.7755102  3.9795918  4.1836735  4.3877551  4.5918367
# [49]  4.7959184  5.0000000

Change this to length.out = 51 and things work as you would expect.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.