RNNs, LSTMs, GRUs, and Sequence Modeling

An on-call assistant sees error logged, then reproduction confirmed, then rollback requested. Reverse those events and the latest evidence changes: an error now arrives after the rollback request. A follow-up event can establish whether the incident is resolved. Event order shapes the next decision, so how can a model track it with one architecture across histories of different lengths?

The previous lesson converted raw class scores, or logits, into a softmax probability distribution. We'll keep that classifier at the tail end, outputting decisions across three actions: ignore, rollback, and escalate. Our job here is to transform a variable-length stream of events into the single, fixed-size vector that the downstream classifier expects.

Order is part of the input

If you shuffle raw event vectors, a bag-of-events sum or average treats the history as identical because those operations are commutative. An order-aware model can distinguish the histories. Compressed representations can still collide, and a model can learn to ignore order that carries no useful signal for its task.

The tiny accumulator below is an arithmetic fixture. Each incoming event updates a running state influenced by earlier events, with older evidence discounted.

Before running it, predict which order leaves the larger final state. The same three values arrive, but the recurrence discounts earlier state by 0.6, so the event in the final position exerts the strongest direct influence:

signals = {
    "error logged": 1.0,
    "reproduction confirmed": 0.4,
    "rollback requested": 0.9,
}

def running_state(events):
    state = 0.0
    for event in events:
        state = 0.6 * state + signals[event]
    return state

incident_order = ["error logged", "reproduction confirmed", "rollback requested"]
reversed_order = list(reversed(incident_order))

print("incident order:", round(running_state(incident_order), 3))
print("reversed order:", round(running_state(reversed_order), 3))
print("same events:", sorted(incident_order) == sorted(reversed_order))

incident order: 1.5
reversed order: 1.564
same events: True

The incident order computes 0.6 * (0.6 * 1.0 + 0.4) + 0.9 = 1.5. Reversing it puts the largest event 1.0 in the final position, producing 0.6 * (0.6 * 0.9 + 0.4) + 1.0 = 1.564. A simple sum would yield 2.3 in both cases; the recurrent update makes sequence order inseparable from the final value.

One cell, reused at every step

A recurrent neural network (RNN) makes this update rule trainable. Its hidden state ($h_t$) is a running vector of numbers carried from one event to the next, serving as working memory. The update rule is called a recurrent cell, and each new event advances the sequence by one timestep.

Start with a scalar hidden state. Multiply the incoming event by 0.8, multiply the previous hidden state by 0.5, add them, and pass the result through hyperbolic tangent (tanh). Like sigmoid, tanh is a smooth S-shaped activation function with output in $(-1, 1)$. It is symmetric around zero: $\tanh(-a)=-\tanh(a)$. A batch's activation mean still depends on its input distribution. With an initial state $h_0 = 0.0$ and a first event of 1.0, the weighted sum is 0.8, giving $h_1 = \tanh(0.8) \approx 0.664$.

For the second event (0.4), combine 0.8 * 0.4 with 0.5 * 0.664. That sum is about 0.652, yielding $h_2 = \tanh(0.652) \approx 0.573$. Run all three steps with these demonstration weights:

import math

events = [1.0, 0.4, 0.9]
w_input, w_state, bias = 0.8, 0.5, 0.0
state = 0.0

for step, event in enumerate(events, start=1):
    state = math.tanh(w_input * event + w_state * state + bias)
    print(f"h_{step} = {state:.3f}")

h_1 = 0.664
h_2 = 0.573
h_3 = 0.764

Three scalar inputs condense into one final scalar summary. When we expand the hidden state into a vector, different coordinates can track distinct patterns in parallel. The fundamental vector RNN recurrence equation is:

Here $x_t \in \mathbb{R}^{d_x}$ is the current input vector, $h_{t-1} \in \mathbb{R}^{d_h}$ is the prior hidden state, and $b_h \in \mathbb{R}^{d_h}$ is a bias vector. The weight matrices $W_{xh} \in \mathbb{R}^{d_h \times d_x}$ and $W_{hh} \in \mathbb{R}^{d_h \times d_h}$ project the new observation and transition the prior memory, respectively. Tanh applies element-wise to squash each hidden coordinate into $(-1, 1)$.

The same parameter matrices $W_{xh}$, $W_{hh}$, and bias $b_h$ execute at step 1, step 2, and step 1,000. Every timestep reuses those shared parameters, so sequence length doesn't add learned copies.^{[1]Reference 1On the Difficulty of Training Recurrent Neural Networkshttps://proceedings.mlr.press/v28/pascanu13.html}

Trace a vector RNN

Represent each incident event with three binary features:

These are one-hot event vectors: exactly one coordinate is 1 and all others are 0. They describe the incoming event at the current timestep, not cumulative flags.

Set up a two-dimensional hidden state ($d_h = 2$) and an initial state $h_0 = [0.0, 0.0]^T$, with biases set to zero for clean manual verification. At the first event $x_1 = [1, 0, 0]^T$, the matrix-vector multiplication $W_{xh} x_1$ extracts the first column of $W_{xh}$, which is $[0.8, 0.1]^T$. Since $h_0 = [0, 0]^T$, the recurrent term $W_{hh} h_0$ contributes nothing. The cell computes $[\tanh(0.8), \tanh(0.1)]^T \approx [0.664, 0.100]^T$.

The tensor dimensions dictate how the calculations scale. Let $d_x$ be input feature width and $d_h$ be hidden width:

Each row of $W_{hh}$ and $W_{xh}$ computes one hidden coordinate. Step through all three events in pure Python:

import math

def matvec(matrix, vector):
    return [sum(weight * value for weight, value in zip(row, vector, strict=True)) for row in matrix]

def add(left, right):
    return [a + b for a, b in zip(left, right, strict=True)]

def tanh_vec(values):
    return [math.tanh(value) for value in values]

events = [
    [1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
    [0.0, 0.0, 1.0],
]
W_xh = [
    [0.8, 0.2, 0.1],
    [0.1, 0.7, 0.4],
]
W_hh = [
    [0.5, 0.1],
    [0.0, 0.6],
]

state = [0.0, 0.0]
for step, event in enumerate(events, start=1):
    state = tanh_vec(add(matvec(W_xh, event), matvec(W_hh, state)))
    print(f"h_{step} = {[round(value, 3) for value in state]}")

print("final width:", len(state))

h_1 = [0.664, 0.1]
h_2 = [0.494, 0.641]
h_3 = [0.39, 0.655]
final width: 2

View larger

Attach a classification head to that final hidden state: a linear projection followed by softmax. This reuses the multiclass head from the previous lesson, with weight matrix $W_{\text{out}} \in \mathbb{R}^{3 \times 2}$ mapping the 2D summary into 3 action logits.

For row 1 (ignore), the pre-activation logit is $0.7 \times 0.390 - 0.3 \times 0.655 = 0.0765$. Rows 2 and 3 produce $0.368$ (rollback) and $0.170$ (escalate). Because $0.368$ is the largest logit, we know rollback will claim the highest softmax probability before running math.exp:

import math

final_state = [0.390, 0.655]
W_out = [
    [0.7, -0.3],   # ignore
    [-0.4, 0.8],   # rollback
    [0.1, 0.2],    # escalate
]
labels = ["ignore", "rollback", "escalate"]

logits = [
    sum(weight * value for weight, value in zip(row, final_state, strict=True))
    for row in W_out
]
shift = max(logits)
weights = [math.exp(logit - shift) for logit in logits]
total = sum(weights)
probabilities = [weight / total for weight in weights]

print("logits:", [round(logit, 3) for logit in logits])
print("probabilities:", [round(probability, 3) for probability in probabilities])
print("prediction:", labels[probabilities.index(max(probabilities))])

logits: [0.076, 0.368, 0.17]
probabilities: [0.291, 0.389, 0.32]
prediction: rollback

With these hand-chosen weights and the rounded final state, the head selects rollback with 38.9% probability. This arithmetic trace provides no accuracy measurement for a trained incident classifier. The head receives the fixed-width state $[0.390, 0.655]^T$; that state can encode some information about the full history and its length.

Count the parameters in a trainable version with one bias vector per affine block: $W_{xh}$ has 6 weights, $W_{hh}$ has 4, and $b_h$ has 2 biases, giving 12 recurrent parameters. A head with 6 weights and 3 trainable biases adds 9, totaling 21. The manual trace used zero biases. These shared weights keep the parameter count fixed as sequence length changes:

input_size, hidden_size, output_size = 3, 2, 3
rnn_parameters = hidden_size * input_size + hidden_size * hidden_size + hidden_size
head_parameters = output_size * hidden_size + output_size

for timesteps in [3, 30, 300]:
    print(f"{timesteps:>3} events -> {rnn_parameters + head_parameters} parameters")

  3 events -> 21 parameters
 30 events -> 21 parameters
300 events -> 21 parameters

The final-state interface is a compression bottleneck: useful evidence from the history must be encoded in a fixed-width vector, and some details can be lost. During training, credit assignment to distant events travels backward through the recurrent channel.

Backpropagation through time: why a distant signal can fade

The forward pass assigned its highest softmax probability to rollback. During training, a target label and cross-entropy loss measure the error of that prediction. Backpropagation through time (BPTT) applies the chain rule backward across the unrolled sequence to compute gradients for the shared weight matrices.^{[1]Reference 1On the Difficulty of Training Recurrent Neural Networkshttps://proceedings.mlr.press/v28/pascanu13.html}

Because $W_{hh}$ acts at every timestep, its gradient accumulates local contributions across time. For a summed loss $L=\sum_{t=1}^T L_t(h_t)$, we can write:

Here $\frac{\partial^+ h_k}{\partial W_{hh}}$ is the immediate local derivative at timestep $k$ with previous state $h_{k-1}$ held constant. The key term governing credit assignment over long sequences is the temporal chain:

To see why this chain causes instability, look at the one-dimensional cell we ran earlier: $h_j = \tanh(0.8 x_j + 0.5 h_{j-1})$. The local derivative from one state to the next is:

Because $\tanh'(z) = 1 - \tanh(z)^2 = 1 - h_j^2$, its mathematical derivative lies in $(0,1]$ for finite inputs. With recurrent weight $0.5$, the local factor is at most $0.5$, attaining that bound at $h_j=0$. Trace the product in our example:

import math

events = [1.0, 0.4, 0.9]
w_input, w_state = 0.8, 0.5
state = 0.0
factor = 1.0

for step, event in enumerate(events, start=1):
    state = math.tanh(w_input * event + w_state * state)
    local = (1.0 - state * state) * w_state
    factor *= local
    print(f"step {step}: h={state:.3f}  local factor={local:.6f}  product so far={factor:.6f}")

step 1: h=0.664  local factor=0.279528  product so far=0.279528
step 2: h=0.573  local factor=0.335820  product so far=0.093871
step 3: h=0.764  local factor=0.207910  product so far=0.019517

The temporal derivative is $\partial h_3 / \partial h_0 \approx 0.0195$. If a final loss supplies derivative $g = \partial L / \partial h_3$, the gradient that reaches $h_0$ is only 0.0195 * g. That's less than 2% of the learning signal surviving after just three recurrent steps.

When every step compounds a factor of $0.6$, the shrinkage is exponentially steep:

for steps in [1, 4, 8, 12]:
    gradient_factor = 0.6 ** steps
    print(f"{steps:>2} recurrent steps: {gradient_factor:.6f}")

 1 recurrent steps: 0.600000
 4 recurrent steps: 0.129600
 8 recurrent steps: 0.016796
12 recurrent steps: 0.002177

A matrix bound, and what it actually proves

In the general multivariate case, the one-step derivative $J_j = \frac{\partial h_j}{\partial h_{j-1}}$ is a Jacobian matrix. Because $h_j = \tanh(W_{hh} h_{j-1} + W_{xh} x_j + b_h)$:

Backpropagating an error gradient column vector from timestep $T$ back to timestep $k$ multiplies by transposed Jacobians in reverse chronological order:^{[2]Reference 2Learning Long-Term Dependencies with Gradient Descent Is Difficulthttps://pubmed.ncbi.nlm.nih.gov/18267787/[1]Reference 1On the Difficulty of Training Recurrent Neural Networkshttps://proceedings.mlr.press/v28/pascanu13.html}

This expression assumes a loss only at the final state. When intermediate positions have losses, their gradient contributions must also be added. Matrix order matters: apply $J_T^{\mathsf{T}}$ first, then work backward. Using the induced Euclidean norm and its submultiplicative property gives:

where $\gamma = \max_j \|\operatorname{diag}(1 - h_j^2)\| \le 1$, and $\sigma_{\max}(W_{hh})$ is the maximum singular value (spectral norm) of the recurrent matrix.

If one uniform bound satisfies $\gamma\sigma_{\max}(W_{hh})<1$ across the span, the temporal Jacobian contracts exponentially. This is a sufficient condition for a vanishing path. The total parameter gradient can still include strong contributions from recent positions.

If the bound exceeds one, growth is possible; actual behavior depends on direction, successive activations, and alignment. For example, a matrix can have a singular value above one while its repeated powers eventually decay. Spectral-radius arguments describe constant linear recurrences; nonlinear, input-dependent RNNs require Jacobian-product analysis.^{[1]Reference 1On the Difficulty of Training Recurrent Neural Networkshttps://proceedings.mlr.press/v28/pascanu13.html}

Float32 decay depends on the path factors and parameter scale. For our illustrative factor, $0.6^{50}\approx8.08\times10^{-12}$ is still representable in float32. A small gradient may have little optimization effect before it underflows; whether a parameter update rounds away also depends on the parameter's scale and learning rate. In an expanding path, large gradients can instead destabilize updates or become nonfinite.

View larger

Truncation and gradient clipping

Two operational techniques manage these dynamics during training:

Truncated BPTT limits how far backward the computational graph extends. Training can split a long sequence into chunks, carrying state values forward while h = h.detach() cuts the graph at each boundary. An LSTM needs both h and c detached. This reduces graph memory but also removes credit assignment across the boundary; the forward state can retain evidence the current loss can no longer train earlier steps to store.

Hold the weights fixed and cut our three-event sequence after its second event. Predict which changes: the final state, the derivative with respect to the first event, or the derivative with respect to the shared recurrent weight?

import torch

def trace(cut_after_second):
    events = torch.tensor([1.0, 0.4, 0.9], dtype=torch.float64,
                          requires_grad=True)
    weight = torch.tensor(0.5, dtype=torch.float64, requires_grad=True)
    state = torch.zeros((), dtype=torch.float64)
    for step, event in enumerate(events, start=1):
        state = torch.tanh(0.8 * event + weight * state)
        if cut_after_second and step == 2:
            state = state.detach()
    # Inspect sensitivity of the final state; a loss multiplies these
    # derivatives by its derivative with respect to that state.
    state.backward()
    return state.item(), events.grad, weight.grad.item()

full_state, full_events, full_weight = trace(False)
cut_state, cut_events, cut_weight = trace(True)
assert full_state == cut_state
assert full_events[0] > 0 and cut_events[0] == 0
print(f"final states: full={full_state:.6f}, truncated={cut_state:.6f}")
print("event derivatives, full:", [round(g, 6) for g in full_events.tolist()])
print("event derivatives, truncated:", [round(g, 6) for g in cut_events.tolist()])
print(f"weight derivative: full={full_weight:.6f}, truncated={cut_weight:.6f}")

final states: full=0.764316, truncated=0.764316
event derivatives, full: [0.031227, 0.111713, 0.332657]
event derivatives, truncated: [0.0, 0.0, 0.332657]
weight derivative: full=0.331003, truncated=0.238277

The boundary preserves the number in the state but severs its history. The final event can still train the shared weight through its local use; it cannot train the weight's earlier uses through this final-state objective. Our fixture performs no optimizer step between chunks. In training, a step between chunks also changes the weights, so the carried state may have been computed with an earlier parameter version.

Gradient norm clipping addresses exploding gradients. If the $L_2$ norm of the gradient vector exceeds a threshold $c$, we rescale the vector:

import math

raw_x = 1.4 ** 12
raw_y = -0.8 * raw_x
raw_norm = math.hypot(raw_x, raw_y)
max_norm = 5.0
scale = min(1.0, max_norm / raw_norm)
clipped = (raw_x * scale, raw_y * scale)

print("raw norm:", round(raw_norm, 3))
print("clipped norm:", round(math.hypot(*clipped), 3))
print("clipped gradient:", [round(value, 3) for value in clipped])

raw norm: 72.604
clipped norm: 5.0
clipped gradient: [3.904, -3.123]

For a finite, nonzero gradient, norm clipping rescales its magnitude while preserving direction in exact arithmetic.^{[1]Reference 1On the Difficulty of Training Recurrent Neural Networkshttps://proceedings.mlr.press/v28/pascanu13.html} It bounds the gradient norm. For plain SGD without momentum or weight decay, the step norm is at most $\eta c$; adaptive optimizers and momentum change that relationship. With $g=0$, leave it unchanged instead of evaluating $c/0$. Detect and diagnose nonfinite gradients before rescaling.

Clipping limits oversized gradients, while gates offer a different memory path. Initialization, regularization, training objectives, and dependency length also affect gradient flow. Long-range learning depends on how these choices interact with the data and task.

LSTM: Additive cell state and gated memory

Gradient clipping limits oversized updates. LSTMs also provide an additive cell-state path that carries stored values across updates without applying a saturating activation to the cell itself.

Long short-term memory (LSTM) addresses this difficulty with two related state vectors:^{[3]Reference 3Long Short-Term Memory.https://doi.org/10.1162/neco.1997.9.8.1735}

• Internal cell state ($c_t \in \mathbb{R}^{d_h}$): A vector the cell updates by addition, so a stored coordinate can move forward without passing through tanh on that step.
• Working hidden state ($h_t \in \mathbb{R}^{d_h}$): The exposed memory vector passed to the next layer and the output classifier.

Information flow is controlled by gates: learned, coordinate-wise affine projections squashed by the sigmoid function ($\sigma(z) \in (0, 1)$). A gate value near 0 acts as a closed valve; a value near 1 leaves the valve fully open. The standard modern formulation includes the forget gate introduced by Gers et al.^{[4]Reference 4Learning to Forget: Continual Prediction with LSTMhttps://doi.org/10.1162/089976600300015015}

Suppose a single coordinate in our cell state holds 0.855 after an error event. A routine health check arrives. We want to retain that earlier error flag while recording almost nothing new from the routine check:

The stored value ($c_t \approx0.840$) differs from the exposed state ($h_t \approx0.549$). Closing the output gate can hide a stored coordinate without immediately erasing it. These states are still coupled: the hidden state depends on the cell, and future gates depend on the hidden state. The example chooses gate values by hand to trace the cell update.

Step through three illustrative transitions in pure Python:

import math

cell = 0.0
steps = [
    ("error recorded",    0.00, 0.95, 0.90, 0.80),
    ("routine check",     0.98, 0.02, 0.10, 0.80),
    ("rollback verified", 0.97, 0.15, 0.50, 0.90),
]

for label, forget, write, candidate, expose in steps:
    cell = forget * cell + write * candidate
    hidden = expose * math.tanh(cell)
    print(f"{label:<20} c={cell:.3f}  h={hidden:.3f}")

error recorded       c=0.855  h=0.555
routine check        c=0.840  h=0.549
rollback verified    c=0.890  h=0.640

The full LSTM equations

In a trainable LSTM, learned affine projections compute these gate vectors from the concatenation of the prior working state $h_{t-1}$ and the current input $x_t$ (denoted $[h_{t-1}, x_t] \in \mathbb{R}^{d_h + d_x}$):

Here $\odot$ denotes the element-wise Hadamard product.

Follow the direct cell-state gradient path

In this standard LSTM without peephole connections, the cell receives $(h_{t-1},c_{t-1})$ as separate state inputs. Holding $h_{t-1}$ and $x_t$ fixed also holds the gates and candidate fixed. The direct cell-state Jacobian is:

The contribution from following only this direct path across several steps is:

Contrast this directly with the vanilla RNN product $\prod \operatorname{diag}(1 - h_j^2) W_{hh}$:

1. Direct diagonal path: No $W_{hh}^{T-k}$ factor acts along this channel.
2. Additive cell update: The cell state passes through no $\tanh$ during the update.
3. Controllable retention: Forget gates near one slow decay. If every gate on a coordinate is 0.999 for 100 transitions, its direct-path factor is about 0.905, retaining 90.5%. At 0.5, the same product is approximately $7.89\times10^{-31}$.

The original LSTM's unit-weight cell self-connection was called the constant error carousel.^{[3]Reference 3Long Short-Term Memory.https://doi.org/10.1162/neco.1997.9.8.1735} A forget-gated LSTM learns retention instead of always copying the cell.^{[4]Reference 4Learning to Forget: Continual Prediction with LSTMhttps://doi.org/10.1162/089976600300015015} The full gradient also travels through hidden states and gate dependencies; the output readout includes $o_t\odot(1-\tanh^2(c_t))$, which can be small. The direct path helps preserve cell-state gradients, while other paths can still vanish or explode.

GRU: Unified hidden state and streamlined gating

LSTMs maintain two state vectors ($c_t$ and $h_t$) and four distinct affine transformations per step. The gated recurrent unit (GRU) collapses the cell state and hidden state into a single vector $h_t$, and replaces the separate forget and input gates with a single update gate ($z_t$).^{[5]Reference 5Learning Phrase Representations using RNN Encoder-Decoder for Statistical Machine Translation.https://arxiv.org/abs/1406.1078}

Instead of independent keep and write decisions, the update gate forms a convex combination:

For an old state of 0.80 and a candidate of -0.20, retaining 95% ($z_t = 0.95$) automatically sets the write weight to 5%: 0.95 * 0.80 + 0.05 * (-0.20) = 0.750. The two mixture fractions always sum to 1.

The second gate is the reset gate ($r_t$), which decides how much of the prior hidden state is visible when computing the proposed candidate $\tilde{h}_t$:

When $r_t \to 0$, the candidate calculation ignores past history and acts like a standard feedforward projection of $x_t$. When $z_t \to 1$, the unit ignores the candidate and carries the past state forward unchanged.

Test the convex blend across three update gate settings:

old_state = 0.80
candidate = -0.20

for keep_gate in [0.95, 0.50, 0.05]:
    new_state = keep_gate * old_state + (1.0 - keep_gate) * candidate
    print(f"z={keep_gate:.2f} -> h={new_state:.3f}")

z=0.95 -> h=0.750
z=0.50 -> h=0.300
z=0.05 -> h=-0.150

View larger

Reset gate placement: Cho et al. versus PyTorch

Here $z$ explicitly means the fraction kept from the old state, matching Cho et al. and PyTorch. Some descriptions call the complementary fraction the "update gate"; compare each definition to see which fraction the name denotes. In Cho et al.'s candidate, the recurrent portion is $U_h(r_t\odot h_{t-1})$: reset scales the old state before its recurrent matrix multiplication.^{[5]Reference 5Learning Phrase Representations using RNN Encoder-Decoder for Statistical Machine Translation.https://arxiv.org/abs/1406.1078} $U_h$ is the recurrent block of the concatenated matrix $W_h$ above.

PyTorch's nn.GRU applies the reset gate after the linear transformation:^{[6]Reference 6GRUhttps://docs.pytorch.org/docs/2.14/generated/torch.nn.GRU.html}

PyTorch documents this deliberate change as an efficiency choice: the recurrent projections can be computed together before applying reset. Actual kernel count depends on backend, shapes, and execution mode.

Matrix multiplication and coordinate-wise scaling generally don't commute: $U_h(r\odot h)\ne r\odot(U_hh)$. Mixing coordinates can therefore change the output, although equal gates or particular inputs can make the two agree:

import numpy as np

old = np.array([0.8, -0.2])
reset = np.array([0.1, 0.9])
recurrent = np.array([[1.0, 2.0], [3.0, 4.0]])

before = recurrent @ (reset * old)
after = reset * (recurrent @ old)
print("reset before matrix:", np.round(before, 3))
print("reset after matrix: ", np.round(after, 3))
print("candidate before:   ", np.round(np.tanh(before), 3))
print("candidate after:    ", np.round(np.tanh(after), 3))
assert not np.allclose(before, after)

reset before matrix: [-0.28 -0.48]
reset after matrix:  [0.04 1.44]
candidate before:    [-0.273 -0.446]
candidate after:     [0.04  0.894]

The two candidates differ significantly (-0.273 versus 0.040). When porting weights between frameworks, always check the exact candidate formulation.

Parameter auditing and tensor geometry

Every gate and candidate vector in an RNN cell emits $d_h$ values. With concatenated input $[h_{t-1}, x_t]$ of width $d_h + d_x$, each affine transformation requires a weight matrix of shape $(d_h, d_h + d_x)$ and a bias vector of shape $(d_h,)$.

Counting parameters per cell reveals the structural scaling:

For these one-layer cells with equal widths and one bias per block, GRU has three quarters as many parameters as LSTM. Its three affine blocks also require less dense projection work than four equal-sized blocks. Wall-clock speed still depends on batching, kernels, hardware, and other operations.

Why does PyTorch's nn.GRU report 42 parameters instead of 36? In production implementations, PyTorch stores separate input-hidden ($b_{in}$) and hidden-hidden ($b_{hn}$) bias vectors for each gate:^{[6]Reference 6GRUhttps://docs.pytorch.org/docs/2.14/generated/torch.nn.GRU.html}

The extra biases don't change hidden width. They aren't all merely redundant storage: PyTorch's candidate recurrent bias is multiplied by reset, $r_t\odot b_{hn}$, so it generally can't be merged into a fixed input bias. Porting a textbook GRU requires checking both reset placement and bias placement.

For a one-layer, unidirectional nn.GRU(input_size=3, hidden_size=2, batch_first=True), an input batch tensor of shape (B, T, 3) outputs a tensor of shape (B, T, 2) and a final hidden state of shape (1, B, 2). Here $B$ counts independent sequences and $T$ counts timesteps. The leading dimension in the final hidden state represents $\text{num\_layers} \times \text{num\_directions}$, not sequence position.

Sequence padding and batch masking

In production workloads, sequences within a batch rarely have the same length. Frameworks pad shorter sequences with dummy values (typically zeros) into a uniform rectangular tensor.

A critical failure mode occurs when developers run an unmasked recurrent loop over padded rows. In an RNN, padding is still processed as an event. Even an all-zero input row triggers another state update:

This generally differs from $h_{t-1}$, though a fixed-point state can remain unchanged. The recurrence determines how a zero row affects memory. Verify the change for our chosen weights:

import math

def matvec(matrix, vector):
    return [sum(weight * value for weight, value in zip(row, vector, strict=True)) for row in matrix]

def add(left, right):
    return [a + b for a, b in zip(left, right, strict=True)]

def tanh_vec(values):
    return [math.tanh(value) for value in values]

W_xh = [
    [0.8, 0.2, 0.1],
    [0.1, 0.7, 0.4],
]
W_hh = [
    [0.5, 0.1],
    [0.0, 0.6],
]

def run(events):
    state = [0.0, 0.0]
    for event in events:
        state = tanh_vec(add(matvec(W_xh, event), matvec(W_hh, state)))
    return [round(value, 3) for value in state]

two_events = [
    [1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
]
print("two real events:      ", run(two_events))
print("plus a zero pad row:  ", run(two_events + [[0.0, 0.0, 0.0]]))
print("plus a junk pad row:  ", run(two_events + [[9.0, -9.0, 9.0]]))

two real events:       [0.494, 0.641]
plus a zero pad row:   [0.302, 0.367]
plus a junk pad row:   [1.0, -0.889]

The state shifts from [0.494, 0.641] to [0.302, 0.367] on a zero row, and to [1.0, -0.889] on junk values. A classifier that reads the post-padding state can therefore base its prediction on filler events.

View larger

PyTorch's pack_padded_sequence packages the valid prefixes and their active batch sizes so the recurrent layer skips padding on CPU or GPU.^{[7]Reference 7torch.nn.utils.rnn.pack_padded_sequencehttps://docs.pytorch.org/docs/2.14/generated/torch.nn.utils.rnn.pack_padded_sequence.html} Lengths must describe right-padded sequences; a tensor of lengths must be on CPU. enforce_sorted=False lets the function sort internally. Each sequence needs a positive length no greater than the available timesteps.

import torch
from torch import nn
from torch.nn.utils.rnn import pack_padded_sequence

torch.manual_seed(7)
gru = nn.GRU(input_size=3, hidden_size=2, batch_first=True)
lengths = torch.tensor([3, 2])
base = torch.tensor([
    [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
    [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]],
])
changed_padding = base.clone()
changed_padding[1, 2] = torch.tensor([9.0, -9.0, 9.0])

def packed_final(batch):
    packed = pack_padded_sequence(
        batch, lengths.cpu(), batch_first=True, enforce_sorted=False
    )
    _, hidden = gru(packed)
    return hidden[-1]

def unpacked_final(batch):
    _, hidden = gru(batch)
    return hidden[-1]

packed_changed = not torch.allclose(packed_final(base), packed_final(changed_padding))
unpacked_changed = not torch.allclose(unpacked_final(base), unpacked_final(changed_padding))

print("batch shape:", tuple(base.shape))
print("packed final shape:", tuple(packed_final(base).shape))
print("packed state changed by padding:", packed_changed)
print("unpacked state changed by padding:", unpacked_changed)
output, hidden = gru(base)
parameter_count = sum(parameter.numel() for parameter in gru.parameters())
print("unpacked output shape:", tuple(output.shape))
print("unpacked hidden shape:", tuple(hidden.shape))
print("GRU parameters:", parameter_count)
assert not packed_changed and unpacked_changed
assert parameter_count == 42
assert output.shape == (2, 3, 2) and hidden.shape == (1, 2, 2)

batch shape: (2, 3, 3)
packed final shape: (2, 2)
packed state changed by padding: False
unpacked state changed by padding: True
unpacked output shape: (2, 3, 2)
unpacked hidden shape: (1, 2, 2)
GRU parameters: 42

Packing ensures the second sequence halts at $t=2$. Its final representation stays unchanged when padded-slot values vary.

Sequence compression and the encoder-decoder bottleneck

The vector RNN we built is an encoder: it takes a variable-length sequence of observations and condenses them into a single fixed-width hidden vector $h_T$.

In 2014, Cho et al.^{[5]Reference 5Learning Phrase Representations using RNN Encoder-Decoder for Statistical Machine Translation.https://arxiv.org/abs/1406.1078} and Sutskever et al.^{[8]Reference 8Sequence to Sequence Learning with Neural Networks.https://arxiv.org/abs/1409.3215} used recurrent encoders and decoders to model a target sequence conditioned on a source sequence. Their interfaces differ: Cho et al. condition the decoder on a fixed context throughout, while Sutskever et al. use the encoder's state to initialize the decoder. Both communicate the source through a fixed-size representation.

Three events already share a width-2 vector. A 100-word paragraph still has to pass through that same final vector, entities and causal links included. The decoder never reads the earlier tokens. It only receives $h_T$. Width is one limit. The cell also has to finish step $t-1$ before it can compute step $t$. What changes if a token can read an earlier token directly?

The sequential bottleneck, and how attention changes it

A token still cannot read an earlier token directly. In a standard nonlinear recurrent cell, the state at step $t$ has to wait for step $t-1$. That wait is a sequential execution bottleneck.

For a standard nonlinear RNN, LSTM, or GRU, each new state depends on the previous state:

Here $s_t=h_t$ for the vanilla RNN or GRU, and $s_t=(h_t,c_t)$ for an LSTM.

The state transition at a position must wait for the previous state. A 4,096-token sequence therefore has 4,096 dependent transitions. Backend fusion and persistent algorithms determine how those transitions map to GPU kernel launches; input projections can be computed across positions in advance and batches run in parallel.^{[9]Reference 9cuDNN Recurrent Neural Network Algorithmshttps://docs.nvidia.com/deeplearning/cudnn/backend/latest/api/cudnn-adv-library.html#cudnnrnnalgo-t}

Accelerators benefit from large matrix multiplications. With a small batch or hidden width, recurrent transitions may expose too little work to use them efficiently. A batch turns the projections into matrix-matrix operations. Utilization and memory traffic depend on shapes and kernels, which may reuse recurrent weights across steps.

Self-attention breaks the serial chain

A Transformer replaces these recurrent transitions with self-attention and position-wise transformations.^{[10]Reference 10Attention Is All You Need.https://arxiv.org/abs/1706.03762} For tokens already present during training:

1. Parallel projections: Queries, keys, and values can be projected across the sequence together (often with a combined projection): $$Q = X W_Q, \quad K = X W_K, \quad V = X W_V$$
2. Direct routing: A position can attend directly to any position permitted by the mask: $$\operatorname{Attention}(Q, K, V) = \operatorname{softmax}\left(\frac{QK^{\mathsf{T}}}{\sqrt{d_k}}+M\right) V$$ A causal mask has $M_{ij}=-\infty$ for future positions $j>i$ and zero for allowed positions. Position 4,096 can receive a weighted value from Position 1 within one layer, while a recurrent model needs 4,095 relays. Attention shortens graph depth while still computing interactions across sequence positions.

Attention also needs information about position to distinguish ordering beyond what the mask supplies. The original Transformer adds positional encodings. During autoregressive generation, future tokens aren't available yet, so generating new tokens still has a sequential dependency.

Fixed-width recurrent state can still require substantial training memory. Full BPTT saves intermediate activations across the sequence; truncation or recomputation changes that cost. The table's parallel-operation counts describe a fixed-depth layer's dependency structure, while GPU kernel counts and latency depend on the implementation.

Full attention's quadratic pairwise work can be expensive for long sequences. FlashAttention tiles the computation to avoid materializing the full $T\times T$ score matrix in GPU HBM; the pairwise arithmetic remains quadratic.^{[11]Reference 11FlashAttention: Fast and Memory-Efficient Exact Attention with IO-Awarenesshttps://arxiv.org/abs/2205.14135}

Some recurrent architectures support a different execution strategy. Mamba's selective state-space recurrence has input-dependent coefficients that can be prepared before its scan, enabling an associative parallel scan during training.^{[12]Reference 12Mamba: Linear-Time Sequence Modeling with Selective State Spaceshttps://arxiv.org/abs/2312.00752} Its recurrent state size and per-token update cost are constant with respect to sequence length for a fixed model. This associative formulation applies to Mamba's scan; standard nonlinear LSTM and GRU recurrences keep their step dependencies, and a parallel scan still has scheduling and reduction costs.

Diagnosing recurrent failure modes

Three mistakes are already visible in the numbers from this chapter. Check each one before you trust a longer training run:

When diagnosing training instability, evaluate gradient norms independently from loss values. A model can show stable forward losses while its early-timestep gradients have completely vanished. Order changes a recurrent state, a repeated factor below one fades, and an LSTM or GRU gate can keep a coordinate instead of forcing it through tanh on every step.