Linear & Polynomial Regression

Five flats sold in the same postcode last month.

A sixth flat, 75 m², comes on the market. What is it worth?

You can eyeball it. Prices go up by roughly £15,000–£50,000 per 10 m², so something like £280,000–£290,000. But "eyeball it" does not scale to 40,000 flats and twelve features, and more importantly, it gives no way to say which of two guesses is better.

Try drawing a straight line through those five points on paper. You will find you can draw many lines that look reasonable. One passes through the first and last points. One splits the difference. One hugs the middle three and ignores the ends. Every one of them gives a different answer for the 75 m² flat, and nothing about the picture says which is right.

That is the actual problem linear regression solves. Not "draw a line" — anyone can draw a line. It is define precisely what makes one line better than another, then find the best one.

Defining "best"

A line is $\hat{y} = \beta_0 + \beta_1 x$: an intercept and a slope. For any choice of those two numbers, each data point has a residual — the vertical gap between what actually happened and what the line predicted:

A good line makes the residuals small. But small in what sense? Three candidate rules, and the choice between them has real consequences:

Least squares picks the third. The objective, usually averaged to give mean squared error, is:

The squaring is a genuine design decision with a genuine cost. Squaring means an error of 20 counts four times as much as an error of 10, and an error of 100 counts a hundred times as much. That makes the fit extremely sensitive to outliers. One mis-recorded sale at £4 million in a dataset of £300,000 flats will drag the whole line towards it, because reducing that one enormous squared error is worth more to the objective than fitting every other point well.

Least squares does not fit the typical point. It fits so as to avoid catastrophic errors, which means a single wrong row can move the entire model.

Solving it exactly

Because the objective is a smooth quadratic in the parameters, it has exactly one minimum, and calculus finds it. For the single-feature case:

Work it through on the five flats. $\bar{x} = 70$, $\bar{y} = 269$.

So $\beta_1 = 3200 / 1000 = 3.2$ and $\beta_0 = 269 - 3.2 \times 70 = 45$.

The fitted line is $\hat{y} = 45 + 3.2x$. The 75 m² flat is predicted at $45 + 3.2 \times 75 = 285$, which is £285,000.

And now the coefficients say something in English: each additional square metre is associated with £3,200 of price, and the intercept of £45,000 represents everything not explained by area — the land, the postcode, the fact that a zero-square-metre flat is not a thing. That intercept is not physically meaningful and does not need to be; it is where the line happens to cross when extended back to zero.

The general case

With many features, stack the data into a matrix $X$ (one row per example, one column per feature, plus a column of ones for the intercept) and solve all coefficients at once:

This is the normal equation. It gives the exact answer in one step, with no iteration and no learning rate to tune. Its limitation is the matrix inversion, which costs roughly $O(p^3)$ in the number of features. At 50 features that is instant; at 50,000 features it is not. It also fails outright when $X^\top X$ cannot be inverted, which happens when two features are perfectly correlated or when you have more features than rows.

Gradient descent, and when you need it

The alternative is to start with arbitrary coefficients and walk downhill. The gradient of the loss with respect to each coefficient is:

and you repeatedly update $\beta_j \leftarrow \beta_j - \alpha \frac{\partial J}{\partial \beta_j}$ with a small step size $\alpha$.

1import numpy as np23def fit_gradient_descent(X, y, lr=0.01, epochs=1000):4    n, p = X.shape5    X_b = np.c_[np.ones(n), X]          # prepend intercept column6    beta = np.zeros(p + 1)7    for _ in range(epochs):8        residual = y - X_b @ beta9        grad = -(2 / n) * X_b.T @ residual10        beta -= lr * grad11    return beta

For plain linear regression on modest data, gradient descent is strictly worse than the normal equation — slower, approximate, and it needs a learning rate. It earns its place when the feature count is huge, when the data does not fit in memory, or when the model is not linear regression at all and no closed form exists.

One practical warning: gradient descent on unscaled features behaves terribly. If area ranges over 50–200 and num_bedrooms over 1–5, the loss surface is a long narrow valley, and the algorithm zig-zags across it instead of running down it. Standardising features first turns the valley into a bowl and can cut the required iterations by an order of magnitude. The normal equation does not care about scaling; gradient descent cares enormously.

Interpreting coefficients, carefully

With several features the model is $\hat{y} = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p$, and the standard reading of $\beta_j$ is "the change in $y$ for a one-unit increase in $x_j$, holding all other features constant".

Two things regularly go wrong with that sentence.

"Holding all else constant" may be physically impossible. If your features are area and num_rooms, you cannot increase area by 1 m² while holding room count fixed in any real building stock — they move together. The coefficient is still mathematically defined, but it describes a comparison between flats that barely exist in your data, so it is estimated from very little information and swings wildly between samples. This is multicollinearity, and its symptom is coefficients with implausible magnitudes or signs. Predictions can remain fine while the individual coefficients become nonsense.

Coefficient size is not importance. A coefficient of 3.2 on square metres and 15,000 on "has parking" are not comparable, because the units differ. Comparing them requires standardising the features first, so every coefficient answers "per one standard deviation".

And the ever-present caveat: these are associations, not causes. A regression on hotel data will happily tell you that rooms with more towels are more expensive. Removing the towels will not raise the price.

What linear regression assumes, and how it breaks

The single most useful diagnostic is the residual plot: predicted value on the x-axis, residual on the y-axis. If the model is right, this is a formless cloud centred on zero. Any structure — a curve, a funnel, a cluster — is the model telling you what it failed to capture.

1import matplotlib.pyplot as plt23resid = y_test - model.predict(X_test)4plt.scatter(model.predict(X_test), resid, alpha=0.5)5plt.axhline(0, color="black", linewidth=1)6plt.xlabel("predicted"); plt.ylabel("residual")

When a straight line is not enough

Suppose the real relationship curves. Sales rise with advertising spend, but with diminishing returns: the first £10,000 buys a lot of new customers, the tenth £10,000 buys very few. A straight line fitted to that data will overpredict at the extremes and underpredict in the middle, and its residual plot will show an obvious arc.

The fix is smaller than it looks. Polynomial regression is still linear regression. You do not change the algorithm at all — you change the features you hand it.

Given one feature $x$, create $x$, $x^2$, $x^3$ and treat those as three separate features. The model

is curved in $x$ but perfectly linear in $\beta$ — and "linear" in "linear regression" refers to the coefficients, not the shape of the curve. All the machinery above applies unchanged, including the exact normal-equation solution.
Python
Copy
1from sklearn.preprocessing import PolynomialFeatures2from sklearn.linear_model import LinearRegression3from sklearn.pipeline import make_pipeline4from sklearn.preprocessing import StandardScaler56model = make_pipeline(7    PolynomialFeatures(degree=3, include_bias=False),8    StandardScaler(),          # essential: x^3 is enormous compared to x9    LinearRegression(),10)11model.fit(X_train, y_train)
The scaler is not decoration. With area around 80, the cubic term is around 512,000. Without scaling, the coefficients for different powers span six orders of magnitude, and any subsequent regularisation would penalise them absurdly unevenly.
The cost: features multiply fast
PolynomialFeatures also creates interaction terms, which is usually what you want but is expensive. Starting from $p$ features at degree $d$, you get $\binom{p+d}{d} - 1$ features:
Original features
Degree 2
Degree 3
Degree 4
3
9
19
34
10
65
285
1,000
50
1,325
23,425
316,250
Fifty features at degree 3 gives 23,425 features. If you have 5,000 rows, you now have far more features than examples, and the model can fit the training data perfectly while learning nothing. This is why polynomial expansion beyond degree 2 or 3 is rare in practice, and why it is nearly always paired with regularisation.
Choosing the degree honestly
Never pick the degree by training error, which falls forever. Use cross-validation.
Python
Copy
1import numpy as np2from sklearn.model_selection import cross_val_score34for d in range(1, 9):5    pipe = make_pipeline(PolynomialFeatures(d, include_bias=False),6                         StandardScaler(), LinearRegression())7    rmse = -cross_val_score(pipe, X_train, y_train, cv=5,8                            scoring="neg_root_mean_squared_error")9    print(f"degree {d}: CV RMSE {rmse.mean():7.2f} +/- {rmse.std():.2f}")
A typical output on genuinely quadratic data with noise:
Degree
Train RMSE
CV RMSE
Read
1
18.40
18.71
Underfits — cannot bend
2
4.92
5.14
Best — matches the truth
3
4.88
5.36
No gain, slight cost
5
4.61
7.90
Gap opening
8
3.95
41.20
Fitting noise
Train RMSE improves at every step and is completely uninformative. CV RMSE has a clear minimum. Prefer the simplest degree within one standard error of the best — here degree 2 wins outright, but when degrees 2 and 3 tie, take 2.
Extrapolation: the failure that bites hardest
A degree-4 polynomial fitted to advertising spend between £0 and £50,000 might fit beautifully. Ask it about £80,000 and it may confidently predict negative sales, because outside the range of the training data a polynomial does whatever its highest-order term tells it to, which is to rocket off to $\pm\infty$.
Linear models extrapolate badly; polynomial models extrapolate catastrophically. The higher the degree, the more violent the behaviour outside the observed range. Guard against it explicitly:
Python
Copy
lo, hi = X_train["spend"].min(), X_train["spend"].max()if not (lo <= new_value <= hi):    raise ValueError(f"{new_value} is outside training range [{lo}, {hi}]")
Refusing to answer is a legitimate and underused option.
Putting it together
Python
WrapCopy
1import pandas as pd2from sklearn.model_selection import train_test_split3from sklearn.metrics import mean_absolute_error, r2_score45df = pd.read_csv("flats.csv")6X = df[["area", "bedrooms", "age_years", "distance_to_station"]]7y = df["price"]89X_train, X_test, y_train, y_test = train_test_split(10    X, y, test_size=0.2, random_state=42)1112linear = make_pipeline(StandardScaler(), LinearRegression()).fit(X_train, y_train)13poly2 = make_pipeline(PolynomialFeatures(2, include_bias=False),14                      StandardScaler(), LinearRegression()).fit(X_train, y_train)1516for name, m in [("linear", linear), ("poly-2", poly2)]:17    pred = m.predict(X_test)18    print(f"{name}: MAE {mean_absolute_error(y_test, pred):,.0f}  "19          f"R2 {r2_score(y_test, pred):.3f}")2021# read the plain linear coefficients, which are on standardised features22coefs = pd.Series(linear[-1].coef_, index=X.columns).sort_values()23print(coefs)
Because those coefficients sit after a scaler, each one reads as "£ per one standard deviation of this feature", which makes them directly comparable. A coefficient of −11,400 on distance_to_station and +32,800 on area says area matters roughly three times as much as station distance across the range these features actually vary over. That statement would be unavailable from raw coefficients.
What this means when you build something
Linear regression earns its place as the first model you fit on any numeric target, not because it usually wins but because of what it tells you cheaply. It trains in milliseconds, gives you a baseline number that any complicated model must beat, and its coefficients and residual plot tell you where the structure it cannot capture is hiding.
Fit it, then look at the residual plot before anything else. A shapeless cloud means a linear model is a defensible final answer and you should think hard before spending a week on gradient boosting for a 2% gain. A clear curve means add polynomial terms — degree 2 first, degree 3 only if cross-validation genuinely rewards it, never degree 8. A widening funnel means model the logarithm of the target instead. Distinct clusters mean there is a categorical feature you have not included.
And when you reach for polynomial features, pair them with regularisation and check the training range before every prediction. Polynomials are the easiest way in all of machine learning to build something that looks perfect on the data you have and produces absurdities on the data you do not.