Lasso regression

What this chapter covers

Lasso is linear regression that fires its own useless features. You take the ordinary least-squares loss and add one term — a penalty on the sum of the absolute values of the weights — and that small change makes the model set some coefficients to exactly zero instead of merely small. The zeros are the point. A lasso fit doesn't just predict; it hands you a shorter list of features and says the rest didn't earn their place.

That's a different job from plain regression, and it's a different job from ridge, the L2 cousin we'll keep comparing against. Ridge shrinks every coefficient smoothly toward zero and stops just short — nothing ever actually lands on zero. Lasso, because of the shape of the L1 penalty, snaps coefficients onto zero one at a time as you turn up the pressure. Same data, same linear model, completely different behavior at the corner. We build it from scratch with coordinate descent and the soft-thresholding operator, watch the coefficients fall to zero on real data, and check the whole thing against scikit-learn's Lasso.

The data is the diabetes dataset: 442 patients, ten baseline measurements each, and a number that says how much the disease progressed a year later. Ten features is small enough to watch every coefficient by name, which is exactly what we want when the story is about which ones survive.

A bit of history

Lasso is young as these things go. Robert Tibshirani introduced it in 1996, in a paper titled "Regression Shrinkage and Selection via the Lasso," and the name is an acronym he coined for what it does: Least Absolute Shrinkage and Selection Operator. Two verbs in one method — it shrinks the coefficients, and it selects among them — and that combination was the new thing. Before lasso, feature selection meant stepwise procedures: add a variable, drop a variable, refit, repeat, on a discrete search that was unstable and hard to reason about. Lasso folded selection into the fit itself. You solve one convex optimization and the irrelevant coefficients come out at zero, for free.

The idea didn't drop out of nowhere. Leo Breiman's non-negative garrote (1995) was circling the same target a year earlier, and the L1 penalty had shown up in signal processing as basis pursuit around the same time. But Tibshirani's framing — an L1 constraint on ordinary regression, motivated as a direct competitor to ridge — is the one that stuck and set off twenty years of follow-on work: least angle regression for computing the whole path, the elastic net for when features travel in correlated groups, the graphical lasso, and on. If you do applied ML today and you've ever wanted a sparse model, this 1996 paper is where the road starts.

The intuition

Every coefficient in a linear model costs you something under lasso. Not a cost that scales with how well the feature predicts — a flat toll, the same per unit of weight whether the weight is already large or nearly gone. That flatness is the whole trick. Ridge charges you the square of the weight, so as a coefficient shrinks the penalty's pull on it fades to nothing, and the coefficient settles at some small nonzero value where the shrinking force and the data's pull balance. Lasso's pull never fades. Right up until a weight hits zero the penalty is still tugging with full strength, so for a feature that isn't pulling its weight in the fit, zero wins outright.

The picture people draw is the constraint region. Bounding the sum of absolute weights makes a diamond (a rotated square in two dimensions, a spiky polytope in more); bounding the sum of squares makes a circle. The best fit is where the loss contours first touch that region, and a diamond has corners that stick out along the axes — the places where some coordinates are exactly zero. Contours tend to hit corners. A circle has no corners, so it gets touched on its side, at a point where every coordinate is small but nonzero. Geometry is the whole difference between "small" and "zero."

That's the promise; here's the payoff on a toy problem where we know the truth. I generated 120 points from twelve features, but only four of them actually drive the target — the other eight have true coefficient zero. Ordinary least squares can't know that, so it fits all twelve, handing every noise feature a small spurious weight. Lasso, at a sensible penalty, sets eight of them to exactly zero and keeps the four that matter. Blue is the truth, and you can see lasso landing on it while OLS smears weight across the noise.

The math

Start from the least-squares loss and stack the data into $X$ with $n$ rows and $d$ columns. Lasso minimizes that loss plus an L1 penalty on the weights:

Here $w \in \mathbb{R}^{d}$ is the weight vector, $b$ is the intercept (never penalized — you don't want to shrink the baseline), and $\lambda \ge 0$ is the knob that sets how hard the penalty pushes. At $\lambda = 0$ this is ordinary least squares. Turn $\lambda$ up and coefficients get pulled toward zero; turn it up far enough and they all reach it. The only change from ridge is the penalty term: ridge uses $\lambda \sum_j w_j^2$, the squared L2 norm, and that one swap is the entire behavioral difference between the two methods.

The absolute value is what makes this interesting and what makes it awkward. It's not differentiable at zero, so there's no normal equation, no one-shot matrix solve. What saves us is that if you freeze every weight but one, the problem in that single weight has a clean closed-form answer. Fix all $w_k$ for $k \ne j$, let $r_j = y - \sum_{k \ne j} w_k x_k - b$ be the residual with feature $j$ left out, and define its correlation with that residual as

Then the weight that minimizes the loss along coordinate $j$ (with features standardized so each column has $\frac{1}{n} x_j^{\top} x_j = 1$) is the soft-thresholding operator applied to $\rho_j$:

Read $S_\lambda$ literally: pull $z$ toward zero by $\lambda$, and if that push would carry it past zero, stop it at zero. Anything with correlation smaller than $\lambda$ in magnitude gets clipped to nothing — that is the snap. Coordinate descent is just this update, one feature at a time, cycled until nothing moves. Because the loss is convex, the order doesn't matter and it converges to the one global minimum.

What it's good at, what it isn't

Lasso earns its keep when you suspect most of your features are dead weight and you'd rather the model tell you which than guess yourself. You get a sparse, readable model out of a single convex fit — no stepwise search, no p-value ritual — and sparsity is worth real money in production. Fewer nonzero coefficients means fewer inputs to collect and serve at prediction time, a model you can actually explain to whoever has to sign off on it, and, when the true signal really is sparse, better generalization than throwing every feature at least squares. It's the same L1 idea that regularizes plenty of larger models; learning it here on ten features is learning it everywhere.

Where it gets unreliable is correlated features, and this is the honest knock on it. When two features carry nearly the same information, lasso tends to pick one almost arbitrarily and zero the other, and which one it picks can flip with a small change in the data or the random seed. The selection is unstable even when the prediction is fine, so don't read "this feature got zeroed" as "this feature doesn't matter" — read it as "one of this correlated group got kept." It also caps out at $n$ selected features when you have more features than samples, and the shrinkage it applies to the survivors biases them low. Those limits are exactly what the elastic net in the next chapter exists to fix, by mixing in a bit of ridge to share credit across correlated groups instead of forcing a winner-take-all.

The data

The diabetes dataset ships with scikit-learn, so there's nothing to download. Each of the 442 rows is a patient with ten baseline features — age, sex, body mass index, average blood pressure, and six blood-serum measurements labeled s1 through s6 — already centered and scaled by the people who packaged it. The target is a quantitative measure of disease progression one year after that baseline, running from 25 to 346. It's a genuine regression problem with modest signal, which makes it a good lasso testbed: enough features that some are plausibly redundant, few enough that we can name every one the model keeps or drops.

Body mass index is the single feature you'd bet on up front, and the data agrees — heavier patients trend toward faster progression. Here's BMI against the target across all 442 patients. The trend is real and loose, the kind of signal a linear model can grab a piece of but not all of.

Build it, one function at a time

Pure NumPy, in the order you'd write it: the model, the score, the objective, the operator that does the work, and the loop that calls it.

Prediction is unchanged from ordinary regression — lasso changes how we choose the weights, not what the model computes with them:

def predict(X, w, b):
    """The linear model itself: yhat = Xw + b.

    Identical to ordinary linear regression — lasso changes how we *choose* w,
    not what the model computes. X is (n, d), w is (d,), b is a scalar, and the
    result is (n,), one predicted number per row.
    """
    return X @ w + b

R² is the number we'll actually report, the fraction of variance the model explains against the baseline of always guessing the mean. The penalty never enters here; R² judges predictions, not the loss we optimized:

def r2_score(X, y, w, b):
    """Coefficient of determination: the fraction of variance we explain.

    1 minus (our squared error / the squared error of always guessing the
    mean). 1.0 is perfect, 0.0 is no better than the mean, negative is worse.
    The penalty never enters here — R^2 judges predictions, not the objective.
    """
    resid = y - predict(X, w, b)
    ss_res = float(np.sum(resid ** 2))
    ss_tot = float(np.sum((y - y.mean()) ** 2))
    return 1.0 - ss_res / ss_tot

Now the objective, the thing coordinate descent is driving down. Mean squared error scaled by a half, plus $\lambda$ times the sum of absolute weights — the data-fit term and the toll, written out:

def lasso_objective(X, y, w, b, lam):
    """The thing we minimize: mean squared error plus the L1 penalty.

    The first term is the usual data-fit, scaled by 1/2 so its gradient comes
    out clean. The second, lambda * sum |w_j|, is the L1 penalty — it charges a
    flat toll per unit of weight, the same toll no matter how big the weight
    already is. That constant slope is what pins small coefficients to exactly
    zero: unlike the squared L2 penalty, the pull toward zero doesn't fade as
    the weight shrinks. lam = 0 gives plain least squares back.
    """
    resid = predict(X, w, b) - y
    mse = float(np.mean(resid ** 2)) / 2.0
    penalty = lam * float(np.sum(np.abs(w)))
    return mse + penalty

Here's the heart of it. The soft-thresholding operator is three lines and it's the only place a coefficient ever becomes exactly zero. Shrink toward zero by $\lambda$, clip at zero, keep the sign:

def soft_threshold(z, lam):
    """The soft-thresholding operator S_lambda(z) — lasso's beating heart.

    It shrinks z toward zero by lam and clips anything smaller than lam to
    exactly zero: S_lambda(z) = sign(z) * max(|z| - lam, 0). This is the
    closed-form minimizer of (1/2)(w - z)^2 + lam*|w| over a single weight w,
    and it's the one place the "snap to zero" happens. Every coordinate update
    below is just this function applied to that coordinate's target value.
    """
    return np.sign(z) * np.maximum(np.abs(z) - lam, 0.0)

And the fit: cycle through the coordinates, and for each one strip its current contribution out of the residual, measure how much that feature still correlates with what's left, and set its weight with the operator above. A running $Xw$ keeps each update cheap. When a full sweep moves nothing, we're at the minimum:

def coordinate_descent(X, y, lam, n_iters=1000, tol=1e-9):
    """Fit lasso by cycling through one weight at a time.

    Assumes X is standardized (each column centered, unit variance) so the
    per-column scale is 1 and the intercept is simply the mean of y — a lasso
    penalty is never applied to the intercept, so we fix it and fit the slopes
    on the centered target. For each coordinate j we strip j's current
    contribution out of the residual, compute rho_j = (1/n) x_j . r_j (how much
    feature j still wants to explain), and set w_j = S_lambda(rho_j) / z_j. We
    keep a running X @ w so each update is O(n), not O(nd), and stop when a full
    sweep moves no weight by more than tol.
    """
    n, d = X.shape
    b = float(y.mean())            # intercept: unpenalized, exact for centered X
    yc = y - b                     # center the target
    z = (X ** 2).sum(axis=0) / n   # per-column scale; == 1 for standardized X
    w = np.zeros(d)
    Xw = np.zeros(n)               # running prediction X @ w, kept in sync
    for _ in range(n_iters):
        max_change = 0.0
        for j in range(d):
            r_j = yc - Xw + w[j] * X[:, j]        # residual without feature j
            rho = float(X[:, j] @ r_j) / n
            w_new = soft_threshold(rho, lam) / z[j]
            if w_new != w[j]:
                Xw += (w_new - w[j]) * X[:, j]    # keep the running fit current
                max_change = max(max_change, abs(w_new - w[j]))
                w[j] = w_new
        if max_change < tol:
            break
    return w, b

One detail that matters more than it looks: the features have to be standardized first. Lasso penalizes every coefficient with the same $\lambda$, so if one feature is measured in thousands and another in fractions, the same penalty means wildly different things to each. Put them on a common scale and the penalty is fair. We fit that scaling on the training data only, so nothing leaks from the test set.

Watch it work

This is the lasso path, and it's the animation to sit with. I fit the model fifty times on the training data, sweeping $\lambda$ from small on the left of the grid to large on the right — from 0.049 up to 49.1, the value where the last coefficient dies. Every frame is a real, fully converged refit, not an interpolation.

The top panel is the ten coefficients, one bar each, by name. The bottom panel counts how many are still nonzero. Watch what happens as $\lambda$ climbs: the bars shrink together, and then one by one a coefficient hits zero and goes gray — that's the snap, the soft-threshold clipping it out of the model. The count in the bottom panel steps down each time another feature drops. Press play, then scrub back and forth over the moment a bar crosses to gray.

Two things to notice. First, the order features leave in is information — the ones that hang on longest (body mass index, s5) are the ones the model leans on hardest, and the ones that drop early were barely contributing. Second, this isn't smooth shrinkage toward zero and then a soft landing; it's shrink, shrink, gone. A coefficient is nonzero and then at some exact $\lambda$ it's zero and stays there. Ridge would give you the shrinking without the gone — every bar would keep getting shorter forever and none would ever turn gray. The gray bars are the entire reason to reach for lasso.

The full implementation

The whole file, no library, top to bottom — predict, R², the objective, the soft-threshold, and coordinate descent, plus the standardize and load helpers. This is the code the animation actually ran:

"""Lasso regression, built from scratch.

Lasso is ordinary least squares with one extra term bolted onto the loss: an
L1 penalty on the weights, lambda * sum |w_j|. That single change is the whole
story. Because the penalty has a corner at zero, the minimizer parks
coefficients at *exactly* zero once lambda gets big enough — the model selects
its own features by switching the useless ones off. Ridge (L2) shrinks weights
smoothly toward zero but never reaches it; lasso snaps them to zero one by one.

You can't get there with the normal equation — the absolute value isn't
differentiable at zero, so there's no clean formula for the minimum. The
workhorse is coordinate descent: hold every weight fixed but one, minimize the
loss over that one weight (which has a closed form, the soft-threshold), move to
the next, and cycle until nothing moves. Pure NumPy — no ML library in this
file. The `# region:` markers are what the chapter's include directives pull in.
"""

import numpy as np
import pandas as pd


# region: predict
def predict(X, w, b):
    """The linear model itself: yhat = Xw + b.

    Identical to ordinary linear regression — lasso changes how we *choose* w,
    not what the model computes. X is (n, d), w is (d,), b is a scalar, and the
    result is (n,), one predicted number per row.
    """
    return X @ w + b
# endregion


# region: r2
def r2_score(X, y, w, b):
    """Coefficient of determination: the fraction of variance we explain.

    1 minus (our squared error / the squared error of always guessing the
    mean). 1.0 is perfect, 0.0 is no better than the mean, negative is worse.
    The penalty never enters here — R^2 judges predictions, not the objective.
    """
    resid = y - predict(X, w, b)
    ss_res = float(np.sum(resid ** 2))
    ss_tot = float(np.sum((y - y.mean()) ** 2))
    return 1.0 - ss_res / ss_tot
# endregion


# region: objective
def lasso_objective(X, y, w, b, lam):
    """The thing we minimize: mean squared error plus the L1 penalty.

    The first term is the usual data-fit, scaled by 1/2 so its gradient comes
    out clean. The second, lambda * sum |w_j|, is the L1 penalty — it charges a
    flat toll per unit of weight, the same toll no matter how big the weight
    already is. That constant slope is what pins small coefficients to exactly
    zero: unlike the squared L2 penalty, the pull toward zero doesn't fade as
    the weight shrinks. lam = 0 gives plain least squares back.
    """
    resid = predict(X, w, b) - y
    mse = float(np.mean(resid ** 2)) / 2.0
    penalty = lam * float(np.sum(np.abs(w)))
    return mse + penalty
# endregion


# region: soft_threshold
def soft_threshold(z, lam):
    """The soft-thresholding operator S_lambda(z) — lasso's beating heart.

    It shrinks z toward zero by lam and clips anything smaller than lam to
    exactly zero: S_lambda(z) = sign(z) * max(|z| - lam, 0). This is the
    closed-form minimizer of (1/2)(w - z)^2 + lam*|w| over a single weight w,
    and it's the one place the "snap to zero" happens. Every coordinate update
    below is just this function applied to that coordinate's target value.
    """
    return np.sign(z) * np.maximum(np.abs(z) - lam, 0.0)
# endregion


# region: coordinate_descent
def coordinate_descent(X, y, lam, n_iters=1000, tol=1e-9):
    """Fit lasso by cycling through one weight at a time.

    Assumes X is standardized (each column centered, unit variance) so the
    per-column scale is 1 and the intercept is simply the mean of y — a lasso
    penalty is never applied to the intercept, so we fix it and fit the slopes
    on the centered target. For each coordinate j we strip j's current
    contribution out of the residual, compute rho_j = (1/n) x_j . r_j (how much
    feature j still wants to explain), and set w_j = S_lambda(rho_j) / z_j. We
    keep a running X @ w so each update is O(n), not O(nd), and stop when a full
    sweep moves no weight by more than tol.
    """
    n, d = X.shape
    b = float(y.mean())            # intercept: unpenalized, exact for centered X
    yc = y - b                     # center the target
    z = (X ** 2).sum(axis=0) / n   # per-column scale; == 1 for standardized X
    w = np.zeros(d)
    Xw = np.zeros(n)               # running prediction X @ w, kept in sync
    for _ in range(n_iters):
        max_change = 0.0
        for j in range(d):
            r_j = yc - Xw + w[j] * X[:, j]        # residual without feature j
            rho = float(X[:, j] @ r_j) / n
            w_new = soft_threshold(rho, lam) / z[j]
            if w_new != w[j]:
                Xw += (w_new - w[j]) * X[:, j]    # keep the running fit current
                max_change = max(max_change, abs(w_new - w[j]))
                w[j] = w_new
        if max_change < tol:
            break
    return w, b
# endregion


def standardize(X_train, X_other=None):
    """Center and scale to unit variance using the TRAIN statistics only.

    Lasso penalizes every coefficient by the same lam, so the features have to
    share a scale or the penalty means something different for each one. We fit
    mu and sigma on train and apply them everywhere, so no test information
    leaks into the transform.
    """
    mu = X_train.mean(axis=0)
    sd = X_train.std(axis=0)
    sd[sd == 0] = 1.0
    Xs = (X_train - mu) / sd
    if X_other is None:
        return Xs, mu, sd
    return Xs, (X_other - mu) / sd, mu, sd


def load_data(path="../data/diabetes.csv"):
    """Diabetes: 442 patients, 10 baseline features, a one-year disease score.

    Returns (X, y, feature_names). The target is a quantitative measure of
    disease progression one year after baseline; the features are age, sex, BMI,
    average blood pressure, and six blood-serum measurements.
    """
    df = pd.read_csv(path)
    target = "target"
    features = [c for c in df.columns if c != target]
    X = df[features].to_numpy(float)
    y = df[target].to_numpy(float)
    return X, y, features

The library version

Nobody hand-rolls coordinate descent in production, and once you've seen it you don't need to. sklearn.linear_model.Lasso minimizes the identical objective — its alpha is our $\lambda$ — and runs the same coordinate-descent-with- soft-thresholding underneath, in compiled Cython. We standardize outside the model and let it fit the intercept, matching our setup exactly:

def sklearn_fit(X_train, y_train, alpha):
    """Fit lasso with sklearn at a given alpha. Returns (weights, intercept) —
    the twin of our coordinate_descent, so the coefficients line up one for one
    and the same features come out at exactly zero."""
    model = Lasso(alpha=alpha, fit_intercept=True, max_iter=100000, tol=1e-12)
    model.fit(X_train, y_train)
    return model.coef_, float(model.intercept_)

The one thing to watch when you first use it is that alpha and our $\lambda$ line up only because both objectives carry the same $\frac{1}{2n}$ in front of the squared error. Libraries differ on that scaling — some drop the half, some divide by $n$, some don't — and if you ever port a penalty between two of them and the sparsity looks wrong, that constant is the first place to check. For scoring we fit on train and report test R² and MSE:

def sklearn_score(X_train, y_train, X_test, y_test, alpha):
    """Fit on train, report test R^2 and test MSE — the face-off numbers."""
    model = Lasso(alpha=alpha, fit_intercept=True, max_iter=100000, tol=1e-12)
    model.fit(X_train, y_train)
    pred = model.predict(X_test)
    r2 = float(model.score(X_test, y_test))
    mse = float(((pred - y_test) ** 2).mean())
    return r2, mse

At the penalty we settle on below, the two fits are indistinguishable — same coefficients, same features zeroed, agreeing to eight decimal places. That's the point of building it: there's one lasso solution for a given $\lambda$, and a correct coordinate descent finds the same one the library does.

Scratch versus library

Which $\lambda$? That's the real question with lasso, and you answer it with data you didn't train on. I split the 442 patients 60/20/20 — 265 to train, 88 to validate, 89 held out for the final test — fit the whole path on the training set, and scored each $\lambda$ on the validation set. The curve below is validation MSE against $\log_{10}\lambda$; the marked point is the minimum, at $\lambda \approx 2.21$.

The curve dips and then climbs, which is the shape you want to see. Too little penalty on the left and the model overfits the training noise; too much on the right and it's zeroed out features it needed and underfits. The bottom of the valley is the trade you want. At that $\lambda$, lasso keeps seven of the ten features and drops three — age, s1, and s4 — setting their coefficients to exactly zero. Body mass index comes out strongest at 29.7, then s5 at 21.3, with blood pressure and s3 next; the dropped three are gone, not small.

Now fit that chosen $\lambda$ three ways on the training data and score on the held-out test set: ordinary least squares with all ten features as a baseline, our lasso, and sklearn's.

Two bars are identical — our lasso and sklearn's both land at test R² 0.362, test MSE 3169, the same coefficients to eight decimals. And here's the part worth underlining: the lasso, using seven features, edges out full least squares using all ten (R² 0.339). It threw away three features and predicted the held-out patients slightly better for it. That's regularization doing its job — the three it dropped were contributing more noise than signal, and cutting them tightened the fit on data the model had never seen. The exact figures live in results.json, regenerated whenever the code changes, so the prose can't drift from what ran.

Sanity-check the survivors with a fitted-versus-actual plot on the test set. A perfect model sits on the diagonal; the spread around it is the disease progression these ten baseline measurements can't account for.

The cloud tracks the diagonal but loosely, which is the honest picture of a dataset where ten measurements explain about a third of the variance. Lasso didn't fix that ceiling — no linear model will — but it got there with three fewer inputs than the full fit needed, which is the trade it exists to make.

Takeaways

Reach for lasso when you have more features than you trust and you want the model to make the cut. It does two jobs in one convex fit — shrink and select — and the select half is the one you can't get from ridge or from plain regression. On the diabetes data it dropped three of ten features and predicted slightly better for it; the win there isn't the fraction of a point of R², it's shipping a model with seven inputs instead of ten and a clean answer to "which features actually matter." That answer is worth more in most rooms than the accuracy.

But trust the sparsity more than the specific selection. The instability under correlated features is real, and I've watched a lasso swap which of two twinned features it keeps between two runs on almost the same data. When your features travel in correlated groups — and real features usually do — lasso's habit of crowning one and zeroing the rest is arbitrary in a way that'll bite you if you over-read it. That's precisely the gap the elastic net closes in the next chapter, by blending the L1 penalty that gives you zeros with the L2 penalty that shares credit across a group. Ridge shrinks, lasso selects, elastic net does both; the soft-threshold you built here is the piece the third one is built on top of. Learn it on ten named features you can watch snap to zero, because the same operator is doing the same job inside far bigger models where you can't.