Neural Networks
from Scratch

01Why neural networks?

Most of the code you've written tells a computer exactly what to do, step by step, branch by branch. Neural networks flip that around. You don't write the rules. You write a flexible function with millions of dials, then turn the dials automatically until the function does what you want: the rules are fitted to examples instead of programmed.

Take an example that's notoriously hard to hand-code: deciding whether a photo contains a cat. You can't write if (pixels look cat-shaped). There's no rule that survives contact with kittens, sphinxes, plush toys, or weird camera angles. But you can collect a few thousand labeled photos, define a function with a few million tunable parameters, and let an optimizer find the parameters that classify the training set correctly. With enough capacity and data, it can generalize to photos it has never seen.

The earliest version of this idea was Frank Rosenblatt's perceptron in 1958^[1]. It worked on toy problems, hit a wall when Minsky and Papert proved in 1969 that a single-layer perceptron can't compute XOR^[3], and neural-network research lost most of its funding and attention through the 1970s. Then in 1986 Rumelhart, Hinton, and Williams published Learning representations by back-propagating errors^[2], the paper most often credited with making backpropagation the standard way to train multi-layer networks.

A footnote on the standard story: backpropagation had already been worked out by Seppo Linnainmaa in 1970 (as reverse-mode automatic differentiation)^[16], and Paul Werbos proposed applying it to neural networks in his 1974 PhD thesis. Rumelhart, Hinton and Williams didn't invent it; their paper is where it reached a wide audience, with experiments showing hidden layers learning useful internal representations.

02The neuron: a tiny decision maker

An artificial neuron does three things, in this order: it takes some numbers in, mixes them together with a weighted sum, and decides whether to "fire" based on whether the sum crosses a threshold.

Concretely: given input numbers x₁, x₂, …, x_n, the neuron has a matching set of weights w₁, w₂, …, w_n and a bias b. The neuron computes a weighted sum:

Take each input, multiply it by its weight, add all those products together, then add the bias. The result z is a single number. Hover any letter above to see what it means.

In vector form, z = w · x + b. That's the dot product of two arrays plus a constant, a few lines of code in any language.

Then the neuron applies an activation function to z to decide its output. Rosenblatt's perceptron used the simplest one possible, a hard threshold:

If the pre-activation z is zero or positive, the neuron fires and outputs 1. If z is negative, it stays silent and outputs 0. That's the original perceptron rule, in one if-statement.

In 2D this has a clean geometric meaning. The equation w₁x₁ + w₂x₂ + b = 0 describes a straight line. Points on one side give z > 0 and the neuron fires; points on the other side give z < 0 and it stays silent. The weights tilt and rotate the line; the bias shifts it.

Play with the controls. Every move changes the line, which changes how the point is classified:

One neuron, one straight line. That is the entire expressive power of a single perceptron. It can separate points that are linearly separable, and nothing else. No choice of weights makes a single perceptron compute XOR^[3]. That limitation is why we have to stack neurons. But before we get there, we need a fix for a more immediate problem.

03Smoothing it out: the sigmoid neuron

A perceptron has a problem: its output jumps. Push z from −0.001 to +0.001 and the output flips from 0 to 1. There's no notion of "almost firing" or "barely firing." That makes the neuron impossible to nudge: a tiny weight change either changes nothing at all, or changes everything. In calculus terms, the step's slope is zero everywhere except at the jump, so a slope-following learning rule gets no signal.

Learning requires gentle nudges: small, predictable changes in weights producing small, predictable changes in output. Rumelhart, Hinton and Williams used a smooth S-shaped curve in place of the hard threshold, the logistic sigmoid:

Plug z into this and you always get a number between 0 and 1. Very large z → answer near 1. Very negative z → answer near 0. Right at z = 0, the answer is exactly 0.5.

The sigmoid takes any real number and squashes it into the range (0, 1). When z is very positive, σ(z) is near 1. When z is very negative, it's near 0. In the middle, around z = 0, the curve is gently sloped, which means nudging z by a little nudges the output by a little. That gentle sloping is what makes learning possible.

The slope itself matters a lot; we'll need it for backpropagation in a few sections. The derivative of the sigmoid has a tidy closed form:

The slope of the sigmoid at any input equals the sigmoid's value at that input, times one minus that value. Handy trick: if we already know σ(z), we get its slope for free.

It's maximal (0.25) right at the inflection point z = 0 and decays toward zero in either direction. Try it: toggle "show derivative" and slide the input.

04Stacking neurons into networks

One neuron can carve the input space with one straight line. Two neurons in parallel can carve it with two lines. Stack a second layer of neurons that takes the first layer's outputs as their inputs, and, given enough neurons, you can carve out arbitrarily complex regions: circles, spirals, anything.

A feedforward neural network is exactly that: layers of neurons, where every neuron in layer l takes inputs from every neuron in layer l − 1. The first layer is the input layer, the last is the output layer, and the ones in between are hidden layers.

The diagram above is a 2-3-2 network: two inputs, three hidden neurons, two outputs. Every line is one weight. Every circle (except the inputs) holds one bias. The whole thing is a function from an input vector to an output vector.

We'll use a compact notation that makes the math much easier to write. For layer l, collect all that layer's weights into a matrix W^(l), all the biases into a vector b^(l), and all the activations into a vector a^(l). The forward pass through one layer is then:

For each layer: (1) multiply the previous layer's outputs by this layer's weight matrix, (2) add this layer's biases, getting the pre-activations z. Then (3) apply the activation function σ to each pre-activation, getting this layer's outputs. The superscript (l) just labels which layer we're talking about.

Two operations per layer: a matrix-vector multiply plus a vector add, then an element-wise activation. The whole forward pass is a few of those in a row, a function that turns an input vector into an output vector with the network's parameters baked in.

05Forward propagation

Forward propagation runs the equations above from input to output: compute z^(l) and a^(l) one layer at a time and pass the activations forward.

Step through it on the network below. Set the inputs, then click STEP to advance one phase at a time, or RUN to animate the whole pass. The output layer applies softmax instead of sigmoid because we want a probability distribution over two classes (more on this in the next section):

In code, forward propagation is a short loop over the layers:

function forward(network, x) {
  let a = x;                              // activations of the previous layer
  const activations = [a];                // keep every layer's activations for backprop

  for (let l = 0; l < network.W.length; l++) {
    const W = network.W[l];               // weight matrix for this layer
    const b = network.B[l];               // bias vector for this layer
    const z = matVecMul(W, a).map((zi, i) => zi + b[i]);
    a = z.map(sigmoid);                  // element-wise activation
    activations.push(a);
  }
  return activations;                      // last entry is the network's output
}

activations[activations.length - 1] is the prediction. We keep the intermediate values around because we'll need them when we compute gradients.

06Measuring wrongness: loss functions

A network with random weights makes random predictions. To improve it, we need a single number telling us how wrong it is. Then we can do calculus on that number to figure out which way to nudge each weight to make it smaller. That single number is called the loss (or cost).

Mean squared error

The simplest loss function, used in Nielsen's book^[4] and most introductory examples, is mean squared error. For an example with true output y and prediction ŷ:

For each training example, compute how far the prediction is from the truth, square that distance, then average across all examples (the ½ is a cosmetic factor that makes the derivative tidier). Bigger loss = more wrong.

Average the squared distance between prediction and target, over all training examples. The ½ is there for convenience: it cancels out when we take the derivative, so we get (ŷ − y) instead of 2(ŷ − y).

Cross-entropy: the better choice for classification

MSE works, but for classification problems there's a better loss called cross-entropy. The intuition: instead of measuring numeric distance, measure how surprised the model is by the correct answer.

First we turn the raw output numbers (called logits) into a probability distribution using softmax:

Raise e to the power of each logit (this makes them positive), then divide each by the total. The result is a list of probabilities that sum to 1 and reflect the relative size of the original logits.

Exponentiate each logit (forces it positive), then divide by the sum so the whole vector adds to 1. Now we have probabilities. The cross-entropy loss is then:

Only one term in the sum is nonzero (the one for the correct class, where y_i = 1). The loss is simply −log(prediction for the correct class). If the model assigns high probability to the right answer, log is near 0 and the loss is small. If it assigns near 0, log goes to −∞ and the loss explodes.

If the true class is class 0 (so y = [1, 0]) and the model predicts ŷ = [0.9, 0.1], the loss is −log(0.9) ≈ 0.105. If the model says ŷ = [0.1, 0.9], the loss is −log(0.1) ≈ 2.303. Confident wrong answers are punished heavily; confident right ones are rewarded. That asymmetry is what we want during training.

07Learning by descent

Now we have a function (the loss) that takes the network's weights and returns a single number. Training means finding the weights that make that number small. With millions of weights, brute-force search is hopeless. The standard strategy is gradient descent.

Imagine you're standing on a hilly landscape in dense fog. You can't see the bottom of the valley, but you can feel the slope under your feet. The gradient is that slope written as a vector: it points uphill, so you step the opposite way, feel the slope again, and step again. Repeat until the ground is flat, which is a low point, though not necessarily the lowest one.

The mathematical version is exactly that. For each weight w, we compute ∂C / ∂w, which is how much the loss changes per unit of weight change. (If the curly-d notation is new to you, expand the box below before reading on.) Then we nudge:

The arrow ← means "assign." This line replaces the old weight with: old weight, minus a small fraction (η) of how much the loss responds to changes in that weight. Do this for every weight, every step, and the loss goes down, provided η is small enough.

The constant η ("eta") is the learning rate, how big a step to take. Too small and training crawls; too big and each step overshoots the valley floor, so the weights oscillate or diverge. It is usually the first setting to tune.

Click anywhere on the loss surface below to drop a starting point, then hit DESCEND; the yellow path traces each step. On the elongated valley the steep direction has curvature 12, so the path bounces between the walls once η passes 1/12 ≈ 0.08 and flies off once it passes 2/12 ≈ 0.17 (the default 0.18 diverges there; lower it and descend again). The other surfaces show a saddle point and two separate minima:

The same landscape, in three dimensions

The contour map above is the usual way to draw it on paper. With the loss as height, the same kind of surface is a 3D landscape. Drag the surface below to rotate it, then hit Descend to see the path slide down the side of the hill:

08Beyond SGD: smarter optimizers

Plain SGD has problems. In a long narrow valley it either zigzags between the walls (large η) or crawls along the nearly flat floor (small η). Near a saddle point it slows almost to a stop, and on a flat region it crawls. Real loss landscapes have all three. Most modern training uses SGD with momentum or an adaptive method such as Adam: variants that pass the gradient through a small running memory before applying it.

Momentum: a ball with mass

Momentum keeps a running average of the recent gradient direction and uses that average to update the weights. The intuition: vanilla SGD is a marble that stops the instant you stop pushing it. Momentum SGD is a ball with mass: it picks up speed in directions where the gradient consistently points the same way, and dampens directions where the gradient keeps flipping sign.

Each step, blend β (~90%) of the previous velocity with a fresh nudge from the current gradient. Then move along that blended velocity. In a long valley, the gradient repeatedly points down the valley floor; momentum accumulates speed in that direction. In zig-zag directions, alternating gradients cancel each other out.

RMSProp: adapt the step size per parameter

Different parameters can have wildly different gradient scales. A weight in the output layer might routinely see gradients of size 1.0; a weight five layers deep might see 0.001. Using the same learning rate for both is a compromise neither one wants. RMSProp^[12] tracks a running average of the squared gradient for each parameter and divides each update by the square root of that average. Big-gradient parameters get smaller steps; small-gradient parameters get bigger ones.

Maintain an exponential moving average of how big each parameter's gradient has been recently. Then scale each update by 1 over the square root of that. The tiny ε (typically 10⁻⁸) prevents division by zero when a parameter's gradients have all been zero so far.

Adam: momentum + RMSProp + bias correction

Adam^[13] combines the previous two ideas. It keeps both a momentum (first-moment) running average and an RMSProp-style (second-moment) running average. Then it applies bias correction so the early-training estimates aren't biased toward zero. Adam and its variant AdamW are the most common default optimizers in deep learning today.

Track both how the gradient has been pointing on average (m) and how big it has been on average (v). Correct for the fact that both averages start at zero (the 1 − β^t terms). Take a momentum-style step, scaled per-parameter by the inverse of the running gradient magnitude.

Racing the four optimizers down the same surface shows the differences. On the Rosenbrock function (a banana-shaped valley that is a standard optimizer stress test; the widget uses a rescaled version) at the default η = 0.05, plain SGD has covered only the first stretch of the curve after 320 steps, Momentum and Adam zigzag across the valley before settling onto its floor, and RMSProp follows the floor from the start and gets closest to the minimum. Toy surfaces like these show how each update rule behaves; they don't rank optimizers on real networks.

09Backpropagation: how we actually get the gradients

Gradient descent needs ∂C/∂w for every weight in the network, potentially millions of them. The naive approach nudges one weight, reruns the whole forward pass, and measures the change in loss, once per weight: millions of forward passes per update. Backpropagation is the algorithm that gets us all the gradients in one efficient sweep, by reusing intermediate values from the forward pass and chaining derivatives layer by layer.

Define δ_j^(l) (pronounced "delta-j-l") as the error at neuron j in layer l: how much the loss changes per unit of the neuron's pre-activation.

Nielsen distills backprop into four equations^[4]:

Output error = (how loss responds to output activations) ⊙ (activation slope at the output). This is the textbook form, written for any output activation + any loss. It's what you'd use for sigmoid + mean-squared error.

For us, the output layer has no activation (the code emits raw logits) and the loss is softmax + cross-entropy. Run those two through the chain rule and the σ′ piece, the softmax Jacobian, and most of the algebra cancel. What you're left with is a one-liner that we already met in §6:

Predicted probabilities minus the one-hot label. No σ′ to evaluate, no Jacobian to assemble. The code in §12 uses exactly this line.

Error at layer l = (next layer's weights transposed × next layer's error), times this layer's activation slope. The transpose "flips" the weights so we can send the error back through them. This is the heart of backprop: we already computed δ for layer l+1, so we can use it to compute δ for layer l. Then repeat all the way to the input.

The bias gradient is just the error itself. No multiplication needed. The bias adds 1 unit per 1 unit of change, so its sensitivity is exactly the error of its neuron.

Weight gradient = (upstream activation) × (downstream error). A weight only matters as much as the input flowing through it (activation of the source neuron) times how off-target its destination neuron is (the error δ at the destination).

The symbol ⊙ is element-wise (Hadamard) product: multiply two equal-length vectors slot by slot. ∇_aC in BP1 is the gradient of the loss with respect to the output post-activations, and depends on which loss you use; for the softmax + cross-entropy combination we use here, the whole BP1 expression collapses to BP1′ above.

Here's what each equation is actually saying:

• BP1. The textbook form, valid for any output activation and loss. With softmax + cross-entropy, BP1 collapses to BP1′: the output error is simply predicted probabilities minus the one-hot target, with no activation-slope factor at all.
• BP2. To get the error at the layer below, take the next layer's error, project it back through the weights (transpose multiplication), and weight by that layer's activation slope.
• BP3. The bias gradient is the error itself, because the bias affects only that one neuron's pre-activation, with derivative 1.
• BP4. The weight gradient is the error of the downstream neuron times the activation of the upstream neuron: a connection gets credit or blame in proportion to the signal that flowed through it.

The four equations are all consequences of the chain rule from multivariable calculus^[4]. The algorithm follows directly from them:

1. Run a forward pass; cache every z^(l) and a^(l).
2. Compute δ^(L) at the output using BP1.
3. For l = L − 1 down to 1: compute δ^(l) using BP2.
4. For each layer, accumulate gradients using BP3 and BP4.
5. Update weights: W^(l) ← W^(l) − η · ∂C/∂W^(l) and similarly for biases.

Per layer, the forward pass does one matrix-vector product and the backward pass does two (one to send δ back through the weights, one outer product for the weight gradient). A full gradient therefore costs about three forward passes, however many weights the network has, instead of one forward pass per weight. That ratio is why backprop scales to billions of parameters.

10Why nobody uses sigmoid anymore

Now that we have backprop in front of us, look back at BP2: every layer multiplies the error by σ′(z). The maximum value of σ′ is 0.25, at the inflection point, and it falls toward zero as |z| grows (0.105 at z = 2, 0.018 at z = 4).

Now imagine a 10-layer network. Each layer multiplies the error by a σ′ of at most 0.25, typically less. After 10 layers those factors alone contribute 0.25¹⁰ ≈ 10⁻⁶: unless the weights happen to be large enough to compensate, the gradient reaching the first layer is around a millionth of the gradient at the output, and those early weights barely move. This is the vanishing gradient problem^[5], and through the 1990s and 2000s it made deep sigmoid and tanh networks very hard to train.

Several fixes arrived around 2006 to 2015 (layer-wise pretraining, better initialization, normalization layers, residual connections). The simplest, popularized around 2010 to 2011, is to replace sigmoid in the hidden layers with the rectified linear unit:

If z is positive, output z. If z is zero or negative, output 0. In code it's one line: return z > 0 ? z : 0.

The derivative is 1 for positive inputs and 0 for negative ones, so the activation contributes no shrinking factor: the gradient passes through unchanged for any neuron that's "on." Combined with better initialization, ReLU made much deeper networks trainable with plain gradient descent.

ReLU has a known failure mode: a neuron whose pre-activation is negative for every training example produces zero output and zero gradient, so it stops learning for good ("dying ReLU"). Leaky ReLU patches this with a small slope below zero (the widget uses 0.1z; 0.01z is also common). GELU multiplies z by the Gaussian cumulative distribution function, giving a smooth curve that behaves like ReLU away from zero. GPT-2 uses it in every feed-forward block^[6]; many recent large language models use gated variants such as SwiGLU instead. Try them in the activation widget above to see the shapes.

Initialization matters too

If all your weights start at zero, every neuron in a layer computes the same thing and learns the same thing. That's a symmetry the gradient can never break, so we initialize randomly. But how big should the random numbers be?

Glorot and Bengio showed in 2010^[7] that for tanh-like activations, keeping the signal's variance roughly constant from layer to layer on the forward pass needs weights with variance 1/n_in, where n_in is the number of inputs to the layer. (Their published "Xavier" rule averages that with the backward-pass condition to get 2/(n_in + n_out); many libraries use the simpler 1/n_in form.) For ReLU networks, He et al. showed in 2015^[8] that you need twice that, variance 2/n_in, because ReLU zeroes the negative half of its inputs and so halves the signal's second moment at every layer. Get the scale wrong on a deep network and the activations grow or shrink geometrically with depth.

The code below uses He initialization (variance 2/n_in) when the activation is ReLU and variance 1/n_in otherwise, drawn from a uniform distribution scaled to hit that variance.

11The training loop, end to end

Everything above bolts together into a short program. In pseudocode:

net = makeNetwork([2, 8, 8, 2])              // shape: 2 in, 2 hidden layers of 8, 2 out
data = loadTrainingData()

for epoch in 1..N:
  shuffle(data)
  for batch in batches(data, size=32):
    batchGradients = zeros_like(net)
    for (x, y) in batch:
      // forward pass (caching every layer's activations), loss, backward pass
      exampleGradients = backprop(net, x, y)
      accumulate(batchGradients, exampleGradients)
    applyUpdate(net, batchGradients / len(batch), lr=η)
  evaluate(net, validation_set)            // loss on held-out data, see §14

Each pass of the outer loop over the whole dataset is one epoch. The same steps (forward pass, loss, gradient via backprop, parameter update, repeat) train everything from LeCun's 1989 USPS digit classifier^[14] to today's large language models; what changes are the layer types, the optimizer, and the scale. Karpathy's Zero-to-Hero series^[15] follows this loop from a single neuron up to a small GPT.

12A working network, in your language

Below is the complete network (initialization, forward pass, backprop, SGD update) in four languages: JavaScript, Python, C++, and Rust. Pick whichever you're most comfortable with. The algorithm is identical in all four; only the syntax changes. Every variable has a descriptive name, and every non-obvious line gets a comment.

// ──────────────────────────────────────────────────────────────
// network.js · feed-forward neural net trained with backprop + SGD
// ──────────────────────────────────────────────────────────────

// Activation functions, paired with their derivatives.
// Note: each derivative is expressed in terms of the *post-activation*
// value `a`, not the pre-activation `z`. For our three choices this
// makes the derivative cheap to compute from values we already have.
const ACTIVATIONS = {
  relu: {
    fn:                       z => z > 0 ? z : 0,
    derivativeFromActivation: a => a > 0 ? 1 : 0,
  },
  tanh: {
    fn:                       z => Math.tanh(z),
    derivativeFromActivation: a => 1 - a * a,
  },
  sigmoid: {
    fn:                       z => 1 / (1 + Math.exp(-z)),
    derivativeFromActivation: a => a * (1 - a),
  },
};

// Build a fresh network. `layerSizes` is an array of integers,
// e.g. [2, 8, 8, 2] = 2 inputs, two hidden layers of 8, 2 outputs.
function makeNetwork(layerSizes, activationName = 'relu') {
  const weights = [];   // weights[layer][neuron][input]
  const biases  = [];   // biases[layer][neuron]

  for (let layer = 1; layer < layerSizes.length; layer++) {
    const fanIn  = layerSizes[layer - 1];   // inputs coming in
    const fanOut = layerSizes[layer];       // neurons in this layer

    // Target weight variance: 2/fanIn for ReLU (He init), 1/fanIn otherwise
    // (Xavier-style). Both keep activations from blowing up or vanishing
    // with depth. A uniform draw on [-limit, limit] has variance limit²/3,
    // so limit = sqrt(3 × targetVariance).
    const targetVariance = activationName === 'relu' ? 2 / fanIn : 1 / fanIn;
    const limit = Math.sqrt(3 * targetVariance);

    const layerWeights = [];
    for (let neuron = 0; neuron < fanOut; neuron++) {
      const weightRow = [];
      for (let input = 0; input < fanIn; input++) {
        weightRow.push((Math.random() * 2 - 1) * limit);
      }
      layerWeights.push(weightRow);
    }
    weights.push(layerWeights);
    biases.push(new Array(fanOut).fill(0));   // start biases at zero
  }

  return { layerSizes, weights, biases, activationName };
}

// Forward pass: turn an input vector into an output prediction.
// Returns both the activations (per layer) and the pre-activations,
// because we'll need both for backprop.
function forward(network, inputVector) {
  const activationsPerLayer    = [inputVector.slice()];
  const preActivationsPerLayer = [];
  const activation = ACTIVATIONS[network.activationName];

  for (let layer = 0; layer < network.weights.length; layer++) {
    const layerWeights        = network.weights[layer];
    const layerBiases         = network.biases[layer];
    const previousActivations = activationsPerLayer[activationsPerLayer.length - 1];

    const preActivations    = [];
    const currentActivations = [];
    const isOutputLayer = (layer === network.weights.length - 1);

    for (let neuron = 0; neuron < layerWeights.length; neuron++) {
      // Weighted sum + bias = pre-activation `z`.
      let z = layerBiases[neuron];
      const weightRow = layerWeights[neuron];
      for (let input = 0; input < weightRow.length; input++) {
        z += weightRow[input] * previousActivations[input];
      }
      preActivations.push(z);

      // Output layer stays as raw logits, we apply softmax separately
      // when we need probabilities. Hidden layers use the activation fn.
      currentActivations.push(isOutputLayer ? z : activation.fn(z));
    }

    preActivationsPerLayer.push(preActivations);
    activationsPerLayer.push(currentActivations);
  }

  return { activationsPerLayer, preActivationsPerLayer };
}

// Softmax: turn the raw output logits into a probability distribution.
// Subtracting the max first is a numerical-stability trick, same answer,
// no overflow.
function softmax(logits) {
  const maxLogit = Math.max(...logits);
  const exponentials = logits.map(z => Math.exp(z - maxLogit));
  const total = exponentials.reduce((sum, value) => sum + value, 0);
  return exponentials.map(value => value / total);
}

// Backpropagation for softmax + cross-entropy loss.
// Returns the gradient of the loss with respect to every weight and
// bias in the network, plus the loss itself.
function backprop(network, inputVector, trueClassIndex) {
  const { activationsPerLayer } = forward(network, inputVector);
  const numLayers = network.weights.length;
  const activation = ACTIVATIONS[network.activationName];

  // Allocate zero-initialized gradient arrays with the same shape
  // as the weights and biases.
  const weightGradients = network.weights.map(layer =>
    layer.map(row => row.map(() => 0))
  );
  const biasGradients = network.biases.map(layer => layer.map(() => 0));

  // BP1: at the output, the error simplifies (for softmax+CE) to
  //      (predicted probability) − (one-hot true label).
  const outputLogits     = activationsPerLayer[numLayers];
  const outputProbabilities = softmax(outputLogits);
  let errorAtCurrentLayer = outputProbabilities.slice();
  errorAtCurrentLayer[trueClassIndex] -= 1;

  // Walk backward through the layers, computing gradients as we go.
  for (let layer = numLayers - 1; layer >= 0; layer--) {
    const previousLayerActivations = activationsPerLayer[layer];

    // BP3 + BP4: gradient of the loss wrt this layer's biases and weights.
    for (let neuron = 0; neuron < network.weights[layer].length; neuron++) {
      biasGradients[layer][neuron] = errorAtCurrentLayer[neuron];
      for (let input = 0; input < previousLayerActivations.length; input++) {
        weightGradients[layer][neuron][input] =
          errorAtCurrentLayer[neuron] * previousLayerActivations[input];
      }
    }

    // BP2: propagate the error backward to the previous layer,
    // scaled by the activation function's derivative.
    if (layer > 0) {
      const errorAtPreviousLayer = new Array(network.weights[layer - 1].length).fill(0);
      for (let source = 0; source < errorAtPreviousLayer.length; source++) {
        let weightedError = 0;
        for (let target = 0; target < network.weights[layer].length; target++) {
          weightedError += network.weights[layer][target][source] * errorAtCurrentLayer[target];
        }
        errorAtPreviousLayer[source] =
          weightedError * activation.derivativeFromActivation(activationsPerLayer[layer][source]);
      }
      errorAtCurrentLayer = errorAtPreviousLayer;
    }
  }

  // Cross-entropy loss for this example.
  const loss = -Math.log(Math.max(1e-9, outputProbabilities[trueClassIndex]));
  return { weightGradients, biasGradients, loss, outputProbabilities };
}

// One SGD step: subtract a small multiple of each gradient.
function applyGradients(network, gradients, learningRate) {
  for (let layer = 0; layer < network.weights.length; layer++) {
    for (let neuron = 0; neuron < network.weights[layer].length; neuron++) {
      network.biases[layer][neuron] -= learningRate * gradients.biasGradients[layer][neuron];
      for (let input = 0; input < network.weights[layer][neuron].length; input++) {
        network.weights[layer][neuron][input] -=
          learningRate * gradients.weightGradients[layer][neuron][input];
      }
    }
  }
}

# ───────────────────────────────────────────────────────────────
# network.py · feed-forward neural net trained with backprop + SGD
#
# Pure stdlib, no NumPy. For real work you would use NumPy or PyTorch,
# but the loops here make every operation explicit and easy to follow.
# ───────────────────────────────────────────────────────────────
import math
import random

# Activation functions paired with their derivatives. Each derivative
# is expressed in terms of the post-activation `a`, not the pre-activation `z`,
# so we can compute it cheaply from values we already have.
def relu(z):                                   return z if z > 0 else 0.0
def relu_derivative(a):                        return 1.0 if a > 0 else 0.0
def tanh(z):                                   return math.tanh(z)
def tanh_derivative(a):                        return 1.0 - a * a
def sigmoid(z):                                return 1.0 / (1.0 + math.exp(-z))
def sigmoid_derivative(a):                     return a * (1.0 - a)

ACTIVATIONS = {
    'relu':    (relu,    relu_derivative),
    'tanh':    (tanh,    tanh_derivative),
    'sigmoid': (sigmoid, sigmoid_derivative),
}

def make_network(layer_sizes, activation_name='relu'):
    """Build a network. layer_sizes = [n_in, n_hidden, ..., n_out]."""
    weights = []   # weights[layer][neuron][input]
    biases  = []   # biases[layer][neuron]

    for layer in range(1, len(layer_sizes)):
        fan_in  = layer_sizes[layer - 1]
        fan_out = layer_sizes[layer]

        # Target weight variance: 2/fan_in for ReLU (He init), 1/fan_in otherwise
        # (Xavier-style). A uniform draw on [-limit, limit] has variance limit²/3,
        # so limit = sqrt(3 × target_variance).
        target_variance = 2 / fan_in if activation_name == 'relu' else 1 / fan_in
        limit = math.sqrt(3 * target_variance)

        layer_weights = [
            [(random.random() * 2 - 1) * limit for _ in range(fan_in)]
            for _ in range(fan_out)
        ]
        weights.append(layer_weights)
        biases.append([0.0] * fan_out)   # biases start at zero

    return {
        'layer_sizes':     layer_sizes,
        'weights':         weights,
        'biases':          biases,
        'activation_name': activation_name,
    }

def forward(network, input_vector):
    """Run one example through the network. Returns activations and pre-activations
    per layer, we need both for backprop."""
    activations_per_layer     = [list(input_vector)]
    pre_activations_per_layer = []
    activation_fn, _ = ACTIVATIONS[network['activation_name']]

    for layer in range(len(network['weights'])):
        layer_weights        = network['weights'][layer]
        layer_biases         = network['biases'][layer]
        previous_activations = activations_per_layer[-1]
        is_output_layer      = (layer == len(network['weights']) - 1)

        pre_activations = []
        new_activations = []
        for neuron_index in range(len(layer_weights)):
            # Weighted sum + bias = pre-activation z.
            z = layer_biases[neuron_index]
            for input_index in range(len(layer_weights[neuron_index])):
                z += layer_weights[neuron_index][input_index] * previous_activations[input_index]
            pre_activations.append(z)

            # Output layer keeps raw logits; hidden layers apply the activation fn.
            new_activations.append(z if is_output_layer else activation_fn(z))

        pre_activations_per_layer.append(pre_activations)
        activations_per_layer.append(new_activations)

    return activations_per_layer, pre_activations_per_layer

def softmax(logits):
    """Turn raw scores into probabilities. Subtracting the max is a stability trick."""
    max_logit    = max(logits)
    exponentials = [math.exp(z - max_logit) for z in logits]
    total        = sum(exponentials)
    return [value / total for value in exponentials]

def backprop(network, input_vector, true_class_index):
    """Compute gradients of the cross-entropy loss with respect to all weights and biases."""
    activations_per_layer, _ = forward(network, input_vector)
    num_layers = len(network['weights'])
    _, activation_derivative = ACTIVATIONS[network['activation_name']]

    # Allocate zero gradients with the same shape as the parameters.
    weight_gradients = [
        [[0.0] * len(row) for row in layer]
        for layer in network['weights']
    ]
    bias_gradients = [[0.0] * len(layer) for layer in network['biases']]

    # BP1: output error for softmax + cross-entropy is just (probs − one_hot).
    output_logits        = activations_per_layer[num_layers]
    output_probabilities = softmax(output_logits)
    error_at_current_layer = list(output_probabilities)
    error_at_current_layer[true_class_index] -= 1.0

    for layer in range(num_layers - 1, -1, -1):
        previous_layer_activations = activations_per_layer[layer]

        # BP3 + BP4: bias and weight gradients for this layer.
        for neuron in range(len(network['weights'][layer])):
            bias_gradients[layer][neuron] = error_at_current_layer[neuron]
            for input_idx in range(len(previous_layer_activations)):
                weight_gradients[layer][neuron][input_idx] = (
                    error_at_current_layer[neuron] * previous_layer_activations[input_idx]
                )

        # BP2: send the error one layer further back, scaled by the activation slope.
        if layer > 0:
            error_at_previous_layer = [0.0] * len(network['weights'][layer - 1])
            for source in range(len(error_at_previous_layer)):
                weighted_error = 0.0
                for target in range(len(network['weights'][layer])):
                    weighted_error += network['weights'][layer][target][source] * error_at_current_layer[target]
                error_at_previous_layer[source] = (
                    weighted_error
                    * activation_derivative(activations_per_layer[layer][source])
                )
            error_at_current_layer = error_at_previous_layer

    loss = -math.log(max(1e-9, output_probabilities[true_class_index]))
    return weight_gradients, bias_gradients, loss, output_probabilities

def apply_gradients(network, weight_gradients, bias_gradients, learning_rate):
    """One SGD step: subtract a small multiple of each gradient."""
    for layer in range(len(network['weights'])):
        for neuron in range(len(network['weights'][layer])):
            network['biases'][layer][neuron] -= learning_rate * bias_gradients[layer][neuron]
            for input_idx in range(len(network['weights'][layer][neuron])):
                network['weights'][layer][neuron][input_idx] -= (
                    learning_rate * weight_gradients[layer][neuron][input_idx]
                )

// ─────────────────────────────────────────────────────────────────
// network.cpp · feed-forward neural net trained with backprop + SGD
// Requires C++11 or later (uses <random>, range-for, auto, etc.).
// Compile: g++ -std=c++11 -O2 network.cpp -o network
// ─────────────────────────────────────────────────────────────────
#include <vector>
#include <cmath>
#include <random>
#include <algorithm>
#include <string>
#include <cstddef>

// Pair an activation function with its derivative-from-activation.
struct Activation {
    double (*apply)(double z);
    double (*derivative_from_activation)(double a);
};

static double relu_fn(double z) { return z > 0.0 ? z : 0.0; }
static double relu_d (double a) { return a > 0.0 ? 1.0 : 0.0; }
static double tanh_fn(double z) { return std::tanh(z); }
static double tanh_d (double a) { return 1.0 - a * a; }
static double sigmoid_fn(double z) { return 1.0 / (1.0 + std::exp(-z)); }
static double sigmoid_d (double a) { return a * (1.0 - a); }

static Activation choose_activation(const std::string& name) {
    if (name == "relu")    return { relu_fn,    relu_d };
    if (name == "tanh")    return { tanh_fn,    tanh_d };
    return { sigmoid_fn, sigmoid_d };
}

struct Network {
    std::vector<int> layer_sizes;
    std::vector<std::vector<std::vector<double>>> weights;  // weights[layer][neuron][input]
    std::vector<std::vector<double>> biases;                 // biases[layer][neuron]
    std::string activation_name;
};

struct ForwardResult {
    std::vector<std::vector<double>> activations;       // includes input as layer 0
    std::vector<std::vector<double>> pre_activations;
};

struct Gradients {
    std::vector<std::vector<std::vector<double>>> weight_gradients;
    std::vector<std::vector<double>> bias_gradients;
    double loss;
    std::vector<double> output_probabilities;
};

Network make_network(const std::vector<int>& layer_sizes,
                     const std::string& activation_name = "relu") {
    Network net;
    net.layer_sizes     = layer_sizes;
    net.activation_name = activation_name;

    std::mt19937 rng(std::random_device{}());
    std::uniform_real_distribution<double> uniform(-1.0, 1.0);

    for (std::size_t layer = 1; layer < layer_sizes.size(); ++layer) {
        int fan_in  = layer_sizes[layer - 1];
        int fan_out = layer_sizes[layer];

        // Target weight variance: 2/fan_in for ReLU (He init), 1/fan_in otherwise
        // (Xavier-style). A uniform draw on [-limit, limit] has variance limit²/3,
        // so limit = sqrt(3 × target_variance).
        double target_variance = (activation_name == "relu") ? 2.0 / fan_in : 1.0 / fan_in;
        double limit = std::sqrt(3.0 * target_variance);

        std::vector<std::vector<double>> layer_weights(fan_out, std::vector<double>(fan_in));
        for (int neuron = 0; neuron < fan_out; ++neuron)
            for (int input = 0; input < fan_in; ++input)
                layer_weights[neuron][input] = uniform(rng) * limit;

        net.weights.push_back(layer_weights);
        net.biases.push_back(std::vector<double>(fan_out, 0.0));
    }
    return net;
}

ForwardResult forward(const Network& net, const std::vector<double>& input_vector) {
    ForwardResult result;
    result.activations.push_back(input_vector);
    Activation activation = choose_activation(net.activation_name);

    for (std::size_t layer = 0; layer < net.weights.size(); ++layer) {
        const auto& layer_weights = net.weights[layer];
        const auto& layer_biases  = net.biases[layer];
        const auto& previous       = result.activations.back();
        const bool is_output_layer = (layer == net.weights.size() - 1);

        std::vector<double> pre_activations(layer_weights.size());
        std::vector<double> activations    (layer_weights.size());

        for (std::size_t neuron = 0; neuron < layer_weights.size(); ++neuron) {
            double z = layer_biases[neuron];
            for (std::size_t input = 0; input < layer_weights[neuron].size(); ++input)
                z += layer_weights[neuron][input] * previous[input];
            pre_activations[neuron] = z;
            activations[neuron]     = is_output_layer ? z : activation.apply(z);
        }

        result.pre_activations.push_back(pre_activations);
        result.activations.push_back(activations);
    }
    return result;
}

std::vector<double> softmax(const std::vector<double>& logits) {
    double max_logit = *std::max_element(logits.begin(), logits.end());
    std::vector<double> exponentials(logits.size());
    double total = 0.0;
    for (std::size_t i = 0; i < logits.size(); ++i) {
        exponentials[i]  = std::exp(logits[i] - max_logit);
        total           += exponentials[i];
    }
    for (double& value : exponentials) value /= total;
    return exponentials;
}

Gradients backprop(const Network& net,
                    const std::vector<double>& input_vector,
                    int true_class_index) {
    ForwardResult fwd = forward(net, input_vector);
    Activation activation = choose_activation(net.activation_name);
    std::size_t num_layers = net.weights.size();

    // Zero-shaped gradient containers.
    Gradients grads;
    grads.weight_gradients.resize(num_layers);
    grads.bias_gradients.resize(num_layers);
    for (std::size_t layer = 0; layer < num_layers; ++layer) {
        grads.weight_gradients[layer].assign(net.weights[layer].size(), {});
        for (std::size_t neuron = 0; neuron < net.weights[layer].size(); ++neuron)
            grads.weight_gradients[layer][neuron].assign(net.weights[layer][neuron].size(), 0.0);
        grads.bias_gradients[layer].assign(net.biases[layer].size(), 0.0);
    }

    // BP1: output error for softmax + cross-entropy is just (probs - one_hot).
    std::vector<double> output_probabilities = softmax(fwd.activations.back());
    std::vector<double> error_at_current_layer = output_probabilities;
    error_at_current_layer[true_class_index] -= 1.0;

    for (int layer = static_cast<int>(num_layers) - 1; layer >= 0; --layer) {
        const auto& previous_layer_activations = fwd.activations[layer];

        // BP3 + BP4: bias and weight gradients for this layer.
        for (std::size_t neuron = 0; neuron < net.weights[layer].size(); ++neuron) {
            grads.bias_gradients[layer][neuron] = error_at_current_layer[neuron];
            for (std::size_t input = 0; input < previous_layer_activations.size(); ++input)
                grads.weight_gradients[layer][neuron][input] =
                    error_at_current_layer[neuron] * previous_layer_activations[input];
        }

        // BP2: propagate error one layer further back.
        if (layer > 0) {
            std::vector<double> error_at_previous_layer(net.weights[layer - 1].size(), 0.0);
            for (std::size_t source = 0; source < error_at_previous_layer.size(); ++source) {
                double weighted_error = 0.0;
                for (std::size_t target = 0; target < net.weights[layer].size(); ++target)
                    weighted_error += net.weights[layer][target][source] * error_at_current_layer[target];
                error_at_previous_layer[source] =
                    weighted_error * activation.derivative_from_activation(fwd.activations[layer][source]);
            }
            error_at_current_layer = error_at_previous_layer;
        }
    }

    grads.loss = -std::log(std::max(1e-9, output_probabilities[true_class_index]));
    grads.output_probabilities = output_probabilities;
    return grads;
}

void apply_gradients(Network& net, const Gradients& grads, double learning_rate) {
    for (std::size_t layer = 0; layer < net.weights.size(); ++layer) {
        for (std::size_t neuron = 0; neuron < net.weights[layer].size(); ++neuron) {
            net.biases[layer][neuron] -= learning_rate * grads.bias_gradients[layer][neuron];
            for (std::size_t input = 0; input < net.weights[layer][neuron].size(); ++input)
                net.weights[layer][neuron][input] -= learning_rate * grads.weight_gradients[layer][neuron][input];
        }
    }
}

// ───────────────────────────────────────────────────────────────────
// network.rs · feed-forward neural net trained with backprop + SGD.
// Uses only `std`. For real work, pull in the `rand` crate and replace
// the tiny PRNG at the bottom of this file.
// ───────────────────────────────────────────────────────────────────

use std::f64;

#[derive(Clone, Copy)]
pub enum Activation { ReLU, Tanh, Sigmoid }

impl Activation {
    /// Apply the activation function to a pre-activation `z`.
    pub fn apply(self, z: f64) -> f64 {
        match self {
            Activation::ReLU    => if z > 0.0 { z } else { 0.0 },
            Activation::Tanh    => z.tanh(),
            Activation::Sigmoid => 1.0 / (1.0 + (-z).exp()),
        }
    }

    /// Derivative expressed in terms of the *post-activation* value `a`.
    /// Cheap to compute from values we already have.
    pub fn derivative_from_activation(self, a: f64) -> f64 {
        match self {
            Activation::ReLU    => if a > 0.0 { 1.0 } else { 0.0 },
            Activation::Tanh    => 1.0 - a * a,
            Activation::Sigmoid => a * (1.0 - a),
        }
    }
}

pub struct Network {
    pub layer_sizes: Vec<usize>,
    pub weights:     Vec<Vec<Vec<f64>>>,  // [layer][neuron][input]
    pub biases:      Vec<Vec<f64>>,
    pub activation:  Activation,
}

pub struct Gradients {
    pub weight_gradients:     Vec<Vec<Vec<f64>>>,
    pub bias_gradients:       Vec<Vec<f64>>,
    pub loss:                 f64,
    pub output_probabilities: Vec<f64>,
}

impl Network {
    /// Build a fresh network with the given shape.
    pub fn new(layer_sizes: &[usize], activation: Activation) -> Self {
        let mut weights: Vec<Vec<Vec<f64>>> = Vec::new();
        let mut biases:  Vec<Vec<f64>>      = Vec::new();

        for layer in 1..layer_sizes.len() {
            let fan_in  = layer_sizes[layer - 1];
            let fan_out = layer_sizes[layer];

            // Target weight variance: 2/fan_in for ReLU (He init), 1/fan_in otherwise
            // (Xavier-style). A uniform draw on [-limit, limit] has variance limit²/3,
            // so limit = sqrt(3 × target_variance).
            let target_variance = match activation {
                Activation::ReLU => 2.0 / fan_in as f64,
                _                  => 1.0 / fan_in as f64,
            };
            let limit = (3.0 * target_variance).sqrt();

            let layer_weights: Vec<Vec<f64>> = (0..fan_out)
                .map(|_| (0..fan_in)
                    .map(|_| (rand_uniform() * 2.0 - 1.0) * limit)
                    .collect())
                .collect();
            weights.push(layer_weights);
            biases.push(vec![0.0; fan_out]);
        }

        Self { layer_sizes: layer_sizes.to_vec(), weights, biases, activation }
    }

    /// Run one input vector through the network. Returns (activations, pre-activations)
    /// per layer; the input is included as `activations[0]`.
    pub fn forward(&self, input_vector: &[f64]) -> (Vec<Vec<f64>>, Vec<Vec<f64>>) {
        let mut activations_per_layer: Vec<Vec<f64>> = vec![input_vector.to_vec()];
        let mut pre_activations_per_layer: Vec<Vec<f64>> = Vec::new();

        for (layer, layer_weights) in self.weights.iter().enumerate() {
            let layer_biases = &self.biases[layer];
            let previous     = activations_per_layer.last().unwrap().clone();
            let is_output_layer = layer == self.weights.len() - 1;

            let mut pre_activations = Vec::with_capacity(layer_weights.len());
            let mut activations     = Vec::with_capacity(layer_weights.len());

            for (neuron, weight_row) in layer_weights.iter().enumerate() {
                let mut z = layer_biases[neuron];
                for input in 0..weight_row.len() {
                    z += weight_row[input] * previous[input];
                }
                pre_activations.push(z);
                activations.push(if is_output_layer { z } else { self.activation.apply(z) });
            }

            pre_activations_per_layer.push(pre_activations);
            activations_per_layer.push(activations);
        }

        (activations_per_layer, pre_activations_per_layer)
    }

    /// Backpropagation for cross-entropy + softmax.
    pub fn backprop(&self, input_vector: &[f64], true_class: usize) -> Gradients {
        let (activations_per_layer, _) = self.forward(input_vector);

        let mut weight_gradients: Vec<Vec<Vec<f64>>> = self.weights.iter()
            .map(|layer| layer.iter().map(|row| vec![0.0; row.len()]).collect())
            .collect();
        let mut bias_gradients: Vec<Vec<f64>> = self.biases.iter()
            .map(|b| vec![0.0; b.len()])
            .collect();

        // BP1: softmax(logits) − one_hot(true_class).
        let output_probabilities = softmax(activations_per_layer.last().unwrap());
        let mut error_at_current_layer: Vec<f64> = output_probabilities.clone();
        error_at_current_layer[true_class] -= 1.0;

        for layer in (0..self.weights.len()).rev() {
            let previous_layer_activations = &activations_per_layer[layer];

            // BP3 + BP4: bias and weight gradients for this layer.
            for neuron in 0..self.weights[layer].len() {
                bias_gradients[layer][neuron] = error_at_current_layer[neuron];
                for input in 0..previous_layer_activations.len() {
                    weight_gradients[layer][neuron][input] =
                        error_at_current_layer[neuron] * previous_layer_activations[input];
                }
            }

            // BP2: send error one layer further back.
            if layer > 0 {
                let mut error_at_previous_layer = vec![0.0; self.weights[layer - 1].len()];
                for source in 0..error_at_previous_layer.len() {
                    let mut weighted_error = 0.0;
                    for target in 0..self.weights[layer].len() {
                        weighted_error += self.weights[layer][target][source] * error_at_current_layer[target];
                    }
                    error_at_previous_layer[source] = weighted_error
                        * self.activation.derivative_from_activation(activations_per_layer[layer][source]);
                }
                error_at_current_layer = error_at_previous_layer;
            }
        }

        Gradients {
            weight_gradients,
            bias_gradients,
            loss: -output_probabilities[true_class].max(1e-9).ln(),
            output_probabilities,
        }
    }

    /// One SGD step: subtract a small multiple of each gradient.
    pub fn apply_gradients(&mut self, grads: &Gradients, learning_rate: f64) {
        for layer in 0..self.weights.len() {
            for neuron in 0..self.weights[layer].len() {
                self.biases[layer][neuron] -= learning_rate * grads.bias_gradients[layer][neuron];
                for input in 0..self.weights[layer][neuron].len() {
                    self.weights[layer][neuron][input] -=
                        learning_rate * grads.weight_gradients[layer][neuron][input];
                }
            }
        }
    }
}

fn softmax(logits: &[f64]) -> Vec<f64> {
    let max_logit = logits.iter().copied().fold(f64::NEG_INFINITY, f64::max);
    let exponentials: Vec<f64> = logits.iter().map(|&z| (z - max_logit).exp()).collect();
    let total: f64 = exponentials.iter().sum();
    exponentials.iter().map(|e| e / total).collect()
}

// Tiny xorshift64 PRNG. Single-threaded only: `static mut` is unsynchronized,
// so concurrent calls would be undefined behavior. For real (multi-threaded
// or production) code, drop this and use the `rand` crate.
fn rand_uniform() -> f64 {
    use std::time::{SystemTime, UNIX_EPOCH};
    static mut STATE: u64 = 0;
    unsafe {
        if STATE == 0 {
            STATE = SystemTime::now()
                .duration_since(UNIX_EPOCH).unwrap()
                .subsec_nanos() as u64 | 1;
        }
        STATE ^= STATE << 13;
        STATE ^= STATE >> 7;
        STATE ^= STATE << 17;
        (STATE as f64) / (u64::MAX as f64)
    }
}

Every line traces back to an equation from earlier sections. The same algorithm, scaled up with matrix libraries, GPUs, attention layers and billions of parameters, is how large language models are trained.

13Train it yourself

The widget below runs a compact copy of the JavaScript above (same algorithm, same initialization, shorter variable names) in your browser. Pick a dataset, tweak the architecture, hit TRAIN. The left panel shows the predicted class probability over the input space, with the white line marking the decision boundary where both classes are equally likely; the chart under the controls plots the average training loss per epoch.

Things to try

• Push the learning rate up to 0.2 on the spiral. The loss curve spikes instead of settling, and the accuracy usually ends well below what η = 0.03 reaches. With ReLU, the oversized steps also tend to push several hidden units permanently negative: the dying-ReLU failure from §10.
• Drop the learning rate to 0.001. The loss falls smoothly, with few upticks, but takes several hundred epochs to get where η = 0.03 gets in about fifty.
• Use only 2 hidden units on the spiral. The network doesn't have the capacity to bend the boundary enough and stalls around 70 to 75% accuracy. At 8 or 16 units it fits the curve.
• Switch from ReLU to sigmoid on the spiral. The loss sits on a plateau for much longer before it starts to fall. With only two hidden layers this is a mild version of the vanishing-gradient problem from §10: every σ′ factor is at most 0.25. Tanh, whose slope reaches 1, trains much closer to ReLU's pace.
• Try XOR with 2 hidden units and sigmoid. The loss usually sits near 0.69 (chance level for two classes) for dozens of epochs, sometimes over a hundred, then drops as the hidden units find the two cuts XOR needs. A single perceptron, with no hidden layer, can never do this: Minsky and Papert's point. With ReLU and only 2 units, runs often stall with dead units.

14Overfitting and how to spot it

There's a problem hiding in everything we've done so far. A network with enough capacity can memorize its training set, noise included. Give 100 noisy points to a network with 1,000 weights and, trained long enough, it will bend its decision boundary around every one of them. On the training data it looks perfect; on new data it does worse than a simpler fit would.

This is overfitting: confusing the noise in the training set for signal in the world. The standard practical fix is to hold out some of your data, never let the network train on it, and measure performance separately on the held-out set during training. The gap between the two numbers tells you whether the model is generalizing.

The widget below trains an over-parameterized network (16 + 16 hidden units, about 350 weights and biases) on 40 noisy points, jittered enough that a few land on the wrong side of the gap between the two classes. The solid green curve is loss on the training set. The dashed red curve is loss on a held-out validation set the network never sees. Validation loss usually bottoms out within the first 50 epochs; after that, training loss keeps falling while validation loss climbs.

What overfitting looks like, in three pictures

The diagnostic pattern is:

• Training loss keeps decreasing. The network is finding patterns in the training data.
• Validation loss decreases for a while, then turns around and starts increasing. Past this point, every additional epoch makes the network worse at the actual problem.
• The decision boundary develops weird bulges around individual training points, trying to capture each one's noisy position.

How to fix it

Four standard cures, often combined:

1. More data. The cleanest fix. With a thousand points, the noise averages out and a network this size no longer has the capacity to memorize individual examples. Often not an option in practice.
2. Smaller network. Fewer parameters means less capacity to memorize.
3. Weight decay (L2 regularization). Penalize the network for having large weights, which pushes it toward smoother functions. Train with the weight-decay slider at 0.01 or above and the validation curve climbs far less than it does at 0.
4. Early stopping. Track validation loss and keep a copy of the weights from its best epoch; once it has stopped improving for a while, stop and use that copy.

There are more techniques: dropout, which randomly zeros out neurons during training; batch normalization; and data augmentation. The four above are the baseline most training setups start from.

15Your turn: the code playground

The playground below runs the same compact JavaScript library as the trainer, exposed as the MPGNet object, so you can build your own network, train it, and predict from the editor. The Python, C++, and Rust ports from §12 do the same thing; copy any of them into your own editor to run there. Hit RUN (or Ctrl+Enter / Cmd+Enter) to execute; output prints on the right.

Edit anything: the data, the architecture, the learning rate, a held-out test set. The library follows the §12 JavaScript with shorter names: makeNetwork(layerSizes, { activation }) returns { W, B } (weights and biases), backprop(net, input, classIndex) returns { dW, dB, loss, probs }, and forward, predict, applyGrads and softmaxVec round it out. It also lives at window.MPGNet for use from the browser console.

16Where to go from here

What we just built is called a multilayer perceptron (MLP) or feedforward neural network. It is the base that later architectures build on, and fully connected MLP blocks still sit inside most of them. On its own it handles small tabular datasets and 2D toys, and struggles with the images, audio and text that made deep learning famous.

Images: convolutions

A 256×256 grayscale image has 65,536 pixels. Connecting every pixel to each of, say, 500 neurons in an MLP's first hidden layer takes about 33 million weights for that layer alone, and the network has to learn the same edge detector separately at every position. Convolutional neural networks (CNNs) exploit the image's structure: a small filter of shared weights slides over the image, producing a feature map that records where its pattern appears. The same algorithm we built (backprop on a chain of differentiable operations) still drives training; only the layer types change.

Sequences: attention

For text, audio, or any data that comes in a sequence, an MLP needs a fixed input size and has no built-in notion of order. The standard sequence models before 2017, recurrent networks, read one token at a time, so relating token i to token i+50 means carrying information through every position in between. Vaswani et al.'s 2017 Attention Is All You Need^[9] replaces recurrence with scaled dot-product attention, which connects every pair of positions directly:

Every token in the sequence has three roles: it asks (Q), it advertises (K), and it has content (V). Compute the similarity between every query and every key, scale down, softmax to get attention weights, and use them to mix the values. The result: each token gets a customized summary of every other token.

Stack blocks of attention followed by a small MLP (with residual connections and normalization around each) and you have a transformer, the architecture behind nearly every major language model since 2018. Karpathy's nanoGPT^[10] is a compact reference implementation of a GPT-style transformer: about 300 lines for the model, 300 for the training loop.

You don't need GPT-4 to do interesting things

Eldan and Li's TinyStories paper^[11] showed that transformers with under 10 million parameters, some with a single transformer block, can write fluent, consistent short stories when trained on a synthetic dataset restricted to words a typical 3- to 4-year-old understands. The authors report that every model in the paper trains on a single V100 GPU in at most 30 hours, so the whole pipeline is within reach of a hobby budget.

Building a tiny transformer takes the same forward/backward/update loop with different layers in the middle. Karpathy's series^[15] and Perivolaropoulos's Train Your Own LLM From Scratch workshop^[19] both walk through one step by step.

The final exam

Five questions covering the whole tutorial. If you can answer all five without scrolling back up, you've got the fundamentals.

17Sources & further reading

Numbered citations refer to the superscripts above. Links are given where a stable free copy exists.

1. Rosenblatt, F. (1958). The perceptron: a probabilistic model for information storage and organization in the brain. Psychological Review, 65(6), 386–408. The original artificial neuron.
2. Rumelhart, D. E., Hinton, G. E., & Williams, R. J. (1986). Learning representations by back-propagating errors. Nature, 323, 533–536. Brought backpropagation into the mainstream.
3. Minsky, M. & Papert, S. (1969). Perceptrons. MIT Press. Proves the limits of single-layer perceptrons, XOR among them; widely cited as a cause of the 1970s decline in neural-network research.
4. Nielsen, M. A. (2015). Neural Networks and Deep Learning. Determination Press. Free online: neuralnetworksanddeeplearning.com. Chapter 2 derives the four backprop equations (BP1–BP4) used in §9; chapter 1 uses the quadratic cost of eq. 6.
5. Hochreiter, S. (1991). Untersuchungen zu dynamischen neuronalen Netzen. Diploma thesis, TU Munich. The first formal analysis of the vanishing-gradient problem.
6. Radford, A., Wu, J., Child, R., Luan, D., Amodei, D., & Sutskever, I. (2019). Language Models are Unsupervised Multitask Learners. OpenAI. PDF. The GPT-2 paper; the released model code (src/model.py) applies GELU in every feed-forward block.
7. Glorot, X. & Bengio, Y. (2010). Understanding the difficulty of training deep feedforward neural networks. AISTATS. The Xavier initialization paper: derives the forward (1/n_in) and backward (1/n_out) variance conditions and proposes 2/(n_in + n_out).
8. He, K., Zhang, X., Ren, S., & Sun, J. (2015). Delving Deep into Rectifiers: Surpassing Human-Level Performance on ImageNet Classification. ICCV. The He / Kaiming initialization paper: variance 2/n_in for ReLU layers, and the 30-layer comparison against Xavier initialization.
9. Vaswani, A., Shazeer, N., Parmar, N. et al. (2017). Attention Is All You Need. NeurIPS. arXiv:1706.03762. The Transformer paper.
10. Karpathy, A. nanoGPT. GitHub: github.com/karpathy/nanoGPT. A compact GPT implementation; its README describes a ~300-line model and a ~300-line training loop.
11. Eldan, R. & Li, Y. (2023). TinyStories: How Small Can Language Models Be and Still Speak Coherent English? arXiv:2305.07759. Coherent story generation below 10 million parameters; all models trained on a single V100 GPU in at most 30 hours.
12. Hinton, G. (2012). Lecture 6e: rmsprop, Divide the gradient by a running average of its recent magnitude. Slide deck from the Coursera Neural Networks course. RMSProp's first public appearance.
13. Kingma, D. P. & Ba, J. (2014). Adam: A Method for Stochastic Optimization. arXiv:1412.6980. Source of the Adam update and its bias correction; Section 3 of the paper notes that without the correction the initial steps are much larger.
14. LeCun, Y., Boser, B., Denker, J. S., Henderson, D., Howard, R. E., Hubbard, W., & Jackel, L. D. (1989). Backpropagation applied to handwritten zip code recognition. Neural Computation, 1(4), 541–551. An early convolutional network trained end to end with backpropagation, on U.S. Postal Service ZIP-code digits; the MNIST work in Gradient-based learning applied to document recognition (LeCun et al., 1998) followed nine years later.
15. Karpathy, A. Neural Networks: Zero to Hero. YouTube series, 2022–2023. karpathy.ai/zero-to-hero.html. Builds from a scalar autograd engine up to a small GPT with the same forward, loss, backward, update loop.
16. Linnainmaa, S. (1970). The representation of the cumulative rounding error of an algorithm as a Taylor expansion of the local rounding errors. M.Sc. thesis, University of Helsinki. Reverse-mode automatic differentiation, which is the algorithm we now call backpropagation, sixteen years before Rumelhart-Hinton-Williams.
17. Cybenko, G. (1989). Approximation by superpositions of a sigmoidal function. Mathematics of Control, Signals, and Systems, 2(4), 303–314. The universal approximation theorem for one hidden layer of sigmoid units.
18. Hornik, K. (1991). Approximation capabilities of multilayer feedforward networks. Neural Networks, 4(2), 251–257. Extends the result to any bounded, nonconstant activation function.
19. Perivolaropoulos, A. (2026). Train Your Own LLM From Scratch. GitHub: github.com/angelos-p/llm-from-scratch. A hands-on workshop that builds a ~10M-parameter GPT which trains on a laptop in under an hour.