LASSO Regression


In-depth Articles

Introduction

The Least Absolute Shrinkage and Selection Operator (LASSO) was introduced by Tibshirani (1996). It is a regularization technique that performs both continuous shrinkage and automatic variable selection, producing sparse, interpretable models.

LASSO is particularly useful when:




The LASSO Objective

The LASSO estimator is defined as the solution to the following penalized least squares problem:

\[ \hat{\beta}^{\text{lasso}} = \arg\min_{\beta} \left\{ \frac{1}{2n} \|y - X\beta\|_2^2 + \lambda \|\beta\|_1 \right\} \]

Expanding:

\[ \hat{\beta}^{\text{lasso}} = \arg\min_{\beta} \left\{ \frac{1}{2n} \sum_{i=1}^{n} \left( y_i - x_i'\beta \right)^2 + \lambda \sum_{j=1}^{p} |\beta_j| \right\} \]

where:

The equivalent constrained formulation is:

\[ \min_{\beta} \frac{1}{2n} \|y - X\beta\|_2^2 \quad \text{subject to} \quad \sum_{j=1}^{p} |\beta_j| \leq t \]

where \(t\) is in one-to-one correspondence with \(\lambda\).




Why L1 Gives Sparsity

Geometric Argument

The \(L_1\) ball \(\{\ \beta : \|\beta\|_1 \leq t\ \}\) is a diamond (in 2D) with corners on the coordinate axes. The contours of the RSS loss function are ellipses centered at the OLS solution. The constrained optimum occurs where the ellipse first touches the constraint region. Because the \(L_1\) ball has corners, this contact point is likely to occur at a corner, where one or more coordinates are exactly zero.

In contrast, the \(L_2\) ball (Ridge) is smooth — contact points generically have all coordinates nonzero.

Subdifferential and Subgradient Condition

The \(L_1\) norm is not differentiable at zero. The subdifferential of \(|\beta_j|\) is:

\[ \partial |\beta_j| = \begin{cases} \{+1\} & \text{if } \beta_j > 0 \\ \{-1\} & \text{if } \beta_j < 0 \\ [-1, +1] & \text{if } \beta_j = 0 \end{cases} \]

The optimality (KKT) condition for the LASSO is:

\[ -\frac{1}{n} x_j'(y - X\beta) + \lambda \cdot s_j = 0, \quad s_j \in \partial|\beta_j| \]

When \(\beta_j = 0\), this requires \(\left| \frac{1}{n} x_j'(y - X\beta) \right| \leq \lambda\), which is satisfiable — hence exact zeros are possible.

Soft-Thresholding Operator

For orthonormal design (\(X'X = nI\)), the LASSO solution has a closed form via the soft-thresholding operator:

\[ \hat{\beta}_j = S_\lambda(\hat{\beta}_j^{\text{OLS}}) = \text{sign}(z_j)(|z_j| - \lambda)_+ \]

where \(z_j = \hat{\beta}_j^{\text{OLS}}\) and \((x)_+ = \max(0, x)\). Coefficients with \(|z_j| \leq \lambda\) are set exactly to zero; others are shrunk toward zero by \(\lambda\).




Coordinate Descent Algorithm

The standard algorithm for solving the LASSO is coordinate descent (Friedman, Hastie, Tibshirani, 2010):

Algorithm:

  1. Initialize \(\beta^{(0)}\) (e.g., zeros or warm start from previous \(\lambda\))
  2. For each coordinate \(j = 1, \ldots, p\):

    \[ \tilde{\beta}_j \leftarrow S_\lambda\left( \frac{1}{n} x_j'(y - X_{-j}\beta_{-j}) \right) \bigg/ \left( \frac{1}{n} x_j'x_j \right) \]

    where \(X_{-j}\beta_{-j}\) is the partial residual excluding predictor \(j\).
  3. Repeat step 2 (cycling through all \(j\)) until convergence: \(\max_j |\beta_j^{(k+1)} - \beta_j^{(k)}| < \epsilon\)

Key properties:




The Regularization Path

The LASSO regularization path traces the solutions \(\hat{\beta}(\lambda)\) as \(\lambda\) varies from \(\lambda_{\max}\) down to 0:

\[ \lambda_{\max} = \frac{1}{n} \|X'y\|_\infty = \max_j \left| \frac{1}{n} x_j'y \right| \]

This is the smallest \(\lambda\) for which all coefficients are zero. As \(\lambda\) decreases:

LARS (Least Angle Regression)

Efron et al. (2004) showed that the entire LASSO path can be computed efficiently using the LARS algorithm:

  • Start with all coefficients at zero
  • Find the predictor most correlated with the residual
  • Move \(\beta_j\) in the direction of correlation until another predictor has equal correlation
  • Move jointly in the equiangular direction
  • If a coefficient hits zero, drop it (LASSO modification of LARS)

The computational cost is similar to a single OLS fit: \(O(np \cdot \min(n,p))\) for the entire path.




Choosing \(\lambda\)

K-Fold Cross-Validation

The standard approach is K-fold CV (typically \(K=10\)):

  1. For each \(\lambda\) in a grid, compute the CV error: \[ \text{CV}(\lambda) = \frac{1}{K} \sum_{k=1}^{K} \frac{1}{|D_k|} \sum_{i \in D_k} (y_i - x_i'\hat{\beta}^{(-k)}(\lambda))^2 \]
  2. Two common choices:
    • \(\lambda_{\min}\): the \(\lambda\) minimizing CV error
    • \(\lambda_{1\text{se}}\): the largest \(\lambda\) within one standard error of the minimum — the most regularized (parsimonious) model whose performance is statistically indistinguishable from the best

Information Criteria

Alternatives to CV:

  • BIC: \(\text{BIC}(\lambda) = n \log(\text{RSS}/n) + \log(n) \cdot \text{df}(\lambda)\)
  • AIC: \(\text{AIC}(\lambda) = n \log(\text{RSS}/n) + 2 \cdot \text{df}(\lambda)\)

where \(\text{df}(\lambda) = |\{j : \hat{\beta}_j(\lambda) \neq 0\}|\) (number of nonzero coefficients).




Properties and Theory

Sign Consistency (Variable Selection)

Under the irrepresentable condition (a requirement on the correlation structure between relevant and irrelevant predictors):

\[ \|X_{\mathcal{S}^c}'X_{\mathcal{S}}(X_{\mathcal{S}}'X_{\mathcal{S}})^{-1} \text{sign}(\beta_{\mathcal{S}}^*)\|_\infty < 1 \]

the LASSO consistently selects the true support \(\mathcal{S} = \{j : \beta_j^* \neq 0\}\) as \(n \to \infty\).

Oracle Property (Adaptive LASSO)

The standard LASSO does not satisfy the oracle property in general. The adaptive LASSO (Zou, 2006) uses data-dependent weights and achieves:

  • Consistent variable selection
  • Asymptotic normality of the nonzero estimates (as if the true support were known a priori)

Prediction Error Bounds

Under restricted eigenvalue conditions, the LASSO achieves the near-minimax prediction error rate:

\[ \frac{1}{n}\|X(\hat{\beta} - \beta^*)\|_2^2 = O_P\left( s \frac{\log p}{n} \right) \]

where \(s = |\mathcal{S}|\) is the true sparsity level.

Degrees of Freedom

For the LASSO with fixed \(\lambda\) (Zou, Hastie, Tibshirani, 2007):

\[ \text{df}(\hat{\beta}^{\text{lasso}}) = E\left[ |\{j : \hat{\beta}_j \neq 0\}| \right] \]

An unbiased estimate of df is simply the number of nonzero coefficients.




Variants

Adaptive LASSO

Uses weighted penalties to achieve the oracle property:

\[ \min_{\beta} \frac{1}{2n}\|y - X\beta\|_2^2 + \lambda \sum_{j=1}^{p} w_j |\beta_j| \]

with weights \(w_j = 1/|\hat{\beta}_j^{\text{init}}|^\gamma\), where \(\hat{\beta}^{\text{init}}\) is a \(\sqrt{n}\)-consistent initial estimator (e.g., OLS or Ridge) and \(\gamma > 0\).

Group LASSO

For grouped predictors (e.g., dummy variables from a factor), groups enter or leave the model together:

\[ \min_{\beta} \frac{1}{2n}\|y - X\beta\|_2^2 + \lambda \sum_{g=1}^{G} \sqrt{p_g} \|\beta_g\|_2 \]

where \(\beta_g\) is the subvector for group \(g\) and \(p_g\) its size.

Fused LASSO

Penalizes differences between adjacent coefficients (useful for ordered/spatial data):

\[ \min_{\beta} \frac{1}{2n}\|y - X\beta\|_2^2 + \lambda_1 \sum_{j=1}^{p} |\beta_j| + \lambda_2 \sum_{j=2}^{p} |\beta_j - \beta_{j-1}| \]

LASSO for GLMs

The \(L_1\) penalty extends naturally to generalized linear models. For example, the logistic LASSO:

\[ \min_{\beta} \left\{ -\frac{1}{n} \sum_{i=1}^{n} \left[ y_i x_i'\beta - \log(1 + e^{x_i'\beta}) \right] + \lambda \|\beta\|_1 \right\} \]

Similarly for Poisson regression, Cox proportional hazards, and other exponential family models. Coordinate descent adapts via iteratively reweighted least squares (IRLS) within each coordinate update.




LASSO vs Ridge

A comparison of the two most common regularization approaches:

LASSO (\(L_1\))Ridge (\(L_2\))
Penalty\(\lambda \sum |\beta_j|\)\(\lambda \sum \beta_j^2\)
SparsityProduces exact zerosShrinks toward zero, never exactly zero
Variable selectionYes (built-in)No
Correlated predictorsTends to select one, discard othersShrinks correlated predictors together
Best whenTrue model is sparse; many irrelevant featuresMany correlated, relevant features
Unique solutionNot guaranteed when \(p > n\)Always unique (\(X'X + \lambda I\) invertible)
ComputationalCoordinate descent / LARSClosed form: \((X'X + \lambda I)^{-1}X'y\)

Elastic Net — when in doubt, combine both penalties:

\[ \min_{\beta} \frac{1}{2n}\|y - X\beta\|_2^2 + \lambda \left[ \alpha \|\beta\|_1 + \frac{1-\alpha}{2} \|\beta\|_2^2 \right] \]

with \(\alpha \in [0,1]\). The Elastic Net inherits sparsity from LASSO and the grouping effect from Ridge, making it robust in practice when the true model structure is unknown.