Linear Regression: The Reasoning, the Math, the Geometry

What this lesson delivers

• A precise statement of the regression problem, in Mitchell's <T, P, E> terms.
• The linear model ŷ = wx + b, its parameters, and its vectorized multi-feature form.
• Why we measure error with mean squared error specifically — and what that choice costs us.
• A hand-derivation of the closed-form solution for the best line — no optimizer required.
• Two geometries: the best-fit line in data space, and least-squares as a projection.
• A from-scratch NumPy implementation, verified against Scikit-learn.

1 · The problem: predicting a number

Regression is supervised learning where the answer is a continuous number: a price, a temperature, a duration, a dosage. Contrast it with classification (Part IV), where the answer is a category. If the quantity you want lives on a number line and "a bit more" or "a bit less" is meaningful, you have a regression problem.

Let's make it concrete with the example we'll carry through the whole lesson. You have data on houses — each one's size in square metres and the price it sold for — and you want to predict the price of a house you haven't seen. In Mitchell's template:

Notation we'll use throughout, borrowed straight from Lesson 9. Each house is a training example (x⁽ⁱ⁾, y⁽ⁱ⁾): x⁽ⁱ⁾ is its feature (size), y⁽ⁱ⁾ is its target (the true price), and i runs from 1 to n, the number of examples. Our prediction for example i is written ŷ⁽ⁱ⁾ ("y-hat") — the hat always means "the model's guess", never the truth.

2 · The model: a straight line, and what its knobs mean

A model is a family of functions, and "choosing a model" means choosing which family before you've seen which member fits. Linear Regression makes the simplest non-trivial choice: assume the target is a straight-line function of the feature.

Two numbers define the whole model — they are the parameters, the things learning will tune:

More than one feature: the vectorized form

Houses have more than a size — bedrooms, age, distance to centre. With d features, each example is a feature vector x = (x₁, …, x_d) and each feature gets its own weight. The model becomes a weighted vote of the features plus a bias — exactly the dot product from Lesson 9:

This is why we spent a lesson on the dot product: it collapses "multiply each feature by its weight and add them up" into a single operation, and stacked over every row of the dataset it becomes the matrix–vector product X w that predicts for all n houses in one line. We'll develop the single-feature case in full because you can see it, then write the code so it already handles the vector case.

3 · The loss: making "best fit" a number we can minimize

There are infinitely many lines. "Best fit" is meaningless until we pick a rule that scores every line with a single number, so that lower is better and there's something concrete to minimize. That rule is the loss function of Lesson 6, and for regression the standard choice is Mean Squared Error.

For one example, the error is the residual — the gap between prediction and truth, ŷ⁽ⁱ⁾ − y⁽ⁱ⁾. We square it (so overshooting and undershooting both count as positive error, and big misses hurt disproportionately), then average over all n examples:

Read this the way it's built. Inside: the residual for one house. Squared: the penalty for that house. Summed and divided by n: the average penalty across the dataset. J is written as a function of (w, b) because that is what we get to change — the data is fixed; the parameters are the variables. Training is now a precise sentence: find the (w, b) that make J as small as possible.

4 · Solving it exactly: the line, by hand

Here is what makes Linear Regression special enough to teach first: we don't have to search for the minimum — we can write it down. The loss J(w, b) is a smooth bowl, so its lowest point is where the ground is flat in every direction: where both partial derivatives are zero. This is the calculus fact from Lesson 9 — at a minimum, the slope is zero — applied in two variables at once.

Set the derivative with respect to the bias to zero first. Differentiating J with respect to b (chain rule: the square's exponent drops to a factor of 2, the inside's derivative w.r.t. b is 1):

Divide out the 2/n and split the sum: it says the average residual is zero, which rearranges to b = ȳ − wẁ, where ẁ and ȳ are the mean feature and mean target. In words: the best line always passes through the centre of mass of the data (ẁ, ȳ). That single fact is worth remembering on its own.

Now the derivative with respect to the slope (inside's derivative w.r.t. w is x⁽ⁱ⁾):

Substitute b = ȳ − wẁ into that and solve for w. The algebra (worth doing once by hand as an exercise) collapses to a strikingly interpretable ratio:

The general case: the normal equation

With many features you can't isolate one weight at a time, but the same idea — set all the partial derivatives to zero simultaneously — produces one compact matrix equation. Fold the bias in by adding a column of 1s to the data matrix X (so b becomes just another weight), and the minimizer is:

This is the normal equation. It gives the exact best weights in one shot, no iteration — we'll use it in the code. But notice the matrix inverse: computing it costs on the order of d³ for d features, which becomes ruinous when d is large (thousands of features), and X^TX can be numerically fragile. That cost is the entire motivation for the next lesson. When you can't afford to solve exactly, you descend iteratively instead — gradient descent — which is why Lesson 11 exists.

5 · The geometry — two ways to see it

Geometry one: the best-fit line (data space)

The picture everyone knows. Plot each example as a point — feature on the horizontal axis, target on the vertical. The model is a line through that cloud. For any point, the residual is the vertical gap from the point up or down to the line. Squaring a residual is, literally, the area of a square built on that gap; minimizing MSE means choosing the line that makes the total area of those squares as small as possible. Hence the classic name, least squares.

price
  |                                    •   ← residual = vertical gap
  |                              ______/|        from point to line
  |                        _____/    • |
  |                  _____/  |
  |            _____/     • | ← another residual
  |      •____/  •
  |  •__/    the fitted line  y_hat = w·x + b
  |_/________________________________________ size
       least squares picks the line that minimizes
       the total squared vertical distance to the points

Two things this picture makes obvious. First, the residuals are measured vertically, not perpendicular to the line — because we're predicting y from x, and error means error in the thing we predict. Second, the line pivots around the data's centre of mass (ẁ, ȳ), exactly as the b = ȳ − wẁ result promised.

6 · From scratch in NumPy

Now we earn the title. Everything above — model, loss, closed-form solution — is a handful of NumPy lines. We'll use the toy house dataset, fit it two ways (the single-feature formula and the general normal equation), and confirm they agree.

Start with the model and the loss, each a direct transcription of its equation:

import numpy as np

# Toy dataset: 6 houses. Feature = size (m2), target = price (thousands $).
x = np.array([ 50,  60,  80,  95, 120, 150], dtype=float)   # x^(i)
y = np.array([150, 175, 210, 260, 300, 360], dtype=float)   # y^(i)

# The model:  y_hat = w * x + b   (a straight line)
def predict(x, w, b):
    return w * x + b                # works element-wise on the whole array

# The loss:  mean squared error, exactly J(w, b) from section 3
def mse(y_hat, y):
    return np.mean((y_hat - y) ** 2)

Next, the single-feature closed-form solution — the covariance-over-variance ratio, then the intercept through the centre of mass:

# --- Fit by the closed-form solution (simple linear regression) ---
x_bar, y_bar = x.mean(), y.mean()

w = np.sum((x - x_bar) * (y - y_bar)) / np.sum((x - x_bar) ** 2)   # cov / var
b = y_bar - w * x_bar                                          # centre of mass

print(w, b)                       # ~ 2.104  47.87   (k$ per m2, and k$ intercept)
print(mse(predict(x, w, b), y))   # ~ 35.1   (training MSE at the optimum; RMSE ~ 5.9 k$)

The slope says each extra square metre is worth about $2,100, and the model's predictions land, on average, within roughly $5,900 of the truth (the root of the MSE, back in the target's own units). Now the same fit via the general normal equation, which will already handle many features when you swap in a wider X:

# --- The general form: fold the bias into the weights ---
# Design matrix: a column of 1s (multiplying b) beside the feature column(s).
X = np.column_stack([np.ones_like(x), x])     # shape (6, 2): [1, size] per row

# Normal equation:  theta = (X^T X)^-1 X^T y  -- but SOLVE, never invert.
theta = np.linalg.solve(X.T @ X, X.T @ y)     # @ is matrix multiply
b, w = theta                                  # -> 47.87, 2.104  (same as above)

Finally, use the fitted model to predict — the whole point. A 100 m² house:

predict(np.array([100.0]), w, b)    # array([258.28])  -> about $258k

The Scikit-learn equivalent

Everything above is what a library does for you. Now that you've built it, here's the three-line version — and it's no longer a black box, because you know precisely what coef_ and intercept_ are and how they were found.

from sklearn.linear_model import LinearRegression

model = LinearRegression()                # fits an intercept by default
model.fit(x.reshape(-1, 1), y)          # sklearn wants X as 2-D: (n_samples, n_features)

model.coef_          # array([2.104...])  <- our w
model.intercept_     # 47.87...           <- our b
model.predict([[100]])   # array([258.28...]) <- identical to our from-scratch answer

7 · What we quietly assumed (and when it breaks)

How this connects to the next lesson

Exercises & reflection

The first two are with pencil; the rest are at the keyboard. If you can do them, Linear Regression is genuinely yours.