Backpropagation From First Principles: What a Neural Network Actually Does

model.fit(x, y) is where my understanding of neural networks used to stop, because that one call hides everything I wanted to know. What is being adjusted? By how much? How does a number computed at the output end up changing a weight five layers back?

It turns out the concepts are simpler than the framework APIs suggest: a network is a pile of multiply-add, its error is a single number, and training is repeatedly asking “which direction should each weight move to make that number smaller?” All of it fits in your head at once.

To force myself to actually hold it there, I built a multi-layer perceptron in plain Java that learns to recognise handwritten digits — no TensorFlow, no DL4J, not even a matrix library. The only dependency is a logging framework: github.com/rustamniraula90/hand-written-number-recognition.

These are my notes from working through it.

---

Concept 0: A network is a function you can shape

Before any learning happens, a network is just a parameterised function. Mine turns 784 numbers (a 28×28 image, flattened) into 10 numbers (confidence per digit) through two intermediate layers:

784  ──[W₁ 784×128]──▶  128  ──[W₂ 128×64]──▶  64  ──[W₃ 64×10]──▶  10
      +bias, ReLU              +bias, ReLU           +bias, softmax

Those weight matrices are the model — 784×128 + 128×64 + 64×10 ≈ 109,000 numbers. Training means finding values for them. Everything else is bookkeeping.

I declared the shape rather than hardcoding it, mostly to keep that idea visible while I worked:

MLPNetworkOptimized network = MLPNetworkOptimized.builder(2)
        .addInputLayer(28 * 28)
        .addOutputLayer(10)
        .addHiddenLayer(128)
        .addHiddenLayer(64)
        .setLearningRate(0.001)
        .setBatchSize(10)
        .setEpoch(20)
        .build();

Why initialization isn’t arbitrary

Start every weight at zero and every neuron in a layer computes the same thing forever — identical gradients, identical updates, no differentiation. Start them too large and activations compound layer over layer until they saturate. So weights start random, but scaled to layer width:

public static double xavier(int in, int out) {
    return random.nextGaussian() * Math.sqrt(2.0 / (in + out));
}

That’s Xavier initialization: variance inversely proportional to fan-in plus fan-out, so signal magnitude stays roughly constant as it propagates.

Design note. I pass the network’s overall input and output sizes for every layer rather than each layer’s own fan-in and fan-out, so all layers share one variance. It trains fine, which taught me something useful about the shape of this parameter: initialization is a conditioning concern and fairly forgiving of imprecision, unlike a gradient sign error, which is fatal immediately. Per-layer fan-in/fan-out is the textbook form and the easy refinement.

Concept 1: The forward pass is multiply, add, squash

Each layer computes output = activation(W·input + b). Written out per neuron, the whole idea is eight lines:

double z = 0.0;
for (int j = 0; j < neuron.weights.length; j++) {
    z += neuron.weights[j] * activations[j];
}
z += neuron.bias;
neuron.output = relu(z);

The activation function is the only nonlinear step, and without it the network collapses. Stack two linear layers and you get another linear layer — depth would buy nothing. ReLU (max(0, x)) is the usual choice: trivially differentiable, and it doesn’t saturate for positive inputs.

I used the bounded variant, ReLU6:

public static double relu(double x) {
    return Math.clamp(x, 0, 6);
}

public static double reluDerivative(double x) {
    return x > 0 ? 1 : 0;
}

Design note. The ceiling of 6 is the trick MobileNet uses to keep activations bounded, and I liked it for stability. The paired derivative here is plain ReLU’s, which returns 1 above the clamp where ReLU6’s true derivative is 0. In this network it makes no observable difference — the small learning rate and scaled init keep pre-activations well below 6 — but writing the two functions next to each other is what made the general rule stick for me: an activation and its derivative are a matched pair, and changing one means revisiting the other. Clamping the derivative at 6 as well is a one-line change and the right call if the ceiling ever binds.

The output layer answers a different question

Hidden layers produce features. The output layer has to produce a probability distribution over 10 classes — ten numbers, non-negative, summing to 1. That’s softmax: exponentiate, then normalise.

public static double[] softmax(double[] x) {
    double sum = 0;
    double[] result = new double[x.length];
    for (int i = 0; i < x.length; i++) {
        result[i] = Math.exp(x[i]);
        sum += result[i];
    }
    for (int i = 0; i < x.length; i++) result[i] /= sum;
    return result;
}

Exponentiating first is what makes it a soft max: it amplifies differences, so a logit slightly ahead becomes a lot more probable, without the hard winner-take-all of picking the maximum.

Design note. Library implementations subtract the maximum logit before exponentiating, because Math.exp overflows around 710. Mine skips that, which is fine for MNIST-scale logits and is the first thing I’d add if inputs got wilder. Finding out why that subtraction exists — by not having it — was more instructive than reading about it.

Concept 2: Error has to be one number

You can’t optimise ten outputs at once. Gradient descent needs a single scalar to reduce, so the ten probabilities and the correct answer collapse into one loss. For classification that’s cross-entropy:

public static double crossEntropy(double[][] actual, double[][] expected) {
    double loss = 0;
    for (int i = 0; i < actual.length; i++) {
        for (int j = 0; j < actual[i].length; j++) {
            loss += expected[i][j] * Math.log(actual[i][j]);
        }
    }
    return -loss;
}

With one-hot targets every term but the correct class multiplies by zero, so the loss is just -log(probability assigned to the right answer). Assign the truth 0.9 and you’re barely penalised; assign it 0.01 and the log punishes you enormously. That asymmetry is the point: confident wrong answers should hurt far more than uncertain ones.

Concept 3: Backpropagation is the chain rule, applied backwards

This is the part frameworks hide, and the part I built this for. We have one number (the loss) and 109,000 weights, and we need ∂loss/∂w for each. Computing them independently would mean 109,000 passes through the network. Backpropagation gets all of them in one backwards pass by reusing work.

The insight that unlocked it for me: for any weight, ∂loss/∂weight = (how much the loss changes with that neuron’s output) × (what that neuron’s input was). The first factor — call it delta — is the only hard part, and each layer’s delta can be computed from the layer after it. So you compute deltas once, right to left, and every weight update reads off a delta you already have.

Step one: the delta at the output layer. This is where the earlier choices pay off. Differentiate cross-entropy composed with softmax and all the exponentials and logarithms cancel:

delta[i][j] = outputs[i][outputs[i].length - 1][j] - expectedOutputs[i][j];

prediction − target. That’s the entire output-layer gradient. Softmax and cross-entropy are paired because their derivatives collapse like this, and seeing it fall out of the algebra was the most satisfying moment of the project.

Step two: update this layer, then push delta backwards. For each weight, gradient = delta × the input that flowed through it, averaged over the batch:

for (int i = weight.length - 1; i >= 0; i--) {          // layers, back to front
    for (int j = 0; j < weight[i].length; j++) {
        for (int k = 0; k < weight[i][j].length; k++) {
            double gradientSum = 0.0;
            for (int l = 0; l < batchSize; l++) {
                gradientSum += delta[l][k] * inputs[l][i][j];
            }
            weight[i][j][k] -= learningRate * (gradientSum / batchSize);
        }
    }
    // ... bias update, same shape ...
    if (i > 0) {
        // delta for the previous layer: weighted sum of this layer's deltas,
        // scaled by the derivative of that layer's activation
        newDelta[j][k] = sum * MathUtil.reluDerivative(outputs[j][i-1][k]);
        delta = newDelta;
    }
}

Three ideas live in that loop:

• weight -= learningRate * gradient is gradient descent. The gradient says which direction increases loss; step the other way, scaled small so you don’t overshoot.
• The previous layer’s delta is this layer’s deltas weighted by the connections between them — error attributed backwards in proportion to contribution.
• Multiplying by the activation derivative closes the chain rule. A neuron whose ReLU output was zero contributed nothing, so it receives no blame.

Concept 4: Mini-batches, and why batch size is a real choice

You could update after every sample (noisy, slow) or after all 10,000 (smooth, expensive, glacial progress). Mini-batching splits the difference: accumulate gradients over N samples, then step once.

I implemented it as buffer-and-fire rather than a nested loop — the backward pass is called after every sample but only acts when the buffer fills:

private void backwardPass(int epoch) {
    if (index == batchSize) {
        index = 0;
        double loss = MathUtil.crossEntropy(finalOutput, this.expectedOutputs);
        tuneParameters();
        this.inputs = new double[batchSize][weight.length][];
        this.outputs = new double[batchSize][weight.length][];
        this.expectedOutputs = new double[batchSize][];
    }
}

It keeps the training loop flat, at the cost of coupling one index to three buffers. Averaging over the batch (gradientSum / batchSize) is what makes the learning rate independent of batch size — otherwise doubling the batch would double every step.

Design note. predict() calls the same forwardPass(), which writes into those training buffers at the current index. Train-then-predict is exactly how the program runs, so this is fine as written; separating inference state from training state is the change I’d make before ever interleaving the two.

The lesson I didn’t expect: data layout decides what you can express

I wrote this network twice, and the comparison taught me more than either version alone.

The first, MLPNetwork, is the design an OO instinct reaches for: a Layer of Neuron objects, each owning its weights, bias, input and output. It reads beautifully. Watch where the weight update sources its input activation:

double input = layerIndex == 0 ? inputs[k][j] : prevLayer.neurons[j].output;
gradientSum += delta[k][i] * input;

For the first layer, inputs[k][j] indexes sample k of the batch. Everywhere else, prevLayer.neurons[j].output is one field on a shared object — it holds whatever the most recent forward pass left there. A neuron object can hold one activation; mini-batch backprop needs one per sample in flight. The structure quietly asserts something the algorithm doesn’t agree with.

The second version, MLPNetworkOptimized, replaced objects with flat arrays and kept per-sample buffers:

private double[][][] inputs;   // [batchIndex][layer][neuron]
private double[][][] outputs;

which is why its update can say inputs[l][i][j] — sample l, layer i, neuron j.

I named it “Optimized” because I was chasing speed: flat double[][][] beats chasing pointers through a hundred thousand objects, and the JIT keeps the inner loop in registers. What I got alongside the speed was a structure that could hold the state the maths actually needs. That reframed the whole exercise for me, and it’s the part that generalises furthest beyond neural networks: choosing the data layout is choosing which behaviours are expressible. Neuron.output looked like good encapsulation and was really a claim about dimensionality.

It also explained something I’d taken for granted. Every framework represents a batch as a tensor with a leading batch dimension, and I’d always read that as a vectorisation detail. It isn’t only that — the batch axis has to exist somewhere, and putting it in the data structure is what makes the gradient definable.

Before any of it: the data is a binary blob

One more concept, because it’s the first wall I hit. MNIST doesn’t ship as CSV — it’s IDX, a big-endian header followed by raw bytes. There’s no library call; you assemble integers from bytes yourself:

private static int readInt(FileInputStream fis) throws IOException {
    return (fis.read() << 24) | (fis.read() << 16) | (fis.read() << 8) | fis.read();
}

Then you assert the magic number (2051 for images, 2049 for labels) so a wrong file fails loudly instead of training on garbage, and you normalise pixels into a small range so the first layer’s weights don’t have to absorb a factor of 255. Mine divides by 254.0, which puts a saturated pixel a hair above 1.0 — immaterial at this scale, and a good reminder that scaling constants deserve the same attention as the algorithm. The most effective verification I found was writing a few parsed samples back out as PNGs and looking at them; a transposed dataset is instantly obvious to an eye and invisible to an assertion.

Design note on the run configuration. main trains on the t10k files (10,000 images) and evaluates on the train files (60,000). Conventionally that’s the other way round, and I kept it because the property that matters held — the evaluation set is strictly disjoint from what it learned on — and it makes for a harsher test: less data to learn from, six times as much to be judged on. Shuffling between epochs is the other knob I left alone; with fixed batch composition repeating 20 times, adding a shuffle is the cheapest available improvement.

What’s worth keeping from this

• A network is multiply-add plus one nonlinearity; depth without the nonlinearity is worthless.
• Softmax and cross-entropy are paired because their gradients collapse to prediction − target.
• Backprop is one reverse pass that reuses each layer’s delta, and that reuse is the entire reason training is tractable.
• An activation and its derivative are one decision, not two.
• Frameworks aren’t hiding the calculus. They’re hiding numerical stability guards, initialization scaling, shuffling and the batch axis — and the clearest way to see the shape of that work is to leave it out and notice what starts mattering.