Statistics
Least Squares: The Line of Best Fit
Least Squares is the recipe that finds the single straight line ŷ = a + bx passing through a cloud of data points so that the total of the squared vertical gaps between the points and the line is as small as it can possibly be. Discovered independently by Legendre (1805) and Gauss (who claimed use from 1795), it is the workhorse of statistics: one closed-form pair of formulas — the normal equations — hands you the optimal slope and intercept exactly, with no searching required.- DiscoveredLegendre 1805; Gauss (claimed 1795)
- Objectiveminimize ∑(yᵢ − a − bxᵢ)²
- Slopeb = Sₓᵧ / Sₓₓ
- Intercepta = ȳ − b·x̄
- Line always passes throughthe centroid (x̄, ȳ)
- OptimalityBLUE — Gauss–Markov theorem
Watch the 60-second explainer
A condensed visual walkthrough — narrated, captioned, under a minute.
The exact rule: what "best fit" means
Given n data pairs (x₁, y₁), …, (xₙ, yₙ), we look for a line ŷ = a + bx. For each point the residual is the vertical miss eᵢ = yᵢ − (a + bxᵢ) — the signed gap between the real value and what the line predicts. Least squares chooses the intercept a and slope b that minimize the sum of squared residuals:
S(a, b) = ∑ᵢ (yᵢ − a − bxᵢ)².
This is not a heuristic — it is a genuine minimization problem with a unique answer whenever the xᵢ are not all identical. The minimizing (a, b) satisfy the two normal equations, obtained by setting ∂S/∂a = 0 and ∂S/∂b = 0:
- ∑ yᵢ = n·a + b·∑ xᵢ
- ∑ xᵢyᵢ = a·∑ xᵢ + b·∑ xᵢ²
Solving them gives the two formulas you actually use. Writing x̄ and ȳ for the means, Sₓₓ = ∑(xᵢ − x̄)² and Sₓᵧ = ∑(xᵢ − x̄)(yᵢ − ȳ):
b = Sₓᵧ / Sₓₓ and a = ȳ − b·x̄.
That second formula says something beautiful: the best-fit line always passes exactly through the centroid (x̄, ȳ) of the data. The line pivots about the center of mass; only its tilt is up for negotiation.
Why squares? The intuition behind the objective
The animation makes the objective literal: hang a small square of area eᵢ² off each point, and least squares is the line that makes the total shaded area as small as possible. But why squares and not, say, absolute values?
Three reasons pile up. First, differentiability: eᵢ² is smooth, so S(a, b) is a friendly upward-opening paraboloid with one bottom — calculus finds it in closed form. |eᵢ| has a kink at zero and needs a linear program instead. Second, the mean falls out: for the flat model ŷ = a (no x at all), minimizing ∑(yᵢ − a)² gives a = ȳ, the arithmetic mean. Squaring is precisely the criterion whose optimum is the mean — a deep consistency with the rest of statistics. Third, geometry: think of y = (y₁, …, yₙ) as one point in ℝⁿ. The predictions ŷ that any line can produce form a 2-dimensional plane (spanned by the all-ones vector and the x vector). Least squares is the orthogonal projection of y onto that plane. Minimizing squared error is minimizing Euclidean distance, and the closest point on a plane is the foot of the perpendicular. The residual vector e = y − ŷ is orthogonal to the plane — which is exactly what the normal equations assert.
The pitfall this exposes: squaring makes a residual of 10 count 100× as much as a residual of 1, so a single wild outlier can drag the whole line toward itself. Least squares chases the mean, and the mean is not robust.
A worked example, carried through
Take five points: (1, 2), (2, 4), (3, 5), (4, 4), (5, 5). We fit ŷ = a + bx by hand.
Step 1 — means. x̄ = (1+2+3+4+5)/5 = 3. ȳ = (2+4+5+4+5)/5 = 20/5 = 4.
Step 2 — deviations and products. With dxᵢ = xᵢ − 3 and dyᵢ = yᵢ − 4:
- (1,2): dx = −2, dy = −2 → product 4, dx² = 4
- (2,4): dx = −1, dy = 0 → product 0, dx² = 1
- (3,5): dx = 0, dy = 1 → product 0, dx² = 0
- (4,4): dx = 1, dy = 0 → product 0, dx² = 1
- (5,5): dx = 2, dy = 1 → product 2, dx² = 4
Step 3 — sums. Sₓᵧ = 4 + 0 + 0 + 0 + 2 = 6. Sₓₓ = 4 + 1 + 0 + 1 + 4 = 10.
Step 4 — slope and intercept. b = Sₓᵧ / Sₓₓ = 6/10 = 0.6. a = ȳ − b·x̄ = 4 − 0.6·3 = 4 − 1.8 = 2.2.
The line of best fit is ŷ = 2.2 + 0.6x. Check the centroid: at x = 3, ŷ = 2.2 + 1.8 = 4 = ȳ. ✓
Step 5 — residuals and fit quality. Predictions are 2.8, 3.4, 4.0, 4.6, 5.2; residuals are −0.8, 0.6, 1.0, −0.6, −0.2. Sum of squared residuals SSE = 0.64 + 0.36 + 1.00 + 0.36 + 0.04 = 2.40. Total variation Sᵧᵧ = ∑(yᵢ − ȳ)² = 4 + 0 + 1 + 0 + 1 = 6. The coefficient of determination R² = 1 − SSE/Sᵧᵧ = 1 − 2.40/6 = 0.60, so the line explains 60% of the variance in y.
What the estimator guarantees: Gauss–Markov
Least squares is not just convenient — under mild assumptions it is provably optimal. Suppose the data truly come from yᵢ = α + βxᵢ + εᵢ where the noise εᵢ has mean 0, constant variance σ² (homoscedasticity), and is uncorrelated across observations. The Gauss–Markov theorem then says the OLS estimates (a, b) are BLUE: the Best Linear Unbiased Estimators. Among all estimators that are linear in the yᵢ and unbiased, none has smaller variance. Notably this requires no normality assumption at all — only mean, variance, and correlation conditions.
Two useful sampling facts follow. The slope estimator is unbiased, E[b] = β, and its variance is Var(b) = σ² / Sₓₓ. Read that formula: to pin down a slope precisely you want your x-values spread far apart (large Sₓₓ). Cramming all your measurements into a narrow x-window is the fastest way to a wobbly slope. If, additionally, the noise is Gaussian, OLS coincides with the maximum-likelihood estimate and t-tests and confidence intervals become exact.
The estimate of the noise scale is s² = SSE/(n − 2), dividing by n − 2 because two parameters (a and b) were fitted — the residuals live in an (n − 2)-dimensional space, not n.
Generalizations and where the idea travels
Least squares is the base case of a much larger structure. Stack the data into a matrix X (a column of 1s plus a column of x-values) and a vector y; the normal equations become the single matrix identity XᵀX β̂ = Xᵀy, solved by β̂ = (XᵀX)⁻¹Xᵀy. This one line fits any linear-in-parameters model: multiple predictors (multiple regression), or polynomials ŷ = a + bx + cx² + … (still "linear" because it is linear in the unknown coefficients). Fitting a parabola to points is least squares with an extra column x².
When the noise variances differ, weighted least squares minimizes ∑ wᵢeᵢ² with wᵢ = 1/σᵢ², down-weighting noisy points. When XᵀX is near-singular (collinear predictors), ridge regression minimizes ∑eᵢ² + λ‖β‖², nudging the solution toward stability — the birth of regularization. And if you drop the vertical-only convention and minimize perpendicular distances, you get total least squares, computed by the singular value decomposition, which treats x and y symmetrically (right when both axes carry measurement error).
A subtle warning that all of these inherit: regression of y on x and regression of x on y give different lines (their slopes multiply to r², the squared correlation). "The" best-fit line is not unique until you declare which variable bears the error.
| Criterion | Minimizes | Closed form? | Sensitivity to outliers |
|---|---|---|---|
| Ordinary Least Squares (OLS) | ∑ residualᵢ² | Yes — normal equations | High (squares punish big gaps) |
| Least Absolute Deviations (LAD / L¹) | ∑ |residualᵢ| | No — linear program | Low (robust, follows the median) |
| Total Least Squares (orthogonal) | ∑ perpendicular distance² | Yes — via SVD/PCA | High, but symmetric in x and y |
Frequently asked questions
Why minimize the squared residuals instead of the absolute residuals?
Squaring gives a smooth objective with a unique closed-form minimum (the normal equations), and its optimum for a constant model is the arithmetic mean — consistent with the rest of classical statistics. Absolute residuals (L¹ / least absolute deviations) instead track the median and are far more robust to outliers, but require solving a linear program and can have non-unique solutions. The trade-off is efficiency vs. robustness.
Does the best-fit line have to pass through the average point?
Yes, always. The intercept formula a = ȳ − b·x̄ rearranges to ȳ = a + b·x̄, meaning the point (x̄, ȳ) — the centroid of the data — lies exactly on the line. The line pivots about the center of mass of the cloud; only its slope is determined by the correlation structure.
What does R² actually measure?
R² = 1 − SSE/Sᵧᵧ is the fraction of the variance in y that the line explains, ranging from 0 (line no better than the flat mean) to 1 (perfect fit). For simple linear regression it equals the square of the Pearson correlation r. In the worked example R² = 0.60, so 60% of y's variation is captured by x. It does not tell you whether a line is the right model — a curved relationship can still show a high R².
When does least squares fail or mislead?
When assumptions break: (1) outliers, since squaring lets one bad point dominate; (2) a genuinely nonlinear relationship, where a line hides the real shape (Anscombe's quartet is the classic demonstration of four datasets with identical regression lines but wildly different pictures); (3) heteroscedastic or correlated errors, which invalidate the standard variance formulas; and (4) near-collinear predictors, which make (XᵀX)⁻¹ unstable.
What are the normal equations and where do they come from?
They are the two equations you get by setting the partial derivatives of S(a,b) = ∑(yᵢ − a − bxᵢ)² to zero: ∑yᵢ = na + b∑xᵢ and ∑xᵢyᵢ = a∑xᵢ + b∑xᵢ². Because S is a convex paraboloid, these stationarity conditions pinpoint the global minimum. Geometrically they state that the residual vector is orthogonal (normal) to the space of fitted values — hence the name.
Why divide by n − 2 when estimating the error variance?
Two parameters, a and b, were estimated from the data, so the residuals are constrained to lie in an (n − 2)-dimensional subspace and are on average smaller than the true errors. Dividing SSE by n − 2 (the degrees of freedom) corrects this bias, giving an unbiased estimate s² of the noise variance σ². For a mean-only model you'd divide by n − 1 instead.