Neural networks from scratch: forward pass, backpropagation and activations

Frameworks like PyTorch let you train neural networks in a few lines, but when the loss stops falling or the gradients vanish, you can't tell where to start unless you understand what happens inside. This article builds multi-layer neural networks entirely in NumPy: the forward pass, backpropagation, three activation functions (sigmoid, tanh, ReLU) and four weight initialization schemes. Every claim comes with real output: gradient checking with finite differences, the XOR problem, a two-moons dataset, and a number-by-number comparison with PyTorch autograd.

This article follows PyTorch basics: tensors, autograd and the training loop, where a single backward() call does all the work; here we open that black box. The math behind the update step comes from Linear regression and gradient descent: implementing it with NumPy, and the matrix multiplication and broadcasting used throughout come from NumPy for machine learning: vectorization and broadcasting. Environment: Apple M4, 16 GB, Python 3.10, NumPy 1.26.0, scikit-learn 1.3.1, PyTorch 2.0.1 (used only to generate data). All data is synthetic.

What are neural networks: stacked linear layers with nonlinearity

A feedforward neural network (MLP) is a chain of layers, and each layer does two things: a linear transform z = a @ W + b, then a nonlinear function a' = f(z). Without the nonlinearity, many stacked linear layers are still just one linear map, so the nonlinearity is what gives neural networks their power. According to the Stanford CS231n notes on neural networks, a network with at least one hidden layer can approximate any continuous function to arbitrary accuracy. In practice this is a weak statement, because it says nothing about how many neurons are needed or how easy the network is to train; deeper networks usually work better even though their representational power is equivalent.

We use a last layer that returns a logit (before the sigmoid) and the binary cross-entropy loss. With this combination, the derivative of the loss with respect to the logit of sample i has a compact form: (sigmoid(z) - y) / n, which is exactly what we use to start backpropagation.

Forward pass and backpropagation in neural networks: the chain rule layer by layer

The CS231n notes on backpropagation describe the core idea: each operation (gate) knows its own local derivative; on the backward pass, we multiply the gradient arriving from downstream by the local derivative and pass it back. For a matrix product D = W·X, the gradients are dW = dD·Xᵀ and dX = Wᵀ·dD; a trick worth remembering is that a gradient always has the same shape as its variable, so checking dimensions helps avoid mistakes. The original paper on this method is Rumelhart, Hinton and Williams, Learning representations by back-propagating errors (Nature, 1986).

Let the delta of layer l be dL/dz for that layer. The formulas we implement come down to three lines:

The whole MLP class for neural networks is only a few dozen lines. The forward pass stores every z and a so the backward pass can reuse them, following the CS231n advice to cache intermediate values from the forward pass:

def forward(self, X):
    self.Z, self.A = [], [X]
    for i, (W, b) in enumerate(zip(self.W, self.b)):
        z = self.A[-1] @ W + b
        self.Z.append(z)
        self.A.append(z if i == len(self.W) - 1 else self.f(z))
    return self.A[-1][:, 0]

def backward(self, y):
    n = len(y)
    z = self.A[-1][:, 0]
    delta = ((1 / (1 + np.exp(-z)) - y) / n)[:, None]
    gW, gb = [None] * len(self.W), [None] * len(self.W)
    for i in reversed(range(len(self.W))):
        gW[i] = self.A[i].T @ delta
        gb[i] = delta.sum(axis=0)
        if i > 0:
            delta = (delta @ self.W[i].T) * self.df(self.Z[i - 1], self.A[i])
    return gW, gb

For neural networks, the activation derivatives are written in terms of the output a so they can reuse values already computed: for sigmoid it is a(1 - a), for tanh it is 1 - a², and for ReLU it is 1 if z is positive and 0 otherwise. CS231n gives the same reason for the sigmoid formula σ(1 - σ).

Gradient checking: trust backprop when it matches finite differences

A bug in backprop is the quietest kind of bug: neural networks still train, just slowly or poorly. The standard check is gradient checking: compare the analytic gradient with the central difference (L(w + e) - L(w - e)) / 2e for every weight. We try it on 3-5-4-1 neural networks with 20 random samples, for all three activations, and report the largest relative error over all weights:

numpy 1.26.0 | scikit-learn 1.3.1

[1. Gradient check: hand-written backprop vs central difference]
  sigmoid  max relative error = 4.82e-08
  tanh     max relative error = 3.62e-09
  relu     max relative error = 1.83e-09

The errors are around 1e-8 to 1e-9, small enough to conclude that backprop is correct. According to the CS231n notes, a relative error above 1e-2 almost certainly means the gradient is wrong, 1e-4 to 1e-2 is a suspicious zone worth a second look, and 1e-7 or below is safe. Note that ReLU is not differentiable at 0, so gradient checking can be off if a point lands near the kink; with continuous random data, the probability of this is very small.

XOR: why neural networks need a hidden layer and nonlinearity

XOR is the classic problem: four points (0,0), (0,1), (1,0), (1,1) with labels 0, 1, 1, 0, and no straight line separates the two classes. We compare a logistic model with no hidden layer against neural networks with one hidden layer of four tanh neurons, both with learning rate 0.5:

[2. XOR]
  no hidden layer (logistic): accuracy = 0.25
  output logits of the no-hidden-layer model: ['1.6e-16', '-6.9e-17', '-7.6e-17', '-3.1e-16'] (near 0 at all 4 points)
  1 hidden layer, tanh, epoch    0: loss = 0.7046
  1 hidden layer, tanh, epoch 1000: loss = 0.0055
  1 hidden layer, tanh, epoch 3000: loss = 0.0013
  1 hidden layer: accuracy = 1.00 | predicted probabilities: [0.0, 0.998, 0.998, 0.002]

The 0.25 needs careful reading: the linear model's logits are near 0 at all four points (around 1e-16), meaning it learned nothing, so the 0.25 carries no information about neural networks. When every logit is 0, the sign of the rounding error around 0 decides the predicted label, and the accuracy is just an artifact of floating-point noise. What matters is that no straight line classifies all four points correctly (at most 3 of 4, or 0.75), while neural networks with a hidden layer reach 1.00 as the loss drops from 0.7046 to 0.0013.

Activations and initialization in neural networks on the two-moons data

The second problem is closer to reality: make_moons generates 1,000 points in two moons with noise 0.25, split into 700 train and 300 test. The majority-class baseline scores 0.510 and logistic regression with no hidden layer scores 0.837. The neural networks have two hidden layers of 16 neurons each, trained for 2,000 epochs of full-batch gradient descent with learning rate 0.5. We cross three activations with three initializations: small (0.01 times a normal distribution), xavier (a normal distribution divided by the square root of n_in) and he (multiplied by the square root of 2/n_in).

[3. Two-moons data: 1,000 samples, noise 0.25]
  majority-class baseline: acc = 0.510
  logistic (no hidden layer): test acc = 0.837
  activation init     loss@100  loss@2000  test acc
  sigmoid   small       0.6931     0.6931     0.490
  sigmoid   xavier      0.3920     0.1828     0.930
  sigmoid   he          0.3482     0.1719     0.940
  tanh      small       0.6929     0.3084     0.837
  tanh      xavier      0.2777     0.1467     0.950
  tanh      he          0.2156     0.1418     0.940
  relu      small       0.6931     0.1810     0.937
  relu      xavier      0.2181     0.1450     0.943
  relu      he          0.2057     0.1447     0.930

Three observations from the table. First, sigmoid with small initialization learns nothing: the loss stays at 0.6931, which is ln 2, the loss of a model that always predicts probability 0.5, and the test accuracy of 0.490 is as bad as guessing. Second, for the same neural networks, switching the initialization from small to xavier lifts tanh from 0.837 to 0.950. Third, with a sensible initialization, the three activations give fairly close results (0.93 to 0.95); with a single run and 300 test samples (each wrong sample is 0.33 percentage points), a 0.01 difference is meaningless. We repeat with 5 initialization seeds (data and split unchanged) to see the spread:

  mean ± standard deviation of test acc over 5 initialization seeds:
  sigmoid  small   0.490 ± 0.000
  sigmoid  xavier  0.916 ± 0.040
  tanh     small   0.835 ± 0.002
  tanh     xavier  0.943 ± 0.004
  relu     he      0.940 ± 0.008

For these neural networks the results are stable: tanh with xavier (0.943) and ReLU with he (0.940) are nearly tied, sigmoid is slower and fluctuates more (standard deviation 0.040), and small initialization breaks sigmoid on every seed.

Uniform initialization in neural networks: symmetry is never broken

The CS231n notes on initialization explain why you must not set all weights equal: if neurons share the same weights, they compute the same output and receive the same gradient, so the symmetry is never broken. There are two cases to separate. In case (a), all zeros is even worse: the tanh activation is 0, so the gradient of every weight in earlier layers is also 0 and the network stays completely frozen. Case (b), setting every weight to a nonzero constant (0.5), illustrates the mechanism properly: the gradients are nonzero and the weights do update, but all 8 hidden neurons still receive the same gradient. We verify both:

[4. Uniform initialization]
  (a) all zeros: loss after 500 steps = 0.6931 (ln 2 = 0.6931) | test acc = 0.490
      max |W| of the two layers after training: 0.0e+00, 0.0e+00 (gradient is 0, network is frozen)
  (b) every W = 0.5: loss after 500 steps = 0.3107 | test acc = 0.837
      hidden-layer W changed: 0.500 -> [0.206, -0.693] | max difference between the 8 columns (8 neurons): 0.0e+00

In (a), the loss never leaves ln 2 and every weight is still 0 after 500 steps, because there is no signal to update anything. In (b), the network does learn (loss 0.3107, test accuracy 0.837, equal to the logistic model) and the weights moved from 0.5 to values like 0.206 and -0.693, but the difference between the 8 weight columns is exactly 0: the eight hidden neurons are still copies of each other, so the network has the representational power of one hidden neuron and cannot beat a linear model. Random initialization breaks the symmetry in neural networks; its scale decides whether gradients survive many layers, as the next section shows.

Vanishing gradients with the depth of neural networks

Why does sigmoid with small initialization stall? Backpropagation multiplies delta by W and by the activation derivative at every layer; if these factors are below 1, the product shrinks toward 0 as it travels backward. We measure the gradient norm of W at each layer in deep neural networks with 8 hidden layers of 32 neurons each, after one backward pass (layer 1 is closest to the input, layer 9 is the output layer):

[5. Gradient norm of W per layer, 8 hidden layers x 32 neurons, after 1 backward pass]
  (layer 1 = closest to input, layer 9 = output layer)
  sigmoid  xavier 1.3e-06 1.2e-05 3.9e-05 1.6e-04 7.4e-04 5.5e-03 3.1e-02 1.2e-01 5.6e-01
  tanh     xavier 4.3e-02 1.4e-01 1.5e-01 1.5e-01 2.0e-01 3.3e-01 3.9e-01 4.0e-01 2.9e-01
  relu     he     1.2e-01 4.4e-01 5.3e-01 4.1e-01 3.3e-01 4.0e-01 2.9e-01 3.5e-01 3.3e-01

With sigmoid, the gradient at the first layer (1.3e-6) is about 4.3e5 times smaller than at the output layer (0.56): the early layers receive almost no learning signal. The sigmoid derivative never exceeds 0.25, so each layer multiplies in another small factor; this is exactly what CS231n warns about when it says sigmoid saturates and makes gradients nearly zero. Tanh shrinks by about 7 times (0.29 down to 0.043), while ReLU with He is almost flat (0.33 down to 0.12). Likewise, measuring the standard deviation of activations per layer shows that initialization decides the forward pass too:

  standard deviation of activations per layer (8 hidden layers):
  tanh  small  0.010 0.001 0.000 0.000 0.000 0.000 0.000 0.000
  tanh  xavier 0.496 0.406 0.345 0.350 0.322 0.283 0.248 0.184
  relu  xavier 0.430 0.319 0.191 0.146 0.091 0.052 0.034 0.025
  relu  he     0.608 0.638 0.540 0.582 0.515 0.414 0.380 0.393

In a neural network using tanh with small weights, the activations collapse to 0 from the second layer onward. With ReLU, xavier-style initialization (dividing by the square root of fan-in, in the spirit of Glorot and Bengio's Understanding the difficulty of training deep feedforward neural networks, AISTATS 2010) shrinks the activations from 0.430 to 0.025 over 8 layers, while He keeps them between 0.38 and 0.64. The reason is that ReLU zeroes about half of its outputs, so it needs twice the variance to compensate; CS231n states that the 1/sqrt(n) derivation assumes a zero mean, which is not true for ReLU, and recommends variance 2/n following Delving Deep into Rectifiers by He et al. (2015).

Comparing neural networks with PyTorch autograd

The final check: build the same neural networks (2-8-8-1, tanh) in PyTorch with exactly the weights of the NumPy version, compute the loss on the same 64 samples, and compare the gradient of every weight matrix:

[6. Comparison with PyTorch autograd]
  torch 2.0.1
  NumPy loss = 0.742936 | PyTorch loss = 0.742936
  max gradient difference between the two: 2.78e-17

The loss matches to 6 decimal places and the gradients differ by only about 1e-17, which is the rounding error of 64-bit floats. Our hand-written backprop for neural networks and autograd perform the same computation: autograd just automates building the graph and applying the chain rule the way we wrote it by hand.

When not to build neural networks yourself, and common mistakes

The NumPy neural networks here are for understanding, not for production: it has no GPU, no optimizers like Adam, no regularization, and full-batch training does not scale to large data. In real projects use PyTorch or an equivalent framework; for tabular data, tree models are often stronger, as seen in Gradient boosting with XGBoost and LightGBM: when it beats deep learning. Common mistakes when writing neural networks:

• Wrong matrix dimensions: always check that a gradient has the same shape as its parameter; mixing up where to transpose is the most common bug.
• Forgetting to divide by n when the loss is a mean: the gradient is n times too large, and the effective learning rate changes with batch size.
• Applying sigmoid before cross-entropy: it overflows easily when |z| is large. The code uses the stable form max(z,0) - z*y + log1p(exp(-|z|)), the same idea as PyTorch's BCEWithLogitsLoss.
• Uniform initialization (all zeros or one constant) makes neurons identical forever, and initialization that is too small in a deep network leaves the loss stuck at ln 2, as measured above.
• Skipping gradient checks: run gradient checking on a small network with a few samples before real training, then turn it off because it is very slow.
• Dead ReLU: CS231n warns that one oversized update can push a neuron into a state where it never activates again; monitor the fraction of dead neurons and lower the learning rate or use Leaky ReLU.
• Drawing conclusions from one seed: as the table above shows, a 0.01 difference is noise; repeat over several seeds.

Neural networks from scratch checklist

• Each layer of neural networks is a linear transform followed by a nonlinearity; without the nonlinearity, many layers equal one.
• The forward pass stores z and a; the backward pass uses the chain rule: dW = a_prevᵀ delta, delta_prev = (delta Wᵀ) times f'(z).
• Always gradient-check with finite differences on a small network; compare with autograd when you can.
• Initialize neural networks randomly with fan-in scaling: xavier for tanh and sigmoid, He for ReLU; never all zeros.
• Monitor neural networks through gradient norms and activation standard deviations per layer to catch vanishing or exploding signals.
• Use a numerically stable loss (logit plus fused cross-entropy) and compare against the majority-class baseline.

The full code is in code/s7-21 (nn_scratch.py, output.txt, requirements.txt), run with python nn_scratch.py. Every number about the neural networks in this article is real output from an Apple M4, 16 GB, with fixed seeds (two runs gave identical output); the data is synthetic.

Next article

Stable deep learning training: optimizers, learning rates and regularization. We have seen that initialization and activation functions decide whether neural networks can learn at all; the next article moves on to the tools that keep training stable as neural networks grow: Adam, learning rate schedules and regularization.