Why does linear regression minimise the sum of squared errors instead of, say, absolute errors?

Ordinary least squares (OLS) minimises squared residuals for both mathematical and statistical reasons. Mathematically, the squared error loss L(β) = Σ(yᵢ - Xᵢβ)² is smooth and differentiable everywhere, so its minimum can be found analytically by setting the gradient to zero — this gives the closed-form normal equation β = (XᵀX)⁻¹Xᵀy. Absolute error |yᵢ - Xᵢβ| has a non-differentiable kink at zero, which prevents a clean closed-form solution and requires iterative optimisation (like linear programming) instead.

Statistically, minimising squared error is equivalent to maximum likelihood estimation under the assumption that residuals are normally distributed with constant variance (homoscedastic). The squared loss heavily penalises large errors — a residual twice as large contributes four times the loss — which makes OLS very sensitive to outliers. This is precisely why Mean Absolute Error (MAE) or Huber loss are preferred when the data contains outliers: they grow linearly rather than quadratically with the error.

import numpy as np
from sklearn.linear_model import LinearRegression

# Closed-form normal equation
X = np.array([[1, 1], [1, 2], [1, 3]])  # design matrix with intercept column
y = np.array([2, 4, 5])
beta = np.linalg.inv(X.T @ X) @ X.T @ y
print(beta)  # [intercept, slope]

# scikit-learn does the same thing under the hood
model = LinearRegression().fit(X[:, 1:], y)
print(model.intercept_, model.coef_)

Geometric intuition: minimising squared error is equivalent to finding the orthogonal projection of y onto the column space of X. This projection is unique and has a clean geometric interpretation — the residual vector is perpendicular to every column of X, which is exactly what XᵀXβ = Xᵀy encodes.