Linear regression is over two hundred years old — Legendre and Gauss used least squares to predict the orbits of comets — and it remains one of the most useful tools in data science. It is also the perfect laboratory: nearly every concept in ML (loss functions, optimisation, regularisation, overfitting, probabilistic interpretation) appears in its simplest form here.
The model#
Given features $\mathbf{x} \in \mathbb{R}^d$, linear regression predicts
where we append a constant feature $x_0 = 1$ so the intercept $w_0$ is absorbed into $\mathbf{w}$. For $n$ examples stacked as rows of the design matrix $\mathbf{X} \in \mathbb{R}^{n \times (d+1)}$, predictions are $\hat{\mathbf{y}} = \mathbf{X}\mathbf{w}$.
"Linear" means linear in the parameters — the features themselves can be non-linear transformations of raw inputs ($x^2$, $\log x$), as we will see next lecture.
The loss: mean squared error#
Why squares? Three reasons: (1) it is smooth and convex, giving a unique closed-form solution; (2) it corresponds to maximum likelihood under Gaussian noise; (3) it penalises large errors heavily. The last reason is also its weakness — outliers dominate the fit.
The closed-form solution#
Setting the gradient to zero:
These are the normal equations. If $\mathbf{X}^\top\mathbf{X}$ is invertible:
Geometric view. $\hat{\mathbf{y}} = \mathbf{X}\hat{\mathbf{w}}$ is the orthogonal projection of $\mathbf{y}$ onto the column space of $\mathbf{X}$; the residual vector is perpendicular to every feature column.
Gradient descent solution#
The closed form costs $O(nd^2 + d^3)$. For very large $d$ or streaming data, use gradient descent:
import numpy as np
rng = np.random.default_rng(0)
n = 200
area = rng.uniform(40, 200, n) # square metres
rooms = rng.integers(1, 6, n)
price = 15 + 0.9 * area + 8 * rooms + rng.normal(0, 10, n) # in lakh taka (synthetic)
X = np.column_stack([np.ones(n), area, rooms])
y = price
# 1) Least squares (stable)
w_ls, *_ = np.linalg.lstsq(X, y, rcond=None)
# 2) Gradient descent on standardised features (for good conditioning)
mu, sd = X[:, 1:].mean(0), X[:, 1:].std(0)
Xs = np.column_stack([np.ones(n), (X[:, 1:] - mu) / sd])
w = np.zeros(3)
for _ in range(2000):
w += 0.1 * 2 / n * Xs.T @ (y - Xs @ w)
# convert back to original units
w_gd = np.r_[w[0] - (w[1:] * mu / sd).sum(), w[1:] / sd]
print("least squares:", w_ls.round(3))
print("grad descent :", w_gd.round(3))Both recover coefficients close to the true values (15, 0.9, 8). Note the standardisation before gradient descent — without it, the very different feature scales create an ill-conditioned problem and gradient descent crawls.
Evaluating a regression model#
- MSE / RMSE — in squared / original units.
- MAE — robust to outliers.
- $R^2$ — fraction of variance explained: $R^2 = 1 - \frac{\sum(y_i - \hat{y}_i)^2}{\sum(y_i - \bar{y})^2}$. $R^2 = 0$ means no better than predicting the mean; it can be negative on test data.
Always plot residuals against predictions and against each feature. Patterns in residuals (curves, funnels) reveal missing non-linearity or non-constant variance (heteroscedasticity).
Interpreting coefficients#
$w_j$ is the expected change in $y$ for a one-unit increase in $x_j$, holding all other features fixed. Cautions:
- Coefficients depend on units; standardise to compare importance.
- With correlated features, individual coefficients can be unstable or have counter-intuitive signs.
- Association is not causation. A coefficient describes a pattern in observational data, not the effect of intervening.
Statistical assumptions#
For valid confidence intervals and hypothesis tests on coefficients, classical theory assumes: linearity, independent errors, constant error variance, and (for small samples) normally distributed errors. The Gauss–Markov theorem states that under the first three, OLS is the best linear unbiased estimator (BLUE). For prediction alone, these assumptions matter less — but checking residuals remains good practice.