The Project | About Us | Contribute | Donations | License
HOME My courses Completed
◄ Previous
◄

Regression Algorithms

Next ►
►

In 05-02 we classified: delayed, yes or no. But Rutalia's head of operations wants more: "don't tell me whether it will be late, tell me how many minutes it will take". Predicting a continuous value is regression, the other half of supervised learning. In this lesson we work through it in depth with the canonical dataset: linear regression with its exact analytical solution (pure algebra with numpy), and then the absolute protagonist, gradient descent — the iterative optimization algorithm left pending in 05-02 which is, literally, the engine of the neural networks in 05-04. We will complete the picture with Ridge/Lasso regularization, polynomial regression (where the bias-variance trade-off from 05-01 gets formalized) and regression's own metrics.

Contents

  1. The problem: predicting delivery_minutes
  2. Simple linear regression: the least-squares line
  3. Multiple linear regression: the analytical solution with numpy
  4. Gradient descent: THE algorithm
  5. Batch, stochastic and mini-batch
  6. Regularization: Ridge and Lasso
  7. Polynomial regression and the bias-variance trade-off
  8. Regression metrics: MSE, RMSE, MAE, R²

The problem: predicting delivery_minutes

Same canonical dataset as 05-01 (seed 42), new target: the continuous column delivery_minutes.

data = generate_dataset()                       # 05-01
X_num = np.column_stack([data["distance_km"], data["weight_kg"],
                         data["departure_hour"], data["weekday"]])
X = np.column_stack([X_num, one_hot(data["zone"], ZONES)])
y = data["delivery_minutes"]                    # ← continuous, no longer binary
X_tr, X_te, y_tr, y_te = train_test_split_manual(X, y)
X_tr_s, X_te_s = standardize(X_tr, X_te)

Remember that the generator's hidden truth is multiplicative (base × congestion × rush hour ×...): a linear model won't be able to capture it entirely. Perfect: we will see precisely how much escapes it and why.

Simple linear regression: the least-squares line

Let's start with a single feature, the most informative one: distance_km. We are looking for the line

minutes ≈ w · distance + b

What does "the best line" mean? The geometric intuition: plot the 2,000 deliveries as points (distance, minutes). Each candidate line leaves a residual at every point: the vertical distance between the actual minutes and the predicted ones. Least squares picks the line that minimizes the sum of squared residuals. Why squared and not the absolute value? Three reasons: it penalizes large errors more, it makes the function differentiable everywhere (key for what's coming) and it has a closed-form solution:

w = cov(x, y) / var(x)        b = ȳ − w·x̄
x = data["distance_km"]
w = np.cov(x, y, bias=True)[0, 1] / x.var()   # slope
b = y.mean() - w * x.mean()                    # intercept
print(f"minutes ≈ {w:.2f}·km + {b:.2f}")

You will get something like minutes ≈ 6.9·km + 6.5: the model has approximately rediscovered the generator's "6 min/km + 5 fixed" — without anyone telling it. The slope is not exactly 6.0 because the multiplicative factors (congestion, rush hour) inflate the average relationship. This is learning from data in its purest form.

Assumptions worth knowing (the linear model assumes them): a genuinely linear relationship, independent residuals with constant variance and no dominant outliers. When they fail, the line still gets computed... but it stops meaning what you think.

Multiple linear regression: the analytical solution with numpy

With the 13 features, we look for the vector w (13 weights) and the bias b:

ŷ = X·w + b

The standard trick is to absorb b by appending a column of ones to X (so b becomes just another weight). Minimizing the sum of squares ‖X·w − y‖² has an analytical solution, the normal equations:

w* = (XᵀX)⁻¹ Xᵀ y

The geometric intuition in one sentence: X·w can only produce vectors in the subspace spanned by the columns of X; the best ŷ is the orthogonal projection of y onto that subspace, and the normal equations say exactly that (the residual must be perpendicular to all the columns: Xᵀ(y − Xw) = 0).

def linear_regression_analytic(X, y):
    """Solves the normal equations. Returns weights (bias included)."""
    Xb = np.column_stack([np.ones(len(X)), X])     # column of ones → bias
    # lstsq is more numerically stable than explicitly inverting XᵀX
    w, *_ = np.linalg.lstsq(Xb, y, rcond=None)
    return w

def predict_linear(w, X):
    return np.column_stack([np.ones(len(X)), X]) @ w

w = linear_regression_analytic(X_tr_s, y_tr)
pred = predict_linear(w, X_te_s)
print(f"Test RMSE: {np.sqrt(np.mean((pred - y_te) ** 2)):.2f} minutes")

Cost: O(n·d²) to form XᵀX plus O(d³) to solve. With d=13 it is instantaneous. So why bother with iterative methods? Because the closed-form solution neither scales nor generalizes: with d = 10⁵ features the O(d³) blows up, and — more importantly — it only exists for linear least squares. For the logistic regression of 05-02 there is no closed formula. For the networks of 05-04, none either. We need a general minimization algorithm.

Gradient descent: THE algorithm

This section is the heart of the lesson — and the bridge to everything left in the module.

The cost function

We define the mean squared error (MSE) as a function of the parameters:

J(w, b) = (1/n) · Σᵢ (ŷᵢ − yᵢ)²     with  ŷᵢ = w·xᵢ + b

Shift your point of view: the data are fixed; the variables are the parameters. J is a surface over parameter space — for least squares, a convex paraboloid with a single minimum (the one the analytical formula finds in one shot).

The idea: descending the mountain blindfolded

Imagine standing on that surface in fog: you only perceive the slope under your feet. Strategy: take a small step downhill and repeat. The gradient ∇J (the vector of partial derivatives) points in the direction of steepest ascent; walking along −∇J is the fastest descent locally.

repeat until convergence:
    w ← w − α · ∂J/∂w
    b ← b − α · ∂J/∂b

α is the learning rate: the step size.

The derivatives

Differentiating J (chain rule over the square; the 2 gets absorbed into the convention):

∂J/∂w = (2/n) · Xᵀ(ŷ − y)          ∂J/∂b = (2/n) · Σ(ŷ − y)

The intuitive reading: the gradient with respect to each weight is the mean error weighted by the corresponding feature. If the model falls short on long deliveries, the distance_km component of the gradient pushes its weight up. The error, quite literally, steers the learning.

Implementation with a descent trace

def gradient_descent(X, y, alpha=0.1, epochs=200):
    """Minimizes the MSE. Returns weights, bias and cost history."""
    n, d = X.shape
    w, b = np.zeros(d), 0.0                    # we start at the origin
    history = []
    for epoch in range(epochs):
        y_hat = X @ w + b                      # 1. current prediction
        error = y_hat - y                      # 2. residuals
        grad_w = (2 / n) * (X.T @ error)       # 3. gradient w.r.t. w
        grad_b = (2 / n) * error.sum()         #    ...and w.r.t. b
        w -= alpha * grad_w                    # 4. step downhill
        b -= alpha * grad_b
        history.append(np.mean(error ** 2))    # cost for this epoch
    return w, b, history

w_gd, b_gd, hist = gradient_descent(X_tr_s, y_tr, alpha=0.1, epochs=200)
for e in (0, 9, 49, 199):
    print(f"epoch {e:3d}: MSE = {hist[e]:8.2f}")

You will watch the MSE plummet in the first epochs and level off: that is convergence. Compare w_gd with the analytical solution: they should match to the third decimal — two radically different algorithms, the same minimum, because the paraboloid only has one.

The learning rate: the delicate piece

α Behavior Symptom in the cost history
Too small (0.0001) Converges, but after a huge number of epochs Decreases at an exasperating crawl
Adequate (≈0.01-0.3 here) Fast, stable convergence Smooth drop to a plateau
Too large (1.5) Each step leaps to the other side of the valley, farther every time The cost grows and diverges to infinity/NaN

Try alpha=1.5 and watch the explosion: it is the most common error in practice. The cost history is your electrocardiogram: look at it always. Note also why we standardize: with features on disparate scales, the surface J is a long, narrow valley and the descent zigzags; with standardized features the valley is round and the same α works for every direction.

Why this algorithm is THE algorithm

Gradient descent only requires the cost to be differentiable. Swap the MSE for the logistic regression cost: same four lines. Chain millions of parameters in a neural network: the same loop, with the gradient computed by backpropagation (05-04). All of modern deep learning is, in essence, this for loop with cleverly computed gradients. That is why we wrote it by hand.

Batch, stochastic and mini-batch

Our version uses all n deliveries at every step (batch). With Rutalia's real millions of records, each step would cost too much. Alternatives:

Variant Data per step Cost per step Trajectory Typical use
Batch All (n) O(n·d) Smooth, deterministic Small datasets
Stochastic (SGD) 1 example O(d) Very noisy, oscillates around the minimum Streaming, online
Mini-batch 32-256 examples O(B·d) Moderate noise The standard (and 05-04's choice)

SGD's noise is not only a defect: on non-convex surfaces (those of neural networks) it helps escape poor local minima. Mini-batch strikes the balance: a reasonably reliable gradient, bounded cost per step and matrix operations that exploit vectorization. Hold on to this concept: it reappears verbatim when training the network in 05-04.

def minibatch_descent(X, y, alpha=0.05, epochs=50, batch=64, seed=42):
    rng = np.random.default_rng(seed)
    n, d = X.shape
    w, b = np.zeros(d), 0.0
    for _ in range(epochs):
        order = rng.permutation(n)                 # shuffle every epoch
        for i in range(0, n, batch):
            chunk = order[i:i + batch]
            err = X[chunk] @ w + b - y[chunk]
            w -= alpha * (2 / len(chunk)) * (X[chunk].T @ err)
            b -= alpha * (2 / len(chunk)) * err.sum()
    return w, b

Regularization: Ridge and Lasso

With 13 features we have data to spare. But picture the real case: hundreds of features (one per zone×hour×day crossing, weather, customer type...). With many features and limited data, the linear model finds enormous weights that exploit spurious correlations in the training set: overfitting, again.

Regularization adds a penalty for large weights to the cost:

Method Cost Effect on the weights
Ridge (L2) MSE + λ·Σwⱼ² Shrinks them all smoothly toward 0; they never reach exactly 0
Lasso (L1) MSE + λ·Σ|wⱼ| Drives irrelevant weights to exactly 0: it selects features

λ controls the strength: λ=0 is ordinary regression; a huge λ crushes everything toward the trivial model. It is chosen with cross-validation (05-01), never with the test set.

Why does L1 produce exact zeros and L2 doesn't? The L2 penalty is a paraboloid (its gradient vanishes smoothly at 0: it pushes proportionally to the weight, less and less), whereas L1 has "corners": it pushes toward 0 with constant force λ, whatever the weight contributes, and weights that contribute less than λ end up pinned at zero. That is why Lasso works as an automatic feature selector: at Rutalia, with hundreds of candidate columns, it would tell you which ones truly matter.

For Ridge, adding it to gradient descent costs one line (grad_w += 2 * lam * w), and it even keeps an analytical solution: w = (XᵀX + λI)⁻¹Xᵀy. In sklearn: Ridge(alpha=...), Lasso(alpha=...) (careful: sklearn calls λ alpha, not the learning rate).

Polynomial regression and the bias-variance trade-off

Rutalia's truth is not linear (the factors multiply). An elegant trick: keep the model linear in the weights but enrich the features with derived terms: distance², distance·rush_hour, etc. The machinery (analytical or gradient) does not change at all: only X grows.

from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline

for degree in (1, 2, 5):
    model = make_pipeline(PolynomialFeatures(degree), LinearRegression())
    model.fit(X_tr_s, y_tr)
    mse_tr = np.mean((model.predict(X_tr_s) - y_tr) ** 2)
    mse_te = np.mean((model.predict(X_te_s) - y_te) ** 2)
    print(f"degree {degree}: MSE train={mse_tr:7.2f}  test={mse_te:7.2f}")

Degree 1 falls short; degree 2 captures the generator's multiplicative interactions and improves the test score; degree 5 nails the training set and worsens the test. It is the U-shaped curve from 05-01, now with a technical name. The generalization error decomposes into:

  • Bias: error from the model's rigidity — the straight line cannot represent the multiplicative truth (degree 1: high bias).
  • Variance: error from sensitivity to the particular sample — the degree-5 polynomial would change a lot if we regenerated the data with a different seed (high variance).
  • Irreducible noise: the generator's rng.normal(0, 4). No model will get below an MSE ≈ 16 (= 4²) on test; if one does on train, it is memorizing noise.

Capacity ↑ ⇒ bias ↓ and variance ↑: the optimum sits at the elbow. Everything we have seen to fight overfitting (k in k-NN, tree pruning, bagging, Ridge/Lasso's λ, polynomial degree) are knobs on the same dial.

A brief mention: regression trees are the trees of 05-02 with leaves that predict the mean of their examples and splits that minimize variance instead of Gini; the regression random forest averages their outputs. They capture nonlinearities without feature engineering and are the strong alternative to the polynomial approach.

Regression metrics: MSE, RMSE, MAE, R²

Metric Formula Units Reading
MSE mean of (ŷ−y)² minutes² The one we optimize; punishes large errors heavily
RMSE √MSE minutes The same information, in interpretable units
MAE mean of |ŷ−y| minutes "Typical" error; robust to outliers
R² 1 − MSE/var(y) dimensionless Fraction of variance explained; 1=perfect, 0=like predicting the mean, <0=worse than the mean 
def regression_metrics(y_true, y_pred):
    mse = np.mean((y_pred - y_true) ** 2)
    return {"MSE": mse, "RMSE": np.sqrt(mse),
            "MAE": np.mean(np.abs(y_pred - y_true)),
            "R2": 1 - mse / np.var(y_true)}

To talk to the business, use RMSE or MAE ("we're off by about 5 minutes per delivery"); to compare models against each other, R². If RMSE ≫ MAE, there are outliers: a few deliveries with huge errors that the square amplifies — investigate them before switching models.

Common Mistakes and Tips

  • Picking the learning rate blindly. If the cost diverges, α is too large; if it drops at a snail's pace, too small. Plot the cost history in every experiment: it is your number-one diagnostic tool.
  • Gradient descent without standardizing. With disparate scales, J's valley is narrow and elongated and a single α doesn't work for every direction: painfully slow convergence or divergence. Always standardize before gradient descent.
  • Regularizing the bias b. The penalty must apply only to the weights w; penalizing b artificially shifts the predictions. (sklearn already handles this correctly; in your manual implementation, mind that detail.)
  • Choosing λ or the polynomial degree by looking at the test set. It is the 05-01 leak in its most frequent disguise. Hyperparameters → cross-validation; test set → once, at the very end.
  • Reporting MSE to the business. "An MSE of 26 squared minutes" means nothing to operations. Translate to RMSE/MAE in minutes.
  • Tip: always keep two references: the trivial model (predicting the mean, R²=0) as the floor above and the irreducible noise as the floor below the error (here RMSE ≈ 4): every honest model lives between the two.

Exercises

  1. The descent's electrocardiogram. Run gradient_descent on the standardized canonical dataset with alpha ∈ {0.001, 0.1, 1.1} and 100 epochs. Print the MSE at epochs 0, 10, 50 and 99 for each α. Classify each behavior (slow / correct / divergent) and verify that with α=0.1 the weights match the analytical solution.

  2. Ridge against polynomial overfitting. Generate degree-5 polynomial features and compare LinearRegression() with Ridge(alpha=10) on train and test (MSE). How much does Ridge close the train-test gap? Also try Lasso(alpha=0.5) and count how many coefficients it leaves exactly at 0 (np.sum(model.coef_ == 0)).

  3. How much room is left? Compute test RMSE and R² for: (a) the trivial model that predicts the train mean, (b) your linear regression, (c) a RandomForestRegressor(n_estimators=200). Knowing that the generator's noise is σ=4, what RMSE is unattainable and how much does each model have left to squeeze?

Solutions

Exercise 1:

for a in (0.001, 0.1, 1.1):
    w_, b_, h = gradient_descent(X_tr_s, y_tr, alpha=a, epochs=100)
    print(f"α={a}: " + "  ".join(f"e{e}={h[e]:.1f}" for e in (0, 10, 50, 99)))

w_analytic = linear_regression_analytic(X_tr_s, y_tr)
w_gd, b_gd, _ = gradient_descent(X_tr_s, y_tr, alpha=0.1, epochs=2000)
print(np.allclose(w_analytic[1:], w_gd, atol=1e-2))   # True

α=0.001: the MSE decreases but at epoch 99 it is still far from the minimum (slow). α=0.1: rapid drop and plateau (correct). α=1.1: the MSE grows epoch after epoch until overflow (divergent: each step leaps past the other side of the valley). With enough epochs, α=0.1 reproduces the analytical solution: two paths, a single minimum, because the linear MSE is convex.

Exercise 2:

from sklearn.linear_model import Ridge, Lasso

pf = PolynomialFeatures(5)
Xp_tr, Xp_te = pf.fit_transform(X_tr_s), pf.transform(X_te_s)

for name, m in [("linear", LinearRegression()),
                ("ridge", Ridge(alpha=10)), ("lasso", Lasso(alpha=0.5))]:
    m.fit(Xp_tr, y_tr)
    tr = np.mean((m.predict(Xp_tr) - y_tr) ** 2)
    te = np.mean((m.predict(Xp_te) - y_te) ** 2)
    print(f"{name}: train={tr:.1f} test={te:.1f}")
print("zeroed coefs in lasso:", np.sum(Lasso(alpha=0.5).fit(Xp_tr, y_tr).coef_ == 0))

The degree-5 linear model shows a substantial train-test gap (overfitting); Ridge shrinks the weights and closes much of the gap while giving up a little train performance. Lasso, moreover, zeroes out the vast majority of the hundreds of polynomial coefficients: the surviving terms point to the interactions that truly matter (distance×hour, distance×congested zone).

Exercise 3:

from sklearn.ensemble import RandomForestRegressor

candidates = {
    "trivial": np.full(len(y_te), y_tr.mean()),
    "linear": predict_linear(linear_regression_analytic(X_tr_s, y_tr), X_te_s),
    "forest": RandomForestRegressor(n_estimators=200, random_state=42)
              .fit(X_tr, y_tr).predict(X_te),
}
for name, pred in candidates.items():
    m = regression_metrics(y_te, pred)
    print(f"{name}: RMSE={m['RMSE']:.2f}  R2={m['R2']:.3f}")

The trivial model marks the ceiling of ignorance (R²≈0, RMSE = standard deviation of y). No model can get below RMSE ≈ 4 (the σ of the generator's noise: information that simply doesn't exist in the features). The linear model lands halfway (bias: it can't represent the multiplicative truth); the forest, which does capture nonlinearities, gets much closer to the floor of 4. The distance from each RMSE to 4 is each model's real room for improvement.

Conclusion

This lesson leaves three layers firmly in place. Linear regression with its dual solution: the analytical one (orthogonal projection, normal equations, lstsq) and the iterative one. Gradient descent developed piece by piece — MSE cost, derivatives, learning rate, convergence history, batch/stochastic/mini-batch variants — and verified against the exact solution. And overfitting control with formal vocabulary: bias versus variance, with Ridge/Lasso regularization and the polynomial degree as knobs, and irreducible noise as the floor no honest model breaks through. But notice the limitation we are dragging along: to capture Rutalia's nonlinear truth we had to manufacture the features ourselves (polynomials, interactions) — knowledge supplied by hand, yet again. What if the model could also learn the transformations? That is exactly what a neural network is: layers of regressions (linear + activations) stacked, trained end to end with the same gradient descent you just implemented. The gradient, though, will have to be propagated backwards through the layers: backpropagation. See you in 05-04.

◄ Previous
◄
Next ►
►
The Project | About Us | Contribute | Donations | License
© Copyright 2026. All rights reserved