Linear regression

The model and the loss

A linear model predicts a number as a weighted sum of features plus an offset:

Absorb the bias by appending a constant feature $x_0 = 1$, so $\hat y = \mathbf{w}^\top\mathbf{x}$ with $\mathbf{w}\in\mathbb{R}^{d+1}$. Stacking the $n$ training inputs as rows of the design matrix $\mathbf{X}\in\mathbb{R}^{n\times(d+1)}$ and the targets as $\mathbf{y}\in\mathbb{R}^n$, all predictions at once are $\hat{\mathbf{y}} = \mathbf{X}\mathbf{w}$ (1.2 Matrices as transformations). The loss is the mean squared error,

a convex quadratic bowl in $\mathbf{w}$: it has exactly one minimum (when $\mathbf{X}$ has full column rank), which is why three different methods must agree.

Way 1 — algebra: the normal equations

Expand the loss: $\|\mathbf{X}\mathbf{w}-\mathbf{y}\|^2 = \mathbf{w}^\top\mathbf{X}^\top\mathbf{X}\mathbf{w} - 2\mathbf{y}^\top\mathbf{X}\mathbf{w} + \mathbf{y}^\top\mathbf{y}$. Using the identities $\partial(\mathbf{w}^\top\mathbf{A}\mathbf{w})/\partial\mathbf{w} = 2\mathbf{A}\mathbf{w}$ for symmetric $\mathbf{A}$ and $\partial(\mathbf{a}^\top\mathbf{w})/\partial\mathbf{w} = \mathbf{a}$ (1.5 Multivariable and matrix calculus):

Setting the gradient to zero gives the normal equations and their solution:

$\mathbf{X}^\top\mathbf{X}$ is a $(d+1)\times(d+1)$ matrix; it is invertible exactly when the columns of $\mathbf{X}$ are linearly independent (no feature is a combination of others, and $n \ge d+1$). Forming and solving costs $O(nd^2 + d^3)$ — instant for $d$ in the hundreds, impossible for $d$ in the millions. In code, never invert: np.linalg.lstsq(X, y) or np.linalg.solve(X.T @ X, X.T @ y) is more stable, and when columns are nearly collinear the SVD-based solver (or ridge, 3.6 Regularization and generalization) is the right tool (1.3 Eigen-things, SVD and low rank).

Way 2 — geometry: projection

$\hat{\mathbf{y}} = \mathbf{X}\mathbf{w}$ ranges over the column space of $\mathbf{X}$ — the set of all vectors that are linear combinations of the feature columns, a $(d+1)$-dimensional flat inside $\mathbb{R}^n$. Least squares asks for the point of that flat closest to $\mathbf{y}$: its orthogonal projection. Closest means the residual $\mathbf{r} = \mathbf{y} - \mathbf{X}\mathbf{w}^*$ is perpendicular to every column of $\mathbf{X}$:

which is the normal equations again. Two consequences you can check on any fit: because the constant column is in $\mathbf{X}$, the residuals sum to exactly zero; and the residuals are uncorrelated with every feature — if a residual plot shows a pattern against a feature, you need a transform of that feature, not the feature itself. The projection is linear in $\mathbf{y}$: $\hat{\mathbf{y}} = \mathbf{H}\mathbf{y}$ with the hat matrix $\mathbf{H} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top$, a projection ($\mathbf{H}^2=\mathbf{H}$) whose diagonal entries measure each point's leverage (1.1 Vectors, dot products and norms for projections).

Worked example

Three points $(x, y)$: $(1,1), (2,2), (3,2)$. With a constant column, $\mathbf{X} = \begin{bmatrix}1&1\\1&2\\1&3\end{bmatrix}$, $\mathbf{y} = (1,2,2)^\top$.

So $\hat y = \tfrac{2}{3} + \tfrac{1}{2}x$. Predictions $(1.167, 1.667, 2.167)$, residuals $(-0.167, 0.333, -0.167)$ — they sum to zero as promised. The sum of squared residuals is $\text{SSE} = 1/6 \approx 0.167$; the total sum of squares around $\bar y = 5/3$ is $\text{SST} = 2/3$; so

the line explains 75% of the variance of $y$.

Way 3 — gradient descent

For large $d$ or streaming data, iterate instead of solving (1.10 Optimization: from gradient descent to Adam):

Because $J$ is a convex quadratic, gradient descent converges to $\mathbf{w}^*$ whenever $\eta < 2/\lambda_{\max}$, where $\lambda_{\max}$ is the largest eigenvalue of the Hessian $\tfrac{2}{n}\mathbf{X}^\top\mathbf{X}$ (1.3 Eigen-things, SVD and low rank). The speed is governed by the ratio $\lambda_{\max}/\lambda_{\min}$: a feature measured in millimetres next to one measured in kilometres makes the bowl a long narrow valley, and the iterates zig-zag. Standardizing features (subtract mean, divide by standard deviation) makes the bowl nearly round and lets you use a much larger step. Mini-batch SGD uses a random subset per step; its gradient is an unbiased estimate of the full gradient, which is all the theory needs (4.5 Optimizers and learning-rate schedules). The interactive demo lets you watch exactly this: add points, start descent, and see the path and the loss curve — then compare with the closed form.

Why squared error? Maximum likelihood

Assume the data were generated by a line plus Gaussian noise: $y_i = \mathbf{w}^\top\mathbf{x}_i + \varepsilon_i$, $\varepsilon_i \sim \mathcal{N}(0, \sigma^2)$ independently. Then $p(y_i\mid\mathbf{x}_i,\mathbf{w}) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\!\big(-\frac{(y_i-\mathbf{w}^\top\mathbf{x}_i)^2}{2\sigma^2}\big)$ and the log-likelihood of the dataset is

The first term does not depend on $\mathbf{w}$, so maximizing the likelihood (1.7 Probability II: Bayes, independence, MLE and MAP) is exactly minimizing the sum of squared errors. This tells you what MSE assumes: symmetric, light-tailed, constant-variance noise. When outliers are common (heavy tails), the Gaussian assumption is wrong and the squared loss lets a few points dominate; the MLE under Laplace noise is the absolute error, and the Huber loss interpolates (4.3 Loss functions catalog). The same recipe — write a likelihood, take the negative log — produces every loss function in this path (1.9 Information theory: entropy, cross-entropy and KL divergence, 3.3 Logistic regression and softmax classification).

Beyond straight lines, and reading a fit

"Linear" means linear in the parameters, not in the inputs. Replace $\mathbf{x}$ by features $\phi(\mathbf{x})$ — powers $x, x^2, x^3$, interactions $x_1x_2$, logs, splines, or the output of a pretrained network — and everything above still holds with $\mathbf{X}$ replaced by $\Phi$. Polynomial degree is the capacity knob from 3.1 What is learning? ERM, generalization and the ML workflow: degree 1 underfits a curve, degree 15 memorizes noise; 3.6 Regularization and generalization shows how to keep high degree and control the damage.

Diagnostics in three plots: predicted vs actual (should hug the diagonal), residuals vs predicted (should be a structureless band — a curve means a missing nonlinearity, a funnel means non-constant variance), and a histogram of residuals (roughly symmetric). Report RMSE in the units of $y$ alongside $R^2$; $R^2$ alone can be high on a useless model with a huge spread of $y$, and it never decreases when you add features, which is why the adjusted version penalizes $d$. Coefficients are comparable only when features are on the same scale.

import numpy as np
from sklearn.linear_model import LinearRegression

Xb = np.c_[np.ones(len(X)), X]                     # add the constant column
w_closed, *_ = np.linalg.lstsq(Xb, y, rcond=None)  # normal equations, stably

w = np.zeros(Xb.shape[1]); lr = 0.1
for _ in range(2000):                              # gradient descent
    w -= lr * (2 / len(y)) * Xb.T @ (Xb @ w - y)

model = LinearRegression().fit(X, y)               # scikit-learn
np.allclose(w_closed[1:], model.coef_), np.allclose(w, w_closed, atol=1e-4)

Where this shows up later

• A linear layer with no activation is this model; stack them with nonlinearities and you have an MLP — 4.1 From perceptron to MLP; universal approximation; the same gradient $\mathbf{X}^\top(\hat{\mathbf{y}}-\mathbf{y})$ reappears for logistic regression — 3.3 Logistic regression and softmax classification.
• Add $\lambda\|\mathbf{w}\|^2$ and the normal equations become ridge — 3.6 Regularization and generalization; the SVD view of the inverse explains why collinearity hurts — 1.3 Eigen-things, SVD and low rank.
• MSE-as-Gaussian-MLE is the template for every loss — 4.3 Loss functions catalog; residual analysis is error analysis — 3.14 Project: end-to-end tabular ML with scikit-learn.
• Linear probes on LLM activations are linear regressions/classifiers — 8.11 Interpretability: what is happening inside; the LoRA update is a low-rank linear map — 8.2 Fine-tuning and PEFT.