Neural Networks from Scratch

Imagine teaching a child to recognise cats. You could write down rules: "furry, four legs, pointed ears, whiskers." That is classical ML — hand-crafted features. Or you could show the child 10,000 photos of cats and non-cats and let them figure out the pattern themselves. They learn features you never named — the specific curve of an ear, the texture of fur, the shape of eyes. That is a neural network.

The child's brain adjusts internal connections after each photo — strengthening what was right, weakening what was wrong. A neural network does exactly this: adjust weights after each prediction based on how wrong it was. That adjustment process is backpropagation. The rule for how much to adjust is gradient descent.

import numpy as np

# ── Activation functions and their derivatives ─────────────────────────
# Derivatives needed for backpropagation

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

def sigmoid_derivative(z):
    s = sigmoid(z)
    return s * (1 - s)   # convenient: derivative in terms of the output

def relu(z):
    return np.maximum(0, z)

def relu_derivative(z):
    return (z > 0).astype(float)   # 1 if z > 0, 0 otherwise

def tanh(z):
    return np.tanh(z)

def tanh_derivative(z):
    return 1 - np.tanh(z) ** 2

def softmax(z):
    # Numerically stable: subtract max before exponentiating
    e = np.exp(z - z.max(axis=1, keepdims=True))
    return e / e.sum(axis=1, keepdims=True)

# ── Verify derivatives numerically ─────────────────────────────────────
# Numerical gradient: (f(z+h) - f(z-h)) / (2h)
def numerical_gradient(f, z, h=1e-5):
    return (f(z + h) - f(z - h)) / (2 * h)

z_test = np.array([[-2.0, -0.5, 0.0, 0.5, 2.0]])

print("Verifying analytical vs numerical derivatives:")
for name, f, df in [
    ('sigmoid', sigmoid, sigmoid_derivative),
    ('relu',    relu,    relu_derivative),
    ('tanh',    tanh,    tanh_derivative),
]:
    analytical = df(z_test)
    numerical  = numerical_gradient(f, z_test)
    max_error  = np.abs(analytical - numerical).max()
    print(f"  {name:<10}: max error = {max_error:.2e}  "
          f"{'✓ correct' if max_error < 1e-4 else '✗ wrong'}")

# ── Dying ReLU demonstration ───────────────────────────────────────────
# Once relu output is 0, its gradient is 0 — the neuron never updates
z_negative = np.array([-5.0, -3.0, -1.0])
print(f"
ReLU on negative inputs: {relu(z_negative)}")
print(f"ReLU derivative:         {relu_derivative(z_negative)}  ← all zero, no learning")

import numpy as np

np.random.seed(42)

# ── Build a 3-layer network for DoorDash delivery time prediction ────────
# Architecture: 4 inputs → 8 hidden → 4 hidden → 1 output
# Task: regression (predict delivery minutes)

def sigmoid(z):    return 1 / (1 + np.exp(-np.clip(z, -500, 500)))
def relu(z):       return np.maximum(0, z)
def relu_deriv(z): return (z > 0).astype(float)

# ── Initialise weights ─────────────────────────────────────────────────
# He initialisation for ReLU: scale by sqrt(2/n_inputs)
# Avoids vanishing/exploding gradients at the start
def init_weights(n_in, n_out, activation='relu'):
    scale = np.sqrt(2.0 / n_in) if activation == 'relu' else np.sqrt(1.0 / n_in)
    return np.random.randn(n_in, n_out) * scale

# Network parameters
W1 = init_weights(4, 8)       # layer 1 weights:  (4, 8)
b1 = np.zeros((1, 8))          # layer 1 bias:     (1, 8)
W2 = init_weights(8, 4)        # layer 2 weights:  (8, 4)
b2 = np.zeros((1, 4))          # layer 2 bias:     (1, 4)
W3 = init_weights(4, 1, 'linear')  # output weights:   (4, 1)
b3 = np.zeros((1, 1))          # output bias:      (1, 1)

total_params = (W1.size + b1.size + W2.size + b2.size + W3.size + b3.size)
print(f"Network architecture: 4 → 8 → 4 → 1")
print(f"Total parameters:     {total_params}")

# ── Forward pass — explicit at every step ─────────────────────────────
def forward_pass(X):
    """
    X: (batch_size, 4) — [distance, traffic, prep_time, order_value]
    Returns prediction and all intermediate values (needed for backprop)
    """
    # Layer 1: linear → ReLU
    Z1 = X  @ W1 + b1    # (batch, 4) @ (4, 8) + (1, 8) = (batch, 8)
    A1 = relu(Z1)          # (batch, 8)

    # Layer 2: linear → ReLU
    Z2 = A1 @ W2 + b2    # (batch, 8) @ (8, 4) + (1, 4) = (batch, 4)
    A2 = relu(Z2)          # (batch, 4)

    # Output layer: linear (no activation for regression)
    Z3 = A2 @ W3 + b3    # (batch, 4) @ (4, 1) + (1, 1) = (batch, 1)
    A3 = Z3                # raw output for regression

    # Cache everything — backprop needs these
    cache = {'X': X, 'Z1': Z1, 'A1': A1, 'Z2': Z2, 'A2': A2, 'Z3': Z3, 'A3': A3}
    return A3, cache

# ── Test on a batch ───────────────────────────────────────────────────
np.random.seed(0)
n_batch = 32
X_batch = np.column_stack([
    np.abs(np.random.normal(4.0, 2.0, n_batch)).clip(0.5, 15),  # distance
    np.random.randint(1, 11, n_batch).astype(float),              # traffic
    np.abs(np.random.normal(15, 5, n_batch)).clip(5, 35),        # prep
    np.abs(np.random.normal(350, 150, n_batch)).clip(50, 1200),  # order value
])
# Standardise
X_mean  = X_batch.mean(axis=0)
X_std   = X_batch.std(axis=0) + 1e-8
X_batch = (X_batch - X_mean) / X_std

pred, cache = forward_pass(X_batch)
print(f"
Forward pass on batch of {n_batch}:")
print(f"  Input shape:  {X_batch.shape}")
print(f"  Z1 shape:     {cache['Z1'].shape}")
print(f"  A1 shape:     {cache['A1'].shape}")
print(f"  Z2 shape:     {cache['Z2'].shape}")
print(f"  Output shape: {pred.shape}")
print(f"  Sample predictions (unnormalised): {pred[:5, 0].round(3)}")

import numpy as np

np.random.seed(42)

def sigmoid(z):     return 1 / (1 + np.exp(-np.clip(z, -500, 500)))
def relu(z):        return np.maximum(0, z)
def relu_deriv(z):  return (z > 0).astype(float)
def init_w(n_in, n_out): return np.random.randn(n_in, n_out) * np.sqrt(2.0 / n_in)

W1 = init_w(4, 8); b1 = np.zeros((1, 8))
W2 = init_w(8, 4); b2 = np.zeros((1, 4))
W3 = init_w(4, 1); b3 = np.zeros((1, 1))

def forward(X):
    Z1 = X  @ W1 + b1; A1 = relu(Z1)
    Z2 = A1 @ W2 + b2; A2 = relu(Z2)
    Z3 = A2 @ W3 + b3; A3 = Z3   # linear output
    return A3, {'X':X,'Z1':Z1,'A1':A1,'Z2':Z2,'A2':A2,'Z3':Z3,'A3':A3}

# ── Backpropagation — every gradient explicit ──────────────────────────
def backward(y_true, cache):
    """
    Compute gradients of MSE loss w.r.t. all weights and biases.
    Every line is one application of the chain rule.
    """
    n   = y_true.shape[0]
    A3  = cache['A3']

    # ── Loss gradient ──────────────────────────────────────────────────
    # MSE: L = mean((A3 - y)^2)
    # dL/dA3 = 2/n * (A3 - y)
    dA3 = 2 / n * (A3 - y_true)          # (batch, 1)

    # ── Output layer gradients ─────────────────────────────────────────
    # Z3 = A2 @ W3 + b3
    # dL/dW3 = A2^T @ dA3
    # dL/db3 = sum(dA3, axis=0)
    # dL/dA2 = dA3 @ W3^T  (pass back to layer 2)
    dW3 = cache['A2'].T @ dA3             # (4, 1)
    db3 = dA3.sum(axis=0, keepdims=True)  # (1, 1)
    dA2 = dA3 @ W3.T                      # (batch, 4)

    # ── Layer 2 gradients — through ReLU ──────────────────────────────
    # A2 = relu(Z2), so dL/dZ2 = dL/dA2 * relu'(Z2)
    dZ2 = dA2 * relu_deriv(cache['Z2'])   # (batch, 4)
    dW2 = cache['A1'].T @ dZ2             # (8, 4)
    db2 = dZ2.sum(axis=0, keepdims=True)  # (1, 4)
    dA1 = dZ2 @ W2.T                      # (batch, 8)

    # ── Layer 1 gradients — through ReLU ──────────────────────────────
    dZ1 = dA1 * relu_deriv(cache['Z1'])   # (batch, 8)
    dW1 = cache['X'].T @ dZ1              # (4, 8)
    db1 = dZ1.sum(axis=0, keepdims=True)  # (1, 8)

    return {'dW1':dW1,'db1':db1,'dW2':dW2,'db2':db2,'dW3':dW3,'db3':db3}

# ── Verify gradients numerically (gradient check) ─────────────────────
# This is the gold standard test for backprop correctness
def mse_loss(y_pred, y_true):
    return np.mean((y_pred - y_true) ** 2)

n_check = 8
X_c = np.random.randn(n_check, 4)
y_c = np.random.randn(n_check, 1) * 10 + 35   # delivery times

pred_c, cache_c = forward(X_c)
grads = backward(y_c, cache_c)

print("Gradient check (analytical vs numerical):")
h = 1e-5
for name, W, dW in [('W1', W1, grads['dW1']), ('W3', W3, grads['dW3'])]:
    # Check a random element
    i, j      = 0, 0
    W_orig    = W[i, j]

    W[i, j]   = W_orig + h
    loss_plus, _ = forward(X_c); loss_plus = mse_loss(loss_plus, y_c)
    W[i, j]   = W_orig - h
    loss_minus, _ = forward(X_c); loss_minus = mse_loss(loss_minus, y_c)
    W[i, j]   = W_orig   # restore

    numerical  = (loss_plus - loss_minus) / (2 * h)
    analytical = dW[i, j]
    error      = abs(numerical - analytical) / (abs(numerical) + abs(analytical) + 1e-10)
    print(f"  {name}[{i},{j}]: analytical={analytical:.6f}  "
          f"numerical={numerical:.6f}  rel_error={error:.2e}  "
          f"{'✓' if error < 1e-4 else '✗'}")

import numpy as np
from sklearn.model_selection import train_test_split

np.random.seed(42)
n = 2000

# DoorDash delivery time dataset
distance  = np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15)
traffic   = np.random.randint(1, 11, n).astype(float)
prep      = np.abs(np.random.normal(15, 5, n)).clip(5, 35)
value     = np.abs(np.random.normal(350, 150, n)).clip(50, 1200)
delivery  = (8.6 + 7.3*distance + 0.8*prep + 1.5*traffic
             + np.random.normal(0, 4, n)).clip(10, 120)

X = np.column_stack([distance, traffic, prep, value])
y = delivery.reshape(-1, 1)

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)

# Standardise
X_mean, X_std = X_tr.mean(0), X_tr.std(0) + 1e-8
y_mean, y_std = y_tr.mean(),  y_tr.std()  + 1e-8
X_tr_sc = (X_tr - X_mean) / X_std
X_te_sc = (X_te - X_mean) / X_std
y_tr_sc = (y_tr - y_mean) / y_std
y_te_sc = (y_te - y_mean) / y_std

# ── Full neural network class ──────────────────────────────────────────
class NeuralNetwork:
    def __init__(self, layer_sizes, learning_rate=0.01):
        self.lr   = learning_rate
        self.W    = []
        self.b    = []
        for i in range(len(layer_sizes) - 1):
            n_in, n_out = layer_sizes[i], layer_sizes[i+1]
            self.W.append(np.random.randn(n_in, n_out) * np.sqrt(2.0 / n_in))
            self.b.append(np.zeros((1, n_out)))

    def relu(self, z):      return np.maximum(0, z)
    def relu_d(self, z):    return (z > 0).astype(float)
    def mse(self, yp, yt):  return np.mean((yp - yt) ** 2)
    def mae(self, yp, yt):  return np.mean(np.abs(yp - yt))

    def forward(self, X):
        self.cache = [X]
        self.Z     = []
        A = X
        for i, (W, b) in enumerate(zip(self.W, self.b)):
            Z = A @ W + b
            self.Z.append(Z)
            A = self.relu(Z) if i < len(self.W) - 1 else Z   # linear output
            self.cache.append(A)
        return A

    def backward(self, y_true):
        n   = y_true.shape[0]
        dA  = 2 / n * (self.cache[-1] - y_true)
        dWs, dbs = [], []
        for i in reversed(range(len(self.W))):
            dZ = dA * self.relu_d(self.Z[i]) if i < len(self.W) - 1 else dA
            dWs.insert(0, self.cache[i].T @ dZ)
            dbs.insert(0, dZ.sum(axis=0, keepdims=True))
            dA = dZ @ self.W[i].T
        # Update weights
        for i in range(len(self.W)):
            self.W[i] -= self.lr * dWs[i]
            self.b[i]  -= self.lr * dbs[i]

    def fit(self, X, y, epochs=200, batch_size=64, verbose=True):
        history = []
        for epoch in range(epochs):
            # Shuffle data each epoch
            idx = np.random.permutation(len(X))
            X, y = X[idx], y[idx]
            # Mini-batch gradient descent
            for start in range(0, len(X), batch_size):
                Xb = X[start:start+batch_size]
                yb = y[start:start+batch_size]
                self.forward(Xb)
                self.backward(yb)
            # Epoch loss
            loss = self.mse(self.forward(X), y)
            history.append(loss)
            if verbose and (epoch + 1) % 50 == 0:
                print(f"  Epoch {epoch+1:3d}: MSE = {loss:.6f}")
        return history

# ── Train ─────────────────────────────────────────────────────────────
nn = NeuralNetwork(layer_sizes=[4, 16, 8, 1], learning_rate=0.005)
print("Training neural network on DoorDash delivery data:")
history = nn.fit(X_tr_sc, y_tr_sc, epochs=200, batch_size=64)

# Evaluate
y_pred_sc = nn.forward(X_te_sc)
y_pred    = y_pred_sc * y_std + y_mean   # un-standardise

from sklearn.metrics import mean_absolute_error
mae = mean_absolute_error(y_te, y_pred)
print(f"
Test MAE: {mae:.4f} minutes")
print(f"(Baseline — always predict mean: "
      f"{mean_absolute_error(y_te, np.full_like(y_te, y_tr.mean())):.4f} min)")

import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)
torch.manual_seed(42)
n = 2000

distance  = np.abs(np.random.normal(4.0, 2.0, n)).clip(0.5, 15)
traffic   = np.random.randint(1, 11, n).astype(float)
prep      = np.abs(np.random.normal(15, 5, n)).clip(5, 35)
value     = np.abs(np.random.normal(350, 150, n)).clip(50, 1200)
delivery  = (8.6 + 7.3*distance + 0.8*prep + 1.5*traffic
             + np.random.normal(0, 4, n)).clip(10, 120)

X = np.column_stack([distance, traffic, prep, value])
y = delivery.reshape(-1, 1)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)

X_mean, X_std = X_tr.mean(0), X_tr.std(0) + 1e-8
y_mean, y_std = y_tr.mean(),  y_tr.std()  + 1e-8
X_tr_sc = (X_tr - X_mean) / X_std
X_te_sc = (X_te - X_mean) / X_std
y_tr_sc = (y_tr - y_mean) / y_std

# Convert to PyTorch tensors
X_train_t = torch.FloatTensor(X_tr_sc)
y_train_t = torch.FloatTensor(y_tr_sc)
X_test_t  = torch.FloatTensor(X_te_sc)

# ── Define network — identical architecture to our from-scratch version ─
model = nn.Sequential(
    nn.Linear(4, 16),    # layer 1: 4 inputs → 16 neurons
    nn.ReLU(),
    nn.Linear(16, 8),    # layer 2: 16 → 8
    nn.ReLU(),
    nn.Linear(8, 1),     # output:  8 → 1 (regression)
)

criterion = nn.MSELoss()
optimizer = optim.SGD(model.parameters(), lr=0.005)

# ── Training loop — forward + backward + update ────────────────────────
print("Training with PyTorch:")
batch_size = 64
dataset    = torch.utils.data.TensorDataset(X_train_t, y_train_t)
dataloader = torch.utils.data.DataLoader(dataset, batch_size=batch_size, shuffle=True)

for epoch in range(200):
    model.train()
    for X_batch, y_batch in dataloader:
        # 1. Forward pass — PyTorch traces all operations
        y_pred = model(X_batch)

        # 2. Compute loss
        loss = criterion(y_pred, y_batch)

        # 3. Backward pass — autograd computes all gradients
        optimizer.zero_grad()   # clear previous gradients
        loss.backward()          # ← replaces our entire backward() function

        # 4. Update weights
        optimizer.step()        # W -= lr * W.grad  for all parameters

    if (epoch + 1) % 50 == 0:
        model.eval()
        with torch.no_grad():
            train_loss = criterion(model(X_train_t), y_train_t).item()
        print(f"  Epoch {epoch+1:3d}: MSE = {train_loss:.6f}")

# ── Evaluate ───────────────────────────────────────────────────────────
model.eval()
with torch.no_grad():
    y_pred_sc = model(X_test_t).numpy()
y_pred = y_pred_sc * y_std + y_mean

mae = mean_absolute_error(y_te, y_pred)
print(f"
Test MAE: {mae:.4f} minutes")

# ── What PyTorch automates ─────────────────────────────────────────────
print(f"
What PyTorch replaces:")
print(f"  loss.backward() → our 30-line backward() function")
print(f"  optimizer.step() → our weight update loop")
print(f"  nn.Linear → our W @ x + b")
print(f"  autograd → the chain rule applied automatically to any graph")

import torch

# ── A custom loss for imbalanced return prediction ─────────────────────
# Production scenario: a retailer's return-prediction model needs a
# custom weighted loss so false negatives (missed returns) cost more
# than false positives. An engineer writes a custom autograd.Function
# for speed instead of relying on autograd to differentiate a Python loop.

class WeightedBCE(torch.autograd.Function):
    @staticmethod
    def forward(ctx, logits, targets, pos_weight):
        probs = torch.sigmoid(logits)
        ctx.save_for_backward(probs, targets, pos_weight)
        loss = -(pos_weight * targets * torch.log(probs + 1e-8)
                  + (1 - targets) * torch.log(1 - probs + 1e-8))
        return loss.mean()

    @staticmethod
    def backward(ctx, grad_output):
        probs, targets, pos_weight = ctx.saved_tensors
        # Gradient of weighted BCE w.r.t. the logits (not the probabilities!)
        # A common bug: differentiating w.r.t. probs and forgetting the
        # extra sigmoid-derivative term needed to get back to the logits.
        n = targets.numel()
        grad_logits = (pos_weight * targets * (probs - 1)
                       + (1 - targets) * probs) / n
        return grad_output * grad_logits, None, None

# ── Exactly the check this module teaches, run on the real op ─────────
torch.manual_seed(0)
logits  = torch.randn(8, dtype=torch.float64, requires_grad=True)
targets = torch.randint(0, 2, (8,), dtype=torch.float64)
weight  = torch.tensor(3.0, dtype=torch.float64)

# gradcheck compares the custom backward() to a numerical finite-difference
# gradient — the exact same idea as this module's numerical_gradient(),
# just wrapped as a standard PyTorch utility teams run before shipping
# any custom autograd.Function.
ok = torch.autograd.gradcheck(
    lambda l: WeightedBCE.apply(l, targets, weight),
    (logits,),
    eps=1e-6, atol=1e-4,
)
print(f"Custom backward matches numerical gradient: {ok}")
print("If this were False, the model would still train — just on a")
print("silently wrong gradient, degrading accuracy with no error raised.")

The weighted-sum-plus-activation computation was loosely inspired by neuroscience in the 1940s and 1950s, but modern neural networks are not modelling how real neurons work — real neurons spike in time, involve thousands of distinct neurotransmitter mechanisms, and do not compute a single differentiable scalar. Treating the biological analogy as literal leads people to expect properties (adaptability, energy efficiency, one-shot learning) that artificial neurons simply do not have. The useful mental model is the one this module actually uses: a neuron is z = Σwᵢxᵢ + b followed by a non-linearity — nothing more, nothing biological required to reason about it correctly.

"Learns its own features" means the network does not need you to hand-craft what a feature looks like — it does not mean every design decision is automatic. You still choose the number of layers, the width of each layer, which activation to use where, how to standardise inputs, how to initialise weights, and the batch size. Get any of these wrong — as the errors section of this module shows with NaN losses and dead networks — and the network learns nothing at all, no matter how good the raw data is. Feature learning replaces manual feature engineering; it does not replace architecture and preprocessing decisions.

Without a bias term, every neuron's pre-activation z = Σwᵢxᵢ is forced to equal zero whenever every input is zero — geometrically, the decision boundary of every layer is forced to pass through the origin. For real data, where the useful separation between classes is rarely centred exactly at the origin, this is a serious restriction, not a cosmetic one. That is why this module initialises biases to zero but weights to a scaled random distribution — the bias needs to be free to shift away from zero during training, it just does not need a random starting point to break symmetry the way weights do.

More capacity only helps if the dataset has enough signal to constrain it. This module's own errors section shows the failure mode directly: a network with far more parameters than training examples will memorise the training set and score much worse on unseen data — the loss keeps dropping on training data while test performance gets worse. Model size has to be matched to dataset size and problem complexity; the fix for "not learning enough" is not always "add more layers," and the fix for overfitting is usually to shrink the network, add regularisation, or get more data — not the reverse.

There is no closed-form solution being solved here. Gradient descent takes small, iterative steps downhill on a loss surface that is almost never convex for a real network — different random initialisations, different mini-batch orderings, and different learning rates will all converge to different sets of weights, often with similar but never identical loss. That is precisely why this module trains with mini-batches over many epochs rather than solving a system of equations once: there is no single correct answer to converge to, only a "good enough" region of weight space reached through repeated, approximate steps.

Starting from input X of shape (batch, n_features): layer 1 computes Z1 = X @ W1 + b1, producing a (batch, n_hidden1) pre-activation, then applies an activation function element-wise to get A1 — the actual output that becomes the input to the next layer. Layer 2 repeats this: Z2 = A1 @ W2 + b2, A2 = activation(Z2). The output layer does one more linear step, Z3 = A2 @ W3 + b3, and either leaves it linear (regression) or applies sigmoid/softmax (classification) depending on the task. Every intermediate value — X, Z1, A1, Z2, A2 — has to be cached during this pass, because backpropagation needs all of them to compute gradients afterward.

Without one, stacking layers is mathematically pointless: a linear function of a linear function is still just a linear function. If every layer only computed Z = A_prev @ W + b with no activation in between, the whole network — regardless of depth — could be collapsed algebraically into a single equivalent linear layer, no more expressive than plain linear regression. The non-linearity inserted after each linear step (ReLU, sigmoid, tanh) is what lets each additional layer actually add representational power, letting the network approximate curved decision boundaries and complex functions instead of only straight lines and hyperplanes.

Batch gradient descent computes the gradient over the entire dataset before taking one step — exact, but far too slow to use on large datasets, and the whole dataset often does not fit in memory. Stochastic gradient descent computes the gradient from a single example per step — very fast per step but extremely noisy, with the gradient direction jumping around from sample to sample. Mini-batch gradient descent, using batches of roughly 32 to 256 samples, is the practical compromise: it is stable enough to converge reliably, small enough to fit in memory and run efficiently on a GPU, and it is what every production deep learning framework defaults to.

If every weight in a layer starts at exactly zero, every neuron in that layer computes the exact same z, applies the exact same activation, and — critically — receives the exact same gradient during backpropagation. Every neuron in the layer updates identically forever, so a layer with 100 neurons behaves like a layer with 1 neuron repeated 100 times; the network never breaks this symmetry on its own. Random initialisation (this module uses He initialisation, scaling by √(2/n_inputs) for ReLU networks) breaks that symmetry so different neurons learn different things, while also keeping the initial activations in a numerically well-behaved range instead of vanishing or exploding.

The NumPy version in this module works, but every time the architecture changes — a new layer, a skip connection, a different activation — the backward() function has to be rewritten by hand, and a single wrong transpose silently produces incorrect gradients with no error thrown. Autograd builds a computational graph automatically during the forward pass and derives the correct backward pass for that exact graph every time, for arbitrary architectures, without anyone re-deriving the chain rule by hand. It also plugs into GPU execution, mixed-precision training, and distributed training — infrastructure that would take far more than "faster NumPy" to reimplement correctly from scratch.