Ridge regression

What this chapter covers

Ridge regression is linear regression with a leash. It fits the same model — a weighted sum of the features plus an intercept — but it changes what "best" means: instead of the weights that minimize squared error alone, you take the weights that minimize squared error plus a penalty on how large those weights get. That one extra term is the whole idea, and it fixes the two places plain least squares embarrasses itself: correlated features that hand each other enormous offsetting coefficients, and a fit that clings so tightly to the training sample it can't generalize.

This chapter builds straight on week 8. If you haven't read the linear regression page, read it first — predict, r2_score, and the normal equation all come back unchanged, and ridge is a two-character edit to that normal equation. We build the penalized objective, the closed-form solution that minimizes it, and the coefficient path: the picture of every weight shrinking toward zero as the penalty grows. Then we run it against scikit-learn's Ridge and confirm they land on the same coefficients to thirteen decimal places.

The real data is the diabetes dataset built into scikit-learn: 442 patients, ten features — age, sex, BMI, blood pressure, and six blood-serum measurements — and a target that measures disease progression a year on. Several of those serum features are strongly correlated, which is exactly the mess ridge exists to clean up.

A bit of history

Ridge regression has two independent origin stories that turned out to be the same math. The name and the statistical framing come from Arthur Hoerl and Robert Kennard, two chemical engineers at DuPont, in a pair of 1970 papers in Technometrics — "Ridge Regression: Biased Estimation for Nonorthogonal Problems." Their problem was practical and industrial: when the predictor variables in a regression are nearly collinear, the least-squares estimates go haywire — huge coefficients, wild signs, numbers that swing violently if you add one more data point. Hoerl and Kennard showed you could trade a little bias for a large cut in variance by adding a small constant to the diagonal of the $X^{\top}X$ matrix before inverting it. They called the plot of coefficients against that constant the "ridge trace," which is where the method got its name.

The same idea had already shown up on the mathematics side. Andrey Tikhonov, working in the Soviet Union through the 1940s and formalized in the 1960s, was solving ill-posed inverse problems — cases where the data doesn't pin down a unique answer — by adding a regularizing term that prefers smaller, smoother solutions. That's why you'll see ridge called "Tikhonov regularization," and why the $\lambda I$ term below is sometimes written with a more general matrix. Statistics found it through unstable regressions; numerical analysis found it through unstable inverses. They're the same fix, and it's the first regularizer most people meet.

The intuition

Start with what goes wrong. Least squares wants the weights that best explain the target, and when two features carry nearly the same information it has no reason to prefer one over the other. It can put a big positive weight on the first and a big negative weight on the second, or the reverse, or any of a thousand near-equivalent combinations — they all fit the training data about equally well. The result is coefficients that are enormous, unstable, and impossible to read. Here are two of the serum features in the diabetes data, s1 and s2. Their correlation is 0.90; they move together almost perfectly.

To a least-squares fit, that near-diagonal line is a trap: it can trade weight between s1 and s2 freely, and it will, landing on whatever offsetting pair happens to squeeze out the last drop of training error. Ridge closes the trap by charging rent. Every unit of coefficient now costs something — specifically, lambda times its square — so a giant positive weight and a giant negative weight that nearly cancel are no longer free. The fit would rather split the difference: two modest weights instead of two enormous ones. That's the whole mechanism. The penalty makes the model prefer small, shared coefficients over large, brittle ones, and small coefficients are exactly what survives contact with new data.

The math

The model is unchanged from linear regression. Stack the data into a matrix $X$ with $n$ rows and $d$ columns, and each row gets a prediction as a weighted sum of its features plus an intercept:

where $w \in \mathbb{R}^{d}$ is the weight vector and $b$ is the scalar intercept. Ordinary least squares picks $w$ and $b$ to minimize the squared error. Ridge minimizes the squared error plus an L2 penalty on the weights:

The first term is the same sum of squared residuals least squares minimizes. The second, $\lambda \lVert w \rVert_2^2 = \lambda \sum_{j=1}^{d} w_j^2$, is the penalty. The scalar $\lambda \ge 0$ is the knob: at $\lambda = 0$ the penalty vanishes and you're back to ordinary least squares; as $\lambda \to \infty$ the penalty dominates and every weight is crushed toward zero. Written with the mean squared error instead of the sum, the same objective is $\text{MSE} + \lambda' \lVert w \rVert_2^2$ — dividing the error term by $n$ just rescales what $\lambda$ means, so the shape is identical. We keep the sum form because it makes $\lambda$ line up exactly with scikit-learn's alpha.

The intercept $b$ is deliberately absent from the penalty. We never punish the model for the overall height of the line — only for the steepness it assigns to each feature. Because the objective is still a convex quadratic, we can set the derivative to zero and solve for the minimum directly, exactly like the normal equation. Folding the intercept into the weights by gluing a column of ones onto $X$, the solution is:

That is the normal equation with one addition: $\lambda I$ on the diagonal (with a zero in the intercept's slot, so it stays unpenalized). That single term does the regularizing. It also does something quietly essential: $X^{\top}X$ can be singular or nearly so when features are collinear, and adding $\lambda I$ lifts every eigenvalue by $\lambda$, guaranteeing the matrix is invertible. Ridge always has an answer, even where least squares divides by zero.

One thing the math assumes: the features are on comparable scales. Since the penalty adds up squared coefficients, a feature measured in small units carries a naturally larger coefficient and gets penalized more, for no good reason. So you standardize first — center every feature to mean zero, scale to unit variance — and then a single $\lambda$ means the same thing across all of them.

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

Ridge earns its place whenever features are correlated or you have more features than you can comfortably support, which on real data is most of the time. It stabilizes coefficients that least squares leaves thrashing, it never fails to invert, and it reliably improves out-of-sample error by trading a bit of bias for a large drop in variance. It's cheap — the same one-shot linear solve as ordinary least squares, plus a diagonal — and it has exactly one hyperparameter to tune. In practice ridge is the default I reach for the moment plain linear regression's coefficients start looking unstable or its test error trails its training error.

What ridge won't do is simplify your model. It shrinks coefficients smoothly toward zero but essentially never sets one exactly to zero, so you keep all your features, just at smaller weights. If you have fifty features and you suspect forty are noise, ridge will quietly shrink all fifty rather than tell you which ten matter — you get a more stable dense model, not a sparse, readable one. That job belongs to the L1 penalty in lasso, which zeroes coefficients outright and performs feature selection, or to elastic net, which mixes L1 and L2 to get some of both. Those are the next two chapters. The rule of thumb I use: ridge when you believe most features carry a little signal and you want them all, calmer; lasso when you believe most features are useless and you want a knife.

The data

The diabetes dataset ships inside scikit-learn, so there's nothing to download. Each of the 442 rows is a patient: ten baseline measurements — age, sex, body mass index, average blood pressure, and six blood-serum measurements labeled s1 through s6 — and a target that quantifies disease progression one year after baseline, ranging from 25 to 346. It's a small, honest regression problem, and it has the property that makes ridge worth teaching on it: the serum features are entangled. We already saw s1 and s2 riding a near-perfect diagonal.

The single clearest signal is body mass index, which tracks progression about as well as any one feature can. Here it is against the target. The relationship is real and clearly linear-ish, and just as clearly noisy — one feature explains part of the story, never all of it.

Build it, one function at a time

Most of ridge is week 8. The model is the same one line — predictions are the features times the weights plus the intercept:

def predict(X, w, b):
    """The linear model: yhat = Xw + b. Same model as ordinary least squares.

    X is (n, d), w is (d,), b is a scalar; the result is (n,). Ridge changes how
    we pick w and b, not what the model computes once we have them.
    """
    return X @ w + b

And R², the number we report, is unchanged: one minus the ratio of our squared error to the error of just guessing the mean.

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. This is
    the number we report, unchanged from linear regression.
    """
    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

The first genuinely new piece is standardization, because ridge needs it. The penalty sums squared coefficients, so it only makes sense if every feature is on the same scale. We center each feature to mean zero and scale it to unit standard deviation, fitting those statistics on the training set and reusing them everywhere so no test information leaks in:

def standardize(X, mu=None, sd=None):
    """Center each feature to mean 0 and scale to unit standard deviation.

    Ridge penalizes the squared length of the weight vector, so it treats every
    coefficient on the same scale. If one feature is in dollars and another in
    kilograms, the raw magnitudes are meaningless and the penalty falls unevenly.
    Standardizing first puts every feature on equal footing so a single lambda
    means the same thing everywhere. Fit mu/sd on train, reuse them everywhere.
    """
    if mu is None:
        mu = X.mean(axis=0)
        sd = X.std(axis=0)
        sd = np.where(sd == 0.0, 1.0, sd)   # guard a constant column
    return (X - mu) / sd, mu, sd

Now the objective — the thing ridge actually minimizes. It's the squared error plus lambda times the squared length of the weight vector, and the intercept stays out of the penalty:

def ridge_objective(X, y, w, b, lam):
    """What ridge minimizes: squared error plus lambda times the squared weights.

    The first term is ordinary least squares — how badly the line misses. The
    second, lam * (w . w), is the L2 penalty: it charges the fit for the length
    of the weight vector, pulling every coefficient toward zero. The intercept b
    is NOT in the penalty — we never punish the model for the height of the line,
    only for the steepness it assigns to the features.
    """
    resid = predict(X, w, b) - y
    return float(resid @ resid + lam * (w @ w))

Minimizing that objective has a closed form, and it's the normal equation with one change. Glue on the column of ones for the intercept, add $\lambda$ to the diagonal — but a zero in the intercept's slot, so we don't penalize it — and solve. That $\lambda I$ is the entire difference between this and ordinary least squares:

def ridge_closed_form(X, y, lam):
    """The exact ridge solution: w = (XtX + lam I)^-1 Xt y, intercept unpenalized.

    Glue a column of ones onto X so the intercept rides along as the first
    weight, then build a penalty matrix that is lam on the diagonal EXCEPT the
    intercept slot, which stays 0. Solving (Xb'Xb + P) theta = Xb'y minimizes the
    ridge objective in one shot. The lam I term is why ridge is numerically kind:
    it lifts XtX away from singular, so the inverse exists even when features are
    collinear and OLS would blow up. `np.linalg.solve` factors instead of
    inverting explicitly — same answer, better conditioned.
    """
    ones = np.ones((X.shape[0], 1))
    Xb = np.hstack([ones, X])                  # (n, d+1)
    P = lam * np.eye(Xb.shape[1])
    P[0, 0] = 0.0                              # do not penalize the intercept
    theta = np.linalg.solve(Xb.T @ Xb + P, Xb.T @ y)
    b = float(theta[0])
    w = theta[1:]
    return w, b

The last piece is the one that makes ridge visual. Refit at a whole range of lambdas and record every weight each time, and you get the coefficient path: the trace of what regularization does to the model as you turn the knob up:

def coefficient_path(X, y, lambdas):
    """Refit ridge at every lambda and record the whole weight vector each time.

    This is the picture that explains ridge. At lambda near 0 you get the OLS
    weights — large, and for correlated features, wild. Crank lambda up and the
    penalty dominates: every weight is pulled smoothly toward 0, the big
    correlated coefficients collapse first, and in the limit the whole vector
    goes flat. L2 shrinks toward zero but almost never exactly to zero — that is
    the contrast with L1 (lasso), which sets coefficients to exactly 0 and
    selects features. Returns a (len(lambdas), d) array of weights.
    """
    weights = np.zeros((len(lambdas), X.shape[1]))
    for i, lam in enumerate(lambdas):
        w, _ = ridge_closed_form(X, y, lam)
        weights[i] = w
    return weights

Run that path over a log-spaced grid of lambdas and plot each coefficient as a line, and you get the classic ridge trace — the picture Hoerl and Kennard put the method's name on. At the left, near zero penalty, the weights are the OLS estimates, spread wide. Sweep right and they all bend smoothly toward zero, the big correlated ones collapsing fastest, none of them ever quite reaching it:

Watch it work

Here's the same sweep as an animation, one real refit per frame. The top panel is every coefficient as a bar; the bottom panel traces the validation error as lambda grows. Watch the left panel first: at small lambda the bars are the full OLS coefficients, and the correlated serum features stick out — s2, s4, s5 all large. As lambda climbs, every bar contracts toward the dashed zero line, and the big brittle ones shrink fastest. By the far right the whole model is nearly flat: the penalty has won, and ridge is predicting almost the same number for everyone.

The bottom panel is why you don't just crank lambda to infinity. Validation error falls at first — a little shrinkage genuinely helps — bottoms out, then climbs again as the model gets too timid to fit anything. The orange line marks the lambda at the bottom of that curve, the one we keep. The caption names the current lambda, the L2 norm of the weight vector, and the validation error at that point.

The minimum of that validation curve sits at lambda ≈ 31.6. That's the value we carry forward: enough penalty to calm the correlated coefficients, not so much that the model goes limp. It's worth noticing this is a genuine U — the left end (no penalty) and the right end (all penalty) are both worse than the middle. Ridge isn't "less is better" or "more is better"; there's a real sweet spot, and finding it is what the validation sweep is for. Left to its own devices, scikit-learn's RidgeCV, which uses leave-one-out cross-validation instead of our single held-out split, lands nearby at alpha ≈ 20 — the same neighborhood, reached a different way.

Here's that validation curve on its own, with the chosen lambda marked, so you can read the trade-off directly:

The full implementation

The whole file, no library, top to bottom — the model, R², standardization, the ridge objective, the closed-form solve, and the coefficient path. This is the code the animation above actually ran:

"""Ridge regression, built from scratch.

Ridge is ordinary least squares with one extra term: a penalty on the size of
the weights. You still fit a linear model — a weighted sum of the features plus
an intercept — but instead of chasing the smallest squared error you chase the
smallest squared error *plus* lambda times the squared length of the weight
vector. That single term is what tames multicollinearity and overfit: it makes
the fit pay for every unit of coefficient it spends, so correlated features stop
handing each other enormous offsetting weights and settle for small, stable
ones instead.

This file builds directly on week 8 (linear regression). `predict` and
`r2_score` are the same functions; the new work is the penalized objective, the
closed-form solution that minimizes it, and the coefficient path that traces
what happens to the weights as lambda grows.

Pure NumPy — no ML library anywhere in this file. Every function shows up in the
chapter one step at a time; the `# region:` markers are what the include
directives pull in.
"""

import numpy as np
import pandas as pd


# region: predict
def predict(X, w, b):
    """The linear model: yhat = Xw + b. Same model as ordinary least squares.

    X is (n, d), w is (d,), b is a scalar; the result is (n,). Ridge changes how
    we pick w and b, not what the model computes once we have them.
    """
    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. This is
    the number we report, unchanged from linear regression.
    """
    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: standardize
def standardize(X, mu=None, sd=None):
    """Center each feature to mean 0 and scale to unit standard deviation.

    Ridge penalizes the squared length of the weight vector, so it treats every
    coefficient on the same scale. If one feature is in dollars and another in
    kilograms, the raw magnitudes are meaningless and the penalty falls unevenly.
    Standardizing first puts every feature on equal footing so a single lambda
    means the same thing everywhere. Fit mu/sd on train, reuse them everywhere.
    """
    if mu is None:
        mu = X.mean(axis=0)
        sd = X.std(axis=0)
        sd = np.where(sd == 0.0, 1.0, sd)   # guard a constant column
    return (X - mu) / sd, mu, sd
# endregion


# region: ridge_objective
def ridge_objective(X, y, w, b, lam):
    """What ridge minimizes: squared error plus lambda times the squared weights.

    The first term is ordinary least squares — how badly the line misses. The
    second, lam * (w . w), is the L2 penalty: it charges the fit for the length
    of the weight vector, pulling every coefficient toward zero. The intercept b
    is NOT in the penalty — we never punish the model for the height of the line,
    only for the steepness it assigns to the features.
    """
    resid = predict(X, w, b) - y
    return float(resid @ resid + lam * (w @ w))
# endregion


# region: ridge_closed_form
def ridge_closed_form(X, y, lam):
    """The exact ridge solution: w = (XtX + lam I)^-1 Xt y, intercept unpenalized.

    Glue a column of ones onto X so the intercept rides along as the first
    weight, then build a penalty matrix that is lam on the diagonal EXCEPT the
    intercept slot, which stays 0. Solving (Xb'Xb + P) theta = Xb'y minimizes the
    ridge objective in one shot. The lam I term is why ridge is numerically kind:
    it lifts XtX away from singular, so the inverse exists even when features are
    collinear and OLS would blow up. `np.linalg.solve` factors instead of
    inverting explicitly — same answer, better conditioned.
    """
    ones = np.ones((X.shape[0], 1))
    Xb = np.hstack([ones, X])                  # (n, d+1)
    P = lam * np.eye(Xb.shape[1])
    P[0, 0] = 0.0                              # do not penalize the intercept
    theta = np.linalg.solve(Xb.T @ Xb + P, Xb.T @ y)
    b = float(theta[0])
    w = theta[1:]
    return w, b
# endregion


# region: coef_path
def coefficient_path(X, y, lambdas):
    """Refit ridge at every lambda and record the whole weight vector each time.

    This is the picture that explains ridge. At lambda near 0 you get the OLS
    weights — large, and for correlated features, wild. Crank lambda up and the
    penalty dominates: every weight is pulled smoothly toward 0, the big
    correlated coefficients collapse first, and in the limit the whole vector
    goes flat. L2 shrinks toward zero but almost never exactly to zero — that is
    the contrast with L1 (lasso), which sets coefficients to exactly 0 and
    selects features. Returns a (len(lambdas), d) array of weights.
    """
    weights = np.zeros((len(lambdas), X.shape[1]))
    for i, lam in enumerate(lambdas):
        w, _ = ridge_closed_form(X, y, lam)
        weights[i] = w
    return weights
# endregion


def load_data(path="../data/diabetes.csv"):
    """Diabetes progression: 442 patients, 10 features, a progression score.

    Returns (X, y, feature_names). y is a quantitative measure of disease
    progression one year after baseline; the features are age, sex, BMI, blood
    pressure, and six blood-serum measurements (s1-s6), several of which are
    strongly correlated with one another.
    """
    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 the regularized normal equation in production. sklearn.linear_model.Ridge is the same model with the same objective — squared error plus an L2 penalty, intercept fitted but not penalized — and its alpha is exactly our lambda, so at the same setting the two produce the same coefficients:

def sklearn_fit(X_train, y_train, alpha):
    """Fit ridge at a fixed alpha. Returns (weights, intercept) — the twin of our
    ridge_closed_form, so the coefficients line up one for one at alpha == lam."""
    model = Ridge(alpha=alpha)            # fit_intercept=True, intercept unpenalized
    model.fit(X_train, y_train)
    return model.coef_, float(model.intercept_)

Under the hood sklearn solves the same regularized system we do; for dense data it factors $X^{\top}X + \lambda I$ directly. The one convenience worth reaching for is RidgeCV, which sweeps a grid of alphas with cross-validation and hands back the best one — the automated version of the validation curve we drew by hand:

def sklearn_ridge_cv(X_train, y_train, alphas):
    """Let sklearn pick alpha by cross-validation over the grid. Returns the
    chosen alpha — sklearn's answer to 'which lambda minimizes the error?'."""
    model = RidgeCV(alphas=alphas)
    model.fit(X_train, y_train)
    return float(model.alpha_)

The thing to actually look at is what the penalty did to the coefficients. Here are the ten weights at lambda 0 — plain least squares — next to the same ten at our chosen lambda. Every ridge bar is shorter than its OLS twin, and the correlated serum features move the most: s2 shrinks from about −11 to −8, s4 from about 11 to 8. The overall L2 norm of the weight vector drops from 45.5 to 40.6. Nothing gets zeroed — that's L2 — but everything gets calmer.

Scratch versus library

Split the 442 patients 60/20/20 — 265 to train on, 88 to pick lambda with, 89 held out to score — standardize on the training set, and fit three ways: our OLS baseline (ridge at lambda ≈ 0), our ridge at the chosen lambda, and sklearn's Ridge at the same alpha. Score on the held-out 89:

Two things to read off this. First, our ridge and sklearn's are the same bar — test R² of 0.363 for both, because their coefficients agree to a maximum absolute difference of about 1.6e-13. It's the same closed form; there's nothing to disagree about. Second, and the point of the whole chapter: ridge beats the unregularized fit, 0.363 against OLS's 0.339. That's a real gain on held-out data, bought purely by refusing to let the coefficients grow. The exact figures live in results.json, regenerated whenever the code changes, so the prose and the chart can't drift from what the code produced.

An R² of 0.363 is modest — ten baseline measurements explain about a third of the variance in how this disease progresses, which is honest for a problem this hard. Here's the fit as predicted versus actual on the held-out patients; a perfect model would put every point on the dashed diagonal:

The cloud tilts along the diagonal — the model has clearly learned something — but it's a wide cloud, and it flattens at the top: the highest-progression patients get predicted too low, the model hedging toward the middle. That hedge is exactly what ridge does. Shrinking the coefficients pulls predictions toward the mean, which is the price you pay for the stability, and on this data it's a price worth paying.

Takeaways

Reach for ridge the moment plain linear regression starts misbehaving. If the coefficients are enormous, or flip sign when you add data, or the test error trails the training error, you have variance to spend, and ridge spends it well: one hyperparameter, the same closed-form solve, a guaranteed-invertible matrix, and reliably better generalization. It is the cheapest insurance in the regression toolbox, and standardizing the features plus a quick validation sweep for lambda is nearly free. On this data it turned a 0.339 into a 0.363 without touching the model class at all.

Know what ridge is and isn't. It shrinks every coefficient smoothly toward zero and keeps all of them — you end up with a stable, dense model, not a short list of the features that matter. If what you want is selection, if you'd rather the model hand you ten coefficients and set the other forty to exactly zero, ridge is the wrong penalty and lasso is the right one. Elastic net splits the difference. All three minimize squared error plus a penalty; they differ only in the shape of that penalty, and that shape is the whole personality of the method. Ridge's L2 keeps everyone at the table, smaller. Lasso's L1 clears the room. The next two chapters are those two penalties, and now you've seen the one they're measured against.