Linear Regression

Linear regression is the oldest algorithm in this course. Legendre published least squares in 1805 to fit comet orbits, in an appendix titled Sur la Méthode des moindres quarrés, and Gauss claimed he had been using it since 1795 — but neither of them called it regression. That word arrived eighty years later, from a study of human height, and it described a phenomenon rather than a method. Two centuries on, the model remains the first thing to try on any continuous target: fast, interpretable, and a baseline that is embarrassingly hard to beat.

This lesson takes one idea and follows it all the way down: fit a line, look at what it got wrong, and use those errors to do better. The errors are the residuals; using them systematically is the gradient.

Where the name comes from

In 1886 Francis Galton published Regression towards Mediocrity in Hereditary Stature. He had collected the heights of 928 adult children and their parents, and he was looking for the strength of heredity. What he found puzzled him.

Very tall parents did have tall children — but on average, the children were closer to the population mean than their parents were. Very short parents had short children who were, again, closer to the middle. Galton first called this "reversion", then "regression towards mediocrity". The name stuck to the technique he used to measure it, which is why an algorithm about fitting lines carries a word that means "moving backwards".

The left panel is the effect itself. The dashed line is identity — what the data would look like if children simply matched their parents. The fitted line is flatter. A mid-parent 3 inches above average produces a child predicted only about 2 inches above average. The red arrow is the gap, and that gap is the "regression".

The slope Galton reported was about \(2/3\), and it is worth seeing where that comes from:

Correlation and slope are not the same number. The slope carries the units; the correlation is what is left after you strip them out. They coincide only when both variables are measured in the same spread — which is exactly the situation the simulator below puts you in, so you can watch the geometry without the unit bookkeeping.

One figure in that calculation deserves a second look: why do children (2.52) vary more than mid-parents (1.80)? Not because parents are a more uniform group. Mid-parent is the average of two people, and averaging shrinks spread: two individuals at \(s \approx 2.5\) average to \(2.5/\sqrt{2} \approx 1.77\), which is essentially the 1.80 Galton reported. The asymmetry is an artefact of how the \(x\) variable was constructed, not a fact about heredity — and noticing that is the same habit of mind the rest of this section is about.

The model

Predict a continuous target as a weighted sum of features:

• \(w_j\) is the change in \(\hat{y}\) per unit change in \(x_j\), holding the other features fixed;
• \(w_0\) (intercept, or bias) is the prediction when every feature is zero.

In matrix form, with a leading column of 1s absorbed into \(X\): \(\hat{y} = Xw\).

Why not just invert \(X\)?

Written as \(\hat{y} = Xw\), the natural reflex is \(w = X^{-1}y\). It does not work, and seeing exactly why is the fastest route into what least squares is actually doing.

Take three observations and two parameters:

First obstacle: \(X\) is not square. It is \(3 \times 2\), and only square matrices have inverses. In general \(X\) is \(n \times d\) with \(n \gg d\) — thousands of rows, a handful of columns.

Second obstacle, and the real one: there is no exact solution anyway. \(y = Xw\) is three equations in two unknowns — an overdetermined system. An exact \(w\) exists only if \(y\) happens to lie in the column space of \(X\), and here it does not:

Appending \(y\) raised the rank, so \(y\) points somewhere the columns of \(X\) cannot reach. No \(w\) whatsoever satisfies \(Xw = y\). That is not bad luck — it is the normal situation, and it is the whole reason regression exists. If an exact solution did exist, every residual would be zero and the line would pass through every point.

What we do instead

Since \(y\) is out of reach, settle for the reachable point closest to it — minimise \(\lVert y - Xw\rVert^2\). Setting the gradient to zero gives the normal equations, and the trick is what happens to the shapes:

\(X^\top X\) is \(2 \times 2\) — square, and invertible whenever \(X\) has full column rank. We never invert \(X\); we invert \(X^\top X\). Two equations, small enough to solve on paper:

The fitted values are \(1.5,\ 4,\ 6.5\) and the residuals \(0.5,\ -1,\ 0.5\). They do not vanish — they cannot — but \(\lVert e \rVert^2 = 1.5\) is the smallest achievable. And note what the normal equations say in this notation: \(X^\top e = 0\), which is exactly the \(\sum e_i = 0\) and \(\sum e_i x_i = 0\) you will meet again below.

Your instinct is the special case

Now drop the third observation, leaving two equations and two unknowns:

Here \(X^{-1}\) exists and \(w = X^{-1}y\) works perfectly — the residual is exactly zero, because a line through two points interpolates them. And the least-squares formula returns the very same answer, as the algebra requires:

So \(w = X^{-1}y\) is not wrong — it is the exactly-determined special case of a more general formula. That general formula, \(X^{+} = (X^\top X)^{-1}X^\top\), is called the Moore–Penrose pseudo-inverse: what "inverse" becomes when the matrix is not square.

Those four words — holding the other features fixed — carry more weight than they look. Read this fitted model for apartment prices:

3200 does not mean "bigger apartments cost R$ 3200 more per m² in this city". It means: among apartments of the same age, one extra m² is associated with R$ 3200 more. If area and age happen to be correlated in your data — say newer buildings are larger — then the raw relationship between area and price mixes both effects, and the coefficient deliberately separates them.

Least squares

We need a rule for choosing \(w\). Ordinary least squares — OLS from here on — picks the line that minimises the sum of squared residuals (SSE, the sum of squared errors): the squared vertical distances between the data and the fitted line.

Why vertical distances, and why squared? Both choices are worth a moment, because both could have gone differently.

Vertical, because the model's job is to predict \(y\) from \(x\). An error is "how wrong was my prediction of \(y\)", which is measured along the \(y\) axis. (Minimising perpendicular distance is a different and perfectly good method — it is what PCA does — but it answers a different question, and it is why the two regression lines in Galton's figure differ.)

Squared, for three reasons that happen to align:

1. it is smooth and differentiable everywhere, so calculus works and there is a closed-form answer;
2. it punishes one large error more than several small ones, which is often what you want;
3. under Gaussian noise it is maximum likelihood — the least-squares line is the most probable line.

Reason 2 is also its weakness: squaring makes the fit hypersensitive to outliers. Minimising \(\sum|e_i|\) instead gives the more robust least absolute deviations, at the price of losing the closed form.

Drag the points below and watch the line chase the minimum. Then drag one point far from the rest and see how much a single observation can move it:

Solving it exactly

Set the gradient to zero, \(\nabla_w J = -2X^\top(y - Xw) = 0\), and you get the normal equations:

For simple regression — one feature — this collapses to two formulas worth memorising:

The second one says something useful: the fitted line always passes through \((\bar{x}, \bar{y})\). Whatever else it does, it pivots about the centre of mass of the data.

Worked by hand

Five points, small enough to check every step:

With \(\bar{x}=3\) and \(\bar{y}=4.2\):

Now the residuals \(e_i = y_i - \hat{y}_i\):

Two things are true of these numbers and are true of every OLS fit with an intercept:

Check them: \(-0.6+0.6+0.8-1.0+0.2 = 0\). ✓

These are not coincidences, they are the normal equations — one per parameter. The residuals are forced to be orthogonal to the intercept column and to every feature column. That has a geometric reading: \(\hat{y}\) is the projection of \(y\) onto the space spanned by the columns of \(X\), and the residual vector is what sticks out perpendicular to that space. Least squares is a right-angle problem in disguise.

from sklearn.linear_model import LinearRegression

model = LinearRegression().fit(X_train, y_train)
model.coef_, model.intercept_
y_pred = model.predict(X_test)

Residuals

Everything the model failed to explain is in the residuals, and reading them is the single most useful diagnostic habit in regression.

The residual plot puts \(e_i\) on the vertical axis against the fitted \(\hat{y}_i\) on the horizontal. A healthy one looks like nothing at all: a structureless band around zero. Any pattern is the model telling you which assumption it just broke.

Try to diagnose these four before reading the verdict — the top panel is what you would normally look at, and the point is that it is not enough:

Outlier and leverage are different things. An outlier has a large residual — the model missed it. A high-leverage point sits far from \(\bar{x}\) along the \(x\) axis, and it can drag the whole line towards itself. The dangerous case is a point with high leverage and a pull: it moves the line so effectively that its own residual ends up small, hiding the damage. That is the fourth case in the simulator. This is why "just drop the big residuals" is bad advice — the worst offenders do not have big residuals.

Measuring quality

Every metric below is a different way of summarising the same vector of residuals \(e_i = y_i - \hat{y}_i\) into a single number. They disagree because they compress differently, so the names are worth reading in full before the formulas.

The acronyms are not arbitrary: each one spells out its own recipe, read right to left. RMSE is Error → Squared → Mean → Root: take each error, square it, average them, then take the square root. Read the name backwards and you have written the code.

Computing all four on the five points

The residuals from the worked example were \(-0.6,\ 0.6,\ 0.8,\ -1.0,\ 0.2\).

Read these back in context: the target values ranged from 2 to 6, so being off by roughly 0.7 on average is a meaningful error, not a rounding detail. MAE and RMSE are quotable to a non-technical audience — "the model misses by about 0.7" — while MSE is 0.48 squared units, which is not a sentence anyone can say out loud. That is the whole reason RMSE exists: it undoes the squaring so the number returns to the scale of the thing you are predicting.

MAE or RMSE?

They are not interchangeable, and the difference is exactly how they treat one big mistake. Two models on five predictions:

Identical MAE, and RMSE more than twice as large for B. MAE says the two models are equally good; RMSE says B is much worse, because it squares before averaging and the single error of 10 contributes 100 to the sum. Neither is right in the abstract — the question is what your problem costs. Delivery time being off by 10 minutes once is usually worse than being off by 2 minutes five times (use RMSE); for a demand forecast the total error may be what matters (use MAE). Choose before you look at the results, not after.

R², and why it is called that

R² compares your model against the laziest possible one. It is built from two sums: SSE, the sum of squared errors your line still makes, and SST, the total sum of squares, which is the error of ignoring \(x\) altogether.

SST, the total sum of squares, is the error you would make by ignoring \(x\) entirely and always predicting the mean \(\bar{y}\). SSE, the error sum of squares, is what your line still gets wrong. So R² answers: of the variation that was there to explain, what fraction did the model remove? Hence "coefficient of determination" — how much of \(y\) is determined by the features.

On the five points: \(R^2 = 1 - 2.40/8.80 = 0.727\), so the line removed about 73% of the original spread. And since \(r_{xy} = 0.853\) and \(0.853^2 = 0.727\), notice that in simple regression R² really is the squared correlation — that is where the symbol comes from. It stops being true the moment you add a second feature.

R² can be negative on test data, which surprises people who learned it as "a percentage". It means the model did worse than the flat line \(\hat{y}=\bar{y}\). Always fit that constant baseline (DummyRegressor) and report it: it is uncomfortable how often an elaborate model barely clears it.

Assumptions behind the inferences

OLS predictions ask for very little. Interpreting coefficients, confidence intervals and p-values leans on the classical assumptions:

1. Linearity — the true relationship is approximately linear in the features;
2. Independence — residuals are uncorrelated with each other (beware time series);
3. Homoscedasticity — residual variance is constant across the range of \(\hat{y}\);
4. Normality of residuals — needed for exact intervals and p-values, not for the fit itself;
5. No severe multicollinearity — highly correlated features make individual coefficients unstable.

Multicollinearity deserves a sentence of its own, because its symptom is counter-intuitive. With two features correlated at \(r = 0.99\), the data can barely tell their coefficients apart: \(X^\top X\) is nearly singular, so if the truth is \((w_1, w_2) = (1, 1)\), the pair \((2, 0)\) fits within about 2% of it and \((3, -1)\) within about 8%. Refit on a fresh sample and the individual estimates move by a whole unit or more — enough to flip a sign — while their sum stays pinned near 2, and the predictions with it. The model is fine; the interpretation of either coefficient alone is worthless.

From the exact answer to the gradient

We have a formula that solves the problem in one shot. So why does the rest of machine learning bother with iteration?

Because \((X^\top X)^{-1}\) stops being available. It costs \(O(d^3)\), so a million features is hopeless. It does not exist when columns are collinear. And it has no equivalent at all for the models in the second half of this course — a neural network has no normal equations. What does generalise is the idea underneath.

Look again at how we got the closed form: we wrote down \(J(w)\), took its gradient, and asked where it vanishes. Gradient descent keeps the first two steps and gives up on the third. Instead of solving \(\nabla J = 0\) algebraically, it walks downhill until it gets there.

For simple regression, write the loss as a mean so the scale does not depend on \(n\):

Differentiate with respect to each parameter, using the chain rule and remembering \(e_i = y_i - w_0 - w_1x_i\):

Two things fall out of this, and both are worth pausing on.

First, set both to zero and you recover exactly the two identities the worked example satisfied: \(\sum e_i = 0\) and \(\sum e_i x_i = 0\). The closed form and the gradient are the same statement, approached from two directions.

Second, look at what the gradient is made of: it is a sum weighted by the residuals. Points the model already predicts well contribute \(e_i \approx 0\) and barely vote. Points it gets badly wrong dominate the step. The update rule

therefore reads, in plain language: move each parameter in the direction that the mistakes are pointing, by an amount proportional to how big the mistakes are. This is the sentence that will still be true, unchanged, when the model is a 100-layer network and \(\nabla J\) comes from backpropagation.

The picture below makes the geometry concrete. On the left is \(J(w_0, w_1)\) as a landscape — every point is a candidate line, the colour is its error, and the white cross is the closed-form optimum. Click to drop a starting guess and step downhill:

Now switch centre x on and run it again from the same kind of starting point. Same data, same algorithm, dramatically different behaviour.

Uncentred, the contours form a long narrow valley: \(w_0\) and \(w_1\) are strongly coupled, because moving the slope of a line whose \(x\) values are all far from zero also swings its height. The gradient points across the valley rather than along it, so the path zigzags and creeps. Centred, the valley becomes a bowl and the same steps go almost straight to the answer.

Choosing \(\eta\), doing this on batches rather than the whole dataset, and adding penalties to \(J\) that shrink the weights, are the subject of the next lesson: Gradient Descent & Regularization.

Class materials

References

• Legendre, A.-M. "Nouvelles méthodes pour la détermination des orbites des comètes." Courcier (1805) — apêndice Sur la Méthode des moindres quarrés**. full text
• Gauss, C. F. "Theoria motus corporum coelestium in sectionibus conicis solem ambientium." Perthes & Besser (1809). full text
• Galton, F. "Regression towards Mediocrity in Hereditary Stature." J. Anthropological Institute 15 (1886). full text
• Hanley, J. A. ""Transmuting" Women into Men: Galton's Family Data on Human Stature." The American Statistician 58 (2004). DOI
• Anscombe, F. J. "Graphs in Statistical Analysis." The American Statistician 27 (1973). DOI
• Cook, R. D. "Detection of Influential Observation in Linear Regression." Technometrics 19 (1977). DOI

The full course bibliography is on the references page.

Quiz