Roadmap › Part III · Supervised Learning — Regression › Module 3 · Regression

Lesson 12 · From scratch

Multiple & Polynomial Regression — and Feature Scaling for Real

Lessons 10 and 11 built the whole machine on the smallest possible problem: one feature, a straight line, a loss you could solve exactly and then descend. This lesson takes that engine out onto real terrain. Real problems have many inputs at once — size, bedrooms, age — and the relationships you care about are rarely straight lines. The remarkable thing you'll see is how little has to change: your vectorized model Xθ already fits any number of features without a single edit, and the same trick that adds a second feature also lets you bend the line into a curve — because a polynomial is still linear in its parameters. Along the way, feature scaling graduates from the quiet prerequisite of Lesson 11 to a first-class, non-negotiable step, and we meet the shadow that will define the next lesson: the more flexible you make the model, the easier it becomes to fit noise instead of signal.

Reading: ~40 min Prerequisites: Lesson 10 (the linear model, MSE, the normal equation) · Lesson 11 (gradient descent, standardization) Code: multiple regression, feature scaling & polynomial features from scratch in NumPy

Why a straight line through one feature was never going to be enough

No house is priced by its floor area alone, and no interesting quantity in the world sits on a perfectly straight line. A useful model has to (1) combine several inputs at once, weighing each by how much it matters, and (2) capture relationships that curve — diminishing returns, acceleration, thresholds. The astonishing part is that both extensions are the same mathematical move you already know: put more columns in the design matrix X. Add a column for bedrooms and you have multiple regression; add a column for size² and you have a curve. The model stays a dot product Xθ, the loss stays a bowl, and every line of NumPy you wrote in Lesson 11 keeps working. This lesson is about seeing that, and about the one new discipline — scaling — that many features make mandatory.

What this lesson delivers

The multiple linear regression model, written once as ŷ = Xθ for any number of features.
Why the normal equation and the gradient-descent loop from Lessons 10–11 fit it with zero changes — only X gets wider.
How to read a multiple-regression coefficient honestly: "holding the others fixed," and why that phrase carries all the weight.
Feature scaling as a first-class step — standardization vs. min–max, why several features on different units force it, and the train-only discipline that keeps it honest.
Polynomial regression: bending the model into curves with x², x³ features while it stays linear in its parameters.
The degree of a polynomial as a dial between underfitting and overfitting — the exact tension Lesson 13 exists to resolve.

1 · From one feature to many

In Lesson 10 a house had a single number attached to it — its size — and the model was a line ŷ = w x + b. A real listing has a whole row of numbers: size, number of bedrooms, age, distance to the centre. Call them x₁, x₂, …, x_d. The natural generalization gives each feature its own weight and adds a single bias:

ŷ = b + w₁x₁ + w₂x₂ + … + w_dx_d

Geometrically, one feature gave a line in a 2-D plane; two features give a plane tilted in 3-D space; d features give a flat hyperplane in (d+1)-dimensional space. You can't picture the high-dimensional version, and you don't need to — the algebra is identical, and that is the whole point of writing it as a dot product.

Fold the bias into the weights exactly as in Lesson 11: prepend a constant 1 to every feature row so that θ = (b, w₁, …, w_d). Then one row's prediction is a dot product, and stacking all n rows into the design matrix X (now with d+1 columns) makes every prediction at once a single matrix–vector product:

ŷ = Xθ where X is n × (d+1), θ is (d+1) × 1

You already wrote the general model

Look back at Lesson 11's vectorized loop. The prediction line was X @ theta, the residual was X @ theta - y, and the gradient was (2/n) * (X.T @ resid). Nothing in any of those lines mentions the number of features. They were never one-feature formulas that happen to work — they are the d-feature formulas, and you tested them on the case d = 1 because it was the one you could check by hand. Widening X from two columns to twenty changes the shape of an array and nothing else. This is the dividend of having vectorized early.

2 · Fitting it changes nothing

Both roads to the parameters carry straight over. The normal equation from Lesson 10 was derived by setting ∇J = 0, and that derivation never assumed a particular number of features — it gives the same expression here:

θ = (X^TX)⁻¹X^Ty

with X^TX now a (d+1) × (d+1) matrix instead of 2 × 2. And the gradient-descent loop is byte-for-byte the one you already ran. Here is multiple regression from scratch — three features, solved exactly — with the one habit worth forming now: never invert a matrix you can solve. Forming (X^TX)⁻¹ explicitly is slower and less numerically stable than asking a solver for θ directly, so we use np.linalg.solve, which computes the same answer the normal equation defines without ever building the inverse:

import numpy as np

# Houses: [size m2, bedrooms, age years]  ->  price (thousands $)
X_raw = np.array([
    [ 50, 1, 20],
    [ 60, 2, 15],
    [ 80, 2, 30],
    [ 95, 3,  8],
    [120, 3, 25],
    [150, 4,  5],
], dtype=float)
y = np.array([150, 175, 210, 260, 300, 360], dtype=float)

# Standardize EACH column to mean 0, std 1 (axis=0 works down the columns)
mu, sigma = X_raw.mean(0), X_raw.std(0)
Xs = (X_raw - mu) / sigma

# Add the bias column of 1s -> design matrix, shape (6, 4)
X = np.column_stack([np.ones(len(y)), Xs])

# Normal equation, solved WITHOUT inverting: theta = (X^T X)^-1 X^T y
theta = np.linalg.solve(X.T @ X, X.T @ y)
print(theta)   # [b, w_size, w_bed, w_age] in standardized space

And the descent loop — identical to Lesson 11, run on the wider X — converges to the very same θ:

theta = np.zeros(X.shape[1])   # now length 4, not 2 -- the only visible difference
alpha, n = 0.1, len(y)

for step in range(2000):
    resid = X @ theta - y            # same line
    grad  = (2 / n) * (X.T @ resid)   # same line
    theta -= alpha * grad            # same line

print(theta)   # matches the solve() answer to many digits

One caveat the extra features introduce: you need enough rows

The normal equation can be solved for a unique θ only when X^TX is invertible, which requires at least as many independent rows as parameters — loosely, n ≥ d+1, with no feature being an exact copy or linear combination of others. With one feature this was almost impossible to violate; with many features it is a real failure mode. Two perfectly correlated columns (say, area in m² and the same area in ft²) make X^TX singular and the inverse undefined — the model can't decide how to split the credit between two identical inputs. This is multicollinearity, and it's one more reason the regularization of Lesson 13 exists: adding a small term to X^TX makes it invertible again.

3 · Reading the coefficients — the two words that matter

A single-feature slope had an easy story: "each extra square metre adds w thousand dollars." With many features the story gains a crucial clause. The coefficient w_j is the change in the prediction for a one-unit increase in feature j, holding all other features fixed. That last phrase is doing enormous work, and forgetting it is the most common way to misread a regression.

Why "holding the others fixed" is not a technicality

Suppose bigger houses also tend to have more bedrooms. If you looked at price against bedrooms alone, the bedroom coefficient would quietly absorb some of the size effect — bedrooms would look more valuable than they are, because they're standing in for the size that comes with them. Multiple regression untangles this: with size also in the model, the bedroom coefficient measures what an extra bedroom adds to two houses of the same size. That is a genuinely different, and usually more honest, quantity. It's also why a coefficient can flip sign when you add a feature — not a bug, but the model re-attributing credit that one feature had been borrowing from another.

Correlation, still not causation

"Holding the others fixed" is a statement about the fitted line, not about the world. A positive age coefficient does not mean renovating a house to be older would raise its price; it means that, among the houses in this sample, and after accounting for the other columns, older ones happened to be associated with higher prices — perhaps because they sit in older, more central neighbourhoods the model can't see. Regression coefficients describe associations in the data you fed it. Turning them into claims about interventions requires assumptions the arithmetic alone can never supply.

4 · Feature scaling, now a first-class step

In Lesson 11 you standardized a single feature to rescue gradient descent from a narrow valley, and I promised the topic a full treatment here. The reason is simple: with one feature, scaling was a convenience; with many features on wildly different units — size in the tens or hundreds, bedrooms in single digits, age in years — it becomes structural. Two distinct things go wrong without it.

Optimization breaks — unequal scales stretch the loss bowl into the long, narrow canyon of Lesson 11, but now in several directions at once. A single learning rate cannot serve a feature ranging over 150 and one ranging over 4; descent crawls or diverges.

Coefficients become incomparable — a raw weight of 2.1 “per m²” and 15 “per bedroom” say nothing about which feature matters more, because a metre and a bedroom aren't the same-sized change. Only after scaling do the weights live on a common footing.

Two standard rescalings cover almost every case. Both are computed per column and are pure, invertible arithmetic:

Standardization vs. min–max normalization

Standardization (z-scoring): z = (x − μ) / σ — subtract the feature's mean, divide by its standard deviation. Each feature ends with mean 0 and spread 1. This is the default for regression and anything gradient-based: it's robust, it doesn't assume bounded inputs, and it centres the data so the bias behaves.

Min–max normalization: x' = (x − x_min) / (x_max − x_min) — squash each feature into the [0, 1] range. Handy when you need bounded inputs or the feature has hard limits, but a single outlier can compress everything else into a sliver, so it's more fragile than standardization.

Neither changes the shape of the relationship or the model's predictions — they change the coordinates the optimizer works in, and you can always convert the learned parameters back to original units, as we did in Lesson 11.

The discipline that prevents a silent bug: fit the scaler on training data only

Compute μ and σ (or the min and max) from the training set, store them, and apply those same numbers to the validation and test sets. Never recompute them on the test data. Doing so lets information about the test set — its mean, its spread — leak into how the model was prepared, which quietly inflates your measured performance and betrays you in production, where tomorrow's data has no mean yet. This is your first concrete instance of data leakage, the subject of Lesson 19, and the reason real pipelines package a scaler as an object that is fit once and then only applied. Our two lines mu, sigma = X_raw.mean(0), X_raw.std(0) are exactly that scaler; the discipline is where you compute them.

5 · Bending the line — polynomial features

Multiple regression combines several inputs, but every one of its predictions is still a flat hyperplane: double any feature and the prediction moves by a fixed amount. Real relationships curve. Price per extra square metre tapers off for mansions; a drug's effect saturates; cooling slows as an object approaches room temperature. How do we fit a curve with a machine that only knows how to fit flats?

The trick is one of the most elegant sleights of hand in classical ML. We don't change the model at all — we change the features. Take a single input x and manufacture new columns from it: x, x², x³. Then fit an ordinary linear model on those columns:

ŷ = b + w₁x + w₂x² + w₃x³

As a function of x, that is a cubic — a genuine curve that can rise, bend, and turn. But look at it as a function of the parameters b, w₁, w₂, w₃: it is a plain weighted sum, exactly the linear form Xθ. The model is linear in its parameters even though it is nonlinear in x, and "linear regression" has always meant the former. So the normal equation, the gradient, the entire Lesson 11 loop — all of it fits a curve without knowing it's doing anything new.

The one idea to keep: "linear" refers to the parameters, not the picture

This is the sentence that unlocks half of classical ML. A model is a linear model if the prediction is a weighted sum of features with the parameters as the weights — regardless of how curved or exotic those features are. Squaring, cubing, taking logs, multiplying two features together (an interaction, like size × bedrooms): each is just another column you compute up front and hand to the same linear machinery. The features can be as nonlinear as you like; as long as you combine them linearly, you're still solving the convex bowl you already know how to solve. Polynomial regression is not a new algorithm — it's feature engineering feeding the old one.

Here it is from scratch. The only new code is a function that builds the powers; everything after it is the fitting you already own. Notice that scaling is now unavoidable: with x up to 150, the x³ column reaches into the millions, so the raw columns span an enormous range of magnitudes — the canyon of Lesson 11, guaranteed.

# One input, but a curved target: y bends as x grows
x = np.array([1, 2, 3, 4, 5, 6, 7, 8], dtype=float)
y = np.array([2.1, 3.9, 6.2, 9.1, 12.8, 17.0, 22.1, 27.8])

def poly_features(x, degree):
    # columns [x, x^2, ..., x^degree]  (bias added separately)
    return np.column_stack([x ** p for p in range(1, degree + 1)])

P = poly_features(x, degree=2)         # shape (8, 2): [x, x^2]

# Standardize the manufactured columns -- ESSENTIAL: x^2 dwarfs x
mu, sigma = P.mean(0), P.std(0)
Ps = (P - mu) / sigma

X = np.column_stack([np.ones(len(y)), Ps])   # [1, x, x^2], standardized
theta = np.linalg.solve(X.T @ X, X.T @ y)    # same fit as always
print(theta)   # the quadratic that best explains the curve

To predict a new point, transform it the same way

A polynomial model isn't just θ — it's θ plus the exact recipe that built the features: the degree, and the μ/σ used to scale each column. To score a new x, you must raise it to the same powers and subtract/divide by the stored training statistics before the dot product. Change the recipe and the parameters are meaningless. This is why, in practice, the transform and the model travel together as a single pipeline object (Lesson 27) — so a stray point can never be scored through the wrong preprocessing.

6 · The degree is a dial — and it can go too far

Polynomial regression hands you a tempting knob: the degree. Degree 1 is the old straight line. Degree 2 gives a single bend, degree 3 an inflection, and so on — each added power lets the curve wiggle once more. It is tempting to think higher is always better, because a higher-degree curve fits the training points more closely. It fits them too closely, and that is the whole problem.

The extreme case that proves the danger

A polynomial of degree n−1 can pass exactly through any n points — zero training error, a perfect fit. With 8 data points, a degree-7 polynomial hits every one of them dead on. And it is almost always useless: to thread all eight points it must contort into wild oscillations between them, swinging far above and below any sensible trend. Feed it a ninth point and it can miss by a mile. It has fit the noise — the random jitter in these samples — and mistaken it for signal. That is overfitting, and the degree dial is the cleanest place in all of ML to watch it happen.

Degree	What the curve does	Failure mode
1 (a line)	Too rigid to follow a real bend; misses the trend systematically.	Underfitting — high bias.
2–3	Follows the genuine shape, ignores the jitter; predicts new points well.	The sweet spot — for this data.
7 (with 8 points)	Threads every training point exactly, oscillating wildly between them.	Overfitting — high variance.

The tension you can now feel, named

Too simple a model can't capture the truth and is wrong in a consistent, structural way — that's bias. Too flexible a model captures the truth and the noise, so its answer swings depending on which random sample it happened to see — that's variance. The degree of a polynomial trades one directly against the other, and there is a middle setting that minimizes the error on data the model hasn't seen. Finding that setting — and the technique of regularization that lets you use a flexible model without letting it overfit — is precisely what Lesson 13 is about. You've just built the perfect apparatus for studying it.

7 · What to keep honest about

The fine print behind a multi-feature, curved fit

More features is not free. Each column added is a parameter to estimate; with too few rows relative to columns, X^TX becomes ill-conditioned or singular, and the fit turns unstable. Flexibility must be paid for in data.
Polynomial features are highly correlated. x, x², and x³ rise together, so they carry overlapping information — a mild multicollinearity that makes coefficients wobbly and is another reason scaling and, later, regularization help.
Polynomials extrapolate disastrously. Outside the range of your training x, a high-degree curve shoots off toward ±∞. Trust these models to interpolate within the data you've seen, never to predict far beyond it.
Scaling is now mandatory, not advisory. The moment features span different magnitudes — and manufactured powers guarantee they do — unscaled gradient descent is the Lesson 11 canyon. Scale first, always.
A perfect training fit is a warning, not a triumph. Near-zero error on the data you trained on tells you almost nothing; only error on held-out data is evidence. Lesson 19 makes that measurement rigorous.

How this connects to the next lesson

Bridge to Lesson 13

You now have a model flexible enough to be dangerous. Adding features and raising the polynomial degree gives you a dial that runs from a rigid line that underfits to a frantic curve that overfits — and you've seen that fitting the training data perfectly is exactly the wrong goal. Lesson 13 — Overfitting, the Bias–Variance Trade-off, and Regularization makes this the main event. We'll define bias and variance precisely, watch the trade-off on a train-vs-test error curve, and meet Ridge and Lasso — two ways to keep a flexible model honest by penalizing large weights, which also cures the multicollinearity and the singular X^TX we flagged along the way. The apparatus is built; next we learn to tune it.

Exercises & reflection

The first is with pencil; the rest extend the code above. If you can do them, multiple and polynomial regression are genuinely yours.

Exercise 1 · Why the model is still “linear”

Explain in your own words why ŷ = b + w₁x + w₂x² counts as a linear model even though its graph against x is a curve. Then decide: is ŷ = b + w₁x + sin(w₂x) a linear model? Why or why not?

Show the answer

The first is linear because it is a weighted sum of fixed columns (1, x, x²) with the parameters as the weights — linearity is a statement about the parameters, and here each parameter appears to the first power and is only multiplied by data. The curve in x is irrelevant. The second is not linear: w₂ sits inside a sine, so the prediction is not a weighted sum of the parameters — you can't write it as Xθ, and the normal equation no longer applies. Feature engineering can only pre-compute columns from the data; a parameter buried inside a nonlinearity leaves the linear world.

Exercise 2 · Fit the three-feature model two ways

Run the section-2 code. Confirm that np.linalg.solve and the gradient-descent loop land on the same θ to several digits. Then convert the standardized weights back toward a sentence: which of size, bedrooms, and age moves the price most, per one standard deviation of change?

Show discussion

The two methods agree because they solve the same convex bowl — one exactly, one by descent. Because the features were standardized, the coefficients are now directly comparable: the one with the largest magnitude is the feature whose typical (one-σ) variation moves the price most. That comparison is only legitimate after scaling — on raw units, a big coefficient might just mean a small-numbered feature. This is exactly the “coefficients become incomparable” point of section 4, made concrete.

Exercise 3 · Feel the unscaled powers explode

Build degree-3 polynomial features for the section-5 data without standardizing, and try the gradient-descent loop at α = 0.1. Watch it overflow. Estimate how small α must get to stay stable, then compare with the standardized version. Why is the problem worse than the single-feature canyon of Lesson 11?

Show discussion

The columns are x (up to 8), x² (up to 64), and x³ (up to 512) — magnitudes spanning nearly three orders even on tiny inputs, so the loss bowl is stretched in several directions at once. The steepest direction (the x³ axis) caps the learning rate at something minuscule while the shallow directions barely move, so descent either diverges or crawls. Standardizing collapses all three columns onto the same spread, rounding the bowl — and the fit converges in a few hundred steps. Manufactured powers make scaling non-negotiable.

Exercise 4 · Overfit on purpose

Fit the 8-point section-5 data at degrees 1, 2, and 7. For each, compute the training error and plot the fitted curve on a fine grid of x between the data points. Describe what the degree-7 curve does between the training points, and predict what would happen at x = 9.

Show discussion

Degree 1 misses the bend (underfits); degree 2 tracks the true shape with small residuals; degree 7 drives the training error to essentially zero but oscillates wildly between the points, and at x = 9 — just outside the data — it shoots off dramatically. This is overfitting and disastrous extrapolation in one picture: perfect on the training points, worthless everywhere else. It is the exact motivation for holding out test data (Lesson 19) and for regularization (Lesson 13).

Exercise 5 · Add an interaction, keep it linear

Back on the house data, add a new column equal to size × bedrooms (an interaction term) alongside the originals, standardize, and refit. Is this still a linear model? What might a nonzero coefficient on the interaction mean in plain English?

Show discussion

Still perfectly linear — the interaction is just one more pre-computed column, a weighted sum with the parameters as weights, so Xθ and the normal equation apply unchanged. A positive interaction coefficient would mean the value of an extra bedroom grows with the size of the house (bedrooms matter more in big homes than small ones) — a relationship a purely additive model can't express. Interactions and polynomials are the same idea: expand the feature set, fit the old linear machine.

One-paragraph summary to lock it in

Multiple regression predicts from many features at once with a weighted sum ŷ = Xθ, and because you vectorized in Lesson 11, the normal equation and the gradient-descent loop fit it with no changes — only X grows wider. Each coefficient reads as the effect of its feature holding the others fixed, a phrase that untangles correlated inputs but never licenses a causal claim. With several features on different units, feature scaling stops being optional: standardization (z = (x−μ)/σ) rounds the loss bowl so descent converges and puts coefficients on a comparable footing — and the scaler must be fit on training data only, your first taste of avoiding data leakage. Polynomial regression bends the model into curves by manufacturing columns x, x², x³ and fitting the same linear machine, because a model is “linear” in its parameters, not its picture — interactions and other nonlinear features work the same way. The polynomial degree is a dial from a rigid, underfitting line to a frantic, overfitting curve that threads every training point and generalizes horribly — the bias–variance tension that Lesson 13 exists to tame.