Supervised Learning: Derivations, Code, and Failure Modes

Supervised learning is the problem of finding a function f: X → Y given a finite sample of (x, y) pairs. Every algorithm in this file is a different answer to the same question: what assumptions should we make about f, and what does “finding” it mean mathematically? The assumptions determine the failure modes. Understanding the assumptions is the difference between applying ML and practicing it.

This file covers five algorithm families in depth. For each: the mathematical foundation, the derivation (not just the formula), a sklearn implementation plus a from-scratch implementation where instructive, when to use it, and where it breaks. The goal is that after working through this, you never treat any of these as a magic box again.

---

1. Linear Regression

Mathematical Foundation

Linear regression assumes the relationship between input X ∈ ℝⁿ and output y ∈ ℝ is:

y = Xθ + ε,  where ε ~ N(0, σ²I)

This Gaussian noise assumption is not decorative. It is the statement: “I believe the true relationship is linear, and all deviation from it is mean-zero iid Gaussian noise.” Maximum likelihood estimation under this assumption is exactly equivalent to minimizing mean squared error. This is why MSE is the “right” loss for linear regression — not because it’s convenient, but because it follows from the probabilistic model.

Derivation: Normal Equation

Minimize L(θ) = ‖Xθ - y‖²:

∂L/∂θ = 2Xᵀ(Xθ - y) = 0
Xᵀ Xθ = Xᵀ y
θ* = (XᵀX)⁻¹ Xᵀ y

This is the normal equation. It gives the exact solution. No iteration. When does it fail?

• When XᵀX is singular (features are linearly dependent → ridge regression adds λI)

• When n >> 1 (inverting an n×n matrix is O(n³) → use gradient descent)

import numpy as np

class LinearRegressionScratch:
    """
    Linear regression via normal equation and gradient descent.
    Both methods should converge to the same θ.
    """
    def fit_normal_equation(self, X, y):
        """O(n^3) - exact solution, fine for n < 10,000 features"""
        X_b = np.c_[np.ones((len(X), 1)), X]  # add bias
        self.theta = np.linalg.lstsq(X_b.T @ X_b, X_b.T @ y, rcond=None)[0]
        return self
    
    def fit_gradient_descent(self, X, y, lr=0.01, n_iters=1000):
        """O(n*d) per iteration - scales to large datasets"""
        X_b = np.c_[np.ones((len(X), 1)), X]
        m, n = X_b.shape
        self.theta = np.zeros(n)
        self.loss_history = []
        
        for i in range(n_iters):
            residuals = X_b @ self.theta - y
            gradient = (2/m) * X_b.T @ residuals
            self.theta -= lr * gradient
            self.loss_history.append(np.mean(residuals**2))
        return self
    
    def predict(self, X):
        X_b = np.c_[np.ones((len(X), 1)), X]
        return X_b @ self.theta

# Verification: both methods should give the same theta
np.random.seed(42)
X = np.random.randn(100, 3)
y = 3*X[:,0] - 2*X[:,1] + 1.5*X[:,2] + 0.5 + np.random.randn(100)*0.1

model = LinearRegressionScratch()
model.fit_normal_equation(X, y)
print("Normal eq theta:", model.theta)

model2 = LinearRegressionScratch()
model2.fit_gradient_descent(X, y, lr=0.1, n_iters=500)
print("GD theta:        ", model2.theta)
# Should be nearly identical

When to Use / Failure Modes

Situation	Linear Regression Behavior
Features linearly correlated (multicollinearity)	Unstable coefficients, high variance → use Ridge
Non-linear relationship	Underfits → use polynomial features or tree models
Outliers in y	MSE is non-robust → use Huber loss or quantile regression
n >> d (tall matrix)	Use gradient descent, not normal equation
d >> n (wide matrix)	Underdetermined → regularize (Lasso/Ridge)

---

2. Logistic Regression

The Probabilistic Frame (The Correct Way to Think About This)

Logistic regression is not a classification algorithm that happens to output probabilities. It is a discriminative probabilistic model that directly estimates P(Y=1|X). The classification decision is secondary.

Model: log(P(Y=1|x) / P(Y=0|x)) = θᵀx

This says: the log-odds of the positive class is a linear function of the input. Solving for P(Y=1|x):

P(Y=1|x) = σ(θᵀx) = 1 / (1 + exp(-θᵀx))

The sigmoid is not chosen arbitrarily. It is the unique function that maps a linear combination of inputs to a valid probability when the log-odds are linear.

Derivation: Maximum Likelihood

Given binary labels y_i ∈ {0,1}, the likelihood of the dataset under the model:

L(θ) = ∏ᵢ P(yᵢ|xᵢ,θ) = ∏ᵢ σ(θᵀxᵢ)^yᵢ · (1-σ(θᵀxᵢ))^(1-yᵢ)

Log-likelihood (easier to work with):

ℓ(θ) = Σᵢ [yᵢ log σ(θᵀxᵢ) + (1-yᵢ) log(1-σ(θᵀxᵢ))]

Maximizing log-likelihood = minimizing binary cross-entropy. This is the derivation that connects Phase 0 (MLE, KL divergence, cross-entropy) to this loss function. It is not a heuristic.

Gradient: ∂ℓ/∂θ = Σᵢ (yᵢ - σ(θᵀxᵢ)) xᵢ

No closed form — use gradient ascent (equivalently, gradient descent on negative log-likelihood).

class LogisticRegressionScratch:
    def __init__(self, lr=0.1, n_iters=1000, tol=1e-6):
        self.lr = lr
        self.n_iters = n_iters
        self.tol = tol
    
    @staticmethod
    def sigmoid(z):
        return 1 / (1 + np.exp(-np.clip(z, -500, 500)))
    
    def fit(self, X, y):
        X_b = np.c_[np.ones((len(X), 1)), X]
        m = len(y)
        self.theta = np.zeros(X_b.shape[1])
        self.loss_history = []
        
        for i in range(self.n_iters):
            p = self.sigmoid(X_b @ self.theta)
            
            # Binary cross-entropy loss
            loss = -np.mean(y * np.log(p + 1e-15) + (1-y) * np.log(1 - p + 1e-15))
            self.loss_history.append(loss)
            
            # Gradient of negative log-likelihood
            gradient = X_b.T @ (p - y) / m
            self.theta -= self.lr * gradient
            
            if i > 0 and abs(self.loss_history[-2] - loss) < self.tol:
                print(f"Converged at iteration {i}")
                break
        return self
    
    def predict_proba(self, X):
        X_b = np.c_[np.ones((len(X), 1)), X]
        return self.sigmoid(X_b @ self.theta)
    
    def predict(self, X, threshold=0.5):
        return (self.predict_proba(X) >= threshold).astype(int)

# Verify against sklearn
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import make_classification

X, y = make_classification(n_samples=500, n_features=5, random_state=42)
model_scratch = LogisticRegressionScratch(lr=0.5, n_iters=2000)
model_scratch.fit(X, y)

model_sk = LogisticRegression(solver='lbfgs', max_iter=2000)
model_sk.fit(X, y)

print(f"Scratch accuracy: {np.mean(model_scratch.predict(X) == y):.4f}")
print(f"sklearn accuracy: {np.mean(model_sk.predict(X) == y):.4f}")
# Should be close. sklearn uses L-BFGS (2nd-order), faster convergence.

Extension: Multinomial (Softmax) Regression

Multiclass case: output K probability scores summing to 1.

P(Y=k|x) = exp(θₖᵀx) / Σⱼ exp(θⱼᵀx)    # softmax

Loss: categorical cross-entropy = -Σₖ yₖ log P(Y=k|x)

This is exactly what a neural network’s final layer does. Logistic regression is a one-layer network with no hidden layers. This is the bridge.

---

3. Decision Trees

Why Decision Trees Matter Beyond Their Direct Use

Decision trees are weak learners on their own. Their value is threefold: (1) they are human-interpretable, (2) they are the building block of Random Forests and Gradient Boosting (the dominant tabular ML methods), (3) understanding their splitting criterion means understanding information theory in action.

Splitting Criterion: Information Gain

A node splits to maximize information gain — the reduction in entropy of the target variable.

IG(S, feature) = H(S) - Σⱼ (|Sⱼ|/|S|) · H(Sⱼ)

Where H(S) = -Σₖ pₖ log₂ pₖ is the Shannon entropy of the label distribution in node S.

Alternatively: Gini impurity = 1 - Σₖ pₖ² is faster to compute, used by sklearn by default. In practice, the two criteria produce nearly identical trees.

from collections import Counter

def entropy(y):
    """Shannon entropy of a label array"""
    if len(y) == 0:
        return 0
    counts = Counter(y)
    probs = np.array(list(counts.values())) / len(y)
    return -np.sum(probs * np.log2(probs + 1e-10))

def information_gain(y_parent, y_left, y_right):
    """Information gain from a binary split"""
    n = len(y_parent)
    return (entropy(y_parent) 
            - (len(y_left)/n) * entropy(y_left) 
            - (len(y_right)/n) * entropy(y_right))

def best_split(X, y):
    """Find the feature and threshold that maximizes information gain"""
    best_ig = -1
    best_feature, best_threshold = None, None
    
    for feature_idx in range(X.shape[1]):
        thresholds = np.unique(X[:, feature_idx])
        
        for threshold in thresholds:
            mask = X[:, feature_idx] <= threshold
            if mask.sum() == 0 or (~mask).sum() == 0:
                continue
            
            ig = information_gain(y, y[mask], y[~mask])
            if ig > best_ig:
                best_ig = ig
                best_feature = feature_idx
                best_threshold = threshold
    
    return best_feature, best_threshold, best_ig

class DecisionTreeNode:
    def __init__(self):
        self.feature = None
        self.threshold = None
        self.left = None
        self.right = None
        self.label = None  # for leaves

class DecisionTreeScratch:
    def __init__(self, max_depth=5, min_samples_split=2):
        self.max_depth = max_depth
        self.min_samples_split = min_samples_split
        self.root = None
    
    def _build(self, X, y, depth):
        node = DecisionTreeNode()
        
        # Stopping criteria
        if (depth >= self.max_depth or 
            len(y) < self.min_samples_split or
            len(np.unique(y)) == 1):
            node.label = Counter(y).most_common(1)[0][0]
            return node
        
        feature, threshold, ig = best_split(X, y)
        
        if ig == 0:
            node.label = Counter(y).most_common(1)[0][0]
            return node
        
        mask = X[:, feature] <= threshold
        node.feature = feature
        node.threshold = threshold
        node.left = self._build(X[mask], y[mask], depth + 1)
        node.right = self._build(X[~mask], y[~mask], depth + 1)
        return node
    
    def fit(self, X, y):
        self.root = self._build(X, y, depth=0)
        return self
    
    def _predict_one(self, x, node):
        if node.label is not None:
            return node.label
        if x[node.feature] <= node.threshold:
            return self._predict_one(x, node.left)
        return self._predict_one(x, node.right)
    
    def predict(self, X):
        return np.array([self._predict_one(x, self.root) for x in X])

When to Use / Failure Modes

Situation	Decision Tree Behavior
No regularization	Grows to memorize training data (depth → n → zero train error, catastrophic test error)
Linear decision boundaries	Uses axis-aligned splits → many levels needed; linear regression dominates
High-cardinality categoricals	Overfits easily; need ordinal encoding or target encoding
Missing values	Sklearn handles via surrogates; implement carefully from scratch

Key insight on max_depth: depth=1 is a “stump” (single split). Stumps are the base learner in many boosting algorithms. depth=∞ is full overfitting. The sweet spot is controlled by max_depth + min_samples_leaf.

---

4. Random Forests

The Core Idea: Variance Reduction via Decorrelated Trees

A single deep decision tree is high variance — it is sensitive to the specific training samples drawn. The bias-variance decomposition for squared loss:

E[(y - f̂(x))²] = Bias(f̂(x))² + Var(f̂(x)) + σ²

Random Forests reduce variance without increasing bias. How?

Bagging (Bootstrap Aggregating):

1. Draw B bootstrap samples (sampling with replacement) from training data

2. Train one tree on each bootstrap sample

3. Aggregate: majority vote (classification) or mean (regression)

If trees were independent with variance σ², the ensemble variance would be σ²/B → 0. But trees are correlated (all trained on similar data), so actual variance is:

Var(ensemble) = ρσ² + (1-ρ)σ²/B

Where ρ is the correlation between trees. As B → ∞, this approaches ρσ² — a floor determined by inter-tree correlation.

The key Random Forest innovation: at each split, consider only max_features randomly selected features (typically √d for classification, d/3 for regression). This decorrelates the trees (reduces ρ), making the ensemble floor lower.

from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt

X, y = load_iris(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Effect of n_estimators on test accuracy
n_trees_range = [1, 5, 10, 25, 50, 100, 200, 500]
train_scores, test_scores = [], []

for n_trees in n_trees_range:
    rf = RandomForestClassifier(n_estimators=n_trees, random_state=42)
    rf.fit(X_train, y_train)
    train_scores.append(rf.score(X_train, y_train))
    test_scores.append(rf.score(X_test, y_test))

# You will see: test accuracy improves and plateaus. Train accuracy stays ~1.0.
# This IS the bias-variance story: more trees reduce variance, bias unchanged.
print(f"Single tree:  test={test_scores[0]:.3f}, train={train_scores[0]:.3f}")
print(f"100 trees:    test={test_scores[6]:.3f}, train={train_scores[6]:.3f}")

# Feature importance: mean decrease in impurity across all trees
for name, imp in zip(['sl','sw','pl','pw'], rf.feature_importances_):
    print(f"  {name}: {imp:.3f}")

When to Use / Failure Modes

Situation	Random Forest Behavior
Tabular data, mixed types	Excellent baseline; hard to beat without boosting
Very high-dimensional sparse data (text)	Not competitive with linear models; prefer TF-IDF + LogReg
Extrapolation beyond training range	Fails — trees can’t extrapolate; they return leaf statistics
Need probability calibration	Out-of-bag probabilities are reasonable but may be overconfident
Interpretability required	Individual trees are interpretable; forest less so → use SHAP

---

5. Gradient Boosting (XGBoost / LightGBM)

What Gradient Boosting Is Actually Doing

Random Forests reduce variance by averaging parallel trees. Gradient Boosting reduces bias by building trees sequentially, each one correcting the errors of the ensemble so far.

The procedure for regression:

1. Initialize: F₀(x) = mean(y)

2. For m = 1..M: a. Compute residuals: rᵢ = yᵢ - Fₘ₋₁(xᵢ) b. Fit a tree hₘ(x) to the residuals c. Update: Fₘ(x) = Fₘ₋₁(x) + η · hₘ(x)

The residuals rᵢ are the negative gradient of MSE loss with respect to the current prediction. This is why it’s called gradient boosting — each tree fits the gradient of the loss, generalizable to any differentiable loss function.

For classification: fit trees to the negative gradient of log-loss. For ranking: fit to the gradient of a ranking loss. The algorithm is the same; only the loss changes.

XGBoost innovation (Chen & Guestrin, 2016): adds regularization directly to the tree-building objective (L1/L2 on leaf weights), second-order Taylor expansion of the loss (uses curvature for faster convergence), and hardware-level optimizations (cache-aware column access, SIMD vectorization). The 2016 paper is worth reading — it’s 15 pages and clearly written.

import xgboost as xgb
from sklearn.datasets import make_classification
from sklearn.model_selection import cross_val_score
import numpy as np

X, y = make_classification(n_samples=5000, n_features=20, n_informative=10, random_state=42)

xgb_model = xgb.XGBClassifier(
    n_estimators=200,
    max_depth=4,           # shallow trees are the norm for boosting
    learning_rate=0.05,    # η — low lr + more trees generally wins
    subsample=0.8,         # row subsampling (like RF's bagging)
    colsample_bytree=0.8,  # feature subsampling per tree
    reg_alpha=0.1,         # L1 regularization on leaf weights
    reg_lambda=1.0,        # L2 regularization on leaf weights
    eval_metric='logloss',
    use_label_encoder=False,
    random_state=42
)

cv_scores = cross_val_score(xgb_model, X, y, cv=5, scoring='roc_auc')
print(f"XGBoost AUC: {cv_scores.mean():.4f} ± {cv_scores.std():.4f}")

# Manual boosting: visualize how error decreases with each tree
xgb_model.fit(X, y, eval_set=[(X, y)], verbose=False)
# evals_result_ shows training loss per tree — you see it decreasing

Key Hyperparameter Intuition

Parameter	Effect	Common Range
n_estimators	More trees = lower bias (until overfitting)	100–5000
learning_rate	Lower = smoother, more trees needed	0.01–0.3
max_depth	Shallow trees → high bias, low variance	3–8
subsample	Row sampling → variance reduction (like RF)	0.6–1.0
colsample_bytree	Feature sampling per tree → decorrelation	0.5–1.0
reg_lambda	L2 regularization → shrinks leaf weights	1.0 (default)

Community-verified tuning truth: n_estimators=1000, learning_rate=0.05 with early stopping consistently outperforms n_estimators=100, learning_rate=0.3. Lower learning rate + more trees wins on almost every tabular benchmark.

---

6. Support Vector Machines (SVMs)

The Geometric Idea

SVMs find the hyperplane that maximizes the margin — the distance to the nearest training points (support vectors) on each side. The margin is 2/‖w‖, so maximizing margin = minimizing ‖w‖.

Primal problem:

minimize  ½‖w‖²
subject to  yᵢ(wᵀxᵢ + b) ≥ 1  for all i

The dual problem (via Lagrangians) is:

maximize  Σᵢ αᵢ - ½ Σᵢ Σⱼ αᵢ αⱼ yᵢ yⱼ xᵢᵀxⱼ
subject to  αᵢ ≥ 0, Σᵢ αᵢ yᵢ = 0

The key: the dual objective only depends on inner products xᵢᵀxⱼ. Replace this with a kernel function K(xᵢ, xⱼ) = φ(xᵢ)ᵀφ(xⱼ) and you can implicitly map to high-dimensional feature spaces. This is the kernel trick.

When to Use SVMs (and When to Stop)

SVMs were the dominant method pre-2012. Today:

• Still competitive: small-to-medium datasets (<100K samples), high-dimensional sparse data (text classification), when data is not abundant

• Outclassed by: XGBoost/LightGBM on tabular data, neural networks on images/text at scale

• Study value: the SVM derivation teaches constrained optimization, duality, kernel methods — concepts that appear in kernel PCA, Gaussian processes, and theoretical ML

from sklearn.svm import SVC
from sklearn.datasets import make_moons
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

# SVMs are sensitive to feature scale — always scale
X, y = make_moons(n_samples=500, noise=0.2, random_state=42)

svm_rbf = Pipeline([
    ('scaler', StandardScaler()),
    ('svm', SVC(kernel='rbf', C=1.0, gamma='scale', probability=True))
])
svm_rbf.fit(X, y)
print(f"RBF SVM accuracy: {svm_rbf.score(X, y):.4f}")

# C = regularization inverse: high C = less regularization = smaller margin
# gamma = RBF kernel bandwidth: high = narrow kernel = more complex boundary

---

Practice Problems

1. Derive Ridge Regression: Add L2 regularization ‖θ‖² to the MSE loss. Derive the closed-form solution. Show that it is equivalent to adding λI to XᵀX before inverting. Explain geometrically why this always makes XᵀX invertible.

2. Decision tree overfitting: Generate a toy 2D dataset with two classes. Train a decision tree with max_depth=None. Plot the decision boundary. Now add max_depth=3. Compare. Write two sentences explaining why the deep tree memorizes training data in terms of entropy and information gain.

3. Bootstrap variance: Write a function that estimates the variance of the sample mean using bootstrapping (1000 bootstrap samples). Compare to the analytical formula Var(X̄) = σ²/n. The gap between these two is the actual bias of the bootstrap estimator of variance.

4. Gradient in gradient boosting: For log-loss L(y, F) = -y·log(σ(F)) - (1-y)·log(1-σ(F)), compute ∂L/∂F. Show that it equals σ(F) - y, the prediction residual. This is why gradient boosting with log-loss fits trees to (y - predicted_probability).

---

Return to README.md · Next: 02_unsupervised_learning.md