Linear regression

1. Intuition

Imagine you sell houses. Bigger houses usually cost more, but not perfectly: two houses of the same size can sell for different prices. You want a rule of thumb like:

price ≈ (price per square foot) × size + (base price)

That rule is a straight line. The slope is the price per unit of size, and the bias (intercept) is where the line crosses zero size. Linear regression picks the slope and bias automatically, using past sales.

How does it judge a line? For each past sale it looks at the residual: the actual price minus the line's prediction. A good line has small residuals. Because residuals can be positive or negative, it squares them (so they cannot cancel, and big misses count extra) and averages them. That average is the mean squared error (MSE). Training means sliding and tilting the line until the MSE is as small as it gets.

With more than one feature (size, bedrooms, age), the line becomes a flat plane (or hyperplane), and each feature gets its own weight. The idea is unchanged.

2. Visualization

Try it (predict first, then check):

1. Press Add outlier and predict which way the best-fit line will tilt before pressing Best fit.
2. Raise the learning rate past 0.04 and press Gradient step a few times. Predict what MSE does first.
3. Switch to L1 and raise λ. At what value does the slope hit exactly zero, and why doesn't L2 do the same?

What to notice:

• Each orange segment is one residual. MSE is the average of their squared lengths.
• A single far-away point produces a long segment whose square dominates the sum. That is why MSE is sensitive to outliers.
• On the right, the loss drops fast at first (the gradient is large when the line is far off) and then slowly (the gradient shrinks near the optimum).

3. The math

Symbols

Prediction

For the whole dataset at once: $\hat{y} = Xw + b$. Check the shapes: $[n, p] \times [p] = [n]$, then $b$ is added to every entry.

Loss: mean squared error

The term $y_i - \hat{y}_i$ is the residual. MSE is the loss the optimizer minimizes. You might report a different metric to humans (RMSE, MAE, R²). Loss and metric are not the same thing.

Gradients

With one feature, $\hat{y}_i = w x_i + b$. Apply the chain rule to each squared term. The derivative of $(y_i - \hat{y}_i)^2$ with respect to $\hat{y}_i$ is $-2(y_i - \hat{y}_i) = 2(\hat{y}_i - y_i)$. Then $\partial \hat{y}_i / \partial w = x_i$ and $\partial \hat{y}_i / \partial b = 1$. Averaging:

Gradient descent moves each parameter a small step against its gradient:

Worked example

Data: $x = [1, 2, 3]$, $y = [2, 4, 5]$. Start from $w = 1$, $b = 0$, and use learning rate $\eta = 0.1$.

Step 1: predictions. $\hat{y}_i = 1 \cdot x_i + 0$

Step 2: residuals $y_i - \hat{y}_i$: $\;2-1 = 1,\;\; 4-2 = 2,\;\; 5-3 = 2$.

Step 3: squared errors: $1^2 = 1,\;\; 2^2 = 4,\;\; 2^2 = 4$.

Step 4: MSE $= (1 + 4 + 4)/3 = 9/3 = 3$.

Step 5: gradient with respect to slope. The errors $\hat{y}_i - y_i$ are $-1, -2, -2$.

Step 6: gradient with respect to bias.

Both gradients are negative. That means increasing $w$ or $b$ would decrease the loss, so the update should increase them.

Step 7: one gradient-descent update.

Check. New predictions are $1.7333 + 0.3333 = 2.067$, then $3.467 + 0.333 = 3.800$, then $5.200 + 0.333 = 5.533$. Residuals are $-0.067, 0.200, -0.533$. New MSE $= (0.0044 + 0.0400 + 0.2844)/3 = 0.1096$. One step took the loss from 3 to about 0.11.

Note

For linear regression you do not need gradient descent. The normal equation $w = (X^\top X)^{-1}X^\top y$ solves it exactly. For this data it gives $w = 1.5$, $b = 0.667$, MSE $= 0.0556$. We use gradient descent here because it is the same algorithm that trains models with no closed form, like neural networks.

4. Implementation

You will fit linear regression two ways from scratch, by gradient descent and by solving the normal equation, on one feature and then on three. Along the way you will see why feature scaling decides whether gradient descent converges in 1,000 steps or crawls.

Build it from scratch

Start with one feature: the worked example as a loop. Each pass predicts, measures the error, and moves w and b against the MSE gradient from section 3.

import numpy as np

def fit_linear_regression(x, y, lr=0.01, steps=2000):
    w, b = 0.0, 0.0
    for _ in range(steps):
        y_hat = w * x + b                      # forward pass: predictions
        error = y_hat - y                      # (ŷ - y), shape [n]
        w -= lr * (2 / len(x)) * np.sum(error * x)
        b -= lr * (2 / len(x)) * np.sum(error)
    return w, b

Predict: on the worked example's three points, the first step (lr = 0.1) moved w from 1.0 to 1.7333 and cut MSE from 3.0 to 0.11. What happens if you set lr = 1.0?

With many features, fold the bias into the weights: add a column of ones to X, so θ = [b, w1, ..., wp] and every prediction is X θ. Then the gradient-descent loop is the same three lines, with a matrix product in place of the sums.

def add_bias(X):
    return np.c_[np.ones(len(X)), X]                       # [n, p] -> [n, p+1], column of 1s for b

def fit_gradient_descent(X, y, lr=0.1, steps=1000):
    Xb = add_bias(X)
    theta = np.zeros(Xb.shape[1])
    n = len(y)
    for _ in range(steps):
        error = Xb @ theta - y                              # [n]
        grad = (2 / n) * Xb.T @ error                       # [p+1], dMSE/dθ
        theta -= lr * grad
    return theta

Setting that gradient to zero gives the normal equation, $X^\top X\,\theta = X^\top y$: one linear solve, no learning rate. Use np.linalg.solve rather than computing the inverse. It is faster and loses less precision.

def fit_normal_equation(X, y):
    """Solve (XᵀX) θ = Xᵀy. np.linalg.solve is stabler than forming the inverse."""
    Xb = add_bias(X)
    return np.linalg.solve(Xb.T @ Xb, Xb.T @ y)             # θ = [b, w1, ..., wp]

Predict: the script fits synthetic houses with size (0.5 to 4 thousand sq ft), bedrooms (1 to 5), and age (0 to 80 years). The true θ is [50, 120, 15, −0.8]. Gradient descent on raw features needs lr = 1e-4 to avoid diverging. After 1,000 steps, will it be close to the normal-equation answer? What changes after standardizing every feature to mean 0 and standard deviation 1?

Use the library

from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(X_train, y_train)          # X_train shape [n, p], y_train shape [n]
print(model.coef_, model.intercept_) # learned w and b
y_pred = model.predict(X_test)

scikit-learn solves the least-squares problem directly with a numerically careful solver (no learning rate to tune), so it matches the normal equation. For very large or streaming data, SGDRegressor is the gradient-descent version, and it needs StandardScaler for the reason you just saw. The normal-equation cost grows with p³, so with hundreds of thousands of features iterative methods win.

Run it

From guide/:

uv run code/03-regression/linear_regression.py
uv run code/03-regression/linear_regression_scratch.py

The first script reproduces the worked example step by step and ends with the one-feature fit:

Synthetic dataset (true w=3, b=8)
  gradient descent  w=2.824  b=8.979
  scikit-learn      w=2.824  b=8.980

The second runs the multi-feature comparison:

normal equation θ = [ 52.288 121.85   13.598  -0.859],  R² = 0.9807
GD raw features, lr=1e-4, 1000 steps: θ = [11.283 47.817 31.452  2.218]
GD standardized, lr=0.1, 1000 steps: θ = [ 52.288 121.85   13.598  -0.859]
condition number of XᵀX: raw 28,185, standardized 1.4
scikit-learn: b = 52.288, w = [121.85   13.598  -0.859]
checks pass: normal equation = sklearn = standardized GD

The estimates differ from the true [50, 120, 15, −0.8] because of the noise in 200 synthetic rows. Scripts: code/03-regression/linear_regression.py (worked example, figure) and code/03-regression/linear_regression_scratch.py (multi-feature, scaling).

5. Engineering

Good use cases. Fast, interpretable baselines for continuous targets. Pricing and demand models where stakeholders want to read the coefficients. Situations with few examples, where a flexible model would overfit.

Assumptions to check.

• Linearity: the target changes roughly linearly with each feature. Plot residuals against predictions; a curve means you are missing a nonlinear term (try polynomial features or a tree model).
• Independent errors: residuals of one example do not predict another's. Time series often break this.
• Constant error variance (homoscedasticity): residual spread does not grow with the prediction. If it fans out, try predicting log(y).
• Low multicollinearity: highly correlated features make individual weights unstable even when predictions are fine. Ridge regression fixes this.

Preprocessing. Scale features (standardize) when using gradient descent or regularization, otherwise the feature with the largest units dominates the step size and the penalty. Encode categorical features (one-hot). Fit every preprocessing step on the training split only.

Cost. Training via the normal equation costs $O(np^2 + p^3)$. Gradient descent costs $O(np)$ per step. Inference is a single dot product, $O(p)$, which is microseconds. Linear models are among the cheapest models to serve.

Production concerns. Coefficients drift when the world changes (for example, price per square foot rises). Monitor the residual distribution in production, not just the error average. Clip or validate inputs, because a linear model extrapolates without limit: a 50,000 sq ft input produces a confident, absurd price.

Watch out

Failure modes. Outliers pull the line toward them, because squared error punishes big misses quadratically (MAE or Huber loss are more robust). Extrapolation far beyond the training range is unreliable. A good fit does not mean a causal relationship: a positive weight on "has a pool" does not mean adding a pool raises the price by exactly that much.

Common mistakes

• Reporting training MSE as if it were test performance.
• Reading coefficient size as feature importance when features are on different scales.
• Forgetting the bias term, which forces the line through the origin.
• Fitting the scaler on all data before splitting (data leakage).
• Using a learning rate so large that the loss grows each step. Watch the loss curve.

6. Knowledge check

Practice exercise. With data $x = [0, 1, 2]$, $y = [1, 3, 5]$, start at $w = 0$, $b = 0$, $\eta = 0.1$. Compute the MSE, both gradients, and the updated parameters.

Implementation challenge. Extend fit_linear_regression to multiple features using matrix operations: y_hat = X @ w + b, grad_w = (2/n) * X.T @ (y_hat - y). Verify shapes ($X$ is $[n,p]$, $X^\top(\hat{y} - y)$ is $[p]$) and compare your weights with LinearRegression on a synthetic dataset with 3 features.

Summary

• Linear regression predicts $\hat{y} = w \cdot x + b$: a weighted sum of features plus a bias.
• It is trained by minimizing MSE, the average squared residual.
• The gradients are $\frac{2}{n}\sum(\hat{y}-y)x$ for the weight and $\frac{2}{n}\sum(\hat{y}-y)$ for the bias. Gradient descent steps against them.
• It is cheap, interpretable, and a mandatory baseline, but it is sensitive to outliers and cannot capture nonlinear patterns without engineered features.

Next: Regularization and regression metrics

Related: Gradient descent · Logistic regression · Model card: linear regression

Parts of this lesson are adapted from AI Engineering from Scratch by Rohit Ghumare (MIT license, notice): Linear Systems, Linear Regression.

Interview angle