Intro to Neural Networks

The simplest model that anyone calls a "neural network" is one neuron. It takes a vector $x \in \mathbb{R}^d$, computes $z = w \cdot x + b$ for some weight vector $w$ and bias $b$, then passes $z$ through a non-linearity $\sigma$ to produce an output $a = \sigma(z)$. With $\sigma$ a threshold and a binary label, this is Rosenblatt's 1958 perceptron. With $\sigma$ a logistic sigmoid, it is logistic regression. Either way, the decision boundary $\{x : w \cdot x + b = 0\}$ is a hyperplane.

A single neuron can solve linearly separable problems and nothing else. The canonical demonstration of this limit is the XOR function, where no hyperplane separates the four points. Minsky and Papert wrote a whole book on it in 1969, which had the side effect of putting connectionist research on hold for a decade. The fix is to compose neurons: feed the outputs of one layer of neurons as the inputs of the next. Each neuron in the next layer sees a non-linear function of the previous layer, and the composition can carve out non-linear boundaries.

The mental picture: a single neuron is a linear classifier. Two neurons in a hidden layer plus a non-linearity is enough to bend the boundary into the XOR shape. By the time you have ten hidden neurons in two layers, you can draw any boundary you want, if you can find the right weights. That last clause is the entire rest of this chapter.

Library path (sklearn, just to make the limit concrete):

from sklearn.linear_model import LogisticRegression
import numpy as np

X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])
y = np.array([0, 1, 1, 0])   # XOR
clf = LogisticRegression().fit(X, y)
print(clf.score(X, y))   # 0.5 — a coin flip. The boundary is linear.

From-scratch path (numpy, one neuron with sigmoid):

import numpy as np

def sigmoid(z: np.ndarray) -> np.ndarray:
    return 1.0 / (1.0 + np.exp(-z))

def predict(x: np.ndarray, w: np.ndarray, b: float) -> np.ndarray:
    return sigmoid(x @ w + b)

# Fit by hand on XOR — try as many w, b as you like, the best you can do is 0.5 accuracy.
w = np.array([1.0, 1.0])
b = -0.5
preds = (predict(np.array([[0, 0], [0, 1], [1, 0], [1, 1]]), w, b) > 0.5).astype(int)
print(preds)   # [0 1 1 1] — solves OR, not XOR.

A multilayer perceptron is what you get by stacking $L$ affine transformations with element-wise non-linearities between them. For one hidden layer with $h$ hidden units, mapping $\mathbb{R}^d \to \mathbb{R}^q$:

Here $X \in \mathbb{R}^{n \times d}$ is a minibatch of $n$ examples, $W^{(1)} \in \mathbb{R}^{d \times h}$ and $W^{(2)} \in \mathbb{R}^{h \times q}$ are weight matrices, and $b^{(1)}, b^{(2)}$ are biases. $\sigma$ is the activation function applied element-wise. $H$ is the hidden representation.

The non-linearity is load-bearing. Drop it, and $O = X W^{(1)} W^{(2)} + b^{(1)} W^{(2)} + b^{(2)} = X W' + b'$ for $W' = W^{(1)} W^{(2)}$, which is a single affine layer. Two boring layers compose to one boring layer. Add $\sigma$ and the composition can express functions no single affine layer can.

The universal approximation theorem (Cybenko 1989, Hornik 1991) says that a one-hidden-layer MLP with enough hidden units can approximate any continuous function on a compact set, to any desired accuracy. This is reassuring and operationally useless: it does not tell you how many units, what initialization, what optimizer, or whether you can find those weights with gradient descent. In practice, deep networks (many narrow layers) outperform wide networks (one fat hidden layer) at the same parameter budget, for reasons that are still partly mysterious and partly explained by inductive bias.

The PyTorch pieces you will see below — nn.Module (the base class every model subclasses), nn.Linear (one affine layer $xW^T + b$), nn.Sequential (a chain of layers run in order), and later loss.backward / opt.step / DataLoader — are used here on faith; another chapter builds the whole API up from the micrograd worldview you implement in this chapter.

model = nn.Sequential(
    nn.Linear(2, 4), nn.ReLU(),
    nn.Linear(4, 4), nn.ReLU(),
    nn.Linear(4, 1),   # last layer linear, no activation
)

class Neuron(Module):
    def __init__(self, n_in, nonlin=True):
        self.w = [Value(random.uniform(-1, 1)) for _ in range(n_in)]
        self.b = Value(0.0)
        self.nonlin = nonlin
    def __call__(self, x):
        act = sum((wi * xi for wi, xi in zip(self.w, x)), start=self.b)
        return act.relu() if self.nonlin else act
    def parameters(self):
        return self.w + [self.b]

class Layer(Module):
    def __init__(self, n_in, n_out, **kwargs):
        self.neurons = [Neuron(n_in, **kwargs) for _ in range(n_out)]
    def __call__(self, x):
        out = [n(x) for n in self.neurons]
        return out[0] if len(out) == 1 else out

class MLP(Module):
    def __init__(self, n_in, n_outs):
        sizes = [n_in] + n_outs
        self.layers = [Layer(sizes[i], sizes[i+1],
                       nonlin=(i < len(n_outs) - 1))
                       for i in range(len(n_outs))]
    def __call__(self, x):
        for layer in self.layers:
            x = layer(x)
        return x

The two implementations compute the same function. The library version hides the parameter tensors behind nn.Linear; the from-scratch version puts them right there as self.W1, self.b1. When you read PyTorch code, nn.Linear(d_in, d_out) stores weight of shape (d_out, d_in) and applies x @ weight.T + bias. The .T is why a from-scratch version that stores W1 as (d_in, d_hidden) differs from nn.Linear's transposed layout. This is the kind of thing nobody tells you and that costs you an hour the first time you try to copy weights between the two.

The non-linearity in the middle of every MLP layer is called the activation function. Five of them matter historically.

Sigmoid ($\sigma(x) = 1/(1 + e^{-x})$) was the original. It is bounded between 0 and 1, smooth, and easy to interpret as a probability. Its gradient is $\sigma(x)(1 - \sigma(x))$, which is at most $0.25$ at $x = 0$ and approaches $0$ for $|x| > 4$. Stack a few sigmoid layers and the gradient at the input is the product of these small numbers. This is the vanishing gradient problem; it is why pre-2010 deep networks did not train.

tanh ($\tanh(x)$) is sigmoid shifted to be zero-centered, in $[-1, 1]$. Better than sigmoid because activations are zero-mean, but still saturates and still vanishes.

ReLU ($\max(0, x)$) was the breakthrough (Nair and Hinton 2010, AlexNet 2012). It is not bounded. Its gradient is exactly 1 for $x > 0$ and exactly 0 for $x < 0$. No saturation on the positive side, no vanishing for active neurons. Cheap to compute. The cost is the dead ReLU problem: a neuron that always outputs 0 has zero gradient and never recovers. Caveat: people overstate this; with decent initialization, dead ReLUs are rare and self-healing during training.

GELU ($x \cdot \Phi(x)$, where $\Phi(x)$ is the standard Gaussian cumulative distribution function — the probability that a standard-normal draw lands below $x$, rising smoothly from 0 to 1) is a smooth approximation to ReLU that does not have a hard zero. It is the default in BERT, GPT-2, and most pre-2023 transformers. Slightly more expensive than ReLU. In practice, marginally better.

SwiGLU ($\text{SiLU}(xW + b) \odot (xV + c)$, where SiLU — also called Swish — is $x \cdot \text{sigmoid}(x)$) is the gated variant used in Llama, PaLM, and most 2023+ transformers. It is two parallel projections multiplied element-wise, with one of them gated by a sigmoid-like function. The math says it has more expressive capacity per parameter than GELU; the empirics agree by a small but consistent margin.

The lesson is not that one is best. It is that you should know what the current default is for the architecture you are building (ReLU for CNNs, GELU for early transformers, SwiGLU for modern LLMs) and not deviate without a reason. The architecture names here — transformers, and the specific models BERT, GPT-2, Llama, PaLM — are the dominant sequence models built on the attention mechanism; you do not need them yet, and another chapter builds one from scratch.

Library path (each one is a one-liner in PyTorch):

import torch.nn.functional as F

x = torch.linspace(-5, 5, 100)
y_sigmoid = torch.sigmoid(x)
y_tanh = torch.tanh(x)
y_relu = F.relu(x)
y_gelu = F.gelu(x)
y_silu = F.silu(x)   # the building block of SwiGLU

From-scratch path (NumPy, with derivatives — you will need these in section 5):

import numpy as np

def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))

def d_sigmoid(x):
    s = sigmoid(x)
    return s * (1 - s)

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

def d_relu(x):
    return (x > 0).astype(x.dtype)

def gelu(x):
    # Hendrycks-Gimpel exact form via erf
    from scipy.special import erf
    return 0.5 * x * (1 + erf(x / np.sqrt(2)))

# d_gelu is messier; in practice you let autograd handle it.

Computing the output of an MLP on a minibatch is one matmul, one element-wise non-linearity, one more matmul, per layer pair. The shapes are the part you have to track.

For a 2-layer MLP, $X \in \mathbb{R}^{n \times d}$, hidden dim $h$, output dim $q$:

The bias $b^{(1)}$ has shape $(h,)$ and broadcasts across the $n$ batch dimension. The softmax in step 4 is across the $q$ class axis: each row of $\hat{Y}$ is a probability distribution over classes. The pre-softmax outputs $Z^{(2)}$ — the raw class scores before normalization, which another chapter called decision-function scores — have a name from here on: logits. They are what cross_entropy expects, because it folds the softmax and the log into one numerically stable step internally.

This is the easy direction. Implementations rarely get the forward pass wrong. The hard direction is the backward pass, and to do that, you need to remember what $Z^{(1)}$, $H$, and $Z^{(2)}$ were. The forward pass caches its intermediate activations precisely so the backward pass can use them.

Library path:

import torch
import torch.nn as nn
import torch.nn.functional as F

class TwoLayerMLP(nn.Module):
    def __init__(self, d, h, q):
        super().__init__()
        self.fc1 = nn.Linear(d, h)
        self.fc2 = nn.Linear(h, q)
    def forward(self, x):
        return self.fc2(F.relu(self.fc1(x)))   # PyTorch caches activations automatically

x = torch.randn(32, 784)
model = TwoLayerMLP(784, 128, 10)
logits = model(x)   # shape (32, 10) — no softmax; cross_entropy expects logits

From-scratch path (with explicit cache for the backward pass):

def forward(X, params):
    W1, b1, W2, b2 = params['W1'], params['b1'], params['W2'], params['b2']
    Z1 = X @ W1 + b1
    H = relu(Z1)
    Z2 = H @ W2 + b2
    cache = {'X': X, 'Z1': Z1, 'H': H, 'Z2': Z2}
    return Z2, cache   # return logits + cache for backward

A neural network is a function $L(\theta) = \ell(f(x; \theta), y)$ where $f$ is the model, $\ell$ is the loss, and $\theta$ are the parameters. Training means minimizing $L$ by gradient descent: $\theta \leftarrow \theta - \eta \nabla_\theta L$. The hard part is computing $\nabla_\theta L$ when $L$ is a 10-layer chain of matrix multiplies and non-linearities.

Backpropagation is the answer. It is the chain rule applied carefully. The careful part is that you have to remember what activations you saw on the way in, and you have to walk the chain in reverse, so each node accumulates gradients from all of its downstream uses.

The most natural way to think about it: a computational graph is a directed acyclic graph where nodes are values (scalars, tensors) and edges are operations (+, *, matmul, relu). The forward pass builds the graph as you compute. The backward pass walks the graph in topological reverse order, calling a _backward function at each node that uses the chain rule to propagate gradients to its inputs.

For a single multiplication node $c = a \cdot b$ with downstream gradient $\partial L / \partial c$:

• $\partial L / \partial a = b \cdot \partial L / \partial c$
• $\partial L / \partial b = a \cdot \partial L / \partial c$

For an addition $c = a + b$:

• $\partial L / \partial a = \partial L / \partial c$
• $\partial L / \partial b = \partial L / \partial c$

For ReLU $c = \max(0, a)$:

• $\partial L / \partial a = \mathbb{1}[a > 0] \cdot \partial L / \partial c$

These are the only local rules you need. Backprop composes them across an entire model.

Karpathy's micrograd is 100 lines of Python that implements reverse-mode autodiff on scalars. Its Value class wraps a single number and tracks the operations that produced it. Every operation produces a new Value and attaches a _backward closure that knows how to propagate gradients through that operation. Calling .backward on the final loss does a topological sort of the graph and runs the closures in reverse order.

The key trick is that closures capture the operation's inputs. When you compute out = a * b, the closure for out._backward keeps references to a and b, so when later asked to propagate out.grad to its inputs, it can do a.grad += b.data * out.grad and b.grad += a.data * out.grad.

Here is the canonical implementation, lightly annotated. Read it once, then again, then close the file and re-implement it.

class Value:
    """Stores a single scalar value and its gradient."""
    def __init__(self, data, _children=(), _op=''):
        self.data = data
        self.grad = 0.0
        self._backward = lambda: None
        self._prev = set(_children)
        self._op = _op

    def __add__(self, other):
        other = other if isinstance(other, Value) else Value(other)
        out = Value(self.data + other.data, (self, other), '+')
        def _backward():
            self.grad += out.grad      # d(a+b)/da = 1
            other.grad += out.grad     # d(a+b)/db = 1
        out._backward = _backward
        return out

    def __mul__(self, other):
        other = other if isinstance(other, Value) else Value(other)
        out = Value(self.data * other.data, (self, other), '*')
        def _backward():
            self.grad += other.data * out.grad   # d(a*b)/da = b
            other.grad += self.data * out.grad   # d(a*b)/db = a
        out._backward = _backward
        return out

    def __pow__(self, other):
        assert isinstance(other, (int, float))
        out = Value(self.data ** other, (self,), f'**{other}')
        def _backward():
            self.grad += (other * self.data ** (other - 1)) * out.grad
        out._backward = _backward
        return out

    def relu(self):
        out = Value(0 if self.data < 0 else self.data, (self,), 'ReLU')
        def _backward():
            self.grad += (out.data > 0) * out.grad
        out._backward = _backward
        return out

    def backward(self):
        # Topological sort: children before parents
        topo, visited = [], set()
        def build_topo(v):
            if v not in visited:
                visited.add(v)
                for child in v._prev:
                    build_topo(child)
                topo.append(v)
        build_topo(self)
        # Seed the gradient of the output as 1 and walk back
        self.grad = 1.0
        for v in reversed(topo):
            v._backward()

That is it. No exotic features, no symbolic differentiation, no compiler. Eighty lines. With this, plus the Neuron / Layer / MLP classes from micrograd.nn (another 50 lines), you can train a 2-hidden-layer MLP on the moons dataset and watch it learn the decision boundary.

Micrograd does autograd one scalar at a time. PyTorch does it on tensors. The shift is mechanical but the bookkeeping is tighter.

For matmul $C = A B$ where $A \in \mathbb{R}^{m \times k}$ and $B \in \mathbb{R}^{k \times n}$:

• $\partial L / \partial A = (\partial L / \partial C) B^T$, shape $(m, k)$
• $\partial L / \partial B = A^T (\partial L / \partial C)$, shape $(k, n)$

For broadcasting (a bias $b \in \mathbb{R}^h$ added to $Z \in \mathbb{R}^{n \times h}$):

• $\partial L / \partial b = \sum_{i} \partial L / \partial Z_{i,:}$, summed over the broadcast axis

For softmax-cross-entropy, the gradient with respect to the pre-softmax logits is famously clean. If $p = \text{softmax}(z)$ and $L = -\sum_i y_i \log p_i$ for one-hot $y$:

• $\partial L / \partial z = p - y$

That last identity is one of the most elegant results in basic deep learning. It is the reason cross-entropy on logits is the universal loss for classification: the gradient is just "predicted minus actual", with no division by tiny numbers or log-of-zero edge cases.

logits = model(x)                  # (n, K), raw scores, NO softmax
loss = F.cross_entropy(logits, y)  # y is int class indices (n,), reduction='mean'
loss.backward()                    # d_loss/d_logits == (softmax(logits) - onehot(y)) / n

def softmax_xent(logits, y):
    z = logits - logits.max(axis=1, keepdims=True)   # stability shift
    exp_z = np.exp(z)
    probs = exp_z / exp_z.sum(axis=1, keepdims=True)
    n = logits.shape[0]
    loss = -np.log(probs[np.arange(n), y] + 1e-12).mean()
    grad = probs.copy()
    grad[np.arange(n), y] -= 1     # probs - one_hot(y)
    grad /= n                      # average over the batch
    return loss, grad

That is the entire algorithm. The from-scratch loss returns both the scalar and the (probs - one_hot(y)) / n gradient; PyTorch's F.cross_entropy(logits, y) followed by loss.backward is doing a generalized version of the same computation for whatever graph the forward pass produced — the full backward through Z2 = H @ W2 + b2, H = relu(Z1), Z1 = X @ W1 + b1 is just the matmul and ReLU gradient rules above composed in reverse.

For about two decades after backprop was popularized (Rumelhart, Hinton, Williams 1986), training deep networks did not work. Networks deeper than two or three hidden layers either failed to learn, learned very slowly, or got stuck in regions where the loss looked flat. Three things were going wrong simultaneously, and they only got sorted out in the late 2000s.

Vanishing gradients. With sigmoid or tanh activations, the gradient at any layer is the product of the local derivatives at all downstream layers. Sigmoid's derivative is bounded by $0.25$. Ten layers in, the gradient at the input is at most $0.25^{10} \approx 10^{-6}$. Learning never propagates back that far. The fix was switching to ReLU, whose derivative is exactly 1 on the active side.

Bad initialization. If weights are drawn from a normal with too-large variance, the pre-activations $z = Wx + b$ blow up, the sigmoid saturates, the gradient vanishes. Too-small variance and the activations collapse to zero. Glorot and Bengio (2010) derived the right scale: $\text{Var}(W) = 2 / (\text{fan\_in} + \text{fan\_out})$, now called Xavier or Glorot initialization. He et al. (2015) refined it for ReLU: $\text{Var}(W) = 2 / \text{fan\_in}$, called Kaiming or He initialization. Modern frameworks do this by default; in 2008 nobody did.

Optimizers and learning rates. Plain SGD with a fixed learning rate gets stuck. Momentum helps. Adaptive methods (RMSprop, Adam) help more. Learning rate warmup was a 2017 invention. In the late 2000s, every paper had to hand-tune all of this, and most papers reported numbers that did not reproduce.

What unblocked deep networks was the combination: ReLU + good init + Adam + lots of data + GPUs. AlexNet (2012) had all five. The chapter on training (another chapter) walks through each in detail; for now, the takeaway is that "deep networks didn't train" was a true statement in 2005 and a false statement in 2015, and the difference is not philosophical, it is mechanical.

MNIST is 70,000 grayscale images of handwritten digits, 28×28 pixels each, labeled 0 through 9. 60,000 are training, 10,000 are test. It has been a benchmark since 1998. Modern models reaanother chapter.7% test accuracy. A 2-layer MLP with 128 hidden units, trained for one epoch — one full pass through the entire training set, minibatch by minibatch — on a free Colab GPU, reaches about 97.5%. This is the "hello world" that every framework demos with.

An MLP wants a flat feature vector, not a grid, so eaanother chapter×28 image is flattened row by row into a single vector of $28 \times 28 = 784$ pixel values (that is the nn.Flatten in the model and the 784 the network's first layer expects). The reason MNIST is the right first model: the input is small enough (784 features) that you can fit a working MLP in memory on any laptop, the classes are clean (digits, not "is this a cat or a dog with a fence in front"), and the dataset is balanced. You can debug your code without confounding it with messy data. Once you can train MNIST, the same code generalizes to Fashion-MNIST, CIFAR-10, and (with convolutions in another chapter) ImageNet.

Library path (the canonical 30 lines, ready to copy into Colab):

import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.utils.data import DataLoader
from torchvision import datasets, transforms

device = "cuda" if torch.cuda.is_available() else "cpu"

tfm = transforms.Compose([transforms.ToTensor(), transforms.Normalize((0.1307,), (0.3081,))])
train_ds = datasets.MNIST("./data", train=True, download=True, transform=tfm)
test_ds = datasets.MNIST("./data", train=False, download=True, transform=tfm)
train_dl = DataLoader(train_ds, batch_size=128, shuffle=True)
test_dl = DataLoader(test_ds, batch_size=512)

class MLP(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(
            nn.Flatten(),
            nn.Linear(784, 128), nn.ReLU(),
            nn.Linear(128, 64), nn.ReLU(),
            nn.Linear(64, 10),
        )
    def forward(self, x): return self.net(x)

model = MLP().to(device)
opt = torch.optim.Adam(model.parameters(), lr=3e-4)

for epoch in range(3):
    model.train()
    for x, y in train_dl:
        x, y = x.to(device), y.to(device)
        loss = F.cross_entropy(model(x), y)
        opt.zero_grad(); loss.backward(); opt.step()

    model.eval()
    correct = 0
    with torch.no_grad():
        for x, y in test_dl:
            x, y = x.to(device), y.to(device)
            correct += (model(x).argmax(1) == y).sum().item()
    print(f"epoch {epoch}: test acc = {correct / len(test_ds):.4f}")

After 3 epochs you should see test accuracy around 0.977. If you see anything below 0.95, something is wrong; check the normalization stats, the loss function (it must be cross_entropy on logits, not on softmaxed probabilities), and that you are calling opt.zero_grad before each backward pass.

The point of doing it from scratch is not to write production code. The point is to convince yourself that the library is not doing anything you cannot do. Here is a complete two-layer MLP trained on MNIST in pure NumPy, with forward, backward, SGD, and accuracy evaluation. It hits roughly 95% test accuracy in a minute on a laptop CPU.

import numpy as np

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

def softmax_xent(logits, y):
    z = logits - logits.max(axis=1, keepdims=True)
    exp_z = np.exp(z)
    probs = exp_z / exp_z.sum(axis=1, keepdims=True)
    n = logits.shape[0]
    loss = -np.log(probs[np.arange(n), y] + 1e-12).mean()
    grad = probs.copy()
    grad[np.arange(n), y] -= 1
    grad /= n
    return loss, grad

class MLP:
    def __init__(self, d_in=784, d_hid=128, d_out=10, seed=0):
        rng = np.random.default_rng(seed)
        # Kaiming init for ReLU
        self.W1 = rng.standard_normal((d_in, d_hid)) * np.sqrt(2 / d_in)
        self.b1 = np.zeros(d_hid)
        self.W2 = rng.standard_normal((d_hid, d_out)) * np.sqrt(2 / d_hid)
        self.b2 = np.zeros(d_out)

    def forward(self, X):
        Z1 = X @ self.W1 + self.b1
        H = relu(Z1)
        Z2 = H @ self.W2 + self.b2
        self._cache = (X, Z1, H)
        return Z2

    def backward(self, grad_logits):
        X, Z1, H = self._cache
        dW2 = H.T @ grad_logits
        db2 = grad_logits.sum(0)
        dH = grad_logits @ self.W2.T
        dZ1 = dH * (Z1 > 0)
        dW1 = X.T @ dZ1
        db1 = dZ1.sum(0)
        return {'W1': dW1, 'b1': db1, 'W2': dW2, 'b2': db2}

    def step(self, grads, lr=0.1):
        self.W1 -= lr * grads['W1']; self.b1 -= lr * grads['b1']
        self.W2 -= lr * grads['W2']; self.b2 -= lr * grads['b2']

def train(model, X_train, y_train, X_test, y_test, epochs=5, bs=128, lr=0.1):
    n = X_train.shape[0]
    for epoch in range(epochs):
        perm = np.random.permutation(n)
        for i in range(0, n, bs):
            idx = perm[i:i+bs]
            logits = model.forward(X_train[idx])
            loss, grad = softmax_xent(logits, y_train[idx])
            grads = model.backward(grad)
            model.step(grads, lr=lr)
        acc = (model.forward(X_test).argmax(1) == y_test).mean()
        print(f"epoch {epoch}: test acc = {acc:.4f}")

That is 52 lines including blank lines. It trains. You can read every operation. No nn.Module, no optimizer.step, no .backward mystery. The only thing it depends on is NumPy.

A neural network that fails to train can fail in many ways. Three account for ≈90% of beginner debugging time.

Dead ReLUs. A neuron whose pre-activation $z = Wx + b$ is always negative outputs zero, has gradient zero, and never updates. With He initialization and reasonable learning rates, this is rare. With a learning rate too high (so the first big step pushes weights into the negative region) or with poor init, half your neurons can die in the first epoch. The diagnostic: after a few hundred steps, log the fraction of hidden units that are non-zero on a fixed batch. If it falls below ~50%, you have dead ReLUs. The fix: lower the learning rate, use Kaiming init, or use Leaky ReLU / GELU.

Vanishing or exploding gradients. Log the gradient norm at each layer for the first 100 steps. If the gradient at the input layer is 1000× smaller than at the output layer, gradients are vanishing. If the loss spikes to NaN, gradients exploded. The fixes: better init, gradient clipping, switching activations, or using LayerNorm (another chapter).

Loss flat at $\log K$, then explodes. A model with $K$ classes that has not learned anything outputs uniform predictions, giving loss $\log K \approx 2.30$ for MNIST. If your loss sits there for hundreds of steps, your gradients are not reaching the early layers, or your learning rate is too small, or your data is wrong (mislabeled, all-zero, normalized incorrectly). Karpathy's recipe: verify loss at init. The very first loss before any training step should be $\log K$ for a balanced classifier with sensible initialization. If it's anything else, debug before you train.

These three diagnostics catch the vast majority of "my model is not training" situations. Build them into every training script you write.