The Least-Squares Method

1 Learning objectives

2 Motivating example

Suppose you plot the eight bike-share observations from Tutorial 1 and, by eye, draw a line through them. You now have a model. But so does the person sitting next to you, and their line is slightly different, and neither of you can say which one is better — you only disagree about which points "look" more important.

Eyeballing is how regression was done before the 1800s, and it has two problems. It is subjective, so two analysts get two answers. And it does not scale: with one feature you can draw a line on paper, but with twenty features you cannot draw anything at all.

What we need is a rule: a single number that scores any candidate line, so that "best" becomes a well-defined question with a computable answer. Adrien-Marie Legendre proposed such a rule in 1805, and Gauss used it to predict the orbit of Ceres from a handful of observations — a prediction that turned out to be accurate enough to relocate the asteroid. The rule is least squares.

Among all straight lines, choose the one that minimises the sum of the squared vertical distances from the data points to the line.

That is the whole idea. The rest of this tutorial is the consequences of it: why squaring is the right choice, what the resulting formula is, and what guarantees come with it.

3 Concept explanation

3.1 The residual

For a single training example i, the model predicts ŷ_i = w x_i + b. The observed value is y_i. Their difference is the residual:

e_i = y_i − ŷ_i

The sign convention matters and is worth fixing in your memory now:

• Positive residual — the model predicted too low. The point sits above the line.
• Negative residual — the model predicted too high. The point sits below the line.
• Zero residual — the point lies exactly on the line.

3.2 The residual sum of squares

You cannot simply add up all the residuals to get a single measure of badness, because positive and negative residuals cancel. In Tutorial 1 you saw this directly: the mean baseline's residuals summed to exactly zero, even though the baseline was a terrible predictor.

The least-squares criterion solves this by squaring each residual before summing:

RSS(w, b) = Σ_i=1ⁿ (y_i − w x_i − b)²

Read this as a function of the parameters. For any pair (w, b) you choose, RSS tells you how badly that particular line fits the data. The least-squares solution is the pair (ŵ, b̂) that makes RSS as small as possible:

(ŵ, b̂) = argmin_{w, b} RSS(w, b)

3.3 Why square?

Squaring is not the only way to combine residuals into a single number. You could sum the absolute values instead — that gives least absolute deviations, which is also a perfectly respectable method. So why does least squares dominate?

Three reasons push least squares ahead:

1. It is differentiable everywhere. The absolute value function has a sharp corner at zero, which makes it awkward to minimise with calculus. The square is smooth, so you can take derivatives and set them to zero — which is exactly how the closed-form solution is derived.
2. It has a closed-form solution. Because the squared criterion is a quadratic in the parameters, its minimum can be written down directly. Least absolute deviations has no such formula and must be solved iteratively.
3. It is the maximum-likelihood solution under Gaussian noise. If you assume the true relationship is linear and the deviations from it are independent draws from a normal distribution with constant variance, then maximising the likelihood of the observed data is exactly the same as minimising RSS. This is a strong justification, and it comes with a whole statistical toolkit for free.

3.4 Deriving the solution — the simplest case

Start with a model that has no intercept, so the line is forced through the origin: ŷ = w x. This case is contrived but the algebra is short, and the result generalises.

The criterion is:

RSS(w) = Σ (y_i − w x_i)²

Expand the square and differentiate with respect to w:

RSS(w)  = Σ (yᵢ² − 2 yᵢ w xᵢ + w² xᵢ²)

dRSS
────  = Σ (−2 yᵢ xᵢ + 2 w xᵢ²)
 dw
      = −2 Σ xᵢ yᵢ  +  2 w Σ xᵢ²

Set the derivative to zero. The factor of 2 cancels, and the condition becomes:

Σ x_i y_i = w Σ x_i²

So the least-squares slope is:

ŵ = Σ x_i y_i / Σ x_i²

The second derivative is 2 Σ x_i², which is positive unless every x_i is zero. So this stationary point is a minimum, not a maximum or a saddle — the RSS curve is a parabola opening upward, and there is exactly one bottom.

3.5 Deriving the solution — with an intercept

Now the real case: ŷ = w x + b. RSS is now a function of two parameters, so we need two partial derivatives, and both must be zero at the minimum.

RSS(w, b) = Σ (y_i − w x_i − b)²

Partial derivative with respect to b:

∂RSS
────  = Σ −2 (yᵢ − w xᵢ − b)
 ∂b
      = −2 ( Σ yᵢ − w Σ xᵢ − n b )   =  0

Dividing by −2 and rearranging:

Σ y_i − w Σ x_i − n b = 0

Dividing through by n and using the notation x̄ = Σx_i/n and ȳ = Σy_i/n:

b̂ = ȳ − ŵ x̄

This is the first normal equation, and it has a memorable consequence: the fitted line always passes through the point (x̄, ȳ), the centroid of the data.

Partial derivative with respect to w:

∂RSS
────  = Σ −2 xᵢ (yᵢ − w xᵢ − b)
 ∂w
      = −2 ( Σ xᵢ yᵢ − w Σ xᵢ² − b Σ xᵢ )   =  0

Dividing by −2:

Σ x_i y_i − w Σ x_i² − b Σ x_i = 0

This is the second normal equation. Now substitute the first into the second to eliminate b:

Σ xᵢ yᵢ  −  w Σ xᵢ²  −  ( ȳ − w x̄ ) Σ xᵢ   =  0

Σ xᵢ yᵢ  −  w Σ xᵢ²  −  n x̄ ȳ  +  w n x̄²   =  0

( Σ xᵢ yᵢ − n x̄ ȳ )   =   w ( Σ xᵢ² − n x̄² )

The two grouped quantities have standard names:

S_xy = Σ x_i y_i − n x̄ ȳ      S_xx = Σ x_i² − n x̄²

S_xy is the sum of cross-products; it measures how much x and y move together. S_xx is the sum of squared deviations of x; it measures how spread out the feature is. Both are directly computable from the raw sums, which is how statistical software has computed regression since long before computers.

The solution is therefore:

3.6 Four properties of the least-squares fit

The normal equations are not just a computational trick — they imply four properties that any least-squares fit must satisfy. They are also the fastest way to check whether a computer has done the arithmetic correctly.

The second property is the one that generalises most usefully. In the matrix form of the next section it becomes Xᵀe = 0 — the residual vector is orthogonal to every column of the design matrix. That single equation is the geometric heart of least squares.

3.7 The matrix form

With more than one feature, writing out the algebra feature by feature becomes unmanageable. The matrix form handles any number of features at once.

First, absorb the intercept into the weight vector by prepending a column of ones to the design matrix. Define:

      ⎡ 1  x₁₁  x₁₂  ⋯  x₁d ⎤              ⎡ w₀ ⎤            ⎡ y₁ ⎤
X  =  ⎢ 1  x₂₁  x₂₂  ⋯  x₂d ⎥  ,   w  =  ⎢ w₁ ⎥  ,   y  =  ⎢ y₂ ⎥
      ⎢ ⋮   ⋮    ⋮   ⋱   ⋮  ⎥              ⎢ ⋮  ⎥            ⎢ ⋮  ⎥
      ⎣ 1  xₙ₁  xₙ₂  ⋯  xₙd ⎦              ⎣ wd ⎦            ⎣ yₙ ⎦

The first column of ones multiplies w₀, which plays the role of the intercept. The model for all n examples at once is ŷ = Xw, and the criterion becomes:

RSS(w) = ‖y − Xw‖² = (y − Xw)ᵀ(y − Xw)

Expanding and differentiating with respect to the vector w (using the standard vector-calculus identities) gives the gradient:

∇RSS = −2 Xᵀ(y − Xw)

Setting this to zero:

Xᵀ ( y − X w )  =  0

Xᵀ y  −  Xᵀ X w  =  0

Xᵀ X w  =  Xᵀ y          ← the normal equations

Multiplying both sides by the inverse of XᵀX:

This is the single most important formula in this tutorial. Everything in the preceding sections — the centroid property, the two normal equations, the residual orthogonality — is a special case of it. And it is what LinearRegression().fit(X, y) is Computing when you call it, though scikit-learn uses a numerically more stable decomposition rather than forming the inverse explicitly.

3.8 The geometric interpretation

There is a way of seeing least squares that makes the orthogonality property obvious rather than algebraic.

Think of the target vector y as a point in n-dimensional space — one dimension per training example. The columns of X span a d+1-dimensional subspace of that space: every possible prediction Xw that the model could ever make lies in that subspace.

Least squares asks: which point in the subspace is closest to y? The answer, from elementary geometry, is the orthogonal projection of y onto the subspace. And the defining property of an orthogonal projection is that the difference between the original vector and its projection is perpendicular to the subspace:

Xᵀ(y − Xŵ) = 0

which is precisely the normal equations. The residual vector is orthogonal to every column of the design matrix. What looked like an algebraic convenience is actually a geometric necessity.

4 Visual illustration

5 Worked example

We return to the bike-share data from Tutorial 1 and fit a least-squares line by hand, then compare it against the mean baseline.

n   = 8
x̄   = 155 / 8   = 19.375
ȳ   = 2700 / 8  = 337.5

Sxy = Σ xᵢ yᵢ  −  n · x̄ · ȳ
    = 60,570   −  8 × 19.375 × 337.5
    = 60,570   −  52,312.5
    = 8,257.5

Sxx = Σ xᵢ²    −  n · x̄²
    = 3,403    −  8 × 19.375²
    = 3,403    −  8 × 375.390625
    = 3,403    −  3,003.125
    = 399.875

ŵ = Sxy / Sxx
   = 8,257.5 / 399.875
   ≈ 20.6502

b̂ = ȳ − ŵ · x̄
   = 337.5 − 20.6502 × 19.375
   = 337.5 − 400.0976
   ≈ −62.5976

RSS  = 831.19
MSE  = 831.19 / 8   = 103.90
RMSE = √103.90      ≈ 10.19 rentals

Baseline (from Tutorial 1):
RSS  = 171,350
RMSE = 146.3 rentals

Improvement in RMSE:  146.3 → 10.19   (a 93% reduction)

R² = 1 − RSS / TSS
   = 1 − 831.19 / 171,350
   ≈ 0.9951

The same computation in code

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

df = pd.DataFrame({
    "temp_c":  [ 8, 12, 15, 18, 21, 24, 27, 30],
    "rentals": [120, 180, 240, 300, 360, 430, 510, 560],
})

x = df["temp_c"].to_numpy()
y = df["rentals"].to_numpy()
n = len(x)

# --- from scratch, using the normal equations ---
Sxy = (x * y).sum() - n * x.mean() * y.mean()
Sxx = (x ** 2).sum() - n * x.mean() ** 2

w = Sxy / Sxx
b = y.mean() - w * x.mean()

print(f"from scratch : w = {w:.4f}, b = {b:.4f}")

# --- the same thing via scikit-learn ---
model = LinearRegression().fit(x.reshape(-1, 1), y)

print(f"scikit-learn : w = {model.coef_[0]:.4f}, b = {model.intercept_:.4f}")

# --- residuals and metrics ---
y_hat = w * x + b
resid = y - y_hat

rss  = (resid ** 2).sum()
rmse = np.sqrt(rss / n)

print(f"residuals sum: {resid.sum():.6f}")
print(f"RSS          : {rss:.2f}")
print(f"RMSE         : {rmse:.2f}")

from scratch : w = 20.6502, b = -62.5976
scikit-learn : w = 20.6502, b = -62.5976
residuals sum: 0.000000
RSS          : 831.19
RMSE         : 10.19

The two approaches agree to four decimal places, as they must. scikit-learn computes the answer via a QR or SVD decomposition rather than forming (XᵀX)⁻¹ explicitly, because the decomposition is numerically more stable — but the mathematical answer is identical.

6 Interactive activity

7 Self-check questions

Answer each question in your own words before revealing the answer.

8 Practice exercises

Work these by hand first, then verify in the notebook. Use the formulas ŵ = S_xy/S_xx and b̂ = ȳ − ŵx̄.

x = [1, 2, 3, 4]
y = [2, 5, 7, 10]

x = [1, 2, 3, 4, 5]
y = [3, 5, 6, 9, 12]

      ⎡ 1  1 ⎤        ⎡  3 ⎤
      ⎢ 1  2 ⎥        ⎢  5 ⎥
X  =  ⎢ 1  3 ⎥ ,  y = ⎢  6 ⎥
      ⎢ 1  4 ⎥        ⎢  9 ⎥
      ⎣ 1  5 ⎦        ⎣ 12 ⎦

              ⎡ n     Σxᵢ  ⎤   ⎡  5  15 ⎤
XᵀX  =        ⎢            ⎥ = ⎢        ⎥
              ⎣ Σxᵢ  Σxᵢ²  ⎦   ⎣ 15  55 ⎦

              ⎡  Σyᵢ  ⎤   ⎡ 35 ⎤
Xᵀy  =        ⎢       ⎥ = ⎢    ⎥
              ⎣ Σxᵢyᵢ ⎦   ⎣127 ⎦

5b  +  15w  =  35
15b +  55w  =  127

9 Problem set

These problems are longer and more open-ended. They feed directly into Assignment 2.

from sklearn.datasets import fetch_california_housing
data = fetch_california_housing(as_frame=True)
df = data.frame

10 Discussion and projects

11 Summary

Key terms

Residual

The difference between an observed value and its fitted value: e_i = y_i − ŷ_i.

Residual sum of squares (RSS)

The sum of the squared residuals across all training examples. The quantity least squares minimises.

Least squares

The fitting principle that chooses parameters to minimise the sum of squared residuals.

Least absolute deviations

An alternative criterion that minimises the sum of absolute residuals. More robust to outliers, but has no closed-form solution.

Sum of cross-products (S_xy)

Σx_iy_i − n x̄ ȳ. Measures how much the feature and target vary together.

Sum of squared deviations (S_xx)

Σx_i² − n x̄². Measures the spread of the feature values.

Normal equations

The system XᵀXw = Xᵀy obtained by setting the gradient of RSS to zero. Their solution is the least-squares estimate.

Design matrix (X)

The n × (d+1) matrix of feature values, with a leading column of ones to absorb the intercept.

Centroid

The point (x̄, ȳ). Every least-squares line with an intercept passes exactly through it.

Orthogonal projection

The closest point in a subspace to a given vector, characterised by the residual being perpendicular to the subspace.

Collinearity

A situation in which one feature is an exact or near-exact linear combination of others, making XᵀX singular or ill-conditioned.

12 Further reading