Linear Regression

Linear regression predicts a continuous target $y$ from a feature vector $\mathbf{x}$. It is simple, fast and easy to interpret, and almost every idea in this course (loss functions, optimisation, regularisation) appears here first.

The model

Given $n$ training examples $\{(\mathbf{x}^{(i)}, y^{(i)})\}_{i=1}^{n}$ with $\mathbf{x}^{(i)} \in \mathbb{R}^d$, we assume

$$\hat{y} = h_{\boldsymbol\theta}(\mathbf{x}) = \theta_0 + \theta_1 x_1 + \dots + \theta_d x_d = \boldsymbol\theta^\top \mathbf{x},$$

where we prepend $x_0 = 1$ to every example so the intercept $\theta_0$ fits into the dot product.

Stacking all examples as rows gives the design matrix $X \in \mathbb{R}^{n \times (d+1)}$ and the prediction vector $\hat{\mathbf{y}} = X\boldsymbol\theta$.

Loss function

We measure how wrong the model is with the mean squared error (MSE):

$$J(\boldsymbol\theta) = \frac{1}{2n} \sum_{i=1}^{n} \left( \boldsymbol\theta^\top \mathbf{x}^{(i)} - y^{(i)} \right)^2 = \frac{1}{2n} \lVert X\boldsymbol\theta - \mathbf{y} \rVert_2^2 .$$

The factor $\tfrac{1}{2}$ is only there to cancel the 2 that appears when we differentiate.

Solving for the parameters

The normal equation

$J$ is a convex quadratic, so its minimum is where the gradient vanishes:

$$\nabla_{\boldsymbol\theta} J = \frac{1}{n} X^\top (X\boldsymbol\theta - \mathbf{y}) = \mathbf{0} \quad\Longrightarrow\quad \boldsymbol\theta^\star = (X^\top X)^{-1} X^\top \mathbf{y}.$$

This is exact but costs $O(d^3)$ to invert $X^\top X$, which becomes slow when $d$ is large.

Gradient descent

Instead of solving in one step, we repeatedly move against the gradient with a learning rate $\alpha$:

$$\boldsymbol\theta \leftarrow \boldsymbol\theta - \alpha \cdot \frac{1}{n} X^\top (X\boldsymbol\theta - \mathbf{y}).$$

Implementation in Python

import numpy as np


def fit_normal_equation(X: np.ndarray, y: np.ndarray) -> np.ndarray:
    """Closed-form solution. X must already include a column of ones."""
    return np.linalg.solve(X.T @ X, X.T @ y)


def fit_gradient_descent(
    X: np.ndarray, y: np.ndarray, lr: float = 0.1, epochs: int = 1_000
) -> np.ndarray:
    n, d = X.shape
    theta = np.zeros(d)
    for _ in range(epochs):
        grad = X.T @ (X @ theta - y) / n
        theta -= lr * grad
    return theta


rng = np.random.default_rng(0)
x = rng.uniform(-1, 1, size=200)
y = 3.0 + 2.0 * x + rng.normal(scale=0.1, size=200)
X = np.column_stack([np.ones_like(x), x])

print(fit_normal_equation(X, y))   # ≈ [3.0, 2.0]
print(fit_gradient_descent(X, y))  # ≈ [3.0, 2.0]

Notice that we use np.linalg.solve instead of computing an explicit inverse — it is faster and numerically more stable.

Evaluating the model

The coefficient of determination tells us what fraction of the variance in $y$ the model explains:

$$R^2 = 1 - \frac{\sum_i (y^{(i)} - \hat{y}^{(i)})^2}{\sum_i (y^{(i)} - \bar{y})^2}.$$

Metric	Formula	Interpretation
MSE	$\frac{1}{n}\sum (y - \hat y)^2$	Average squared error, in squared units of $y$
RMSE	$\sqrt{\text{MSE}}$	Same units as $y$
$R^2$	see above	1 is perfect, 0 is no better than predicting $\bar y$

Summary

• Linear regression models $\hat y = \boldsymbol\theta^\top\mathbf{x}$ and minimises the MSE.
• The normal equation gives an exact solution; gradient descent scales to large $d$.
• Always scale features and evaluate on held-out data.

Exercise. Show that $X^\top X$ is invertible if and only if the columns of $X$ are linearly independent. What happens to the normal equation when two features are identical?