Linear regression

From 1941 to 2025 Almaty had 31,046 days, and for each the mean temperature is known. Is the city getting warmer, and if so, how fast? By eye nothing can be made out in this cloud: January differs from July by almost thirty degrees, while the quantity sought is tenths of a degree per decade. Linear regression answers with one system of six equations. Below it is solved right in the browser: switch the model’s columns on and off, drag the edges of the window of years and watch the answer change.

Switch off the yearly wave. The model “mean plus trend” explains a thousandth of the spread, R² = 0.001, and its typical miss is the same ten degrees as with the mean alone. Switch the wave back on: R² jumps to 0.85, and the miss falls to 3.9 degrees. Yet the trend barely changes: +0.13 and +0.12 degrees per ten years. The seasons and the warming do not get in each other’s way, and regression finds both at once. Now press “from 1976”: the trend becomes four times steeper. Why, and which number to believe, we will work out below.

The idea

The model is a weighted sum of columns: the prediction for row i is β₀ + β₁x_i1 + … + β_d−1x_i,d−1. The weights are chosen so that the sum of squared misses over all rows is smallest. At that point the misses are not linearly related to any column, and from this comes a system of d equations in d unknowns, however many rows the data has.

Let us start with a line y ≈ a + bx. The miss at a point is the residual e_i = y_i − a − bx_i. Residuals have different signs, so their squares are added: such a sum is smooth, it has derivatives, and its minimum is found by a formula. This is the method of least squares. Setting the derivatives with respect to a and b to zero, we get

The line passes through the centre of mass of the cloud (x̄, ȳ), and its slope is how consistently x and y deviate from their means, divided by the spread of x. A derivation through projections and a widget where the squared residuals are drawn literally are in “Mathematics, the Queen of the Sciences”, in the experiment on the best line. Squares have a downside too, also visible there: a miss of ten units weighs as much as a hundred misses of one, and a single gross error in the data turns the line.

Any number of columns

Put the data into a table X: n rows, d columns, the first column all ones (it gives the intercept β₀). The predictions for all rows are the vector Xβ, the residuals y − Xβ. The derivative of the sum of squared residuals with respect to the weight β_j is −2Σ_i x_ij(y_i − ŷ_i). At the minimum all d derivatives are zero, that is, the vector of residuals is perpendicular to every column:

These are the normal equations; “normal” comes from that very perpendicular. The matrix X^TX is square, d × d: its entry (j, k) is the sum of x_ijx_ik over all rows. A million rows or thirty-one thousand, the system is still six by six.

If the columns are dependent, say the temperature is recorded twice, in degrees Celsius and Fahrenheit, there are infinitely many minima: one column can be replaced by the other. Gaussian elimination will then hit a zero pivot. NumPy in this case returns the shortest of all the solutions, but it is wiser not to feed the same thing into the model twice.

“Linear” in the coefficients

Temperature over a year is not a line but a wave. How can it be fitted with linear regression? Any sinusoid with a period of a year splits into a sine and a cosine: A sin(ωt + φ) = A cos φ · sin ωt + A sin φ · cos ωt. The amplitude and shift enter the formula nonlinearly, but the coefficients of sin ωt and cos ωt are ordinary numbers that multiply columns. So we put two columns into the table, sin ωt and cos ωt with ω = 2π / 365.2425 days, and least squares will find both the size of the wave and the day it peaks. Two more columns at double frequency, a half-year wave, let the yearly course be asymmetric: warming up more slowly than cooling down, or the other way round.

The word “linear” in the name refers to the coefficients: the prediction is a linear combination of columns. The columns themselves can be any functions of the data: sines, squares, logarithms. Here is what the models from the scene give, from one column to six:

columns                      trend °C/10y  wave °C      R²  miss °C
mean                                    -        -  0.0000    10.13
mean + trend                       +0.130        -  0.0010    10.12
mean + wave                             -    13.21  0.8509     3.91
mean + trend + wave                +0.121    13.21  0.8518     3.90
+ half-year wave                   +0.122    13.21  0.8553     3.85

Here R² is the coefficient of determination: the share of the spread around the mean that the model explains, and the miss is the root of the mean squared residual. The yearly wave in Almaty is ±13.2 degrees around a mean of 8.76. With the half-year correction the coldest day of the model year is 8 January, the warmest 24 July, with 26.6 degrees between them. The model temperature rises for 197 days and falls for 168: the city cools down faster than it warms up.

Why did the trend hardly notice the wave? Over a whole number of years the sine and cosine add up to almost zero, and the trend column is almost unrelated to them: X^TX is close to block form, and the equations for the trend and for the wave hardly engage each other. Were the columns pairwise perpendicular, each coefficient would be computed separately, as if the other columns did not exist. For dependent columns this is not so: add one and all the others change.

Code

The first listing is everything needed in general form: a line by the formula with means, Gaussian elimination with pivoting and a fit with any number of columns through the normal equations. We check them on four sets of points that will come up again: the formula and the system give the same thing.

The second listing reads all 31,046 days of Almaty. Note what it keeps: not the rows but, for each year, the matrix X^TX, the vector X^Ty, the sum of squares of y and the number of days. That is enough for everything else. A window of years is the sum of the matrices of its years. A model with some of the columns is a submatrix: dropping a column from X means dropping a row and a column from X^TX. Even R² needs no second pass: at the minimum the sum of squared residuals equals Σy² − β^TX^Ty; expand |y − Xβ|² and substitute the normal equations.

The scene at the start of the article does exactly this: one pass over the days when the page loads, and after that every movement of a window edge is the addition of a few dozen 6 × 6 matrices and the solution of one system.

What it costs

Each row adds all d² pairwise products of its numbers to X^TX (by symmetry half would do), and d more to X^Ty. The pass over the data costs O(n·d²), solving the system by Gaussian elimination O(d³), and the memory is O(d²): the rows can be read from a file or the network and forgotten at once. For Almaty that is about a million multiplications and a six-by-six system.

On logarithmic axes the lines are almost parallel, with a slope close to one: the time grows in proportion to the number of rows. The distance between the lines is the price of columns: with three times as many columns there are almost nine times as many multiplications per row, while the measured time grew about sevenfold; part of it goes on the loop over rows itself.

Since the sums add up, they can be stored as prefix sums: then the slope over any window takes O(d²) operations, without a single pass over the data. The article’s second problem is built on this.

Where does the formula stop working? If there are hundreds of thousands of columns, words in texts, pixels, features of ad impressions, the matrix X^TX does not fit in memory, and d³ takes forever. Then the minimum is sought by gradient descent: step by step, on random batches of rows, without forming the system at all. Neural networks are trained this way too; the first training run in the course chapter on machine learning is the same line, found by descent.

How far to trust the slope

Here is the same trend over different windows of years. The first number is over the days, with the model with all six columns, as in the scene. The second is a line through the 85 (or fewer) yearly means and its standard error:

1941-2025: +0.12 °C per 10 years; through the yearly means +0.12 ± 0.04
1941-1975: -0.72 °C per 10 years; through the yearly means -0.72 ± 0.14
1976-2025: +0.49 °C per 10 years; through the yearly means +0.49 ± 0.07
2000-2025: +0.98 °C per 10 years; through the yearly means +0.98 ± 0.17

The standard error of the slope, √(Σe² / (m − 2) / Σ(x − x̄)²) for m points, says how much the slope would wander if the same years were lived again with different weather. The formula assumes the residuals of different points are independent. For years this is close to the truth, for days it is not: a cold week is seven residuals that resemble each other, not seven independent ones. Computed over the days, the error would come out many times smaller and would mislead. So it is computed over the yearly means.

Four windows, four different answers, and each is honest for its window. Since 1976 the trend is +0.49 ± 0.07 degrees per ten years, since 2000 about a degree. And for 1941–1975 it is minus 0.72: in this data the warmest year of the whole series is 1941, at 11.72 °C, warmer than any of the recent ones. Points at the edges of a window pull the line hardest: in the slope formula each point enters with the factor x_i − x̄, and at the edge it is largest. That is how the warm 1940s bring the trend over all 85 years down to +0.12.

Whether the warm 1940s are right, regression does not know; it draws a line through what it was given. The data here is the ERA5 reanalysis: a model of the atmosphere into which observations are mixed. According to ECMWF, in 1940 there were about 17 thousand of them a day for the whole planet, in 2022 25 million, and the early decades of the reanalysis are therefore less reliable than the later ones. Hence the rule: a number of “so many degrees per decade” without the window and the data source means little, and the window is best chosen before looking at the answer.

And one last thing. Point the magnifier at the years beyond the edge of the window: there the model curve is dashed. That is already a forecast, the formula continued where nobody has checked it. A line that describes 1976–2025 well promises nothing about 2100: a linear model is a description of the data in the window, not a law of nature.

What a single number hides

In 1973 the statistician Frank Anscombe published four sets of eleven points. All four have the same means and spreads of x and y, the same line y = 3 + 0.5x and R² = 0.67; the first listing prints exactly that. Take a look at them.

The first set is an ordinary cloud around a line. The second lies on a smooth arc, and a line is the wrong model for it; switch to residuals and they line up in an arch. The cure is the same trick as with the seasons: add a column x², and R² becomes 0.999999 (the last line of the listing). In the third set ten points lie exactly on a line, and the eleventh has skewed it: click on it, and the line falls into place. In the fourth the whole line is set by one point: the ten others all have the same x = 8. Remove it and the slope is not defined at all: the formula’s denominator becomes zero.

Hence a habit cheaper than any theory: after fitting, look at the residuals. If a pattern is left in them, an arc, a wave, a funnel, the model has not caught it. In the magnifier of the first scene switch to “residuals”: the yearly wave is no longer in them, but something else is visible: in winter the points are scattered more widely than in summer. Winter weather in Almaty is more changeable than summer weather, and a model with a constant error does not take that into account.

History

The method of least squares was first published by Adrien-Marie Legendre in 1805, in an appendix to a book on computing the orbits of comets. In 1809 Carl Friedrich Gauss set it out in more detail, linked it with the law of distribution of errors and remarked that he had been using it since 1795; the dispute over priority was never settled. “Mathematics, the Queen of the Sciences” tells this story in detail.

The word “regression” came later and from another science. In 1886 Francis Galton published the paper “Regression towards Mediocrity in Hereditary Stature”. He collected the heights of 928 adult children and their parents, 205 couples; to compare women’s heights with men’s, he multiplied them by 1.08. It turned out that a child’s deviation from the average height is on average about two thirds of the deviation of the “mid-parent”: children of tall parents are on average tall too, but shorter than their parents, and children of short parents taller. Galton called this regression towards mediocrity.

Today it is clear that this is a property not of heredity but of any two quantities that are imperfectly related. The slope of the line “children on parents” equals the correlation coefficient times the ratio of the spread of children’s heights to that of parents’ heights, and with equal spreads it is less than one unless the relation is perfect. Fit the line “parents on children” and it too comes out shallow: the parents of very tall children are on average shorter than their children. The word stuck and now means any fitting of a relationship to data, and linear regression the kind that is linear in its coefficients.

Problems

In the first problem the model is a line through the origin: the formulas from the article do not fit it, but the derivation is the same. In the second there are thousands of windows, and going over the data for each is too slow.