S085 · Linear Regression from Scratch (numpy)

Module M11: Classical ML · Session 85 of 130 · Track: ML 📈

What you'll be able to do after this session

• Explain Linear Regression from Scratch (numpy) to a friend in one minute, out loud, no notes.
• Recognise it when you see it in real code / systems / papers.
• Complete the hands-on exercise at the end without looking up help.

Prerequisites

If you can't answer these in 30 seconds each, redo the linked session first:

• S020 — Calculus II — Gradients & Gradient Descent from Scratch

Pre-read (skim before the session)

• 🎥 Intuition (5–15 min) · Linear Regression, Clearly Explained — StatQuest — Josh Starmer with a napkin and a fitted line.
• 🎥 Deep dive (30–60 min) · Gradient Descent, Step-by-Step — StatQuest — the algorithm that trains basically everything.
• 🎥 Hands-on demo (10–20 min) · Linear Regression from scratch in Python — Sentdex — write the class yourself, no sklearn.
• 📖 Canonical article · scikit-learn — Linear Regression — the reference implementation & assumptions.
• 📖 Book chapter / tutorial · ISLR ch.3 — Linear Regression (free PDF) — the gold-standard beginner treatment.
• 📖 Engineering blog · Uber Engineering — Forecasting at Uber — where simple linear models still power billion-dollar decisions.

(a) Intuition · 5 min

The one-paragraph mental model. Why this concept exists, what problem it solves, and the analogy that makes it stick.

In one sentence: The 'hello world' of ML.

Linear regression is the "hello world" of machine learning because it's the simplest possible model that still captures the whole ML workflow. The idea: draw the straight line (or hyperplane, in higher dimensions) that best fits your data. "Best" here has a precise meaning — the line that minimises the sum of squared vertical distances between the line and each point. That's it.

Why squared? Two reasons. First, squaring makes over- and under-predictions both count as positive errors (no cancellation). Second, and more magically, squaring gives us a smooth convex loss function — one big bowl with a single bottom. Convex means gradient descent can never get lost in a local minimum; it always finds the global best. That's a huge deal.

The equation is ŷ = w₀ + w₁·x₁ + w₂·x₂ + … + wₙ·xₙ. Each wᵢ is a weight the model learns. If x₁ is square feet and w₁ = 5000, then every extra sqft adds ₹5000 to the predicted price. Coefficients are interpretable — you can look at them and know what the model thinks. That interpretability is why linear regression is still deployed in banks, insurance, healthcare — anywhere a regulator will ask "why did the model say this?".

There are two ways to solve for w: the normal equation (w = (XᵀX)⁻¹Xᵀy, closed-form, exact, O(n·p²+p³)) and gradient descent (iterative, scales to millions of features). Small dataset → normal equation. Big dataset → gradient descent. Same answer.

Gotcha to carry forever: Linear regression assumes the relationship is linear in the parameters, not in the features. You can add x², log(x), cross-terms — as long as they're combined linearly with weights, it's still linear regression.

(b) Visual walkthrough · 15 min

Diagrams, worked examples, step-by-step visuals.

One data point, one line, one loss surface — the picture that explains 90% of ML:

Worked example — solve by hand. Data: (1, 2), (2, 4), (3, 5). Fit ŷ = w·x + b.

Mean of x = 2, mean of y = 11/3 ≈ 3.67.

w = Σ(xᵢ-x̄)(yᵢ-ȳ) / Σ(xᵢ-x̄)²

w = 3/2 = 1.5. b = ȳ - w·x̄ = 3.67 - 1.5·2 = 0.67.

Prediction at x=4: ŷ = 1.5·4 + 0.67 = 6.67.

Worked example 2 — the loss surface. For any (w, b), L(w, b) is the sum of squared errors. Plot it: it's a paraboloid — one bowl. The point at the bottom of the bowl is (w=1.5, b=0.67). Gradient descent walks downhill from wherever you start. In 2D you'd draw contour lines and arrows pointing to the centre.

Assumptions to remember: linearity, independence of errors, homoscedasticity (constant error variance), normality of errors (for confidence intervals). Violate them and predictions still work; the p-values stop working.

(c) Hands-on · 20 min

Code you type and run. Not pseudo-code — real, runnable, copy-pasteable.

Build linear regression from scratch (no sklearn for the fit), then verify with sklearn.

# pip install numpy scikit-learn matplotlib
import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error
 
rng = np.random.default_rng(42)
 
# ---------- 1. Load + scale (GD needs it) ----------
data = fetch_california_housing()
X, y = data.data, data.target
X = StandardScaler().fit_transform(X)
X = np.hstack([np.ones((len(X), 1)), X])   # add bias column
n, p = X.shape
print("X:", X.shape, "y:", y.shape)
 
# ---------- 2. Closed-form (Normal Equation) ----------
w_ne = np.linalg.solve(X.T @ X, X.T @ y)     # (XᵀX)w = Xᵀy   — never invert
rmse_ne = np.sqrt(mean_squared_error(y, X @ w_ne))
print(f"NORMAL EQN  rmse={rmse_ne:.4f}  w[:3]={w_ne[:3]}")
 
# ---------- 3. Batch gradient descent from scratch ----------
def gd(X, y, lr=0.05, epochs=200):
    w = rng.normal(size=X.shape[1]) * 0.01
    for ep in range(epochs):
        y_hat = X @ w
        error = y_hat - y                     # (n,)
        grad = (2 / len(y)) * X.T @ error     # (p,)
        w -= lr * grad
        if ep % 40 == 0:
            loss = np.mean(error ** 2)
            print(f"  ep={ep:3d}  loss={loss:.4f}  |w|={np.linalg.norm(w):.3f}")
    return w
 
w_gd = gd(X, y)
rmse_gd = np.sqrt(mean_squared_error(y, X @ w_gd))
print(f"GRADIENT DES rmse={rmse_gd:.4f}  w[:3]={w_gd[:3]}")
 
# ---------- 4. Verify against sklearn ----------
lr = LinearRegression(fit_intercept=False).fit(X, y)   # bias already in X
rmse_sk = np.sqrt(mean_squared_error(y, lr.predict(X)))
print(f"SKLEARN     rmse={rmse_sk:.4f}  w[:3]={lr.coef_[:3]}")
 
# ---------- 5. Interpret the coefficients ----------
names = ["bias"] + list(data.feature_names)
for name, coef in sorted(zip(names, w_ne), key=lambda t: -abs(t[1]))[:5]:
    print(f"  {name:12s}  weight = {coef:+.3f}")

Observe when you run it:

1. X.shape = (20640, 9) — 8 features + 1 bias column.
2. Normal equation and sklearn produce identical w to 4 decimals — same math.
3. Gradient descent converges to almost the same w after 200 epochs.
4. Loss drops smoothly — that's the convexity working. No wobble.
5. MedInc (median income) has the largest positive coefficient — richer areas, pricier houses. Makes sense.

Try this modification: change lr=0.05 to lr=2.0. The loss will explode (diverges). Then try lr=0.001 — the loss barely moves in 200 epochs (too slow). Learning-rate tuning is the single most important knob in gradient descent.

(d) Production reality · 10 min

How this shows up in real systems, gotchas, war stories, common bugs.

Where linear regression still eats. Fraud scoring in finance (regulators demand interpretability), pricing models in insurance, A/B test lift analysis, forecasting demand (Uber, DoorDash), attribution modelling in marketing, and every "what's the effect of X on Y controlling for Z" causal-inference question. If you can get away with linear regression, you should — it trains in milliseconds, generalises well, and you can explain it in a meeting.

Feature scaling matters for gradient descent, not for normal equation. GD converges way faster if all features are on similar scales. Sklearn's StandardScaler (zero mean, unit variance) is the default.

Multicollinearity — two features that are perfectly correlated (e.g., "sqft" and "sqm"). XᵀX becomes singular; normal equation blows up; GD converges to garbage weights that flip sign every retraining. Detect with the VIF (variance inflation factor) or a correlation matrix. Fix by dropping one feature or by regularising (S087 — Ridge/L2).

Outliers dominate MSE loss because errors are squared. One data-entry error saying a house is 100 million sqft can move all your weights. Fix: robust regression (Huber loss, RANSAC) or clean the data.

"Ecological fallacy" in production. Your model was trained on individual houses; someone plugs in a neighbourhood average and asks for a prediction. Predictions look reasonable but the interpretation is wrong. This is why data contracts matter.

When linear regression breaks. Non-linear relationships (add polynomial features or switch to tree-based models); high-dimensional sparse data (LR + L1 regularisation, aka Lasso); millions of examples with billions of features (Vowpal Wabbit or LinkedIn's Photon-ML).

Ops story. Uber's ETA models started as linear regressions per city and stayed that way for years before eventually moving to XGBoost and then deep learning. Sometimes the boring answer wins.

(e) Quiz + exercise · 10 min

1. Write the linear regression prediction equation and the MSE loss function.
2. What are the two ways to fit linear regression, and when would you use each?
3. Why must features be scaled before gradient descent, but not before the normal equation?
4. What is multicollinearity, how do you detect it, and what does it break?
5. Why is MSE loss dominated by outliers, and what's one alternative loss that isn't?

Stretch (connects to S020 — Gradient Descent): For a dataset with n=1M and p=100 features, estimate the memory and time cost of one epoch of batch gradient descent vs mini-batch (batch size 512). Which one does industry actually use, and why?

Explain-out-loud test

If you can't teach these 3 things to a friend without notes, redo the session:

1. What is Linear Regression from Scratch (numpy)? (one sentence, no jargon)
2. When would you use it? (one concrete example)
3. When would you NOT use it? (one anti-pattern)

What comes next

• Next session (S086): Logistic Regression — Sigmoid, Cross-Entropy, from Scratch
• Previous (S084): The ML Mental Model — Features, Labels, Train/Val/Test
• Hub: The 6-Month Learning Plan

Part of a 130-session evergreen learning series. Session structure: (a) intuition · (b) visual walkthrough · (c) hands-on · (d) production reality · (e) quiz + exercise. Duration: 60 minutes.