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 an easy implementation bug to make 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 sequence unrolls into 10,000 cell applications. Keeping the states costs O(T·n) memory, in addition to other saved activations. The ordinary nonlinear recurrence introduces a sequential dependency across time, although operations within a step and across examples can run in parallel. Full BPTT may be feasible, but truncated BPTT reduces memory and compute by processing chunks of k steps. Values such as 25 or 35 are illustrative choices, not universal defaults. Carry the value of h across each chunk boundary while cutting its connection to the earlier gradient graph.

That distinction is the whole trick, and it is worth stating twice. The forward state is continuous — it can retain information from earlier chunks, if the model has learned to preserve it. 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):

Here, dropping one of three timestep contributions changes this gradient by about 7.6%; that is not a measurement of a 33% end-to-end runtime saving. Truncation is a biased approximation. When the omitted terms are small, it can save work with little gradient error. When distant credit matters, it discards useful information. Signed terms can also cancel, so the error need not improve monotonically with every added step. Choose k by measuring the task, memory budget, and learning behavior; gated cells can also be limited by too short a window.

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: a product of time-varying Jacobians. The shared recurrent matrix appears in every factor, while each timestep contributes its own diagonal tanh slope. A repeated matrix power is a useful approximation, but carefully chosen orthogonal or unit-gain cases can preserve a signal longer than the generic attenuating or amplifying cases.

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⁻⁷. Float32 can represent both this value and 4.6 × 10⁻¹⁸. Its precision is relative, not an absolute cutoff at 10⁻⁷. A small gradient contribution can be lost when added to much larger terms or become too small to change a parameter. That depends on the surrounding scales and optimizer; attenuation is not the same as underflow.

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

A common safeguard against unusually large gradients is gradient norm 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.

For a finite nonzero gradient, this preserves its direction and caps its norm. With plain SGD, the learning rate then sets the update norm; adaptive optimizers also transform the gradient. Clipping alone does not guarantee a loss decrease or repair values that already overflowed. Thresholds such as 1 to 5 are common starting points, but the useful scale depends on the setup.

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']   <- strongly attenuated 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 you can use when large gradients are a problem. 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 direct keep path passes through untouched. Its Jacobian contribution is the identity, while the full cell Jacobian also includes derivatives through the gate and candidate.
• 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 — about six orders of magnitude larger than this vanilla-cell example. At z = 0.5 the direct path is about eight times larger, but both still attenuate strongly. At z = 0.9 this direct path is smaller. These products isolate the keep path; they are not the full GRU Jacobian. 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. The original analysis says that a reset gate near zero lets a unit behave as if it were reading the first symbol of a sequence: the candidate ignores the previously computed state, while the final blend may still retain it. The broader lesson is that the gates let different units adapt to different time scales. Nobody assigns those time scales. They emerge.

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 can reduce kernel-launch overhead. The actual speedup depends on hardware, batch size, and the implementation:

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 + b_in + r ⊙ (U_nh + b_hn)), rather than before it. PyTorch also uses h = z⊙h_old + (1−z)⊙n, so its update gate is the complement of this lesson’s write gate. Matrix multiplication and elementwise reset do not generally commute; copying the same numbers into both implementations does not guarantee matching or even small numerical differences. Compare the complete formulas and parameter conventions first.

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 has exactly 25% fewer matrix-and-bias parameters 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

For autoregressive sequence generation, training commonly feeds the ground-truth token from step t − 1 while predicting the token at t. This is teacher forcing. Feeding model samples instead is another possible training policy; teacher forcing is common, but it is not used by every kind of sequence model.

Ground-truth prefixes give a supervised next-token loss at every position even when the model is initially wrong. They do not decouple the recurrence: h_t still depends on earlier states and shared parameters, and BPTT still couples their gradients.

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. One decoding control 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. Greedy or beam search, top-k, nucleus sampling, and repetition controls are other choices; temperature alone does not guarantee coherence or creativity.

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:                      # contiguous streams, stable batch ordering
    logits, h = model(xb, h)                            # otherwise reset h = None per batch
    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 exactly uniform reference, ln(V), which also initialises the displayed moving average. A lower training loss is improvement on these chunks relative to uniform; it does not measure held-out learning.
• Generated text — sampled from the current weights at the chosen temperature, regenerated as training proceeds. Compare samples across steps; no fixed progression is guaranteed.
• 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; dark means the direct keep path is strong. Persistent dark rows are candidates for further analysis, not proof of a semantic role.
• Reset gate r — same layout. Bright means more prior state enters the candidate; dark means less. Boundary-aligned patterns can suggest a hypothesis, which requires controlled tests.

Reading the loss curve, in numbers

The vertical axis is nats per character. The corpus has V = 24 distinct characters, so an exactly uniform predictor scores ln(24) = 3.178 nats, the dashed line. Divide by ln 2 = 0.693 to get 4.59 bits per character. Random initial logits are often near uniform, but an arbitrary nonuniform predictor can score above or below this reference.

Before attributing a loss change to memory, compute context-free baselines. For this 273-character corpus, the empirical unigram distribution has entropy 2.653 nats: for example, e appears 51 times and space 50 times, while k appears 3 times and v once. A frequency-only predictor evaluated on that same distribution attains this value, but the training curve does not say which parameters learn frequency information or when.

Context-conditioned predictions can exploit repeated phrases such as the gate and the state. The input weights, recurrent state, gates, and readout work jointly, so a loss drop does not isolate bigrams or persistent memory. Test those hypotheses with controlled ablations or interventions.

Loss values do not identify which statistics the network uses: unigram, short-context, and longer-state solutions can overlap in average cross-entropy. This tiny repeated corpus can eventually be fitted very closely, but exact loss-at-step values are properties of this implementation and seed. A plateau says only that the measured objective has stopped improving; test optimisation, width, truncation, data, and implementation hypotheses with controlled runs.

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 recurrent affine work. With H = 32 and V = 24, one forward timestep uses three blocks times (H·V + H²) = 3 × (768 + 1024) = 5,376 multiply-accumulates. Treating backward as roughly twice that gives about 403,200 recurrent multiply-accumulates for a 25-step training chunk. This estimate excludes the readout, nonlinearities, sampling, and browser overhead; wall time must be measured on the target device.

At H = 64, the same recurrent estimate is 16,896 multiply-accumulates per forward timestep, about 3.14 times the H = 32 value, because the H² term grows faster than H·V. The lab can show local timing differences, but browser scheduling and hardware determine whether they affect frame rate.

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.

Stacking lets higher recurrent layers transform the sequence representations produced below them. It does not by itself guarantee a longer usable temporal horizon or a particular abstraction in each layer. Deeper stacks add optimisation and compute costs; useful depth and choices such as residual connections depend on the task and implementation.

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 often an advantage 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.

For strictly causal generation it is unavailable. 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 fixed-context bottleneck and 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.

A 256-dimensional state has the same shape at every source length, but dividing dimensions by tokens is not an information measure. Cho and colleagues reported that their 2014 systems performed worse on longer sentences and on sentences with more unknown words. The often-cited change near 20 words is a feature of that experiment, not a universal encoder limit.

Bahdanau, Cho and Bengio called the fixed-length vector a conjectured bottleneck and proposed an attention-based alternative: retain encoder states and compute a weighted source summary for each decoder step. Their model learned soft alignments and improved translation results; the paper does not establish a universal length threshold or instant removal of every length effect.

Measure the bottleneck

A fixed-width vector must support a harder learned compression task as source length and complexity grow. There is no justified conversion from float dimensions to a fixed number of recoverable token bits: the representation is task-dependent, lossy, distributed, and learned jointly with the decoder. The evidence here is empirical translation quality, not a dimensionality proof.

The lineage: what each successor kept and changed

Recurrent work continued while attention and convolution drew more attention, then gated and memory-efficient variants brought the fixed-state trade-off back into focus. The architectures below make different choices about how the state is updated and how much parallel training they can recover.

The Transformer replaced recurrent sequence mixing with attention, enabling parallel processing across training positions. Recent recurrent and state-space architectures explore fixed or compressed decode state together with scan-friendly or otherwise parallel training forms. A standard GRU update is not an associative linear recurrence, so these systems make different compromises rather than forming one uniform family or a single line of descent.

See it: the family tree

Click any node to read a compact comparison on two useful axes: how much decode state is fixed-size (vertical) and whether training parallelises over time (horizontal). Positions and connecting lines are conceptual, not benchmark scores or claims of direct genealogy.

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 — non-causal context, not a direct left-to-right generator
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), one learned source summary
        o, _ = self.dec(self.emb(tgt_in), c) # this simple decoder receives the source through c
        return self.out(o)
# Source lengths 5 and 500 produce the same-shaped c.
# Attention instead lets each decoder step read across encoder states.

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:

The 50,000-way readout has about 9.3 times as many parameters as the recurrent stack in this example. Direct input/output weight tying requires matching embedding and hidden widths; with d = 256 and n = 512, it needs an added projection or a different width choice.

The recurrent stack costs about 2.75 million multiply-accumulates per token, but the dense 50,000-way readout adds about 25.6 million. Actual latency depends on hardware and kernels. The recurrence also makes token steps sequential.

The two-layer GRU carries 2 × 512 float32 values, or 4,096 bytes, independent of prefix length. For a two-layer standard full-attention decoder whose key and value widths are each 512, an uncompressed float32 KV cache adds 2 × 2 × 512 × 4 = 8,192 bytes per token, or about 82 MB at 10,000 tokens. Grouped-query attention, lower-precision caches, sliding windows, and compression change this comparison, so deployment requires the actual architecture and cache policy.

Five failure patterns to investigate

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 — 2024 LSTM variants with exponential gating and distinct scalar- and matrix-memory structures.
• Linear Attention & RWKV — time-decay linear attention with parallel and recurrent forms.
• SSM / Mamba — selective state-space models and a hardware-aware scan.
• Griffin & Hawk — gated linear recurrence, with local attention interleaved in Griffin.

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.