Neural Networks from Scratch

The neuron

An artificial neuron takes inputs x1…xd, multiplies each by a weight, adds a bias, and passes the total through a nonlinear activation function g:

a = g(w1x1 + … + wdxd + b) = g(w · x + b)

With g = sigmoid this is exactly logistic regression from M55. The name comes from a loose analogy with biological neurons that "fire" when their inputs are strong enough; the analogy stops there. The power comes from layers: the outputs of many neurons become the inputs of the next layer. A network with an input layer, one or more hidden layers, and an output layer is a multi-layer perceptron (MLP), also called a fully connected or dense network.

Why bother? A single neuron's decision boundary is a straight line. The classic counterexample is XOR: output 1 when exactly one of two inputs is 1. No line separates {(0,1), (1,0)} from {(0,0), (1,1)}. A hidden layer can bend space so that a line works afterwards. Explore a 2-4-1 network below: pick an input, watch values flow forward, and switch to training mode to see it learn XOR.

Remember: neuron = weighted sum + bias + nonlinearity. Layers of neurons can represent boundaries no single line can.

Activation functions

Without a nonlinearity, stacking layers is pointless: W2(W1x) is just (W2W1)x, another linear map. The activation is what lets depth add power. Three you must know:

• Sigmoid σ(z) = 1/(1 + e−z): squashes to (0, 1). Derivative σ(z)(1 − σ(z)), at most 0.25. Used for binary outputs; rarely in hidden layers now.
• tanh tanh(z): squashes to (−1, 1), zero-centered. Derivative 1 − tanh²(z), at most 1. Better than sigmoid for hidden layers, still saturates.
• ReLU (rectified linear unit) max(0, z): derivative 1 for positive inputs, 0 for negative. Cheap, does not saturate for positive values, the default in modern hidden layers. A neuron whose input is always negative stops learning ("dead ReLU"); variants such as Leaky ReLU and GELU soften that.

Remember: the activation supplies the nonlinearity; sigmoid/tanh saturate (tiny derivatives for large |z|), ReLU does not for positive inputs, which is why it is the hidden-layer default.

The forward pass in matrix form

Process a whole batch of n examples at once. Put them in rows of X (shape (n, d)). A dense layer with h neurons has a weight matrix W of shape (d, h) (column j holds neuron j's weights) and a bias vector b of shape (h,):

Z = X W + b (shape (n, h)),   A = g(Z)

NumPy broadcasting adds b to every row. A two-layer network for binary classification is then

Z1 = X @ W1 + b1        # (n, d) @ (d, h) -> (n, h)
A1 = np.tanh(Z1)         # hidden activations
Z2 = A1 @ W2 + b2        # (n, h) @ (h, 1) -> (n, 1)
P  = sigmoid(Z2)         # probability of class 1

For k classes the output layer has k neurons and the softmax turns their scores (logits) into probabilities that sum to 1: softmax(z)j = ezj / Σm ezm. Subtract the row maximum before exponentiating; it does not change the result and prevents overflow. The loss is cross-entropy: minus the log of the probability given to the true class.

Remember: each layer is g(XW + b); W is (inputs, outputs); softmax + cross-entropy for multi-class outputs; print shapes.

Backpropagation and the computational graph

Training needs the gradient of the loss with respect to every weight. Backpropagation is the chain rule applied efficiently: run the forward pass, keep the intermediate values, then walk backwards from the loss, multiplying local derivatives. It is easiest to see on a computational graph, where each node is one simple operation.

Each node only needs to know its own local derivative. Multiplication u = w·x sends back x to w (and w to x); addition passes gradients through unchanged; sigmoid multiplies by σ(z)(1 − σ(z)). With binary cross-entropy, the product ∂L/∂p · σ'(z) simplifies to just p − y, which is why output-layer gradients look so clean.

For the two-layer network, backprop in matrix form (averaging over the batch of n) is:

dZ2 = (P - y) / n                 # (n, 1)   output error
dW2 = A1.T @ dZ2                  # (h, 1)   "error times input" again
db2 = dZ2.sum(axis=0)             # (1,)
dA1 = dZ2 @ W2.T                  # (n, h)   send the error back through W2
dZ1 = dA1 * (1 - A1 ** 2)         # (n, h)   times tanh's local derivative
dW1 = X.T @ dZ1                   # (d, h)
db1 = dZ1.sum(axis=0)             # (h,)

Every gradient has the same shape as the parameter it belongs to: a quick sanity check. The authoritative check is a numerical gradient: nudge one weight by a tiny ε both ways and measure the change in loss, (L(w + ε) − L(w − ε)) / 2ε. If it matches your analytic gradient to about 6 decimal places, your backprop is right.

Remember: forward computes and caches values; backward multiplies local derivatives from the loss back to each weight; gradients match parameter shapes; verify with a numerical check.

Training a two-layer network: XOR and moons

Put the pieces together: initialize, then repeat forward, loss, backward, update. First, XOR, the problem a single neuron cannot solve.

Now a slightly bigger problem with a curved boundary, plotted:

The neural-net visualizer can train XOR too; set it to train mode and compare its loss curve with the cell's.

Remember: the training loop is forward → loss → backward → update; neural losses have local minima, so seeds and learning rates matter.

Mini-batches and stochastic gradient descent

So far every step used the whole dataset (full-batch gradient descent). With a million images that is far too slow per step. Stochastic gradient descent (SGD) estimates the gradient from a small random mini-batch (typically 32–512 examples), takes a step, and moves on. One pass through all the data is an epoch. The gradient estimates are noisy, but there are many more steps per epoch, and the noise even helps escape poor regions. Shuffle the data every epoch so batches differ.

Plain SGD has popular refinements: momentum (keep a running average of gradients, like a ball rolling downhill) and Adam (per-parameter adaptive step sizes). Try momentum on the Rosenbrock valley below, then compare with plain descent.

Remember: SGD uses mini-batches for cheap, noisy, frequent steps; shuffle each epoch; momentum and Adam are the usual upgrades.

Initialization and vanishing gradients

If all weights start at zero, every hidden neuron computes the same thing and receives the same gradient, forever: the symmetry never breaks. So weights start random. But the scale matters. Each layer multiplies the signal by its weights; too large and activations explode, too small and they shrink toward zero layer after layer. The fix is to scale the random weights by the number of inputs: Xavier/Glorot initialization (variance about 1/fan_in) for tanh and sigmoid, He initialization (variance 2/fan_in) for ReLU.

Backpropagation multiplies one local derivative per layer. Sigmoid's derivative is at most 0.25, so after 10 sigmoid layers the gradient reaching the first layer can be smaller than 0.25¹⁰ ≈ 10⁻⁶: early layers barely learn. That is the vanishing gradient problem, and it held deep learning back for years. ReLU (derivative 1 when active), careful initialization, normalization layers, and residual (skip) connections are the cures.

Remember: random, fan-in-scaled initialization (Xavier for tanh/sigmoid, He for ReLU); sigmoid stacks vanish gradients; ReLU and skip connections keep deep networks trainable.

Overfitting and regularization

A network with thousands of weights can memorize a small dataset. Everything from M55 applies: watch training vs validation loss, and regularize:

• Weight decay (L2): add λ · Σw² to the loss, which adds 2λw to each weight's gradient and shrinks weights every step.
• Early stopping: track validation loss each epoch and keep the weights from the best epoch.
• Dropout: during training, randomly set each hidden activation to zero with probability p (say 0.5) and scale the survivors by 1/(1 − p). Neurons cannot rely on specific partners, so the network learns redundant, robust features; it behaves like averaging many thinned networks. At prediction time dropout is off.
• More data and data augmentation (flipped or cropped images, for example).

Remember: weight decay, early stopping, and dropout fight overfitting; dropout only runs during training, with inverted scaling.

PyTorch: the same network, automated

Writing backprop by hand teaches you what is happening. In practice you use a framework that computes gradients automatically. PyTorch is the most widely used for research and much of industry. Its core ideas map one-to-one onto what you just built:

Tensors and autograd

An nn.Module and a training loop

The GPU idea: a graphics card runs thousands of simple arithmetic operations in parallel, and neural networks are mostly big matrix multiplications, exactly that workload. Moving the model and the tensors to the same device (.to("cuda"), or "mps" on Apple silicon) can make training 10–100× faster for large models. For the tiny networks in this module, the CPU is plenty.

Remember: tensors + autograd + nn.Module + optimizer = the loop you wrote, with zero_grad → forward → loss → backward → step.

Where deep learning shines

On small tabular datasets, gradient-boosted trees (M56) usually beat neural networks with far less effort. Deep learning wins when the input has rich structure and there is a lot of data, because architectures can build that structure in:

Remember: trees for tables, deep learning for images, audio, and text; start from pretrained models when you can.

Practice

NumPy is available. The tests check values, shapes, and gradients numerically.

import numpy as np

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

def relu_grad(z):
    return (np.asarray(z) > 0).astype(float)

import numpy as np

def softmax(Z):
    Z = np.asarray(Z, dtype=float)
    E = np.exp(Z - Z.max(axis=1, keepdims=True))
    return E / E.sum(axis=1, keepdims=True)

import numpy as np

def forward(X, layers):
    A = np.asarray(X, dtype=float)
    for i, (W, b) in enumerate(layers):
        A = A @ W + b
        if i < len(layers) - 1:
            A = np.maximum(0.0, A)
    return A

import numpy as np

def numerical_gradient(f, w, eps=1e-6):
    w = np.array(w, dtype=float, copy=True)
    grad = np.zeros_like(w)
    for idx in np.ndindex(w.shape):
        old = w[idx]
        w[idx] = old + eps
        up = f(w)
        w[idx] = old - eps
        down = f(w)
        w[idx] = old
        grad[idx] = (up - down) / (2 * eps)
    return grad

import numpy as np

def backprop(X, y, W1, b1, W2, b2):
    n = len(X)
    A1 = np.tanh(X @ W1 + b1)
    P = 1 / (1 + np.exp(-(A1 @ W2 + b2)))
    loss = float(-np.mean(y * np.log(P) + (1 - y) * np.log(1 - P)))
    dZ2 = (P - y) / n
    dZ1 = dZ2 @ W2.T * (1 - A1 ** 2)
    grads = {"W1": X.T @ dZ1, "b1": dZ1.sum(axis=0), "W2": A1.T @ dZ2, "b2": dZ2.sum(axis=0)}
    return loss, grads

Mini-project

Quiz

Recap

• A neuron is g(w·x + b); layers of neurons with nonlinear activations can represent curved boundaries such as XOR's.
• Sigmoid and tanh saturate; ReLU is the hidden-layer default.
• Forward pass: A = g(XW + b) per layer; softmax + cross-entropy for classification.
• Backprop = chain rule on the computational graph; gradients have parameter shapes; verify with numerical gradients.
• Train with forward → loss → backward → update, using shuffled mini-batches (SGD), possibly with momentum or Adam.
• Initialize randomly with fan-in scaling; deep sigmoid stacks suffer vanishing gradients.
• Regularize with weight decay, early stopping, and dropout.
• PyTorch automates it: tensors, autograd, nn.Module, optimizers, devices. Install with pixi add pytorch.
• CNNs for images, transformers for sequences and language; trees still win on most small tables.

Going further

• PyTorch: Learn the Basics
• PyTorch autograd mechanics
• Pixi documentation (search for "system requirements" and "PyTorch" before setting up a GPU environment)
• scikit-learn: MLPClassifier (a ready-made MLP for small problems)
• CS231n notes: backpropagation, intuitions
• Next: M59 · ML in Production.