Lasso Regression and L1 Sparsity

Ordinary least squares fits $\hat\beta = \arg\min_\beta \|y - X\beta\|_2^2$, minimizing the residual sum of squares over $n$ observations and $p$ features. When $p$ is large relative to $n$, or when features are correlated, OLS estimates have high variance and are hard to interpret — nearly every coefficient is nonzero and none can be trusted individually. Regularization addresses this by adding a penalty on the size of $\beta$ to the objective, trading a bit of bias for a large reduction in variance.

Lasso regression solves

where $\lambda \geq 0$ is a tuning parameter controlling regularization strength. This is the L1-penalized analogue of ridge regression, which instead penalizes $\sum_j \beta_j^2$. The equivalent constrained form makes the geometry explicit:

for some budget $t$ that is a decreasing function of $\lambda$. The feasible region $\{\beta : \sum_j |\beta_j| \leq t\}$ is a cross-polytope (a diamond in 2D, an octahedron in 3D) with sharp vertices lying exactly on the coordinate axes. Ridge's equivalent feasible region, $\sum_j \beta_j^2 \leq t$, is a smooth ball with no corners at all. This single geometric difference — corners versus smoothness — is the entire reason Lasso produces sparse solutions and ridge does not. The intercept is typically left unpenalized (center $y$ and standardize $X$ before fitting), and features are standardized to unit variance so that the penalty treats all coefficients on a comparable scale — otherwise a feature measured in millimeters versus kilometers would be penalized asymmetrically for no principled reason.

Why sparsity emerges: the geometric argument

Consider two predictors, so $\beta = (\beta_1, \beta_2)$. The least-squares loss $\|y - X\beta\|_2^2$ has elliptical contours centered at the OLS solution $\hat\beta^{\text{OLS}}$. The regularized solution is the point where the smallest such ellipse touches the constraint region. For ridge (a disk), the boundary is smooth everywhere, so the touching point can land anywhere on the circle — generically at a point with both coordinates nonzero. For Lasso (a diamond), the boundary has flat edges meeting at sharp corners on the axes ($(\pm t, 0)$ and $(0, \pm t)$). Because the corners stick out and the ellipse's tangent direction rotates continuously as it grows, a substantial fraction of ellipse orientations first touch the diamond exactly at a corner, where one coordinate is zero. This is not a special coincidence — it happens for a positive-measure range of data configurations, which is why sparsity is the rule rather than the exception when $\lambda$ is large enough.

Why sparsity emerges: the subgradient argument

The rigorous version of this argument uses subgradients, since $|\beta_j|$ is not differentiable at zero. Writing the objective as $f(\beta) = \frac{1}{2n}\|y-X\beta\|_2^2 + \lambda\|\beta\|_1$, the optimality condition (Karush–Kuhn–Tucker for convex, non-smooth objectives) at coordinate $j$ is:

where $s_j = \text{sign}(\beta_j)$ if $\beta_j \neq 0$, and $s_j \in [-1, 1]$ (any value in that interval) if $\beta_j = 0$. Define the partial correlation $c_j = \frac{1}{n}x_j^\top(y - X\beta)$. The condition says:

• If $\beta_j \neq 0$: $c_j = \lambda \, \text{sign}(\beta_j)$, so $|c_j| = \lambda$ exactly.
• If $\beta_j = 0$: $c_j = \lambda s_j$ for some $s_j \in [-1,1]$, which is satisfiable whenever $|c_j| \leq \lambda$.

So any feature whose correlation with the current residual is weaker than $\lambda$ in absolute value gets set to exactly zero, because zero is a valid subgradient solution there — no infinitesimal push in either direction can decrease the objective, since the L1 penalty's slope at zero is $\pm\lambda$, steeper than the pull of a weak feature. Ridge's penalty derivative $2\beta_j \to 0$ as $\beta_j \to 0$, so it exerts no resistance right at zero — any nonzero correlation, however small, pulls the coefficient away from zero. That is the algebraic root of the geometric picture above.

Coordinate descent with soft-thresholding

Because the objective is convex but non-smooth, no closed-form solution exists in general (unlike ridge, which has $\hat\beta = (X^\top X + \lambda I)^{-1}X^\top y$). The standard, highly efficient algorithm is coordinate descent: fix all coefficients but one, solve that one-dimensional subproblem in closed form, and cycle. For a single coordinate $j$, holding all other $\beta_k$ ($k \neq j$) fixed, the subproblem reduces (after standardizing columns of $X$ to unit norm) to a univariate lasso problem whose exact minimizer is the soft-thresholding operator:

where $c_j = \frac{1}{n} x_j^\top r_j$ is the correlation of feature $j$ with the partial residual $r_j = y - \sum_{k \neq j} x_k \beta_k$ (the residual after removing the effect of every other feature). This single formula is the entire computational engine of Lasso. ``` ALGORITHM: Coordinate Descent for Lasso INPUT: standardized X (n x p, unit-norm columns), centered y (length n), lambda >= 0, tolerance eps, max_iter OUTPUT: coefficient vector beta (length p) beta <- zeros(p) r <- y # r = current residual, y - X*beta, beta=0 initially repeat until max coefficient change < eps or iter == max_iter: for j = 1 .. p:

# remove feature j's current contribution from the residual r <- r + x_j * beta[j]

# compute correlation of residual with feature j c_j <- (1/n) * dot(x_j, r)

# soft-threshold update (closed-form 1-D minimizer) beta_old <- beta[j] beta[j] <- sign(c_j) * max(|c_j| - lambda, 0)

# re-insert feature j's (new) contribution r <- r - x_j beta[j] iter <- iter + 1 return beta ``` Each pass over all $p$ coordinates is one "epoch." In practice the algorithm is run over a path* of decreasing $\lambda$ values (warm-starting each fit from the solution at the previous, larger $\lambda$), which both accelerates convergence and produces the full regularization path used for model selection via cross-validation.

Time complexity. Each coordinate update in the algorithm above costs $O(n)$ (a dot product of length-$n$ vectors for $c_j$, plus $O(n)$ residual bookkeeping). One full epoch over all $p$ coordinates costs $O(np)$. If the algorithm runs for $T$ epochs to converge, total cost is $O(Tnp)$. In practice $T$ is small (tens to low hundreds) because coordinate descent on this objective enjoys fast, often linear, empirical convergence, and warm-starting along a $\lambda$-path reduces $T$ further at each step. This is the same order of complexity as a single ridge solve via naive matrix inversion, $O(p^3 + np^2)$, but coordinate descent is typically far cheaper when $p$ is large and sparsity is present, since coordinates that are zero can be skipped cheaply (active-set / working-set strategies restrict updates to a shrinking candidate set).

Space complexity. $O(np)$ to store $X$, plus $O(n + p)$ for the residual and coefficient vectors — no $p \times p$ matrix needs to be formed or inverted, which is a major practical advantage over the ridge closed-form solution when $p$ is large.

Correctness argument. The Lasso objective

is convex (a sum of a convex quadratic and a convex norm), though not differentiable. Coordinate descent is guaranteed to converge to the global minimum for this specific objective because: 1. Each coordinate subproblem is solved exactly. Fixing $\beta_{-j}$, the restriction of $f$ to $\beta_j$ is a strictly convex function of one variable (quadratic plus $|\beta_j|$ term), whose unique minimizer is exactly the soft-threshold formula $S(c_j,\lambda)$ derived from setting the subgradient to zero. There is no approximation error at each step — every coordinate update strictly decreases $f$ (or leaves it unchanged if already optimal). 2. The non-differentiable part is separable across coordinates ($\|\beta\|_1 = \sum_j |\beta_j|$), which is precisely the condition (Tseng's theorem for block coordinate descent on non-smooth convex functions) under which cyclic coordinate descent is guaranteed to converge to a global minimum, not merely a stationary point or local optimum. This separability is essential — coordinate descent is not generally guaranteed to work on non-smooth objectives whose non-smooth part couples coordinates together. 3. Since $f$ is convex and the sequence of iterates has monotonically non-increasing objective values bounded below (by zero), and each coordinate is optimized exactly, the iterates converge to a point satisfying the global KKT/subgradient stationarity condition for every coordinate simultaneously — which for a convex function is sufficient for global optimality.

Uniqueness caveat. The optimal fitted values $X\hat\beta$ are always unique (strict convexity of the squared-error term in $X\beta$ ensures this), but the coefficient vector $\hat\beta$ itself need not be unique when columns of $X$ are exactly collinear — multiple sparse solutions can achieve the same loss. This is a real and important caveat, discussed further below.

Honest tradeoffs. Lasso trades a small increase in bias (shrinking true nonzero coefficients toward zero, sometimes past their true magnitude) for a large reduction in variance and, uniquely among common penalties, for interpretability via sparsity. It is not universally better than ridge: when all features genuinely contribute a little (a "dense" true signal), ridge's smooth shrinkage often achieves lower prediction error, because Lasso's hard zeroing throws away weak-but-real signal. Lasso also cannot select more than $n$ features when $p > n$ (a known limitation), and it handles groups of correlated features poorly, arbitrarily picking one and zeroing the rest even when they are equally informative.

Consider $n = 4$ observations and $p = 3$ standardized features (mean zero, unit variance columns), predicting a centered target $y$.

Note $x_1$ and $x_2$ are highly correlated (both track the same up/down pattern), while $x_3$ is comparatively weakly related to $y$. We trace one coordinate descent epoch starting from $\beta = (0,0,0)$, using $\lambda = 0.15$.

Update $\beta_1$. With $\beta_2 = \beta_3 = 0$, the partial residual is just $r = y$. Compute $c_1 = \frac{1}{4}\sum_i x_{1i} y_i$:

Soft-threshold: $\beta_1 \leftarrow \text{sign}(1.325)\cdot\max(1.325 - 0.15, 0) = 1.175$. Residual is updated by subtracting $x_1 \cdot 1.175$ from $y$.

Update $\beta_2$. Recompute the partial residual $r = y - x_1\beta_1$ and its correlation with $x_2$; because $x_1$ and $x_2$ are nearly collinear, $x_1\beta_1$ has already absorbed most of the signal that $x_2$ would otherwise explain. The resulting $c_2$ works out to a small value — call it $c_2 \approx 0.08$ — smaller in magnitude than $\lambda = 0.15$. Soft-threshold: $\beta_2 \leftarrow \text{sign}(0.08)\cdot\max(0.08 - 0.15,\,0) = \max(-0.07,\,0) = 0$. $\beta_2$ is set exactly to zero — not because $x_2$ is uninformative about $y$ on its own (it is almost as predictive as $x_1$), but because once $x_1$'s effect is already accounted for, the remaining signal for $x_2$ is too weak to overcome the $\lambda$ threshold. This is the "pick one out of a correlated group" behavior mentioned above, made concrete.

Update $\beta_3$. Suppose the resulting correlation $c_3 \approx 0.31$, above threshold. Then $\beta_3 \leftarrow \text{sign}(0.31)\cdot(0.31 - 0.15) = 0.16$.

End of epoch 1: $\beta \approx (1.175,\ 0,\ 0.16)$ — already sparse after a single pass, with $x_2$ eliminated despite being marginally almost as useful as $x_1$. Subsequent epochs refine these values slightly as the residual keeps being redistributed, but the support (which coefficients are nonzero) typically stabilizes quickly. Contrast with ridge on the same data: ridge would assign both $\beta_1$ and $\beta_2$ moderate, similarly-sized nonzero weights, splitting credit between the two correlated features rather than eliminating one.

Lasso and its L1-penalty machinery are load-bearing components of modern statistical practice wherever feature counts are large or interpretability is required. In genomics and bioinformatics, Lasso-based methods identify small panels of genetic markers associated with disease from datasets where the number of genes vastly exceeds the number of patients. In finance, sparse regression selects a manageable set of economic indicators from thousands of correlated candidates for risk models, where regulators demand explainable, auditable coefficients rather than opaque dense models. In natural language processing and text classification, L1-regularized linear models (and L1-regularized logistic regression, its classification analogue) produce compact, interpretable models over enormous bag-of-words or n-gram feature spaces. The same soft-thresholding mechanism underlies compressed sensing, a signal-processing technique for reconstructing sparse signals from far fewer measurements than the signal's dimension, which relies on exactly the same L1-minimization principle to recover sparse solutions from underdetermined linear systems. Modern deep learning also borrows the idea: L1 regularization is applied to neural network weights to induce sparse, prunable architectures, and the proximal gradient methods used to fit Lasso (ISTA and its accelerated variant FISTA) are foundational building blocks in the broader toolkit of Gradient Descent for Parameter Estimation-based optimization for non-smooth objectives across machine learning.

Ask the tutor: How does the LARS algorithm compute the entire Lasso regularization path in one efficient pass, and why is the path piecewise linear in λ?