Neural Networks Introduction, in depth: neurons, layers as matrix multiplies, why nonlinearity matters, backprop by hand and a network from scratch

A neuron is a weighted sum and a bend

A single artificial neuron takes inputs x1 to xn, multiplies each by a weight, adds a bias and passes the result through an activation function: output = f(w1 x1 + ... + wn xn + b). The weights say how much each input matters and in which direction, the bias shifts the threshold, and f bends the result.

Common choices for f: ReLU, max(0, z), which is cheap and the default for hidden layers; tanh, which squashes to the range -1 to 1 and is smooth; and sigmoid, 1 / (1 + e^-z), which squashes to 0 to 1 and is mostly used at the output to turn a score into a probability. Despite the biological name, it is just a linear function followed by a fixed nonlinear one.

Layers are matrix multiplications

A layer is many neurons reading the same inputs. Put each neuron's weights in a column and the whole layer becomes one matrix multiply: for a batch of inputs X with shape (batch, inputs), weights W with shape (inputs, outputs) and bias b with shape (outputs), the layer computes f(X W + b), shape (batch, outputs). Getting shapes right is most of the practical work, so it is worth doing once by hand.

Take a classic digit classifier: 28 by 28 pixel images flattened to 784 inputs, one hidden layer of 128 units, 10 outputs. The hidden layer has 784 x 128 = 100,352 weights plus 128 biases, 100,480 parameters. The output layer has 128 x 10 + 10 = 1,290. The total is 101,770 parameters. For a batch of 32 images the shapes flow (32, 784) to (32, 128) to (32, 10). Each output row is ten scores, called logits, one per digit.

Depth means stacking layers; width means more units per layer. Neither helps without the bend.

Why the bend matters: XOR

Without activations, two layers collapse into one: (X W1) W2 = X (W1 W2), and any stack of linear layers is a single linear map. A linear model can only separate classes with a straight line, or a flat plane in more dimensions. XOR, where the output is 1 when exactly one of two bits is 1, cannot be separated that way. Fit the best linear function to the four XOR points by least squares and you get weights of zero on both inputs and a bias of 0.5: it predicts 0.5 everywhere, which is the honest answer of a model that cannot see the structure.

Now add a hidden layer with ReLU. Set h1 = ReLU(x1 + x2) and h2 = ReLU(x1 + x2 - 1), and output y = h1 - 2 h2. For inputs (0,0), (0,1), (1,0), (1,1) the hidden values are (0,0), (1,0), (1,0) and (2,1), and the outputs are 0, 1, 1, 0. That is exactly XOR. The second unit only switches on when both bits are set, and the output uses it to cancel the first. The bend let the network carve the input space into regions and treat them differently, which is the entire trick, repeated thousands of times in real models.

Training finds such weights for us. The rest of the article is about how.

The loss: telling the network what good means

To learn, the network needs a single number saying how wrong it is. For regression, mean squared error between prediction and target is the usual choice. For binary classification, use binary cross-entropy: with predicted probability p and label y, the loss is -(y log p + (1 - y) log(1 - p)), averaged over the batch. It punishes confident wrong answers heavily, which is what you want. For multi-class problems, apply softmax to the logits and use cross-entropy against the correct class.

A useful sanity number: a classifier that knows nothing outputs p = 0.5 for a binary task, giving a loss of ln 2, about 0.693. If your untrained binary model starts far from 0.693, check the output layer and the loss wiring before anything else. The run below starts at 0.708, right where it should.

Backpropagation for a two-layer network

Training means adjusting every weight in the direction that lowers the loss, which requires the gradient: how the loss changes with each weight. Backpropagation is the chain rule applied layer by layer from the output back, reusing the activations cached during the forward pass so the whole gradient costs roughly two forward passes.

For our network, sigmoid output plus binary cross-entropy has a famously clean derivative: dL/dz2 = (p - y) / N, the prediction error divided by the batch size. From there each step is one line. The output weights' gradient is the hidden activations transposed times that error: dW2 = a1^T dz2. The error flows back through W2 to the hidden layer, dz2 W2^T, and is multiplied by the derivative of tanh, which is 1 - a1^2. That gives dz1, and dW1 = X^T dz1. Each bias gradient is its layer's error summed over the batch.

Notice the pattern: every weight gradient is input-to-that-layer transposed times error-at-that-layer. That is why frameworks can automate it. They record the forward operations and apply each one's local derivative in reverse. You should still do it by hand once, because when training breaks, the explanation is almost always in this chain: a derivative that vanishes, a shape that broadcasts silently, or activations that saturate.

The whole thing in NumPy

import numpy as np

X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]], float)
y = np.array([[0], [1], [1], [0]], float)

def train(hidden, seed, steps=2000, lr=0.5):
    rng = np.random.default_rng(seed)
    W1 = rng.normal(0, 1 / np.sqrt(2), (2, hidden)); b1 = np.zeros(hidden)
    W2 = rng.normal(0, 1 / np.sqrt(hidden), (hidden, 1)); b2 = np.zeros(1)
    for step in range(steps + 1):
        # forward
        z1 = X @ W1 + b1
        a1 = np.tanh(z1)
        z2 = a1 @ W2 + b2
        p = 1 / (1 + np.exp(-z2))
        loss = -np.mean(y * np.log(p + 1e-12) + (1 - y) * np.log(1 - p + 1e-12))
        # backward: chain rule, layer by layer, reusing cached activations
        dz2 = (p - y) / len(X)
        dW2 = a1.T @ dz2;  db2 = dz2.sum(0)
        dz1 = (dz2 @ W2.T) * (1 - a1 ** 2)      # tanh'(z) = 1 - tanh(z)^2
        dW1 = X.T @ dz1;   db1 = dz1.sum(0)
        # update
        W1 -= lr * dW1; b1 -= lr * db1; W2 -= lr * dW2; b2 -= lr * db2
        if step in (0, 100, 500, 1000, 2000):
            print(step, round(float(loss), 4))
    return p.ravel()

print(train(hidden=4, seed=0).round(3))

Run with four hidden units and seed 0, the loss goes 0.7081 at step 0, 0.5539 at step 100, 0.0148 at step 500, 0.0055 at step 1000 and 0.0022 at step 2000. Final probabilities are 0.000, 0.997, 0.997 and 0.003: XOR, learned.

Then a more instructive experiment. Over 20 random seeds, with two hidden units, the minimum the hand-made solution needed, 12 of 20 runs classified all four points correctly after 2,000 steps; the other 8 did not. With four or eight hidden units, it solved all 20. Two units are enough to represent the answer but give the optimiser few routes to it; extra width makes the loss surface friendlier. That is a small version of a big practical lesson: networks are usually made larger than strictly necessary because it makes them easier to train, and regularisation and data, not tiny size, are the tools against overfitting.

Initialisation matters too. The code scales random weights by one over the square root of the input count, so the variance of each layer's output stays roughly constant. Start all weights at zero and every hidden unit computes the same thing and receives the same gradient forever. The reasons and the variants are covered in Weight Initialization.

The same network in PyTorch

import torch, torch.nn as nn

X = torch.tensor([[0., 0.], [0., 1.], [1., 0.], [1., 1.]])
y = torch.tensor([[0.], [1.], [1.], [0.]])

model = nn.Sequential(nn.Linear(2, 4), nn.Tanh(), nn.Linear(4, 1))   # logits out
loss_fn = nn.BCEWithLogitsLoss()        # sigmoid + BCE, numerically stable
opt = torch.optim.SGD(model.parameters(), lr=0.5)

for step in range(2000):
    opt.zero_grad()
    loss = loss_fn(model(X), y)
    loss.backward()                     # autograd does what the NumPy backward did
    opt.step()

print(torch.sigmoid(model(X)).detach().squeeze())

Three things changed. Layers are objects that own their weights and initialise them sensibly. The output layer emits logits, and BCEWithLogitsLoss applies the sigmoid inside the loss, which avoids taking the log of a probability that has rounded to exactly zero. And loss.backward() replaces the hand-written backward pass using autograd, which is the same chain rule computed from a recorded graph. The optimiser call opt.step() applies the update. For real work you would add a data loader, a held-out split as described in Train, Validation and Test Splits and usually the Adam optimiser instead of plain SGD.

What the hardware does with it

Almost all the work in a neural network is matrix multiplication: the X W in every layer, forward and backward. A matrix multiply of shapes (m, k) and (k, n) costs about 2 m k n floating-point operations, and the backward pass needs two more multiplies of similar size per layer, one for the weight gradient and one to pass the error back. That is why training costs roughly three times the forward pass, and why GPUs, which are built to run thousands of multiply-adds in parallel with dedicated matrix units, dominate this field. Activations like tanh and ReLU are cheap elementwise operations by comparison; their cost is mostly memory traffic.

Memory is the other constraint. Each parameter needs its value and its gradient; the Adam optimiser adds two more running averages. In 32-bit floats that is 16 bytes per parameter before activations, so our 101,770-parameter digit model needs about 1.6 MB, and a one-billion-parameter model needs 16 GB for weights, gradients and optimiser state alone. Activations, the cached a1 and z1 values backprop needs, grow with batch size and depth and often exceed that. The transformer-scale version of this arithmetic is worked through in Backprop: chain rule, memory and kernels.

Universal approximation, stated honestly

You will read that neural networks can approximate any function. The universal approximation theorems do say that a network with one hidden layer and a suitable nonlinearity can approximate any continuous function on a bounded region as closely as you like, given enough hidden units. They say nothing about how many units that takes, which can be astronomically many, nothing about whether gradient descent will find the weights, and nothing about generalising from finite data to inputs you have not seen. They tell you the representation is not the bottleneck. Training, data and architecture are, which is where the practical effort goes.

Failure modes and how to debug them

• Loss does not move. Learning rate too small, gradients not reaching the weights (forgot zero_grad or step, or detached a tensor), or dead ReLUs from a bad initialisation.
• Loss becomes NaN. Learning rate too large, log of zero in a hand-written loss, or unnormalised inputs in the thousands. Normalise inputs to roughly zero mean and unit variance.
• Stuck at chance. Labels shuffled relative to inputs, or the wrong loss for the output, such as softmax output fed to a loss that applies softmax again.
• Training loss falls, validation loss rises. Overfitting. Get more data, add augmentation or weight decay, or stop earlier.
• Silent shape bugs. A (batch,) target broadcast against a (batch, 1) prediction makes a (batch, batch) loss that still trains, badly. Assert shapes.

The first thing to try when anything is wrong is to see whether the model can memorise one batch:

# The single most useful debugging step: can the model memorise ONE small batch?
xb, yb = next(iter(train_loader))
for step in range(500):
    opt.zero_grad()
    loss = loss_fn(model(xb), yb)
    loss.backward()
    opt.step()
print(loss.item())   # should approach zero; if it does not, the bug is in code, not data

If a model cannot drive the loss on 32 examples to nearly zero, it has a bug, and no amount of extra data or tuning will fix it.

What to do next

1. Type in the NumPy XOR network and run it; then change the hidden size to 2 and count how many of 20 seeds succeed.
2. Replace tanh with ReLU in the NumPy code, update the derivative in the backward pass, and confirm the loss still falls.
3. Check one gradient numerically: nudge a single weight by 1e-5, recompute the loss and compare the slope to your dW entry.
4. Port the network to PyTorch and train the 784-128-10 digit classifier; verify the 101,770 parameter count with a one-line sum.
5. Make the overfit-one-batch check the first step of every new model you build.
6. Read about the update rule in depth, then about initialization, before moving on to convolutional networks or transformers.