Linear regression

What this chapter covers

Linear regression is the model you reach for first when the thing you're predicting is a number instead of a label. It fits a straight line — a plane, really, once you have more than one feature — through your data and reads predictions off it. That sounds too simple to matter, and it's exactly why it matters: it's the workhorse baseline for every regression problem, the model whose coefficients you can actually read and argue about, and the thing a fancier model has to beat before it's worth the trouble.

We build it twice. First the closed-form solution — the normal equation, one linear solve that hands you the best-fit line in a single shot. Then gradient descent, which crawls to the same answer one downhill step at a time. The closed form is what you'd ship; the gradient version is what you'd watch, because it's the same optimization loop that trains everything later in the course, and here it's simple enough to see. We run both against scikit-learn and confirm all three land on the same line.

The real data is California housing: 20,640 census block groups, eight measurements each — median income, house age, location, a few more — and a target of median house value. We'll also keep a tiny synthetic set on the side, forty points in one dimension, so you can actually watch the line move.

A bit of history

Least squares is one of the oldest ideas in applied mathematics, and it comes with a genuine priority fight. Adrien-Marie Legendre published the method in 1805, in an appendix on the orbits of comets, and gave it the name we still use — moindres carrés, least squares. Carl Friedrich Gauss published his own version in 1809 and claimed he'd been using it since 1795, to predict the position of the dwarf planet Ceres after it disappeared behind the sun. The astronomers had lost Ceres; Gauss fit a curve to the handful of observations they had, told them where to point their telescopes, and it was there. That's the founding story of the whole field, and it's a regression.

Gauss went further than Legendre — he tied least squares to the normal distribution and showed it was the optimal estimator under Gaussian noise, the result we now call the Gauss-Markov theorem. Francis Galton gave it the name "regression" almost a century later, in the 1880s, studying how the heights of children regressed toward the mean of the population relative to their parents. The word stuck even though it describes a quirk of one dataset rather than the method. Two hundred years on, the normal equation Gauss and Legendre argued over is the same one in the code below, unchanged.

The intuition

You have points, and you want a line through them. Not any line — the one that sits as close to all the points as it can, so that when a new x shows up you can read off a sensible y. "As close as it can" is the whole game, and least squares makes it precise: measure closeness by the vertical gap between each point and the line, square those gaps so they can't cancel out, and pick the line that makes the total smallest.

Squaring is the choice that defines the method, and it's worth pausing on. It means a point twice as far off the line contributes four times the penalty, so the fit bends hard to avoid large misses and shrugs at small ones. That's what makes the line settle onto the bulk of the data — and it's also why a single wild outlier can drag the whole line toward it. Every property of linear regression, good and bad, traces back to that squared gap.

Here's the idea on forty synthetic points. The line is the least-squares fit: the one line, out of all possible lines, where the squared vertical gaps add up to the smallest total.

The math

Stack your data into a matrix $X$ with $n$ rows and $d$ columns — one row per example, one column per feature. The model gives each row a prediction as a weighted sum of its features plus an intercept:

Here $w \in \mathbb{R}^{d}$ is the weight vector — one coefficient per feature — and $b$ is the scalar intercept, the prediction when every feature is zero. Fitting means choosing $w$ and $b$ to minimize the mean squared error over the $n$ examples:

This surface is a bowl — convex, one minimum — so "the best line" is well-defined and there's no local minimum to get stuck in. Set the derivatives to zero and you can solve for the minimum directly. Fold the intercept into the weights by gluing a column of ones onto $X$, and the solution is the normal equation:

That's the closed form: one matrix inverse, one line of code, exact. When $X^{\top}X$ is expensive or badly conditioned you don't want to form the inverse — you solve the linear system instead, or descend the bowl by gradient. The gradient of the loss is where the residual $r = Xw + b - y$ falls out front:

Gradient descent just walks against it, scaled by a learning rate $\eta$, until the steps stop moving anything:

Same destination as the normal equation, reached on foot. We'll do both and check they agree.

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

The case for linear regression is interpretability and speed. Every coefficient is a sentence: hold everything else fixed, move this feature by one unit, and the prediction moves by exactly this much. That's a claim a domain expert can confirm or reject, which is worth more than a fraction of a point of accuracy in most rooms I've worked in. It fits in closed form, so training is a single linear solve — no tuning, no epochs, no random seed. And with the model this constrained, there's very little to overfit; it's the high-bias, low-variance end of the dial, the safe default.

The case against it is that the world is often not linear, and linear regression can only draw one flat surface. If the true relationship curves, or two features only matter in combination, a plane can't capture it and no amount of data will fix that — you'd have to add the curve or the interaction as a feature yourself. It also leans on assumptions that quietly fail: it expects the noise to have roughly constant spread across the range (homoscedasticity) and the errors to be independent, and when those break the coefficients still come out but the confidence you'd put on them shouldn't. Correlated features (multicollinearity) make the coefficients unstable — you'll see the signs go funny in a minute — and a single outlier, thanks to that squared loss, can lever the whole line. None of this makes it useless. It makes it a baseline you understand, which is the most useful thing a model can be.

The data

California housing is the standard first real regression dataset, and it's built into scikit-learn so there's nothing to download. Each of the 20,640 rows is a census block group — a few hundred to a few thousand people — with eight features: median income, median house age, average rooms and bedrooms, population, average occupancy, and latitude/longitude. The target is the median house value in the block, in units of $100,000. One quirk to know going in: the target is capped at 5.0, so a wall of blocks sit at exactly $500,000, which will come back to bite the fit at the top end.

The single strongest signal is median income, which makes sense — richer areas have pricier houses. Here's median income against median house value on a sample of blocks. The relationship is clearly linear-ish and clearly noisy, and you can see the cap: that flat ceiling of points at 5.0.

Build it, one function at a time

Six short functions, pure NumPy, in the order you'd actually write them: smallest testable piece first, then the two ways to fit.

The model is one line. Predictions are the features times the weights plus the intercept — a matrix-vector product and an add:

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

    X is (n, d) — n rows, d features. w is (d,), b is a scalar. The result is
    (n,), one predicted number per row. This one line is the whole model;
    everything else in the file is about choosing w and b.
    """
    return X @ w + b

Because that's a plain matrix product, the same function predicts one row or a whole dataset, one feature or eight. Next, the thing we're trying to make small — the mean squared error, the average of the squared residuals:

def mse_loss(X, y, w, b):
    """Mean squared error: the average squared gap between guess and truth.

    This is the thing we minimize. Squaring makes every miss positive and
    punishes big misses far more than small ones, which is what pins the fit
    to the bulk of the data (and what makes outliers hurt).
    """
    resid = predict(X, w, b) - y
    return float(np.mean(resid ** 2))

The loss tells you how wrong a line is; R² tells you whether that's any good, by comparing your error against the error of just guessing the mean every time. It's the number people actually report, so we implement it directly:

def r2_score(X, y, w, b):
    """Coefficient of determination: fraction of the 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.
    """
    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 fit. The closed form solves for the minimum in one shot. We glue a column of ones onto the features so the intercept comes along as an ordinary weight, then solve the normal equations. Note we solve the system rather than literally inverting $X^{\top}X$ — same answer, better numerics:

def normal_equation(X, y):
    """The closed-form OLS solution: w = (XtX)^-1 Xt y, intercept and all.

    We glue a column of ones onto X so the intercept rides along as the first
    weight, then solve the normal equations. `np.linalg.solve` factors the
    matrix instead of forming the inverse explicitly — same answer, better
    conditioned. Splits the fitted vector back into (w, b) at the end.
    """
    ones = np.ones((X.shape[0], 1))
    Xb = np.hstack([ones, X])                 # (n, d+1)
    theta = np.linalg.solve(Xb.T @ Xb, Xb.T @ y)
    b = float(theta[0])
    w = theta[1:]
    return w, b

That's linear regression, done. Everything from here is the second route to the same answer, the one we can watch. Gradient descent needs the slope of the loss at the current guess, so it knows which way is downhill. Differentiate the MSE and the residual falls out front:

def gradient(X, y, w, b):
    """The slope of the MSE at the current (w, b): which way is downhill.

    Differentiate the mean squared error. The residual falls out front, so
    the gradient is just the residual projected back onto each feature (and
    summed, for the intercept), scaled by 2/n.
    """
    n = X.shape[0]
    resid = predict(X, w, b) - y              # (n,)
    grad_w = (2.0 / n) * (X.T @ resid)        # (d,)
    grad_b = (2.0 / n) * float(resid.sum())
    return grad_w, grad_b

A step is just moving against that gradient by the learning rate:

def gd_step(w, b, grad_w, grad_b, lr):
    """One downhill step: move against the gradient by the learning rate."""
    return w - lr * grad_w, b - lr * grad_b

And the fit is starting from a flat line at zero and taking those steps until they stop changing anything. Because the loss is a bowl, this can't get stuck — it rolls to the same bottom the normal equation jumps to:

def gradient_descent(X, y, lr, n_iters, record=False):
    """Fit by repeatedly stepping downhill from w = 0, b = 0.

    With a small enough learning rate this walks to the same (w, b) the normal
    equation hands you in one shot — the MSE surface is a bowl, so there's a
    single minimum to fall into. When `record` is on we log (w, b, loss) BEFORE
    each step so the chapter can animate the line converging.
    """
    w = np.zeros(X.shape[1])
    b = 0.0
    history = []
    for _ in range(n_iters):
        if record:
            history.append((w.copy(), b, mse_loss(X, y, w, b)))
        grad_w, grad_b = gradient(X, y, w, b)
        w, b = gd_step(w, b, grad_w, grad_b, lr)
    if record:
        history.append((w.copy(), b, mse_loss(X, y, w, b)))
    return w, b, history

Watch it work

This is gradient descent fitting the synthetic line, one real step per frame. The top panel is the data with the current line drawn over it; the line starts flat at zero — slope 0, intercept 0 — and rotates and lifts into place as the steps run. The bottom panel traces the MSE on a log axis, so you can watch the loss fall off a cliff early and then flatten as the line settles. The caption names the iteration, the current slope and intercept, and the current MSE.

Press play. Every frame is one call to the gradient-and-step loop from impl.py, at learning rate 0.12. Reset and run it as many times as you like.

Two things to notice. The line moves fastest at the start, when the residuals are largest and the gradient is steep, and it barely moves at the end — that's the loss curve flattening, the model creeping the last little bit into the minimum. And the place it lands is not the line we generated the data from. We built the points with slope 2.6 and intercept 5.0, but the least-squares fit comes out at slope 3.00 and intercept 5.10, because forty noisy points don't pin the line exactly. After 60 iterations gradient descent has matched the closed-form solution to two decimals; both agree the best line through these points has slope 3.00, not the 2.6 we started from. That gap between the true generating line and the best fit to a finite sample is not a bug — it's the thing statistics is about.

The full implementation

The whole file, no library, top to bottom — the model, the two metrics, the normal equation, and gradient descent. This is the code the animation above actually ran:

"""Ordinary least squares linear regression, built from scratch.

A linear model predicts a number as a weighted sum of the features plus an
intercept. "Fitting" means choosing the weights that make the mean squared
error smallest. There are two ways to get there and this file has both: the
closed-form normal equation (one linear solve, exact) and gradient descent
(take a step downhill, repeat) so we can watch the fit converge.

Pure NumPy — no ML library anywhere in this file. Every function below shows
up in the chapter one step at a time; the `# region:` markers are what the
book'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.

    X is (n, d) — n rows, d features. w is (d,), b is a scalar. The result is
    (n,), one predicted number per row. This one line is the whole model;
    everything else in the file is about choosing w and b.
    """
    return X @ w + b
# endregion


# region: mse_loss
def mse_loss(X, y, w, b):
    """Mean squared error: the average squared gap between guess and truth.

    This is the thing we minimize. Squaring makes every miss positive and
    punishes big misses far more than small ones, which is what pins the fit
    to the bulk of the data (and what makes outliers hurt).
    """
    resid = predict(X, w, b) - y
    return float(np.mean(resid ** 2))
# endregion


# region: r2
def r2_score(X, y, w, b):
    """Coefficient of determination: fraction of the 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.
    """
    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: normal_equation
def normal_equation(X, y):
    """The closed-form OLS solution: w = (XtX)^-1 Xt y, intercept and all.

    We glue a column of ones onto X so the intercept rides along as the first
    weight, then solve the normal equations. `np.linalg.solve` factors the
    matrix instead of forming the inverse explicitly — same answer, better
    conditioned. Splits the fitted vector back into (w, b) at the end.
    """
    ones = np.ones((X.shape[0], 1))
    Xb = np.hstack([ones, X])                 # (n, d+1)
    theta = np.linalg.solve(Xb.T @ Xb, Xb.T @ y)
    b = float(theta[0])
    w = theta[1:]
    return w, b
# endregion


# region: gradient
def gradient(X, y, w, b):
    """The slope of the MSE at the current (w, b): which way is downhill.

    Differentiate the mean squared error. The residual falls out front, so
    the gradient is just the residual projected back onto each feature (and
    summed, for the intercept), scaled by 2/n.
    """
    n = X.shape[0]
    resid = predict(X, w, b) - y              # (n,)
    grad_w = (2.0 / n) * (X.T @ resid)        # (d,)
    grad_b = (2.0 / n) * float(resid.sum())
    return grad_w, grad_b
# endregion


# region: gd_step
def gd_step(w, b, grad_w, grad_b, lr):
    """One downhill step: move against the gradient by the learning rate."""
    return w - lr * grad_w, b - lr * grad_b
# endregion


# region: gradient_descent
def gradient_descent(X, y, lr, n_iters, record=False):
    """Fit by repeatedly stepping downhill from w = 0, b = 0.

    With a small enough learning rate this walks to the same (w, b) the normal
    equation hands you in one shot — the MSE surface is a bowl, so there's a
    single minimum to fall into. When `record` is on we log (w, b, loss) BEFORE
    each step so the chapter can animate the line converging.
    """
    w = np.zeros(X.shape[1])
    b = 0.0
    history = []
    for _ in range(n_iters):
        if record:
            history.append((w.copy(), b, mse_loss(X, y, w, b)))
        grad_w, grad_b = gradient(X, y, w, b)
        w, b = gd_step(w, b, grad_w, grad_b, lr)
    if record:
        history.append((w.copy(), b, mse_loss(X, y, w, b)))
    return w, b, history
# endregion


def load_data(path="../data/california_housing.csv"):
    """California housing: 20,640 block groups, 8 features, median house value.

    Returns (X, y, feature_names). y is the median house value in units of
    $100,000, the column the model learns to predict.
    """
    df = pd.read_csv(path)
    target = "MedHouseVal"
    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 normal equation in production, and once you understand it you shouldn't either. sklearn.linear_model.LinearRegression is the same model with the same objective — ordinary least squares — and it fits the intercept by default:

def sklearn_fit(X_train, y_train):
    """Fit OLS with sklearn. Returns (weights, intercept) — the twin of our
    normal_equation, so the coefficients line up one for one."""
    model = LinearRegression()            # fit_intercept=True by default
    model.fit(X_train, y_train)
    return model.coef_, float(model.intercept_)

Under the hood sklearn doesn't form $(X^{\top}X)^{-1}$ any more than we do. It solves the least-squares system with an SVD, which is more robust than the normal equation when features are nearly collinear — and California housing has some of that, with average rooms and average bedrooms moving together and latitude and longitude tracing the coastline. On well-conditioned data the two routes give identical coefficients; on ugly data the SVD degrades more gracefully. For scoring we let it fit and report test R² and MSE:

def sklearn_score(X_train, y_train, X_test, y_test):
    """Fit on train, report test R^2 and test MSE — the face-off numbers."""
    model = LinearRegression()
    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

The other thing worth reading off the fitted model is the coefficients, because this is where linear regression earns its keep. Median income comes out around 0.44 — one extra unit of median income (that's $10,000) predicts about $43,900 more in median house value, holding the rest fixed. Latitude and longitude are both strongly negative, around -0.42 and -0.44, which is the model discovering the map: prices fall as you move north and inland from the expensive southern coast. And average bedrooms and average rooms come out with opposite signs (+0.65 and -0.11) even though both are "more house" — that's multicollinearity, the two correlated features splitting the credit between them in a way that makes each coefficient hard to read on its own. Interpretability is the selling point, but it comes with fine print.

Scratch versus library

Split the data 80/20 — 16,512 blocks to train on, 4,128 held out — fit on the training half three ways, and score on the held-out half. Our normal equation, our gradient descent (on standardized features, since the raw scales differ by orders of magnitude), and sklearn:

All three bars are the same length: test R² of 0.585, test MSE of 0.528. They're not close, they're identical — our normal equation matches sklearn's coefficients to twelve decimal places, and gradient descent, given enough steps, lands on the same predictions. There's only one least-squares line, and every correct method finds it. 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.585 means the model explains about 59% of the variance in house values from these eight features — decent for a straight-out-of-the-box linear fit, and nowhere near the ceiling. To see where it leaks, plot the residuals against the fitted values. In a clean fit this is a shapeless cloud centered on zero; structure here means the model is missing something:

The cloud isn't shapeless, and that's the lesson. There's a diagonal streak of points running down the right side — those are the capped blocks, the ones truly worth $500,000 or more that the model, having no idea about the cap, predicts as high as $717,000 and then eats a big negative residual on. And the spread of the residuals widens as the fitted value grows, a fan shape: that's heteroscedasticity, the constant-variance assumption failing in plain sight. The model is more wrong, and wrong less predictably, on expensive houses than on cheap ones. Linear regression will happily hand you coefficients over data like this; the residual plot is where it tells you what it couldn't fit.

Takeaways

Fit the line first. Same advice as the threshold, one rung up: before the gradient boosting and the neural net, run ordinary least squares, read the R², and look at the coefficients. That number is the bar your fancier model has to clear, and those coefficients are a free sanity check on whether the features even point the right way. If income comes out negative, you have a bug or a data problem, and you want to find that out from a model you can read, not from a black box.

The two builds in this chapter are the two halves of the rest of the course. The normal equation is the last time you'll get the answer in closed form — it's a gift that linear regression's loss is a convex bowl with a formula for the bottom. Everything harder gives that up, and when it does, gradient descent is the tool that still works: the exact loop you watched fit this line is what trains logistic regression, and boosting, and every neural network later on, just on a loss surface with no formula and, eventually, no guarantee of a single minimum. Learn it here, on the one problem where you can check its work against an exact answer.

And keep the residual plot. The coefficients tell you what the model learned; the residuals tell you what it couldn't. The fan and the diagonal streak in that last chart are the model being honest about its own assumptions — linearity and constant variance — right where they break. Reach for something heavier when the residuals show structure the line can't absorb, which on real data is most of the time. But start here, because a linear fit you understand is worth more than a complicated one you don't.