RNN & GRU
Recurrent Architectures

Chapter 1: The Vanilla RNN Cell

Chapter 0 ended with a scalar notebook. One number is not much of a memory — it can encode "how excited am I" and nothing else. Promote it to a vector and the cell becomes the real thing: the Elman RNN, named for Jeffrey Elman’s 1990 paper Finding Structure in Time, and still the cell everyone means by "vanilla RNN".

Every symbol, in words, before we touch it:

• x_t — the input at time t, a vector in R^d. For characters it is a one-hot: all zeros except a 1 at the index of the character. d is the vocabulary size.
• h_t — the hidden state after reading x_t, a vector in Rⁿ. n is yours to pick: 128, 256, 512. It is the entire memory budget.
• W_x — the input-to-hidden matrix, shape n × d. It answers "what does this new symbol contribute to memory?"
• W_h — the hidden-to-hidden matrix, shape n × n. It answers "how should the old memory be rewritten before the new symbol lands on it?" This is the matrix that will cause every problem in Chapter 3.
• b — a bias in Rⁿ, a per-unit resting level.
• tanh — applied elementwise, squashing every component into (−1, 1). Without it the whole recurrence collapses into one big matrix product and the network can only represent linear functions of the input sequence.

Two design decisions are already baked in and both matter. First, the same W_x, W_h, b are used at every timestep — this is weight sharing across time, and it is why a 3-token sequence and a 3,000-token sequence cost the same number of parameters. Second, h_t is the only channel from the past to the future. Everything the model will ever know about characters 1 through t must survive inside those n numbers. That is a brutal bottleneck, and later chapters are largely about widening it or defending it.

The output head: state is not prediction

h_t is a summary, not an answer. To get a prediction you add a readout layer — one more matrix, mapping the state to one score per vocabulary item, then a softmax to turn scores into probabilities:

W_y has shape d × n (vocabulary rows, state columns). The vector y_t is called the logits — unnormalised scores, one per possible next character, free to be any real number. Softmax exponentiates them and divides by the sum, producing a proper distribution. Keeping the readout separate from the recurrence is a real engineering choice: it lets you swap vocabularies without retraining the cell, and it lets you tie W_y to the input embedding to save parameters, a trick standard in language models.

Do it by hand: three steps, every number

Vocabulary {a, b, c}, so d = 3. State size n = 2. Concrete weights, chosen to be readable rather than trained:

Input sequence: a, b, c. Because each x is one-hot, W_xx_t is simply a column of W_x — column 1 for a, column 2 for b, column 3 for c. That is worth internalising: for one-hot inputs the input matrix is a lookup table, which is exactly why real implementations use an embedding lookup instead of a matrix multiply.

Step 1 — read a, so W_xx₁ = column 1 = [1.0, 0.0]. The state is still zero, so W_hh₀ = [0, 0].

Step 2 — read b, so W_xx₂ = column 2 = [−0.5, 0.8]. Now the recurrent term is live. Multiply row by row:

Step 3 — read c, so W_xx₃ = column 3 = [0.2, −0.3].

Five decimal places through the intermediates is not fussiness: round each step to four and the third state drifts by one in the last digit, so the printed library output would no longer match. Rounding error compounds along a recurrence exactly the way gradients do — a small preview of Chapter 3.

Three observations from those nine multiplications. Unit 1 went 0.80 → 0.02 → 0.18: it is being driven mostly by the input, since W_x row 1 has the large entries. Unit 2 went −0.10 → 0.62 → −0.03: it holds on longer, because W_h[2][2] = 0.6 is the biggest recurrent weight and feeds unit 2 back into itself. And every pre-activation landed inside roughly (−1.2, 1.2), where tanh is steep and informative. Chapter 3 is about what happens when they do not.

See it: step the cell, watch the state

The widget runs exactly the arithmetic above. Press Step to advance one character; the diagram shows the input column, the recurrent contribution, the pre-activation, and the resulting state as signed bars. The recurrent-weight slider scales all of W_h at once: turn it to 0 and the cell becomes memoryless (each state depends only on the current character); turn it up past 1.5 and watch the pre-activations saturate against the flat tails of tanh, which is a preview of Chapter 3’s disaster.

From scratch, then the one-liner

Here is the cell in NumPy with the exact weights from the hand calculation, so you can check the printout against the numbers above. This is the whole forward pass — there is nothing hidden.

python — Elman RNN forward, from scratch
import numpy as np

Wx = np.array([[ 1.0, -0.5,  0.2],
               [ 0.0,  0.8, -0.3]])      # (n=2, d=3)
Wh = np.array([[ 0.5, -0.2],
               [ 0.1,  0.6]])              # (n=2, n=2)
b  = np.array([ 0.1, -0.1])                # (n=2,)

def rnn_step(x, h):
    # x: (d,) one-hot   h: (n,) previous state  ->  (n,) new state
    return np.tanh(Wx @ x + Wh @ h + b)

def rnn_forward(idxs, d=3, n=2):
    h  = np.zeros(n)
    hs = []
    for i in idxs:
        x = np.zeros(d); x[i] = 1.0       # one-hot for this character
        h = rnn_step(x, h)
        hs.append(h.copy())                # keep EVERY state — backprop needs them
    return np.array(hs)                    # (T, n)

print(np.round(rnn_forward([0, 1, 2]), 4))
# [[ 0.8005 -0.0997]
#  [ 0.0202  0.6171]
#  [ 0.1845 -0.0277]]   <-- matches the hand calculation exactly

That hs.append(h.copy()) is not bookkeeping for the demo — it is mandatory. The backward pass in Chapter 2 needs h_t at every t to compute 1 − h_t², so the forward pass must retain the whole state trajectory. Memory during training therefore scales as O(T × n), which is precisely why truncated backpropagation exists.

And the library equivalent. PyTorch calls the two matrices weight_ih_l0 and weight_hh_l0, and carries two biases instead of one (mathematically redundant, kept for cuDNN compatibility):

python — the same cell via torch.nn.RNN
import torch, torch.nn as nn

rnn = nn.RNN(input_size=3, hidden_size=2, batch_first=True)
with torch.no_grad():                        # paste our weights in
    rnn.weight_ih_l0.copy_(torch.tensor(Wx))
    rnn.weight_hh_l0.copy_(torch.tensor(Wh))
    rnn.bias_ih_l0.copy_(torch.tensor(b))
    rnn.bias_hh_l0.zero_()                    # fold both biases into one

x = torch.eye(3)[[0,1,2]].unsqueeze(0).double()   # (1, 3, 3)
out, hT = rnn.double()(x)
print(out.squeeze(0).round(decimals=4))
# tensor([[ 0.8005, -0.0997], [ 0.0202,  0.6171], [ 0.1845, -0.0277]])

Same three rows. The library is the loop, compiled. Everything you will meet later — gates, clipping, stacking — is a modification of the six lines above, not a replacement for them.

Chapter 2: Unrolling and BPTT

We have a cell that runs forward. To train it we need gradients, and gradients need a graph. The problem is that the recurrence has a cycle: h feeds into itself. Backpropagation is defined on directed acyclic graphs. So the first move is a bookkeeping trick that makes the cycle disappear.

Unrolling: a loop is a deep network in disguise

Unrolling means writing out one copy of the cell per timestep, laid left to right, with h_t−1 wired into the copy at t. A 3-step sequence becomes a 3-layer feed-forward network. A 100-step sequence becomes a 100-layer one. The copies are not independent — they share W_x, W_h, b, the same physical numbers used three times or a hundred times.

Once unrolled, it is an ordinary DAG and ordinary backpropagation applies. Running backprop over the unrolled graph has a name: backpropagation through time, or BPTT. It is not a new algorithm. It is the chain rule on a graph that happens to be very deep and to reuse the same weights at every depth.

That weight sharing has one consequence you must never forget, because it is the single most common implementation bug in hand-written RNNs. Each parameter receives a gradient contribution from every timestep it participated in, and those contributions add up:

The chain rule, derived step by step

Take the smallest interesting case: a scalar RNN, h_t = tanh(w·x_t + u·h_t−1), three timesteps, and a loss on the final state only, L = ½(h₃ − y)². Scalar so every quantity is a number you can check.

Start at the output and walk backwards. The derivative of the loss with respect to the last state is immediate:

Next, how does h₃ depend on h₂? Write a₃ for the pre-activation w·x₃ + u·h₂. Then h₃ = tanh(a₃), and the derivative of tanh is 1 − tanh², so:

That factor — squash-derivative times recurrent weight — is the atom of the whole subject. Every step backwards multiplies by one of them. Chapter 3 is nothing but the study of what happens when you multiply many of them together.

Now the parameter u. It enters the graph three times, once per step, so its gradient is a sum of three terms. The term from step t is "the gradient arriving at h_t, times the local derivative of h_t with respect to u holding h_t−1 fixed", and that local derivative is (1 − h_t²)·h_t−1:

and the gradient arriving at an earlier state is the gradient at the later one dragged back through the atoms:

Those two lines are the entire BPTT algorithm for a scalar cell. In the vector case the products become matrix products and u becomes W_h^T, but the structure is identical.

Now do it with actual numbers

Set w = 0.8, u = 0.9, inputs x = [1, 1, 1], target y = 0.5, h₀ = 0. Forward first:

The three squash derivatives, which we will reuse constantly:

Notice they are already shrinking as we move forward in time: 0.559, 0.217, 0.151. The state is drifting toward the flat part of tanh, and the flat part of tanh has a small derivative. Backwards:

Now assemble ∂L/∂w, which collects one term per timestep. Since x_t = 1 at every step, the local factor for w at step t is just (1 − h_t²):

And for u, where the local factor at step t is (1 − h_t²)·h_t−1. The t = 1 term vanishes because h₀ = 0 — a small but real fact: the very first step teaches the recurrent weight nothing, because there was no previous state to multiply.

Truncated BPTT: cutting the tail off honestly

A 10,000-character document unrolls into a 10,000-layer network. Storing every state costs O(T·n) memory, and the backward pass costs O(T) sequential steps you cannot parallelise. Nobody does this. Instead, truncated BPTT processes the sequence in chunks of k steps (k = 25 and k = 35 are the classic choices), backpropagates only within a chunk, and carries the value of h across the chunk boundary while cutting the gradient there.

That distinction is the whole trick, and it is worth stating twice. The forward state is continuous — the model’s memory really does span the whole document. Only the credit assignment is truncated: the optimiser never asks "how should I change W to have made the state 300 steps ago better?" In PyTorch this is one line, h = h.detach(), at the chunk boundary.

What does truncation cost? Use the numbers we just computed. Full BPTT over three steps gave a w-gradient sum of 0.19604. Truncate to k = 2 (drop the t = 1 term):

A 7.6% error for a 33% saving. And it gets better, not worse, as k grows: the terms you are dropping are the ones that were already tiny. That is the honest justification for truncation — not that long-range gradients do not matter, but that in a plain RNN they are already negligible, so discarding them costs almost nothing. In a gated cell the same terms are much larger, which is exactly why gated cells benefit from bigger k.

See it: drag the truncation window

Each bar below is the magnitude of one timestep’s contribution to the weight gradient, computed by actually running the recurrence over twelve steps. (Magnitudes, because with an alternating input the signed terms partly cancel and would hide the decay we are here to see.) Drag the window slider to keep only the last k steps and watch how much of the total survives — the readout gives you the exact percentage. Then raise the recurrent weight and watch the far bars grow: the truncation you were getting away with at u = 0.3 starts throwing away real gradient by u = 0.9.

One deliberate difference from the hand calculation: this widget drives the cell with a small alternating input rather than a constant one. That keeps the state near zero, where the tanh derivative is close to 1 and the per-step factor really is about u. With a large constant drive the state saturates near ±0.95, the tanh derivative collapses to 0.1, and the per-step factor stays far below 1 no matter how large you make u — which is itself worth knowing, and is exactly why Chapter 3’s slider has a separate control for the average tanh derivative.

BPTT from scratch, then autograd

The from-scratch version below computes the same 0.0826 and 0.0647. Read the loop backwards: it starts with the gradient at the final state and drags it back, accumulating into dw and du at every step.

python — BPTT by hand, scalar cell
import numpy as np

w, u, y = 0.8, 0.9, 0.5
xs = [1.0, 1.0, 1.0]

# ---- forward, keeping every state ----
hs = [0.0]
for x in xs:
    hs.append(np.tanh(w*x + u*hs[-1]))
loss = 0.5 * (hs[-1] - y)**2

# ---- backward through time ----
dw = du = 0.0
dh = hs[-1] - y                      # dL/dh_T
for t in range(len(xs), 0, -1):        # t = 3, 2, 1
    da  = dh * (1 - hs[t]**2)          # through the tanh
    dw += da * xs[t-1]                 # += NOT = : shared weight
    du += da * hs[t-1]                  # += NOT = : shared weight
    dh  = da * u                          # hand the gradient one step further back

print(round(loss,4), round(dw,4), round(du,4))
# 0.0887 0.0826 0.0647   <-- exactly the hand calculation

Autograd does the same thing, and confirming that the two agree is the standard way to trust a hand-written backward pass. If you ever write your own cell, run this comparison before anything else:

python — the same gradients from torch autograd
import torch

w = torch.tensor(0.8, requires_grad=True)
u = torch.tensor(0.9, requires_grad=True)

h = torch.tensor(0.0)
for _ in range(3):
    h = torch.tanh(w*1.0 + u*h)          # the unrolled graph is built here
(0.5*(h - 0.5)**2).backward()          # BPTT = backward() on that graph
print(w.grad.item(), u.grad.item())   # 0.0826  0.0647

# truncation in real code is exactly one line, at the chunk boundary:
#     h = h.detach()   # keep the VALUE, cut the gradient path

Chapter 3: Vanishing and Exploding

Chapter 2 left a clue in plain sight: the gradient contribution from three steps back was ten times smaller than the one from the last step. Now we find out why, and discover that the answer forbids plain RNNs from doing the one thing we built them for.

The product of Jacobians

To push a gradient from time T back to time 0, backpropagation multiplies together the step-to-step derivatives. For a vector cell, each of those is a matrix — a Jacobian, the matrix of all partial derivatives of one vector with respect to another. Differentiate h_t = tanh(W_xx_t + W_hh_t−1 + b) with respect to h_t−1:

D_t is a diagonal matrix holding the tanh derivative for each unit, and W_h is the same recurrent matrix at every step. Chain them over T steps and the gradient that reaches the beginning of the sequence is:

There it is: the same matrix, raised to (roughly) the T-th power. Powers of a matrix do exactly two things — they blow up or they collapse — and which one depends on a single number.

Deriving the bound

Take norms. The norm of a product is at most the product of the norms, so:

Bound each factor separately. For the diagonal: tanh has derivative 1 − tanh², which is 1 at the origin and shrinks toward 0 in the tails, so ‖D_t‖ ≤ γ = 1 always, and in practice sits well below 1 — a unit sitting at h = 0.9 contributes only 1 − 0.81 = 0.19. For the weight matrix, the relevant norm is the largest singular value, written σ_max — intuitively, the largest factor by which W_h can stretch any vector in a single step. It is often loosely called the spectral radius, though strictly the spectral radius is the largest absolute eigenvalue, and that is the one that governs repeated multiplication; the measurement script at the end of this chapter makes the difference concrete. Put them together:

The exponent is the sequence length. That is the whole result, and it is Pascanu, Mikolov and Bengio’s 2013 analysis in one line: if γ·σ_max < 1 the gradient is guaranteed to vanish geometrically; if it exceeds 1 the gradient is free to explode geometrically. Since γ ≤ 1, the sufficient condition for vanishing is σ_max < 1 — which is where random initialisations usually land, and where weight decay actively pushes you.

What geometric decay actually looks like

"Geometric" hides how violent this is, so put numbers on it. Take a well-behaved network: σ_max = 0.9, and states typically around h ≈ 0.7 so the average tanh derivative is 1 − 0.49 ≈ 0.5. The per-step factor is 0.9 × 0.5 = 0.45.

Check the T = 20 row by hand so it is not magic: ln(0.45) = −0.7985, times 20 is −15.97, and e^−15.97 = 1.16 × 10⁻⁷. A float32 has about seven significant decimal digits, so at T = 20 the long-range gradient is already competing with rounding error, and at T = 50 it is numerically zero. The model does not learn the dependency slowly. It does not learn it at all.

Now the other direction. Set σ_max = 1.5 and suppose the states sit near zero where tanh’ ≈ 1, so the per-step factor is 1.5:

A gradient scaled by 640 million, hit with a learning rate of 0.01, moves every weight by millions. One optimiser step and the weights are inf; the step after that the loss is NaN and the run is dead. This is the classic RNN training experience: loss falling nicely, then a single vertical spike, then nothing.

See it: turn the dial through the phase change

The plot below is a log-scale view of the bound (γ·σ_max)^T as T grows, evaluated at whatever per-step factor your two sliders produce. Slide the spectral radius: below the critical line the curve dives, above it the curve climbs, and right where the product hits 1 it walks flat. Then lower the average tanh derivative — a proxy for "the states are saturated" — and notice that the safe region shrinks. Saturation makes vanishing worse, which is why weight initialisation for RNNs is unusually fussy.

The cheap fix: gradient clipping

Explosion has a blunt, effective, universally used cure, called gradient clipping: a rescaling applied to the gradient after it is computed but before the optimiser uses it. Compute the global gradient norm across every parameter; if it exceeds a threshold c, rescale the entire gradient to have norm exactly c:

Worked example. Suppose a batch produces a gradient with global norm ‖g‖ = 47.3 and you clip at c = 5.0. The scale factor is 5.0 / 47.3 = 0.1057, and every component is multiplied by it: a component that was 12.0 becomes 12.0 × 0.1057 = 1.268; one that was −0.4 becomes −0.0423. Check the result: the new norm is 47.3 × 0.1057 = 5.0, as required.

Crucially this preserves the direction of the gradient and only shortens it. You still step downhill; you just refuse to take a mile-long stride off a cliff. Typical thresholds are c = 1 to 5 for language models. The cost is essentially zero and it is standard practice in every RNN codebase.

Measure it yourself

Do not take the plot on faith. This script builds a random recurrent matrix with a chosen spectral radius and measures the actual product of Jacobians at every distance. Run it, watch the numbers fall off a cliff — and then look carefully at the ρ = 1.5 row, which does something the bound alone would not have told you.

python — measure the decay directly
import numpy as np

def make_Wh(n, rho, seed=0):
    # Scale so the SPECTRAL RADIUS (largest |eigenvalue|) is rho. Repeated
    # products are governed by the spectral radius; sigma_max bounds only ONE step.
    rng = np.random.default_rng(seed)
    W = rng.standard_normal((n, n))
    ev = np.max(np.abs(np.linalg.eigvals(W)))
    return W * (rho / ev)

def log10_grad_norm(rho, T=60, n=64):
    Wh = make_Wh(n, rho)
    rng = np.random.default_rng(1)
    h = np.zeros(n); J = np.eye(n); acc = 0.0; out = []
    for t in range(T):
        x = rng.standard_normal(n) * 0.5              # stand-in for Wx @ x_t
        h = np.tanh(Wh @ h + x)
        J = (1 - h**2)[:, None] * (Wh @ J)          # diag(1-h^2) @ Wh @ J
        s = np.linalg.norm(J, 2)
        acc += np.log10(s); J /= s                  # renormalise — never under/overflows
        out.append(acc)
    return out                                      # log10 of the accumulated norm

for rho in (0.5, 0.9, 1.5, 2.5):
    g = log10_grad_norm(rho)
    print(rho, [f"1e{g[t]:+.1f}" for t in (4, 9, 19, 49)])

# 0.5 ['1e-1.4', '1e-3.2', '1e-7.1', '1e-18.8']   <- dead well before step 20
# 0.9 ['1e-0.2', '1e-0.9', '1e-2.4', '1e-7.5']   <- at the noise floor by step 50
# 1.5 ['1e+0.8', '1e+0.8', '1e+0.6', '1e-0.0']   <- flat: saturation cancels the growth
# 2.5 ['1e+1.6', '1e+2.0', '1e+2.8', '1e+5.3']   <- exploding

And the two-line library form of clipping, which belongs in every RNN training loop you ever write. Note that it goes after backward() and before step() — it operates on gradients that exist but have not yet been applied:

python — clipping in the training loop
loss.backward()
total = torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=5.0)
optimizer.step()
optimizer.zero_grad()

# clip_grad_norm_ RETURNS the pre-clip norm — log it. A histogram of that
# number is the single most useful diagnostic for an RNN that is about to die.

Chapter 4: Gates — the GRU Cell

Chapter 3 handed us a precise diagnosis. Every backward step multiplies by D_tW_h, and any matrix repeatedly multiplied by itself either dies or explodes. So do not ask for a better W_h. Ask for a different multiplier — one the network can set to 1 when it wants a memory to survive.

Derive the fix from the requirement

Write down what we need and let the architecture fall out. We want an update where the network can choose, per unit and per timestep, between two behaviours:

• Keep: h_t = h_t−1. The state passes through untouched. The Jacobian ∂h_t/∂h_t−1 is the identity, so the gradient passes backwards untouched too. No decay, no explosion.
• Replace: h_t = h̃_t, some freshly computed content. The old value is overwritten, as an RNN cell must sometimes do.

A hard switch between the two is not differentiable, and we need gradients. So take the standard trick: replace the discrete choice with a continuous blend, and let a learned number control the blend. Call that number z:

where ⊙ is elementwise multiplication. At z = 0 this is exactly "keep". At z = 1 it is exactly "replace". In between it interpolates. This is a convex combination: the two coefficients are non-negative and sum to 1, so h_t is always a weighted average of the old state and the new candidate and can never blow up in magnitude. That property is free, and it is worth a lot.

z has to come from somewhere, and it has to be in [0, 1]. Compute it the way we compute everything else — a linear map of the input and the previous state — and squash it with the logistic sigmoid σ(a) = 1/(1 + e^−a), whose range is exactly (0, 1):

This is the update gate, and it is a vector of length n, not a scalar. Unit 7 can be keeping while unit 8 replaces, on the very same timestep. That per-unit independence is what lets one GRU hold a long-range fact in a few units while the rest of the state churns with local detail.

The highway, in numbers

That claim deserves arithmetic rather than assertion, so put the same 20-step distance from Chapter 3 next to three different gate settings. The GRU’s per-step backward multiplier is (1 − z_t); the vanilla cell’s was γ·σ_max ≈ 0.45.

Three things fall out of those four numbers. At z = 0.1 the gradient arrives at 12% strength after twenty steps — six orders of magnitude better than the vanilla cell, and comfortably above the float32 noise floor. At z = 0.5 the GRU is no better than the vanilla cell; a gate stuck at one half buys you nothing. And at z = 0.9 it is dramatically worse. The gate does not make gradients survive. It makes survival a choice the network can make per unit, and it will only choose it where the unit is carrying something the loss actually rewards remembering.

This is also why a trained GRU develops the striped structure you will see in Chapter 7: a handful of units sit near z = 0 for long stretches, and those are the ones with a live gradient path back hundreds of steps. The rest churn at z near 1 and remember nothing beyond a few characters, which is correct — most of what a character model needs is local.

The second gate: permission to ignore

One gate is not quite enough. The candidate h̃_t is the "fresh content", and the obvious way to compute it is exactly the vanilla RNN cell: tanh(W_hx_t + U_hh_t−1 + b_h). But that forces the candidate to be a function of the whole previous state, always. Consider reaching a sentence boundary in a language model: the right move is to propose content that ignores the previous clause entirely. With only an update gate, the cell can overwrite the state, but the replacement it computes is still contaminated by what it is replacing.

So add a second gate, applied to h_t−1 inside the candidate:

This is the reset gate. At r = 1 the candidate reads the full previous state, exactly like a vanilla cell. At r = 0 the candidate is computed from the input alone — the cell proposes a completely fresh start. Cho and colleagues, introducing the GRU in 2014, described exactly this behaviour: units with frequently-active reset gates learn short-lived local features, while units that keep r high learn long-range ones.

The complete cell, four lines

Shapes, because they are the thing people get wrong: every W is n × d, every U is n × n, every b is length n, and z, r, h̃, h are all length n. Three of each, one set per gate plus one for the candidate. Data flow for one timestep, concretely, with d = 128 and n = 256: a (128,) input and a (256,) previous state go in; three (256,) pre-activations are computed; two are squashed by sigmoid into gates and one by tanh into a candidate; they combine into a (256,) new state. Nothing else crosses the timestep boundary.

One step, by hand, every number

n = 2, d = 3 with vocabulary {a, b, c}. Previous state h_t−1 = [0.5, −0.2], input character b (so every W·x picks column 2). Weights:

The update gate. First the recurrent contribution, row by row:

The reset gate. Same shape of calculation:

The gated state, then the candidate. Apply r elementwise before the U_h multiply — this is the step people skip:

The blend. Unit by unit, spelled out:

Read the story in those numbers. Unit 1 had z = 0.67, so it took two thirds of its new value from the candidate — it is writing. Unit 2 had z = 0.38, so it kept 62% of its old value; the candidate pulled it upward from −0.2 toward +0.25, but only partway, landing at −0.03 — it is mostly remembering. And unit 1’s reset gate was 0.52, so its candidate was computed from a half-attenuated view of the past, while unit 2 at r = 0.67 let two thirds of the past through. Four independent decisions in one timestep, all learned, all differentiable.

See it: drive the gates yourself

Set the gates by hand and watch the arithmetic redo itself. There are four sliders, not two, because z and r are vectors: each hidden unit gets its own pair. The widget opens on exactly the values computed above — z = [0.6748, 0.3775], r = [0.5200, 0.6660] — so the readout should match the hand calculation digit for digit, and “Back to the hand example” returns you there.

Each number line shows h_t−1 and h̃_t as anchors with h_t sliding between them — the convex combination made visual. Pin a unit’s z to 0 and that unit freezes no matter what its candidate says; pin it to 1 and its history is erased in one step; sweep r and the candidate itself moves, because a different amount of past feeds into it. Try u1 writes, u2 holds: one unit overwriting while the other keeps, on the same timestep, is the per-unit independence the prose above claims — and the two gradient meters at the bottom show what it costs, with (1 − z) and its twentieth power for each unit side by side.

From scratch, then nn.GRUCell

python — GRU forward, from scratch, matching the hand numbers
import numpy as np
sigmoid = lambda a: 1 / (1 + np.exp(-a))

Wz = np.array([[0.0, 0.6, 0.0], [0.0, -0.4, 0.0]])   # only column 'b' matters here
Uz = np.array([[0.3, 0.1], [-0.2, 0.5]]);  bz = np.array([0.0, 0.1])
Wr = np.array([[0.0, -0.3, 0.0], [0.0, 0.7, 0.0]])
Ur = np.array([[0.2, -0.4], [0.1, 0.3]]);  br = np.array([0.2, 0.0])
Wh = np.array([[0.0, 0.5, 0.0], [0.0, 0.4, 0.0]])
Uh = np.array([[0.7, 0.2], [-0.1, 0.9]]);  bh = np.array([0.0, 0.0])

def gru_step(x, h):
    z     = sigmoid(Wz @ x + Uz @ h + bz)          # how much to WRITE   (n,)
    r     = sigmoid(Wr @ x + Ur @ h + br)          # how much to READ    (n,)
    h_bar = np.tanh(Wh @ x + Uh @ (r * h) + bh)    # the candidate       (n,)
    return (1 - z) * h + z * h_bar, z, r, h_bar   # convex blend

x = np.array([0., 1., 0.])                          # one-hot for 'b'
h, z, r, hb = gru_step(x, np.array([0.5, -0.2]))
print(np.round(z,4), np.round(r,4), np.round(hb,4), np.round(h,4))
# [0.6748 0.3775] [0.52 0.666] [0.5753 0.2488] [ 0.5508 -0.0306]
# every number matches the hand calculation above

PyTorch packs the three gates into one fused matrix for speed — weight_ih is (3n, d) and weight_hh is (3n, n), stacked in the order r, z, n (PyTorch calls the candidate n). One matmul instead of three is a real kernel-level win, and it is why library GRUs are several times faster than a loop of three separate multiplies:

python — the fused library cell
import torch, torch.nn as nn

cell = nn.GRUCell(input_size=3, hidden_size=2)
print(cell.weight_ih.shape, cell.weight_hh.shape)   # (6, 3) (6, 2)  == (3n, d), (3n, n)

# whole sequences at once, cuDNN-fused:
gru = nn.GRU(input_size=128, hidden_size=256, num_layers=1, batch_first=True)
x = torch.randn(8, 50, 128)          # (batch, time, d)
out, hT = gru(x)
print(out.shape, hT.shape)         # (8, 50, 256)  (1, 8, 256)
print(sum(p.numel() for p in gru.parameters()))   # 296,448 — see Chapter 5

One caution when you compare your implementation against PyTorch: their candidate applies the reset gate after the hidden matmul, as tanh(W_nx + r ⊙ (U_nh + b_hn)), rather than before it. It is a different but equally valid variant — cuDNN-friendly, because the hidden matmul does not depend on r and can be issued earlier. Expect small numerical differences, not a bug.

Chapter 5: GRU vs LSTM

The GRU was not the first gated cell. Hochreiter and Schmidhuber solved the vanishing-gradient problem in 1997 with the LSTM — Long Short-Term Memory — seventeen years before Cho and colleagues proposed the GRU. The GRU is best understood as an answer to the question "how much of the LSTM is load-bearing?" So: put them side by side, count what each costs, and see what the literature actually found.

The LSTM, in the same notation

Two structural differences jump out, and everything else follows from them.

First: the LSTM has two state vectors, the GRU has one. The LSTM keeps an internal cell state c — the actual long-term memory, never shown to the outside — and a separate hidden state h, which is a gated view of c produced by the output gate. The GRU has no such separation: its single h is both the memory and the output. So an LSTM can hold a fact in c while refusing to expose it (o ≈ 0), which is genuinely useful for latching a value that should influence a decision only much later. A GRU that wants to hide something must encode it in a direction the readout ignores — possible, but it has to learn that rather than getting it architecturally.

Second: the LSTM’s keep and write decisions are independent; the GRU’s are tied. In the GRU the coefficients are (1 − z) and z, so writing more necessarily means keeping less. In the LSTM, f and i are separate sigmoids. Set f = 1 and i = 1 and the cell state accumulates — the LSTM can add to memory without erasing any of it, which the GRU structurally cannot. That extra freedom is also the LSTM’s liability: c has no bound, can drift to large magnitudes, and needs that final tanh(c) to keep h in a sane range.

Count the parameters by hand

Take d = 128 (input width) and n = 256 (state width). Every gate and every candidate has the same three pieces: an input matrix W of shape n × d, a recurrent matrix U of shape n × n, and a bias b of length n. Count one block:

Now multiply by the number of blocks. The GRU has three (z, r, candidate); the LSTM has four (f, i, o, candidate):

So a GRU is exactly 25% smaller than an LSTM at the same width, and the general formulas are GRU = 3(nd + n² + n) and LSTM = 4(nd + n² + n). If you want to spend the saving on width instead, solve 3(nd + n²) = 394,240 for n at d = 128: n² + 128n − 131,413 = 0, giving n = (−128 + √(16,384 + 525,652))/2 = (−128 + 736.2)/2 ≈ 304. A 304-wide GRU costs about what a 256-wide LSTM costs — a genuinely useful trade when memory, not depth, is the binding constraint.

The gradient path, compared with numbers

Both cells beat the vanilla RNN for the same structural reason — an additive path whose backward multiplier the network controls — but the two multipliers are different quantities and it is worth seeing them side by side.

In the GRU, the state update is a convex blend, so the per-step Jacobian is approximately diag(1 − z_t). In the LSTM, the cell update c_t = f⊙c_t−1 + i⊙g has per-step Jacobian approximately diag(f_t). Same shape of result, different gate:

Identical for the two gated cells, and nine orders of magnitude better than the plain cell. Push the gate further: at z = 0.01 (or f = 0.99) the 30-step factor is 0.99³⁰ = 0.740 — the gradient arrives essentially intact. That is the whole benefit, and it is available to both architectures, which is why the literature keeps finding them within noise of each other.

The practical difference is where the default sits. An untuned LSTM initialises b_f at zero, so f starts at σ(0) = 0.5, and 0.5³⁰ = 9.3 × 10⁻¹⁰ — a cell that forgets fast before it has learned anything. Setting b_f = 1 moves the start to σ(1) = 0.731, and 0.731³⁰ = 8.3 × 10⁻⁵, which is nearly five orders of magnitude more gradient at initialisation. That single line of initialisation is most of the reported LSTM-versus-GRU gap in the papers that tuned it.

Compute cost, not just parameter cost

Parameters are memory; multiply-accumulates (MACs — one multiply plus one add, the unit of arithmetic work) are time. Per timestep, each block does one n×d matrix-vector product and one n×n product, so the GRU costs 3(nd + n²) MACs and the LSTM 4(nd + n²). At our sizes that is 295,000 versus 393,000 MACs per token — again the 3:4 ratio. For a 1,000-token sequence, roughly 0.3 versus 0.4 GMAC per sequence per layer.

But the number that actually decides wall-clock time on a GPU is neither of those. Both cells are sequential: step t cannot begin until step t − 1 finishes. A 1,000-token sequence means 1,000 dependent kernel launches, each doing a tiny matrix-vector product that leaves an A100 almost entirely idle. That is why the fused cuDNN kernel exists, why batch size matters far more for RNNs than for Transformers, and ultimately why the field moved on. Chapter 8 returns to this.

See it: the same sequence, both cells

Both cells are run below on an identical input sequence with identically-scaled random weights. Toggle between them and watch the state trajectories: the LSTM panel shows both c and h, so you can see the output gate hiding part of the memory; the GRU shows its single state. Slide the input drive to push the cells harder — the LSTM’s cell state grows past ±1 (it is unbounded) while the GRU’s state cannot leave [−1, 1], because a convex combination of things in [−1, 1] stays in [−1, 1]. The parameter panel recomputes both counts live as you change d and n.

When each one actually wins

Both cells, from scratch, side by side

python — the two cells, same interface
import numpy as np
sig = lambda a: 1/(1+np.exp(-a))

def gru_step(x, h, P):
    z = sig(P['Wz']@x + P['Uz']@h + P['bz'])
    r = sig(P['Wr']@x + P['Ur']@h + P['br'])
    c = np.tanh(P['Wh']@x + P['Uh']@(r*h) + P['bh'])
    return (1-z)*h + z*c                    # ONE state out

def lstm_step(x, h, c, P):
    f = sig(P['Wf']@x + P['Uf']@h + P['bf'])   # keep how much of c?
    i = sig(P['Wi']@x + P['Ui']@h + P['bi'])   # write how much of g? (INDEPENDENT of f)
    o = sig(P['Wo']@x + P['Uo']@h + P['bo'])   # reveal how much of c?
    g = np.tanh(P['Wg']@x + P['Ug']@h + P['bg'])
    c = f*c + i*g                            # unbounded accumulator
    return o*np.tanh(c), c                   # TWO states out

def n_params(d, n, gates):
    return gates * (n*d + n*n + n)
print(n_params(128,256,3), n_params(128,256,4))   # 295680 394240

python — the library swap is one word
import torch.nn as nn
gru  = nn.GRU (128, 256, batch_first=True)   # out, h    = gru(x)
lstm = nn.LSTM(128, 256, batch_first=True)   # out,(h,c) = lstm(x)   <- note the tuple

# the LSTM trick everyone forgets: bias the forget gate toward remembering.
# PyTorch packs biases as [i, f, g, o], so the forget slice is n : 2n.
for name, p in lstm.named_parameters():
    if 'bias' in name:
        n = p.shape[0] // 4
        p.data[n:2*n].fill_(1.0)      # f-gate bias = 1 -> sigmoid(1) = 0.73 keep by default

Chapter 6: Training a Character-Level GRU

We have a cell that remembers and a gradient that survives. Now assemble an actual trainable model: read a stream of characters, predict the next one at every position, and learn from the mistakes. This is the character-level language model, the standard proving ground for recurrent architectures, and it is exactly what Chapter 7 will run live in this page.

The loss: cross-entropy, computed by hand

At each timestep the readout produces d logits — one score per vocabulary item — and softmax turns them into probabilities. The cross-entropy loss is the negative log of the probability assigned to the character that actually came next:

Concretely, with a vocabulary of 5 and logits y = [2.0, 0.5, −1.0, 0.3, 1.2], target index 0. Exponentiate each:

Sanity-check that against the baselines. A uniform guesser scores −ln(1/5) = ln 5 = 1.6094, so this model is already well ahead of chance. A perfect model scores 0. The unit is nats; divide by ln 2 = 0.6931 to get bits per character, which is the number people quote for character models — 0.6444 / 0.6931 = 0.93 bits per character here.

The gradient of cross-entropy with respect to the logits is famously clean. All the softmax algebra cancels and what remains is "predicted minus actual":

Read the sign structure: the correct class gets a negative gradient (push its logit up), every wrong class gets a positive one (push theirs down), and the magnitudes sum to zero so the total "probability mass" is conserved. This single vector is where every gradient in the network originates.

Teacher forcing: what goes in at step t

There are two things you could feed the model at step t during training: the character the model predicted at step t − 1, or the character that actually appeared. Feeding the truth is called teacher forcing, and it is what essentially every sequence model is trained with.

The reason is conditioning, not laziness. Feed back the model’s own early-training predictions and they are noise, so the state is garbage, so the loss carries no useful signal about the cell’s parameters — the model has to get step 1 right before step 2 can teach it anything. Teacher forcing decouples the steps: every position gets a clean, correct history and learns its own conditional independently. Training is faster and vastly more stable.

The training loop, assembled

Sampling: temperature, with numbers

To generate, feed a seed character, take the output distribution, draw a character from it, feed that back, repeat. The only knob is temperature T, which divides the logits before softmax:

Take the same logits [2.0, 0.5, −1.0, 0.3, 1.2] and watch the top probability move:

Check the T = 0.5 row by hand: exponentiate [4.0, 1.0, −2.0, 0.6, 2.4] to get [54.598, 2.718, 0.135, 1.822, 11.023], sum 70.297, and 54.598/70.297 = 0.7767. Halving the temperature took the top character from 52% to 78%. In the limit T → 0 sampling becomes argmax; as T → ∞ it becomes uniform. T is not a property of the model — it is a decoding choice you make after training, and you can change it per generation.

The complete GRU, forward and backward, from scratch

This is the entire model — no autograd, no library cell, every derivative written out. It is the same algorithm the Chapter 7 lab runs in JavaScript. Read the backward pass slowly: each block undoes exactly one line of the forward pass, in reverse order.

python — char-level GRU, forward + backward, complete
import numpy as np
sig = lambda a: 1/(1+np.exp(-a))

def init(V, H, seed=0):
    g = np.random.default_rng(seed);  s = 0.08
    P = {}
    for k in ('z','r','h'):                 # three blocks: update, reset, candidate
        P['W'+k] = g.standard_normal((H,V))*s   # (H,V) input->hidden
        P['U'+k] = g.standard_normal((H,H))*s   # (H,H) hidden->hidden
        P['b'+k] = np.zeros(H)
    P['Wy'] = g.standard_normal((V,H))*s        # (V,H) readout
    P['by'] = np.zeros(V)
    return P

def loss_and_grads(P, inputs, targets, h0, V, H):
    """inputs/targets: lists of int char ids, len T. h0: (H,) carried state."""
    T = len(inputs)
    hs = {-1: h0.copy()}; zs, rs, cs, ps = {}, {}, {}, {}
    loss = 0.0

    # ---------- FORWARD ----------
    for t in range(T):
        i = inputs[t]                                  # one-hot => COLUMN LOOKUP
        hp = hs[t-1]
        zs[t] = sig(P['Wz'][:,i] + P['Uz']@hp + P['bz'])
        rs[t] = sig(P['Wr'][:,i] + P['Ur']@hp + P['br'])
        cs[t] = np.tanh(P['Wh'][:,i] + P['Uh']@(rs[t]*hp) + P['bh'])
        hs[t] = (1-zs[t])*hp + zs[t]*cs[t]
        y     = P['Wy']@hs[t] + P['by']
        y    -= y.max()                                # stability: shift before exp
        e     = np.exp(y); ps[t] = e/e.sum()
        loss += -np.log(ps[t][targets[t]] + 1e-12)

    # ---------- BACKWARD ----------
    G  = {k: np.zeros_like(v) for k, v in P.items()}
    dh_next = np.zeros(H)
    for t in reversed(range(T)):
        i, hp = inputs[t], hs[t-1]
        dy = ps[t].copy(); dy[targets[t]] -= 1          # p - onehot
        G['Wy'] += np.outer(dy, hs[t]); G['by'] += dy
        dh = P['Wy'].T @ dy + dh_next               # output path + recurrent path

        # h = (1-z)*hp + z*c
        dz  = dh * (cs[t] - hp)
        dc  = dh * zs[t]
        dhp = dh * (1 - zs[t])                        # <-- the identity highway

        # c = tanh(Wh[:,i] + Uh@(r*hp) + bh)
        da  = dc * (1 - cs[t]**2)
        G['Wh'][:,i] += da; G['bh'] += da
        G['Uh'] += np.outer(da, rs[t]*hp)
        drh  = P['Uh'].T @ da                        # grad wrt (r*hp)
        dr   = drh * hp
        dhp += drh * rs[t]

        # z = sigmoid(...)  and  r = sigmoid(...)
        daz = dz * zs[t] * (1 - zs[t])
        G['Wz'][:,i] += daz; G['bz'] += daz
        G['Uz'] += np.outer(daz, hp); dhp += P['Uz'].T @ daz

        dar = dr * rs[t] * (1 - rs[t])
        G['Wr'][:,i] += dar; G['br'] += dar
        G['Ur'] += np.outer(dar, hp); dhp += P['Ur'].T @ dar

        dh_next = dhp                                  # hand it one step further back

    for k in G: np.clip(G[k], -5, 5, out=G[k])       # explosion insurance
    return loss/T, G, hs[T-1]

Four details in there are the difference between working and not. dh = Wy.T@dy + dh_next — gradient arrives at h_t from two paths and both must be summed. Every parameter gradient uses += because of weight sharing. G['Wh'][:,i] += da touches only one column, because a one-hot input activates one column. And dhp = dh*(1-z) is the identity highway from Chapter 4, visible in the code: when z is near 0, the gradient passes back multiplied by ~1.

Wrap it in Adagrad — per-parameter learning rates that shrink for parameters receiving large gradients — and you have a complete trainer:

python — the training loop, truncated BPTT + Adagrad + sampling
text = open('input.txt').read()
chars = sorted(set(text)); V = len(chars); H, K, LR = 128, 25, 0.1
stoi = {c:i for i,c in enumerate(chars)}
data = [stoi[c] for c in text]

P  = init(V, H)
M  = {k: np.zeros_like(v) for k,v in P.items()}       # Adagrad accumulators
h  = np.zeros(H); pos = 0; smooth = -np.log(1.0/V)      # start at chance

for step in range(20000):
    if pos + K + 1 >= len(data):
        pos, h = 0, np.zeros(H)                        # wrap: only here do we reset memory
    xs, ys = data[pos:pos+K], data[pos+1:pos+K+1]      # targets = inputs shifted by one
    L, G, h = loss_and_grads(P, xs, ys, h, V, H)        # h carries; its GRAD does not
    for k in P:
        M[k] += G[k]*G[k]
        P[k] -= LR * G[k] / (np.sqrt(M[k]) + 1e-8)
    smooth = 0.999*smooth + 0.001*L
    pos += K

def sample(P, h, seed_id, n, temp=1.0):
    out, i = [], seed_id
    for _ in range(n):
        z = sig(P['Wz'][:,i] + P['Uz']@h + P['bz'])
        r = sig(P['Wr'][:,i] + P['Ur']@h + P['br'])
        c = np.tanh(P['Wh'][:,i] + P['Uh']@(r*h) + P['bh'])
        h = (1-z)*h + z*c
        y = (P['Wy']@h + P['by']) / temp        # temperature acts on LOGITS
        e = np.exp(y - y.max()); p = e/e.sum()
        i = np.random.choice(len(p), p=p)          # feed the SAMPLE back — exposure bias lives here
        out.append(i)
    return out

And the same model in PyTorch, for the version you would actually ship. Note how detach() is the one line implementing truncation, and view(-1, V) flattens time into the batch dimension so a single cross-entropy call covers all 25 positions:

python — the same model, library edition
import torch, torch.nn as nn

class CharGRU(nn.Module):
    def __init__(self, V, H):
        super().__init__()
        self.emb = nn.Embedding(V, H)      # the one-hot column lookup, made explicit
        self.gru = nn.GRU(H, H, batch_first=True)
        self.out = nn.Linear(H, V)
    def forward(self, x, h=None):
        e, h = self.emb(x), h
        o, h = self.gru(e, h)              # o: (B, T, H)
        return self.out(o), h             # logits: (B, T, V)

model = CharGRU(V, 128); opt = torch.optim.Adam(model.parameters(), lr=2e-3)
h = None
for xb, yb in loader:                      # xb, yb: (B, 25) int64
    logits, h = model(xb, h)
    h = h.detach()                        # TRUNCATION — value kept, graph cut
    loss = nn.functional.cross_entropy(logits.view(-1, V), yb.view(-1))
    opt.zero_grad(); loss.backward()
    nn.utils.clip_grad_norm_(model.parameters(), 5.0)
    opt.step()

Chapter 7: The Live Lab

Everything from Chapters 1 through 6 is now running in this page. Below is a real character-level GRU — the exact cell from Chapter 4, the exact backward pass from Chapter 6, translated to JavaScript and training in your browser on a small corpus, with no server and no library. Press Train and watch the loss fall.

This is not an animation of training. Every number you see is computed: the loss curve is the running cross-entropy, the sampled text comes from the current weights, and the gate strips are the actual z and r values from the most recent forward pass. If you set the learning rate to 5 and the model diverges, it diverged.

What each panel is telling you

• Loss curve — smoothed cross-entropy in nats per character. The dashed line is the chance baseline, ln(V), where an untrained model starts. Distance below that line is how much the model has actually learned.
• Generated text — sampled from the current weights at the chosen temperature, regenerated as training proceeds. Early on it is noise; then spaces appear in plausible places; then fragments of real words.
• Update gate z — one column per recent timestep, one row per hidden unit (the first 32 of them, once you push the slider past 32). Bright means z near 1 (this unit is overwriting itself); dark means z near 0 (this unit is holding a value). Rows that stay dark across many columns are your long-term memory units.
• Reset gate r — same layout. Bright means the candidate reads the past; dark means it ignores it. Watch for dark columns at word boundaries — the cell learning to start fresh.

Reading the loss curve, in numbers

The vertical axis is nats per character, and nats are only meaningful against a baseline. The corpus here has V = 24 distinct characters, so a model that has learned nothing scores ln(24) = 3.178 nats, which is the dashed line. Divide by ln 2 = 0.693 to convert: 3.178 / 0.693 = 4.59 bits per character. That is the price of naming one of 24 symbols with no knowledge at all.

Now watch what each stage of learning buys. The first thing any character model finds is the unigram distribution — how often each character appears, ignoring context. In this corpus space and e and t are far more common than k or v, and exploiting that alone typically drops the loss to around 2.7–2.9 nats within the first few dozen steps. No memory is involved; that is what the readout bias b_y learns, on its own, almost immediately.

The next drop is the bigram structure: after t comes h, after a space comes a word-initial letter, after q comes u. That needs the input matrix W_h but barely any state, and it takes the loss into the low twos. Everything below that is genuinely recurrent — it requires the state to be carrying information from more than one character back, and it is the part the gates make possible.

So read the curve in three regimes: still-high and slowly falling means letter frequencies; the middle band means pairs and triples; below that is real context. At the default H = 32 the model blows through all three in roughly a hundred steps — it is under 0.25 by step 200 — because 270 characters is small enough to compress almost perfectly. To watch the regimes separately, slow it down: at learning rate 0.001 the curve parks at 2.63 and stays there, which is the “frequencies and nothing else” plateau made visible. A curve that flattens there is telling you about capacity or step size, not about the data.

One caveat on the number itself: the loss displayed is an exponential moving average with a decay of 0.97, so it lags the instantaneous loss by roughly 1/(1 − 0.97) ≈ 33 steps. That smoothing is why the curve looks clean, and also why a divergence you cause with a huge learning rate takes a moment to show up on the plot even though the sampled text degrades instantly.

Experiments to run right now

The point of a lab is to break things on purpose. Each of these produces a distinct, diagnosable failure — run them and match the symptom to the cause you already know.

Why the browser can do this at all

Count the work. With H = 32 hidden units and V = 24 vocabulary items, one timestep costs three blocks × (H·V + H²) = 3 × (768 + 1024) = 5,376 multiply-accumulates forward, and roughly twice that backward. A 25-step chunk is therefore around 25 × 5,376 × 3 ≈ 400,000 operations per training step — a few milliseconds of plain JavaScript. Eight steps per animation frame is comfortable, even on a phone.

Push H to 64 and the H² term quadruples while the H·V term merely doubles, so cost grows faster than width — the same n² scaling that made Chapter 5’s parameter count dominated by the recurrent matrices. You can feel the frame rate change. That is the recurrence tax, measured on your own hardware.

The lab is the real test: if you can predict what each slider will do before you move it, and explain the result afterwards using z, r and the loss curve, you understand the GRU. One question to check that, then we follow the lineage forward.

Chapter 8: Beyond the Cell

A single GRU layer reading left to right is the atom. Real systems compose it three ways — stack it, run it in both directions, and split it into an encoder and a decoder — and the third of those runs into a wall that changed the whole field. Then we follow the lineage forward to the architectures that inherited from it.

Stacking: depth in the feature direction

A stacked (or deep) RNN feeds the output sequence of layer 1 into layer 2 as its input sequence. Layer 1 sees characters and produces a sequence of 256-dim states; layer 2 sees that sequence and produces its own. Each layer keeps its own hidden state and its own weights.

Count the cost, because it is not what people guess. With d = 128 and n = 256, layer 1 costs 3(256×128 + 256² + 256) = 295,680. Layer 2’s input is now 256-wide, not 128, so it costs 3(256×256 + 256² + 256) = 3 × 131,328 = 393,984. Two layers total 689,664 — not double the first layer but 2.33×, because the second layer’s input matrix grew. Every layer after the first costs the same 393,984.

Depth here buys abstraction, not reach. Layer 1 learns character-level regularities; layer 2 learns patterns over layer 1’s features. The temporal horizon does not improve just by stacking — each layer still faces the same recurrence. Two to four layers is the practical range; past that, gradients through depth and time become their own problem and you need residual connections between layers.

Bidirectional: reading the future

A bidirectional RNN runs two independent cells — one left-to-right, one right-to-left — and concatenates their states at each position. Output width doubles to 2n; parameters double to 591,360 for our example layer.

This is a strict upgrade for tagging tasks — named-entity recognition, part-of-speech, acoustic frame classification — because the label at position t genuinely depends on words after t. "Bank" is disambiguated by what follows it as much as by what precedes it.

And it is completely unusable for generation. To run the backward pass the model needs the whole sequence up front, but a language model generating token t by definition has not produced tokens t + 1 onward. The information is not merely unavailable, it does not exist yet. This is the same constraint that forces causal masking in decoder-only Transformers, and it is worth recognising as one constraint wearing two costumes.

Encoder–decoder: the bottleneck that created attention

For translation you need a variable-length input to become a variable-length output. The 2014 sequence-to-sequence answer: an encoder RNN reads the source sentence and its final hidden state becomes the context vector c; a decoder RNN is initialised with c and generates the target one token at a time.

Put a number on the squeeze. A 256-dim state is 256 floats regardless of input length. Encoding a 5-token sentence gives roughly 51 dimensions per token; a 50-token sentence gives 5.1; a 100-token paragraph gives 2.6. The budget is constant while the demand grows linearly. Cho and colleagues measured exactly this in 2014: encoder–decoder BLEU — an n-gram overlap score for translation quality, running from 0 to 100 — held up to about 20 source words and then fell off sharply, and the fall tracked sentence length rather than sentence difficulty.

Bahdanau, Cho and Bengio’s fix in 2015 was not a better encoder. It was to stop insisting on one vector: keep every encoder state h₁…h_T, and let the decoder compute a fresh weighted average of them at every output step, with weights it chooses based on what it is generating right now. That weighted read is attention. The length curve flattened out immediately.

Put a number on the squeeze

The information argument is worth making precisely, because it is the reason attention exists and people usually wave at it. Suppose a source vocabulary of 30,000 word types. Naming one word takes at most log₂(30,000) = 14.87 bits. A 50-word sentence therefore carries at most 50 × 14.87 = 743 bits of identity information; a 500-word document carries 7,436 bits.

Against that, a 256-dimensional float32 context vector is nominally 256 × 32 = 8,192 bits — but that is a wild overestimate of what it can actually convey, because the values are continuous, noisy, and the decoder has to read them through a learned map that is far from a lossless decoder. A realistic effective capacity is a few bits per dimension at best, so call it several hundred bits. The 50-word sentence sits comfortably inside that. The 500-word document does not, and no training run changes the arithmetic.

This is why the failure is gradual and length-dependent rather than abrupt. Cho and colleagues plotted BLEU against source length in 2014 and saw exactly that shape: flat to about 20 words, then a steady decline. The model was not getting worse at translation; it was running out of channel.

The lineage: what each successor kept and changed

The recurrent idea did not die; it went quiet for six years and came back once people wanted constant-memory inference again. Every architecture below keeps the fixed-size state that Chapter 1 introduced. What they change is how the state is updated — and almost all of them change it for the same reason: to make training parallel over time.

Notice the pattern in the "Changed" column. Between 2014 and 2017 the field traded the fixed state away for parallel training. From 2023 onward it has been trying to buy the fixed state back without giving the parallelism up — by finding recurrences whose update is linear enough in the state to be computed with an associative scan. The GRU’s nonlinear tanh candidate is precisely what blocks that, which is why none of the modern cells keep it.

See it: the family tree

Click any node to read what it kept, what it changed, and how it sits on the two axes that matter: how much of the state is fixed-size (vertical) and whether training parallelises over time (horizontal). The GRU sits in the top-left — fully fixed state, fully sequential training — and every arrow out of it is somebody trying to move right without falling down.

All three compositions, in code

python — stacked, bidirectional, and encoder-decoder
import torch, torch.nn as nn

# 1) STACKED — depth in features, not in time
deep = nn.GRU(128, 256, num_layers=2, dropout=0.2, batch_first=True)
out, h = deep(torch.randn(4, 50, 128))
print(out.shape, h.shape)     # (4,50,256)  (2,4,256)  <- one state PER LAYER
print(sum(p.numel() for p in deep.parameters()))   # 691,200 (with torch's 2 biases)

# 2) BIDIRECTIONAL — tagging only, never generation
bi = nn.GRU(128, 256, bidirectional=True, batch_first=True)
out, h = bi(torch.randn(4, 50, 128))
print(out.shape, h.shape)     # (4,50,512)  (2,4,256)  <- width DOUBLES

# 3) ENCODER-DECODER — and the bottleneck, made visible
class Seq2Seq(nn.Module):
    def __init__(self, V, H):
        super().__init__()
        self.enc = nn.GRU(H, H, batch_first=True)
        self.dec = nn.GRU(H, H, batch_first=True)
        self.emb = nn.Embedding(V, H); self.out = nn.Linear(H, V)
    def forward(self, src, tgt_in):
        _, c = self.enc(self.emb(src))       # c: (1,B,H) — EVERYTHING the source became
        o, _ = self.dec(self.emb(tgt_in), c) # the decoder only ever sees c
        return self.out(o)
# src of length 5 and src of length 500 produce the SAME size c.
# That single fact is what attention was invented to remove.

Chapter 9: Cheat Sheet & Connections

Everything in one place, with every symbol spelled out in words. Nothing new is introduced here — this is the page to come back to.

Every equation, every symbol

Which architecture, and why

Self-test: one GRU step, start to finish

Before the failure table, one worked example to check yourself against. Same cell as Chapter 4 — n = 2, previous state h_t−1 = [0.5, −0.2], the same U matrices and biases — but now the input character is a, so the W columns are W_z = [−0.2, 0.5], W_r = [0.4, −0.1], W_h = [0.3, −0.6]. Cover the answers and compute it; the recurrent contributions U_zh = [0.13, −0.20] and U_rh = [0.18, −0.01] are unchanged from Chapter 4, because h_t−1 is the same.

Compare with Chapter 4, where the same state met the character b and landed at [0.5508, −0.0306]. One character of difference moved unit 2 from −0.03 to −0.45, because a drives its candidate strongly negative and opens its update gate wider (z₂ = 0.60 here versus 0.38 there). That is the cell distinguishing two inputs from the same state — exactly the capability Chapter 0 said an MLP could not have.

Sizing a real one: the deployment budget

The other calculation worth being able to do from memory is the cost of a model you might actually ship. Take a two-layer GRU language model, input width d = 256, state width n = 512, vocabulary 50,000. One block costs n·d + n² + n:

Look at those last two lines together. The readout is nine times larger than the entire recurrent stack. That is the ordinary situation for character- and word-level language models at modest state widths, and it is why weight tying (sharing the readout matrix with the input embedding) is such a common trick: it removes more parameters than any change you could make to the cell.

Compute per token is 3(nd + n²) for layer 1 plus 3(2n²) for layer 2 = 1,179,648 + 1,572,864 = 2.75 million multiply-accumulates, which any phone does in well under a millisecond. The binding constraint is not arithmetic; it is that those steps are sequential.

And the payoff line, the one that keeps recurrence alive. Decode-time state is 2 layers × 512 floats × 4 bytes = 4,096 bytes, constant forever. A Transformer of the same width caches a key and a value per layer per token: 2 × 2 × 512 × 4 = 8,192 bytes per token, so at 10,000 generated tokens it is holding 82 MB and still growing. Four kilobytes versus eighty-two megabytes is the whole argument of Chapter 8’s bottom half, in two numbers.

The five things that actually break

Where to go next

If you want more foundations on recurrence:

• Recurrent Networks — RNN and LSTM built from absolute zero, with the memory intuition developed at a slower pace than here.
• RNN Sequences (CS231n) — the one-to-many, many-to-one and many-to-many patterns, image captioning, and Karpathy’s char-RNN in its original setting.
• CS224N Lecture 4: RNN Language Models — the Stanford treatment: n-grams, perplexity, and RNN language modelling as a course lecture.

If you want to follow the lineage forward:

• Transformer — where the attention patch from Chapter 8 became the whole architecture. Self-attention, multi-head, KV cache.
• xLSTM — the LSTM’s inventor fixing its two flaws in 2024: exponential gating and a parallelisable matrix memory.
• Linear Attention & RWKV — the family that writes attention as a recurrence, so it trains in parallel and decodes in constant memory.
• SSM / Mamba — state space models, the selective scan, and input-dependent recurrence at scale.
• Griffin & Hawk — a gated linear recurrence interleaved with local attention; the hybrid answer.

If you want to build with it:

• microGPT — every line of a working GPT, the natural next step once you have trained a character model here.
• All Gleams — the full catalogue.