Ordinary Least Squares

Why Does This Exist?

Linear regression hands you a formula for the weights and moves on. Two questions get skipped: why that formula, and why no serious library computes it that way.

Both answers come out of one picture. A greenhouse logs its water use for three days alongside the sunlight hours each day got, so two columns go into the fit: a column of ones for the intercept, and sunlight hours.

Now the shift that makes everything follow. A list of three numbers is a point in three-dimensional space, so your water readings are one specific point out there, call it $y$. Your two columns are two other points, which you may scale and add, and every combination sweeps out a flat plane through the origin. That plane is the column space of $X$: every prediction your features can ever make, everything else unreachable (vector spaces and rank generalises it).

So your data is a point, your reachable predictions are a plane, and $y$ is not on it. Find the closest point on the plane.

Think of It Like This

How It Actually Works

Write $\hat{y} = Xw$ for the fitted values and $r = y - \hat{y}$ for the residual. Whatever $w$ you pick, $\hat{y}$ lands on the plane. Making $r$ as short as possible means making it perpendicular to the plane, and a vector is perpendicular to a plane exactly when it is perpendicular to everything spanning it, so to every column of $X$:

Those are the normal equations, where normal means perpendicular, nothing to do with the normal distribution. The condition is a list of dot products, one per column, each required to be zero.

Worked example: three days, slope 1.5

Take sunlight hours $x = [0, 1, 2]$ with water use $y = [1, 2, 4]$. Sums are $n = 3$, $\sum x = 3$, $\sum y = 7$, $\sum x^2 = 5$, $\sum xy = 0 + 2 + 8 = 10$. The slope is:

The intercept is $b_0 = (7 - 1.5 \times 3) / 3 = 2.5 / 3 \approx 0.83$. Predictions are $0.83$, $2.33$ and $3.83$. Residuals are $0.17$, $-0.33$ and $0.17$. Add them: $0.17 - 0.33 + 0.17 = 0$. Dot them with $x$: $0 \times 0.17 + 1 \times (-0.33) + 2 \times 0.17 = 0.01 \approx 0$, which is zero up to rounding. Both columns sit square to the leftover, exactly as the diagram shows.

One consequence you can check by hand

The intercept column is all ones, and its dot product with the residual is the sum of the residuals. So residuals sum to exactly zero whenever the model has an intercept. Duplicate columns fall out of the same picture: two identical columns buy no second dimension, so the plane stays a line and $X^\top X$ cannot be inverted. The nearest point is still one point, but the weights reaching it are not unique, decided by matrix rank. Correlated features generalise this headache in multicollinearity.

Squared error is an assumption about noise

Squared error looks like the default. It is a claim.

Suppose each day's water use is the true linear value plus independent noise from a bell curve of the same width every day. The probability of your three readings, as a function of $w$, is a product of those densities, and its logarithm is a sum of squared gaps between observed and predicted. Maximising that log is minimising squared error. Same problem. So squared loss is maximum likelihood under Gaussian noise of constant variance, and heavy tails would make absolute error the likelihood-maximising choice instead. You were never neutral. You were assuming.

Why nobody forms the normal equations

The condition number of a matrix measures how much it amplifies small relative errors, so a large one means digits get eaten during the solve. Forming $X^\top X$ squares it: a design matrix at a condition number of a million becomes a system at a trillion, and double precision carries roughly sixteen digits, so you have spent twelve before the solver starts. QR and SVD work on $X$ directly and never form that product, paying the condition number once rather than twice. That is why library solvers factor instead of using the formula everyone learns, and why ridge regression adds a penalty when the columns nearly coincide.

Code

Both claims at once with plain arithmetic: slope 1.5, residuals summing to zero, residual square to the sunlight column.

1def ols_line(x: list[float], y: list[float]) -> tuple[float, float]:2    n: int = len(x)3    sx: float = sum(x)4    sy: float = sum(y)5    sxx: float = sum(v * v for v in x)6    sxy: float = sum(a * b for a, b in zip(x, y))7    slope: float = (n * sxy - sx * sy) / (n * sxx - sx * sx)8    intercept: float = (sy - slope * sx) / n9    return intercept, slope10
11
12xs: list[float] = [0.0, 1.0, 2.0]13ys: list[float] = [1.0, 2.0, 4.0]14b0, b1 = ols_line(xs, ys)15pred: list[float] = [b0 + b1 * v for v in xs]16resid: list[float] = [b - a for a, b in zip(pred, ys)]17print(f"{b1:.1f} {b0:.2f}")  # -> 1.5 0.8318print(f"{sum(resid):.1f}")  # -> -0.019print(f"{sum(a * b for a, b in zip(xs, resid)):.1f}")  # -> -0.0

The printout matches the hand computation: slope 1.5 with intercept 0.83, residual sum negative zero, dot product with the feature negative zero. Change the last reading from 4 to 10 and rerun: the slope jumps, which shows how one point drags a squared-error fit.

Watch Out For

The Quick Version

• Ordinary least squares projects the target onto the column space of the features, and the residual leaves at a right angle.
• The normal equations say one thing: make the residual perpendicular to every column.
• On readings 1, 2 and 4 at inputs 0, 1 and 2, the slope is 1.5 with intercept about 0.83.
• Those residuals are 0.17, minus 0.33 and 0.17, summing to zero and orthogonal to the feature.
• Minimising squared error assumes Gaussian noise of constant variance.
• Forming the normal equations squares the condition number, so libraries factor instead.
• Gauss-Markov gives best linear unbiased only, which evaporates under unequal or correlated errors.

What to Read Next

• Linear Regression
• Ridge Regression
• Multicollinearity