closed/form

import numpy as np

y = np.array([2.0, 3.5, 4.0, 5.5])   # targets

# Model: predict one constant c for every input.
# Loss:  MSE(c) = mean((y - c)^2)
# Set the derivative to zero:
#   d/dc MSE = -2 * mean(y - c) = 0   =>   c = mean(y)

c = y.mean()          # the closed-form answer
print(c)              # 3.75 -- no training loop needed

import torch

y = torch.tensor([2.0, 3.5, 4.0, 5.5])
c = torch.zeros(1, requires_grad=True)      # start from a guess
opt = torch.optim.SGD([c], lr=0.1)

for _ in range(200):                        # walk downhill instead
    loss = ((y - c) ** 2).mean()
    opt.zero_grad(); loss.backward(); opt.step()

print(c.item())       # ~3.75 -- converges to the mean

import numpy as np

X_train = np.array([[1200., 3], [1500., 4], [np.nan, 2], [900., 2]])
X_test  = np.array([[1100., 3]])

# 1) Impute with the TRAIN median
med = np.nanmedian(X_train[:, 0])
X_train[np.isnan(X_train[:, 0]), 0] = med

# 2) Standardize -- fit on TRAIN ONLY, apply everywhere
mu, sigma = X_train.mean(axis=0), X_train.std(axis=0)
X_train_z = (X_train - mu) / sigma
X_test_z  = (X_test  - mu) / sigma      # reuse train stats. No leakage.

# 3) One-hot encode a categorical column
city   = np.array([0, 2, 1, 0])         # ids for ["SEA", "NYC", "LA"]
onehot = np.eye(3)[city]                 # rows of the identity matrix

# 4) Log-transform a skewed feature
price     = np.array([300_000., 5_000_000., 450_000.])
price_log = np.log1p(price)              # log(1 + x): safe at zero

import torch
import torch.nn.functional as F

X_train = torch.tensor([[1200., 3], [1500., 4], [float('nan'), 2], [900., 2]])
X_test  = torch.tensor([[1100., 3]])

# 1) Impute with the TRAIN median
med = X_train[:, 0].nanmedian()
X_train[:, 0] = torch.nan_to_num(X_train[:, 0], nan=med.item())

# 2) Standardize with TRAIN stats
mu, sigma = X_train.mean(0), X_train.std(0)
X_train_z = (X_train - mu) / sigma
X_test_z  = (X_test  - mu) / sigma       # same rule: train stats only

# 3) One-hot -- this exact op feeds embedding layers later in DL
city   = torch.tensor([0, 2, 1, 0])
onehot = F.one_hot(city, num_classes=3).float()

# 4) Tame skew
price_log = torch.log1p(torch.tensor([3e5, 5e6, 4.5e5]))

# DL note: BatchNorm/LayerNorm are this chapter's idea,
# applied continuously INSIDE the network.

import numpy as np

def train_test_split(X, y, test_frac=0.2, seed=0):
    idx = np.random.default_rng(seed).permutation(len(y))
    cut = int(len(y) * (1 - test_frac))
    tr, te = idx[:cut], idx[cut:]
    return X[tr], X[te], y[tr], y[te]

def k_fold_indices(n, k=5, seed=0):
    idx = np.random.default_rng(seed).permutation(n)
    return np.array_split(idx, k)        # k disjoint validation folds

def accuracy(y, p):
    return (y == p).mean()

def precision_recall_f1(y, p):           # binary: 1 = positive
    tp = ((p == 1) & (y == 1)).sum()
    fp = ((p == 1) & (y == 0)).sum()
    fn = ((p == 0) & (y == 1)).sum()
    prec = tp / (tp + fp + 1e-12)
    rec  = tp / (tp + fn + 1e-12)
    f1   = 2 * prec * rec / (prec + rec + 1e-12)
    return prec, rec, f1

def r2_score(y, yhat):
    ss_res = ((y - yhat) ** 2).sum()
    ss_tot = ((y - y.mean()) ** 2).sum()
    return 1 - ss_res / ss_tot

import torch

def accuracy(y, p):
    return (y == p).float().mean()

def precision_recall_f1(y, p):
    tp = ((p == 1) & (y == 1)).sum().float()
    fp = ((p == 1) & (y == 0)).sum().float()
    fn = ((p == 0) & (y == 1)).sum().float()
    prec = tp / (tp + fp + 1e-12)
    rec  = tp / (tp + fn + 1e-12)
    f1   = 2 * prec * rec / (prec + rec + 1e-12)
    return prec, rec, f1

# In real projects: the `torchmetrics` package implements all of
# these (plus ROC-AUC) with GPU + batching support:
#
#   from torchmetrics.classification import BinaryF1Score
#   f1 = BinaryF1Score(); f1(preds, target)
#
# But you should be able to write each one from memory --
# interviewers love asking, and bugs here poison every experiment.

import numpy as np

rng = np.random.default_rng(0)
X = rng.uniform(0, 10, size=(100, 1))
y = 3.0 * X[:, 0] + 4.0 + rng.normal(0, 1.5, 100)   # truth: w=3, b=4

Xb = np.c_[X, np.ones(len(X))]        # append a bias column of 1s

# ---- Route A: closed form ----
w = np.linalg.lstsq(Xb, y, rcond=None)[0]   # solves min ||Xb w - y||^2
print(w)                                    # ~[3.0, 4.0], instantly

# ---- Route B: gradient descent on the same bowl ----
w, lr = np.zeros(2), 0.01
for _ in range(2000):
    err  = Xb @ w - y                  # (n,)  residuals
    grad = 2 * Xb.T @ err / len(y)     # "errors, pushed back through X"
    w   -= lr * grad
print(w)                               # same answer, iteratively

import torch, torch.nn as nn

X = torch.rand(100, 1) * 10
y = 3.0 * X + 4.0 + torch.randn(100, 1) * 1.5

model   = nn.Linear(1, 1)              # w and b live inside the module
loss_fn = nn.MSELoss()
opt     = torch.optim.SGD(model.parameters(), lr=0.01)

for _ in range(2000):
    loss = loss_fn(model(X), y)        # forward
    opt.zero_grad()
    loss.backward()                    # autograd computes X^T(err) for you
    opt.step()

print(model.weight.item(), model.bias.item())   # ~3.0, ~4.0

# Realize what this means: nn.Linear + MSELoss + SGD
# *is* linear regression. Your first "neural network"
# was this chapter all along.

import numpy as np

# ---- Ridge: closed form survives; just add lambda*I ----
def ridge_fit(X, y, lam=1.0):
    d = X.shape[1]
    return np.linalg.solve(X.T @ X + lam * np.eye(d), X.T @ y)

# ---- Lasso: coordinate descent + soft-thresholding ----
def soft(z, t):                        # shrink toward 0, clamp at 0
    return np.sign(z) * np.maximum(np.abs(z) - t, 0.0)

def lasso_fit(X, y, lam=0.1, iters=300):
    n, d = X.shape
    w = np.zeros(d)
    col_sq = (X ** 2).sum(axis=0)
    for _ in range(iters):
        for j in range(d):             # optimize one weight at a time
            r_j  = y - X @ w + X[:, j] * w[j]   # residual w/o feature j
            rho  = X[:, j] @ r_j
            w[j] = soft(rho, lam * n / 2) / col_sq[j]
    return w

# The soft() clamp is what produces EXACT zeros:
# if a feature's correlation with the residual is below the
# threshold, its weight snaps to 0 -- feature deleted.

import torch, torch.nn as nn

d = 10
model = nn.Linear(d, 1)
mse   = nn.MSELoss()

# ---- Ridge is hiding in an argument you've already used ----
opt = torch.optim.SGD(model.parameters(), lr=0.01,
                      weight_decay=1e-3)   # weight_decay == L2 penalty

# ---- L1 has no flag; add it to the loss yourself ----
lam = 1e-3
for _ in range(2000):
    loss = mse(model(X), y) + lam * model.weight.abs().sum()
    opt.zero_grad(); loss.backward(); opt.step()

# DL bridge: dropout, early stopping, and data augmentation are
# all this chapter's idea -- "make fitting the training set harder" --
# expressed without a penalty term. (Caveat for later: with Adam,
# true weight decay and L2 diverge; that's why AdamW exists.)

import numpy as np

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

def logreg_fit(X, y, lr=0.1, iters=3000):
    Xb = np.c_[X, np.ones(len(X))]      # bias trick again
    w  = np.zeros(Xb.shape[1])
    for _ in range(iters):
        p    = sigmoid(Xb @ w)          # predicted P(y=1), shape (n,)
        grad = Xb.T @ (p - y) / len(y)  # the tidy BCE gradient
        w   -= lr * grad
    return w

def predict_proba(X, w):
    return sigmoid(np.c_[X, np.ones(len(X))] @ w)

def predict(X, w, thresh=0.5):          # thresh is a business decision!
    return (predict_proba(X, w) >= thresh).astype(int)

# Multiclass? Replace sigmoid with softmax and BCE with
# cross-entropy; the gradient keeps the exact same X^T(p - y) shape.

import torch, torch.nn as nn

X = torch.randn(200, 2)
y = (X[:, 0] + X[:, 1] > 0).float().unsqueeze(1)   # a linear truth

model   = nn.Linear(2, 1)
loss_fn = nn.BCEWithLogitsLoss()   # sigmoid + BCE fused: numerically
opt     = torch.optim.SGD(model.parameters(), lr=0.1)   # stable in one op

for _ in range(3000):
    logits = model(X)               # NOTE: no sigmoid in the forward pass
    loss   = loss_fn(logits, y)
    opt.zero_grad(); loss.backward(); opt.step()

probs = torch.sigmoid(model(X))     # sigmoid only at inference time

# Multiclass version -- the pattern you already know from DL:
#   model = nn.Linear(d, num_classes)
#   loss  = nn.CrossEntropyLoss()   (takes raw logits, fuses softmax)
# Softmax regression IS the final layer of every classifier net.

import numpy as np

def knn_predict(X_train, y_train, X_new, k=5):
    # pairwise squared distances via broadcasting:
    # (m,1,d) - (1,n,d) -> (m,n,d) -> sum over d -> (m,n)
    d2 = ((X_new[:, None, :] - X_train[None, :, :]) ** 2).sum(-1)

    nn_idx = np.argsort(d2, axis=1)[:, :k]   # k closest per query row
    votes  = y_train[nn_idx]                 # (m, k) neighbor labels

    # majority vote per row
    return np.array([np.bincount(v).argmax() for v in votes])

# Regression variant: return y_train[nn_idx].mean(axis=1)
# Weighted variant: weight votes by 1/d -- closer neighbors count more.

# That (m,n) distance matrix is O(m*n*d) time and O(m*n) memory:
# fine for thousands of points, brutal for millions. Real systems
# use KD-trees / ball trees, or approximate NN (FAISS, HNSW) --
# the same tech behind vector databases and RAG retrieval.

import torch

def knn_predict(X_train, y_train, X_new, k=5):
    d   = torch.cdist(X_new, X_train)         # (m, n) pairwise distances
    idx = d.topk(k, largest=False).indices    # k nearest, no full sort
    votes = y_train[idx]                       # (m, k)
    return votes.mode(dim=1).values            # majority label per row

# Three habits worth stealing here:
#  * torch.cdist -- optimized pairwise distances, GPU-ready
#  * topk instead of full argsort -- O(n log k), not O(n log n)
#  * this exact pattern (embed -> cdist -> topk) is literally
#    how retrieval / kNN-augmented models work in modern DL.

import numpy as np

class GaussianNB:
    def fit(self, X, y):
        self.classes = np.unique(y)
        self.mu  = np.array([X[y == c].mean(0) for c in self.classes])
        self.var = np.array([X[y == c].var(0)  for c in self.classes])
        self.var += 1e-9                        # never divide by zero
        self.log_prior = np.log(np.bincount(y) / len(y))
        return self                             # "training" is done.

    def predict(self, X):
        # log N(x | mu, var) per class/sample/feature, then the naive
        # step: SUM over features (a product, in log space)
        ll = -0.5 * (np.log(2 * np.pi * self.var[:, None, :])
             + (X[None] - self.mu[:, None]) ** 2 / self.var[:, None, :])
        scores = ll.sum(-1) + self.log_prior[:, None]   # (C, n)
        return self.classes[scores.argmax(0)]

# Multinomial NB (spam): P(word|class) = word counts per class,
# plus Laplace smoothing (add 1) so unseen words don't zero out
# the whole product.

import torch
from torch.distributions import Normal

def gnb_fit(X, y, C):
    mu = torch.stack([X[y == c].mean(0) for c in range(C)])
    sd = torch.stack([X[y == c].std(0) + 1e-6 for c in range(C)])
    log_prior = torch.log(torch.bincount(y).float() / len(y))
    return mu, sd, log_prior

def gnb_predict(X, mu, sd, log_prior):
    # Normal(...).log_prob broadcasts to (C, n, d);
    # summing over features is the "naive" independence step
    ll = Normal(mu[:, None], sd[:, None]).log_prob(X[None]).sum(-1)
    return (ll + log_prior[:, None]).argmax(0)

# torch.distributions is a hidden gem: log_prob, sampling, KL --
# the same API powers VAEs and policy-gradient RL later.
# Also note: log-space arithmetic is the same trick behind
# log_softmax + NLLLoss in deep learning.

import numpy as np

def gini(y):
    p = np.bincount(y) / len(y)
    return 1 - (p ** 2).sum()

def best_split(X, y):
    best = (None, None, np.inf)             # (feature, threshold, score)
    for j in range(X.shape[1]):
        for t in np.unique(X[:, j]):
            L = X[:, j] <= t
            if L.all() or not L.any():
                continue
            score = (L.mean() * gini(y[L])          # weighted child
                     + (~L).mean() * gini(y[~L]))   # impurity
            if score < best[2]:
                best = (j, t, score)
    return best

def grow(X, y, depth=0, max_depth=4, min_leaf=3):
    if depth == max_depth or len(y) < 2 * min_leaf or gini(y) == 0:
        return {"leaf": np.bincount(y).argmax()}
    j, t, _ = best_split(X, y)
    if j is None:
        return {"leaf": np.bincount(y).argmax()}
    L = X[:, j] <= t
    return {"feat": j, "thresh": t,
            "left":  grow(X[L],  y[L],  depth+1, max_depth, min_leaf),
            "right": grow(X[~L], y[~L], depth+1, max_depth, min_leaf)}

def predict_one(node, x):
    while "leaf" not in node:
        go_left = x[node["feat"]] <= node["thresh"]
        node = node["left"] if go_left else node["right"]
    return node["leaf"]

# An honest note instead of a fake idiom:
#
# Trees do NOT train by gradient descent. Choosing a split is a
# discrete, non-differentiable search -- there's no gradient to
# follow, so there is no natural "PyTorch tree". In practice you
# reach for scikit-learn (single trees, forests) or
# XGBoost / LightGBM (boosted trees, Chapter 10).
#
# What DOES translate is the impurity math, on tensors:

import torch

def gini(y, num_classes):
    p = torch.bincount(y, minlength=num_classes).float() / len(y)
    return 1 - (p ** 2).sum()

def entropy(y, num_classes):
    p = torch.bincount(y, minlength=num_classes).float() / len(y)
    p = p[p > 0]                       # 0*log(0) := 0
    return -(p * p.log2()).sum()

# Deeper connection: entropy here IS the entropy in cross-entropy
# loss. Information gain = mutual information between a split
# and the label -- the same quantity, one field over.

import numpy as np
# uses grow() / predict_one() from Chapter 9

# ---- Bagging: many deep trees on bootstrap samples ----
def bagging_fit(X, y, n_trees=25, seed=0):
    rng = np.random.default_rng(seed)
    forest = []
    for _ in range(n_trees):
        idx = rng.integers(0, len(y), size=len(y))  # with replacement
        forest.append(grow(X[idx], y[idx], max_depth=8))
    return forest

def bagging_predict(forest, x):
    votes = np.array([predict_one(t, x) for t in forest])
    return np.bincount(votes).argmax()
# Random forest = the same loop + random feature subsets in best_split.

# ---- Gradient boosting (regression) with stumps ----
def fit_stump(X, r):                    # depth-1 regression tree
    best, best_err = None, np.inf
    for j in range(X.shape[1]):
        for t in np.unique(X[:, j]):
            L = X[:, j] <= t
            if L.all() or not L.any():
                continue
            lv, rv = r[L].mean(), r[~L].mean()
            err = ((r[L]-lv)**2).sum() + ((r[~L]-rv)**2).sum()
            if err < best_err:
                best_err, best = err, (j, t, lv, rv)
    return best

def stump_predict(s, X):
    j, t, lv, rv = s
    return np.where(X[:, j] <= t, lv, rv)

def boost_fit(X, y, rounds=100, lr=0.1):
    pred, stumps = np.full(len(y), y.mean()), []
    for _ in range(rounds):
        residual = y - pred              # what we still get wrong
        s = fit_stump(X, residual)       # small tree on the errors
        pred += lr * stump_predict(s, X) # humble corrective step
        stumps.append(s)
    return y.mean(), stumps

# The deep-learning bridge, made concrete: a "residual" IS the
# negative gradient of MSE with respect to current predictions.
import torch

y    = torch.tensor([3.0, 5.0, 8.0])
pred = torch.full((3,), y.mean().item(), requires_grad=True)

loss = 0.5 * ((y - pred) ** 2).sum()
loss.backward()

print(-pred.grad)     # == y - pred: the residuals, via autograd
print(y - pred)       # identical.

# So gradient boosting = gradient descent in FUNCTION space:
#   step direction: a tree fit to -dL/dF   (not -dL/dw!)
#   step size:      the learning rate nu
# Swap in BCE's gradient and the same loop does classification.
# This viewpoint is why it's called *gradient* boosting.

# Ensembling shows up in DL too, same variance logic:
#   * deep ensembles (train k nets, average) -- SOTA uncertainty
#   * dropout ~= an implicit ensemble of subnetworks
#   * model souping / EMA weights ~= averaging in weight space

import numpy as np

def svm_fit(X, y, C=1.0, lr=0.01, epochs=1000):
    """Linear soft-margin SVM by subgradient descent. y in {-1,+1}."""
    n, d = X.shape
    w, b = np.zeros(d), 0.0
    for _ in range(epochs):
        margins = y * (X @ w + b)
        viol = margins < 1            # inside the street (or wrong side)
        # subgradient of  0.5||w||^2 + C * sum(hinge)
        grad_w = w - C * (y[viol, None] * X[viol]).sum(0)
        grad_b =    - C * y[viol].sum()
        w -= lr * grad_w / n
        b -= lr * grad_b / n
    return w, b

def rbf_kernel(A, B, gamma=0.5):
    d2 = ((A[:, None] - B[None]) ** 2).sum(-1)
    return np.exp(-gamma * d2)

# Kernel SVMs are solved in the dual (SMO algorithm) -- that's what
# libsvm / sklearn.svm.SVC do. The prediction becomes a weighted
# vote over support vectors:
#   f(x) = sum_i alpha_i * y_i * K(x_i, x) + b

import torch

d = X.shape[1]
w = torch.zeros(d, requires_grad=True)
b = torch.zeros(1, requires_grad=True)
opt = torch.optim.SGD([w, b], lr=0.01)

lam = 1e-3                              # plays the role of 1/C
for _ in range(1000):
    margins = y * (X @ w + b)           # y in {-1, +1}
    hinge = torch.clamp(1 - margins, min=0).mean()
    loss  = lam * 0.5 * (w ** 2).sum() + hinge
    opt.zero_grad(); loss.backward(); opt.step()

# Notice what this is: a linear layer + a different loss function.
# swap hinge -> BCE and you have logistic regression;
# swap hinge -> MSE and you have linear regression.
# "Model architecture" and "loss" are independent dials.

# The kernel trick echoes in DL as: fixed kernel similarity vs.
# learned embeddings + dot products (attention!). Attention scores
# q . k are inner products in a *learned* feature space.

import numpy as np

def kmeans(X, k=3, iters=100, seed=0):
    rng = np.random.default_rng(seed)
    centers = X[rng.choice(len(X), k, replace=False)]
    for _ in range(iters):
        # assign: (n,1,d)-(1,k,d) -> (n,k) squared distances
        d2 = ((X[:, None] - centers[None]) ** 2).sum(-1)
        labels = d2.argmin(1)
        # update: each center -> mean of its members
        new = np.array([X[labels == j].mean(0) if (labels == j).any()
                        else centers[j]              # empty cluster guard
                        for j in range(k)])
        if np.allclose(new, centers):
            break                                    # converged
        centers = new
    return centers, labels

def inertia(X, centers, labels):
    return ((X - centers[labels]) ** 2).sum()        # J, for the elbow plot

# k-means++ seeding: pick the first center at random, then pick
# each next center with probability proportional to its squared
# distance from the nearest already-chosen center.

import torch

def kmeans(X, k=3, iters=100, seed=0):
    g = torch.Generator().manual_seed(seed)
    centers = X[torch.randperm(len(X), generator=g)[:k]]
    for _ in range(iters):
        labels = torch.cdist(X, centers).argmin(1)      # assign
        new = torch.stack([
            X[labels == j].mean(0) if (labels == j).any() else centers[j]
            for j in range(k)
        ])
        if torch.allclose(new, centers):
            break
        centers = new                                    # update
    return centers, labels

# Why bother with a torch version? Scale. cdist + argmin on GPU
# clusters millions of points fast -- this exact loop is used to
# build vector-DB indexes (IVF) and to quantize embeddings.
#
# GMMs in practice: sklearn.mixture.GaussianMixture. In torch,
# EM's E-step is Normal(...).log_prob + softmax over components --
# the same pieces you met in Chapter 8.

import numpy as np

def dbscan(X, eps=0.5, min_pts=5):
    """Teaching version: O(n^2) memory, honest logic."""
    n = len(X)
    d = np.sqrt(((X[:, None] - X[None]) ** 2).sum(-1))
    neighbors = [np.flatnonzero(d[i] <= eps) for i in range(n)]
    labels, cluster = np.full(n, -1), 0        # -1 = noise / unvisited

    for i in range(n):
        if labels[i] != -1 or len(neighbors[i]) < min_pts:
            continue                            # claimed, or not core
        labels[i] = cluster                     # new cluster seed
        queue = list(neighbors[i])              # flood-fill (BFS)
        while queue:
            j = queue.pop()
            if labels[j] == -1:
                labels[j] = cluster             # border or core joins
                if len(neighbors[j]) >= min_pts:
                    queue.extend(neighbors[j])  # core: keep expanding
        cluster += 1
    return labels                               # -1s that remain = noise

def single_linkage(X, k=2):
    """Naive agglomerative clustering down to k clusters."""
    clusters = [[i] for i in range(len(X))]
    d = np.sqrt(((X[:, None] - X[None]) ** 2).sum(-1))
    while len(clusters) > k:
        best, pair = np.inf, None
        for a in range(len(clusters)):
            for b in range(a + 1, len(clusters)):
                dist = min(d[i, j] for i in clusters[a]
                                   for j in clusters[b])
                if dist < best:
                    best, pair = dist, (a, b)
        a, b = pair
        clusters[a] += clusters.pop(b)          # merge closest pair
    return clusters
# Production: scipy.cluster.hierarchy.linkage + dendrogram plots.

import torch

# Both algorithms are built on one object you already know how
# to make fast: the pairwise distance matrix.
D = torch.cdist(X, X)                       # (n, n), GPU-ready

# DBSCAN's expensive part, vectorized:
eps, min_pts = 0.5, 5
neighbor_mask  = D <= eps                   # (n, n) boolean
neighbor_count = neighbor_mask.sum(1)       # per-point density
is_core        = neighbor_count >= min_pts  # core points in one line

# The flood-fill itself stays plain Python (it's inherently
# sequential), but with the mask precomputed on GPU it's cheap.

# Practical stack for clustering at scale:
#   sklearn.cluster.DBSCAN / AgglomerativeClustering  (CPU, solid)
#   HDBSCAN -- density clustering w/ varying densities; the modern
#              default for clustering *embeddings* from deep models
# Typical modern pipeline: encoder -> embeddings -> UMAP -> HDBSCAN.
# Classical clustering didn't die; it moved into latent space.

import numpy as np

def pca(X, n_components=2):
    Xc = X - X.mean(0)                    # center. ALWAYS center.
    U, S, Vt = np.linalg.svd(Xc, full_matrices=False)

    components = Vt[:n_components]         # principal directions (m, d)
    var  = S ** 2 / (len(X) - 1)           # variance along each PC
    ratio = var / var.sum()                # the "receipt"

    Z = Xc @ components.T                  # projected coords (n, m)
    return Z, components, ratio[:n_components]

def reconstruct(Z, components, mean):
    return Z @ components + mean           # back to d dims (lossy)

# Choosing m: cumulative ratio -- keep enough PCs to explain
# e.g. 95% of variance:
#   m = np.searchsorted(np.cumsum(ratio), 0.95) + 1

import torch

Xc = X - X.mean(0)

# Route 1: exact SVD (same math as the NumPy tab)
U, S, Vh = torch.linalg.svd(Xc, full_matrices=False)
Z = Xc @ Vh[:2].T

# Route 2: randomized -- fast for big matrices
U, S, V = torch.pca_lowrank(Xc, q=2)
Z = Xc @ V[:, :2]

# The DL bridge is exact, not an analogy:
# a LINEAR autoencoder with tied weights,
#     encode: z = W x      decode: x_hat = W^T z
# trained on reconstruction MSE, learns the PCA subspace.
# Nonlinear autoencoders (add ReLUs) generalize PCA to curved
# manifolds -- and embeddings from deep encoders are the modern
# answer to "reduce my dimensions". PCA is their linear ancestor,
# and still the first thing to try.

import numpy as np

def perceptron_fit(X, y, epochs=20):
    """Rosenblatt, 1958. y in {-1, +1}."""
    w, b = np.zeros(X.shape[1]), 0.0
    for _ in range(epochs):
        for xi, yi in zip(X, y):
            if yi * (xi @ w + b) <= 0:      # misclassified
                w += yi * xi                 # nudge toward correct side
                b += yi
    return w, b

# Guaranteed to converge IF the data is linearly separable.
# If not, it thrashes forever -- and that limitation (XOR!)
# froze neural-net research for a decade. The fix was depth:
# hidden layers make XOR trivially separable in learned space.

# You've already written the rest of the bridge yourself:
# Karpathy's micrograd (Zero to Hero, episode 1) is exactly
# "Chapter 4's gradient, applied recursively via the chain rule."

import torch.nn as nn

# Logistic regression IS a one-layer neural network:
clf = nn.Sequential(
    nn.Linear(d, 1),                # Chapter 6, verbatim
)                                    # + BCEWithLogitsLoss

# Add ONE hidden layer, and it stops being "classical":
mlp = nn.Sequential(
    nn.Linear(d, 32), nn.ReLU(),    # learned feature engineering
    nn.Linear(32, 1),               # ...then Chapter 6 again
)

# That's the whole difference. The first Linear+ReLU is doing the
# job you did by hand in Chapter 2 and via kernels in Chapter 11.

# When classical still wins (knowing this is a superpower):
#  * tabular data -- boosted trees beat nets shockingly often
#  * small data   -- nets overfit; ridge/NB/forests shine
#  * interpretability, latency, or a baseline in an afternoon
# Rule of thumb: always build the classical baseline FIRST.
# If your net can't beat ridge regression, the net is lying to you.