Backpropagation Through Time: RNN Training Algorithm

Training a feedforward network is conceptually straightforward: run a forward pass through a fixed-depth computation graph, then run backpropagation to compute gradients with respect to every parameter. The graph has a clear start and end, and each layer gets visited exactly once in each direction. Every parameter is used at most once per forward pass, so its gradient is computed from a single path through the network.

Recurrent neural networks break that assumption. As we explored in the RNN Architecture chapter, an RNN processes a sequence one step at a time, feeding its hidden state forward from step to step. The same weight matrix $\mathbf{W}_{hh}$ gets applied at every timestep, which means a single parameter participates in every step of the computation. At step 1, it shapes the first hidden state. At step 2, it shapes the second. At step 50, it shapes the fiftieth. When you ask: "How should I adjust $\mathbf{W}_{hh}$ to reduce the loss?", the answer requires aggregating gradient contributions from all fifty steps simultaneously.

This is the problem that backpropagation through time (BPTT) solves. It is not a fundamentally new algorithm. BPTT is the same chain rule that powers all neural network training, applied to a computation graph that extends across time. The key insight is that by "unrolling" the RNN, you convert the looping recurrent structure into a deep but regular feedforward graph, and then you apply backpropagation to that unrolled graph exactly as you would to any other network.

The result is an elegant algorithm that correctly computes gradients for all shared parameters across all timesteps. But the unrolling operation exposes a fundamental tension: the deeper the unrolled graph (the longer the sequence), the longer the gradient paths, and the more severe the numerical problems that arise from multiplying many matrices together. BPTT's structure creates challenges that feedforward backprop never faces, chief among them the vanishing and exploding gradient problems we will examine closely in this chapter and the next.

This chapter explains how BPTT works from first principles, derives the gradient equations that flow backward through time, covers truncated BPTT for practical training on long sequences, analyzes the computational and memory costs in detail, and shows how to implement BPTT correctly in PyTorch. By the end, you will understand how the algorithm works, why it behaves the way it does on long sequences, and what practitioners have learned to do about it.

To understand BPTT, you first need a clear picture of what "unrolling" means, why it is conceptually useful, and what it costs in terms of computation and memory.

A standard RNN update at timestep $t$ is:

where:

• $\mathbf{h}_t \in \mathbb{R}^d$: the hidden state at timestep $t$, a $d$-dimensional vector that summarizes sequence history up to that point
• $\mathbf{h}_{t-1} \in \mathbb{R}^d$: the hidden state at the previous timestep, carrying information forward
• $\mathbf{W}_{hh} \in \mathbb{R}^{d \times d}$: the recurrent weight matrix, applied to the previous hidden state at every step
• $\mathbf{W}_{xh} \in \mathbb{R}^{d \times k}$: the input weight matrix, mapping inputs to the hidden space
• $\mathbf{x}_t \in \mathbb{R}^k$: the input vector at timestep $t$
• $\mathbf{b}_h \in \mathbb{R}^d$: the bias vector
• $\tanh$: the hyperbolic tangent activation applied elementwise, squashing values to $(-1, 1)$

This equation describes a computation that is deeply self-referential: to compute $\mathbf{h}_t$, you need $\mathbf{h}_{t-1}$. To compute $\mathbf{h}_{t-1}$, you need $\mathbf{h}_{t-2}$. The chain extends all the way back to $\mathbf{h}_0$, the initial hidden state (typically a zero vector). This is why RNNs are powerful: they can, in principle, let information from the very beginning of a sequence influence predictions at the very end.

When you express an RNN as a recurrence relation, you are describing a loop: the same operation executed repeatedly, each time consuming the hidden state from the previous iteration. This compact representation is convenient for describing the model, but it is not the right representation for computing gradients.

Gradient computation requires a directed acyclic graph (DAG). The chain rule can only traverse a graph that does not have cycles, because cycles would create circular dependencies in the gradient equations. An RNN's loop is inherently cyclic: $\mathbf{h}_t$ depends on $\mathbf{h}_{t-1}$ which depends on $\mathbf{h}_{t-2}$, and so on. The way to resolve this is to expand the loop along the time dimension and treat each step as a distinct node in the graph.

Unrolling makes this expansion explicit. Instead of showing the RNN as a loop that re-uses a single node, you draw out a separate copy of the computation for each timestep. For a sequence of length $T$, you get a graph with $T$ hidden state nodes connected in a chain: $\mathbf{h}_1 \to \mathbf{h}_2 \to \cdots \to \mathbf{h}_T$. Each node receives the previous hidden state and the current input and produces the next hidden state. The weight matrices $\mathbf{W}_{hh}$ and $\mathbf{W}_{xh}$ are shared across all copies, meaning they are the same physical parameters, but the intermediate values they produce at each step are all distinct and stored separately.

This is an important subtlety. Unrolling does not create new parameters; it creates new nodes in the computation graph that reference the same parameters. The implication for gradient computation is that each weight participates in $T$ computation paths, and its total gradient must account for all of them.

An output $\hat{\mathbf{y}}_t$ is computed from $\mathbf{h}_t$ at each timestep (or at the final timestep only, depending on the task):

where:

• $\mathbf{W}_{hy} \in \mathbb{R}^{c \times d}$: the output weight matrix, mapping hidden state to output space (e.g., a vocabulary in language modeling)
• $\mathbf{b}_y \in \mathbb{R}^c$: the output bias
• $c$: the output dimensionality (e.g., vocabulary size)

In many applications, an output is produced at every timestep. Language modeling is the canonical example: at each position in a sentence, the model predicts the next word. In sequence classification, outputs are produced only at the final step. In sequence-to-sequence models, outputs appear only during the decoding phase. BPTT applies in all of these cases, though the gradient structure differs slightly depending on where losses originate.

The total loss over a sequence is the sum of per-timestep losses:

where:

• $\mathcal{L}$: total sequence loss
• $\mathcal{L}_t$: loss at timestep $t$, measuring prediction error at that step
• $T$: sequence length

The additivity of the loss is important: it means the gradient of the total loss is the sum of gradients from each individual timestep. Each $\mathcal{L}_t$ only directly depends on $\hat{\mathbf{y}}_t$, which depends on $\mathbf{h}_t$. But $\mathbf{h}_t$ is a function of all earlier hidden states and therefore of all earlier inputs and all applications of $\mathbf{W}_{hh}$. So even though $\mathcal{L}_t$ does not directly reference steps before $t$, its gradient with respect to $\mathbf{W}_{hh}$ reaches all the way back through the entire history.

With the unrolled computation graph established, standard backpropagation applies. The challenge is that the loss at time $t$ depends on $\mathbf{h}_t$, which depends on $\mathbf{h}_{t-1}$, which depends on $\mathbf{h}_{t-2}$, and so on. Computing the gradient of the total loss with respect to a parameter like $\mathbf{W}_{hh}$ requires following these chains back through time and aggregating every contribution.

The gradient of the total loss with respect to $\mathbf{W}_{hh}$ decomposes into per-timestep contributions:

where:

• $\frac{\partial \mathcal{L}}{\partial \mathbf{W}_{hh}}$: total gradient for the recurrent weight matrix
• $\frac{\partial \mathcal{L}_t}{\partial \mathbf{W}_{hh}}$: the gradient contribution from the loss at timestep $t$

Each per-timestep gradient $\frac{\partial \mathcal{L}_t}{\partial \mathbf{W}_{hh}}$ must account for the fact that $\mathcal{L}_t$ depends on $\mathbf{h}_t$, which itself depends on $\mathbf{W}_{hh}$ through all previous timesteps. Applying the chain rule:

where:

• $\frac{\partial \mathcal{L}_t}{\partial \mathbf{h}_t}$: gradient of the loss at time $t$ with respect to the hidden state at time $t$, obtained directly from the output layer
• $\frac{\partial \mathbf{h}_t}{\partial \mathbf{h}_k}$: Jacobian capturing how the hidden state at time $t$ depends on the hidden state at an earlier time $k$, computed as a product of local Jacobians
• $\frac{\partial \mathbf{h}_k}{\partial \mathbf{W}_{hh}}$: how the hidden state at time $k$ directly depends on the weight matrix (via the term $\mathbf{W}_{hh}\mathbf{h}_{k-1}$ in the RNN update)

The last term, $\frac{\partial \mathbf{h}_k}{\partial \mathbf{W}_{hh}}$, captures the direct effect: at timestep $k$, the operation $\mathbf{W}_{hh}\mathbf{h}_{k-1}$ uses $\mathbf{W}_{hh}$ explicitly, so differentiating this term with respect to $\mathbf{W}_{hh}$ gives $\mathbf{h}_{k-1}^T$ (as an outer product structure). The middle term carries this local gradient signal all the way from timestep $k$ forward to timestep $t$, where the loss occurs.

The middle term $\frac{\partial \mathbf{h}_t}{\partial \mathbf{h}_k}$ is the key quantity in understanding both how BPTT works and why it is numerically challenging. It is computed via the chain rule through the sequence of hidden state transitions from $k$ to $t$:

where:

• $\prod_{j=k+1}^{t}$: a product of $t-k$ Jacobian matrices, one for each timestep between $k$ and $t$
• $\frac{\partial \mathbf{h}_j}{\partial \mathbf{h}_{j-1}}$: the Jacobian of hidden state $j$ with respect to hidden state $j-1$

Each local Jacobian captures two effects: the tanh nonlinearity applied elementwise at step $j$, and the recurrent weight matrix $\mathbf{W}_{hh}$ that mixes hidden state dimensions. Working through the differentiation:

where:

• $\text{diag}(1 - \mathbf{h}_j^2)$: a diagonal matrix whose entries are the tanh derivative at each hidden unit at timestep $j$. The tanh derivative is $\frac{d}{dz}\tanh(z) = 1 - \tanh^2(z)$, so the diagonal entries are $1 - h_{j,i}^2$ for each dimension $i$.
• $\mathbf{W}_{hh}^T$: the transpose of the recurrent weight matrix, arising because we differentiate $\mathbf{W}_{hh}\mathbf{h}_{j-1}$ with respect to $\mathbf{h}_{j-1}$

The tanh derivative $(1 - \tanh^2(z))$ is a number between 0 and 1 for any input $z$. When $z$ is large in magnitude (the tanh is saturated), the derivative is near 0. When $z$ is near zero, the derivative is near 1. This saturation behavior is the root cause of the vanishing gradient problem: at timesteps where hidden units are saturated, the local Jacobian multiplies the incoming gradient by something close to zero, annihilating it.

The product of these Jacobians across $t - k$ timesteps is what makes BPTT expensive and unstable. When $t - k$ is large (long sequences), you are multiplying many matrices together. If the eigenvalues of the resulting product are less than 1 in magnitude, the gradient vanishes exponentially. If greater than 1, it explodes exponentially. This is the vanishing and exploding gradient problem, analyzed in depth in the next chapter.

Because the same weight matrix $\mathbf{W}_{hh}$ is used at every timestep, its gradient is the sum of contributions from every timestep in the sequence. The complete gradient combines contributions from all future losses and all past hidden state paths:

where:

• The outer sum $\sum_{t=1}^{T}$ iterates over each timestep's loss
• The inner sum $\sum_{k=1}^{t}$ iterates over all past timesteps that contributed to $\mathbf{h}_t$ through $\mathbf{W}_{hh}$
• The product $\prod_{j=k+1}^{t}$ is the chain of Jacobians connecting time $k$ to time $t$

This double sum has $O(T^2)$ terms in total, though in practice the backward pass is structured to avoid redundant work through the delta recurrence described below. The same accumulation applies symmetrically to $\mathbf{W}_{xh}$, with $\frac{\partial \mathbf{h}_k}{\partial \mathbf{W}_{xh}}$ replacing the terminal term.

To see what this equation computes, think about a specific loss, say $\mathcal{L}_{10}$ at timestep 10 of a 20-step sequence. This loss propagates a gradient signal backward through steps 10, 9, 8, ..., 1. This contributes to $\mathbf{W}_{hh}$'s gradient via ten different paths of different lengths. Meanwhile, $\mathcal{L}_{20}$ propagates through all 20 steps, and $\mathcal{L}_1$ propagates through just 1 step. The total gradient for $\mathbf{W}_{hh}$ is the aggregate of all these signals, each attenuated by the product of Jacobians along its path.

In practice, BPTT works by defining a per-timestep error signal $\boldsymbol{\delta}_t$ that captures the total gradient flowing into hidden state $\mathbf{h}_t$ from all future losses. This quantity is evaluated right to left, from $t = T$ down to $t = 1$:

where:

• $\boldsymbol{\delta}_t \in \mathbb{R}^d$: the total error signal at timestep $t$, combining local and future gradient contributions
• $\frac{\partial \mathcal{L}_t}{\partial \mathbf{h}_t}$: the local gradient from the loss at time $t$ only (obtained from the output layer's backward pass)
• $\left(\frac{\partial \mathbf{h}_{t+1}}{\partial \mathbf{h}_t}\right)^T \boldsymbol{\delta}_{t+1}$: the gradient flowing back from the next timestep, multiplied by the transposed Jacobian
• $\boldsymbol{\delta}_{T+1} = \mathbf{0}$ by convention (no loss beyond the last timestep)

The base case initializes $\boldsymbol{\delta}_T = \frac{\partial \mathcal{L}_T}{\partial \mathbf{h}_T}$, the gradient from the final timestep's local loss. Then the recurrence propagates this signal backward, accumulating contributions at each step. At each intermediate timestep $t$, the error signal has two components: one from the local loss $\mathcal{L}_t$ at that step, and one from all the future error signals that have been folded together through the backward recurrence. The second component is what allows gradients from the end of the sequence to reach parameters that contributed to early-sequence hidden states.

Once all $\boldsymbol{\delta}_t$ values have been computed, the weight gradients follow directly from outer products summed across all timesteps:

where:

• $\boldsymbol{\delta}_t \odot \mathbf{g}_t$: the elementwise product of the error signal with the tanh derivative vector $\mathbf{g}_t = (1 - \mathbf{h}_t^2)$, which applies the local activation gradient before accumulating
• $\mathbf{h}_{t-1}^T$: the previous hidden state as a row vector, so the outer product gives a $d \times d$ matrix for $\mathbf{W}_{hh}$
• $\mathbf{x}_t^T$: the input vector at time $t$ as a row vector, giving a $d \times k$ matrix for $\mathbf{W}_{xh}$

This is the same outer-product structure as a feedforward network, just summed over all timesteps rather than evaluated once. Each outer product contributes a rank-1 update to the gradient matrix, and the total gradient is the sum of $T$ such rank-1 updates. This structure also reveals why increasing hidden size $d$ increases the parameter count and gradient computation quadratically: the recurrent weight matrix $\mathbf{W}_{hh}$ is $d \times d$, and computing outer products at each of the $T$ timesteps scales as $O(T \cdot d^2)$.

It is worth stepping back to appreciate why gradient flow through time is necessary rather than a mathematical artifact. Consider training an RNN to model language. To predict the word "cat" at position 20 in a sentence that begins "The fluffy orange ___", the model needs to recognize that the adjectives "fluffy" and "orange" are relevant. Gradient signal flowing backward from the prediction error at position 20 to the words at positions 2, 3, and 4 is what teaches the model to pay attention to those early adjectives. If that gradient signal cannot reach those positions (because it has been attenuated to near-zero), the model cannot adjust its weights to learn those long-range associations. BPTT is the mechanism by which learning happens; gradient decay is the mechanism by which it fails.

Full BPTT requires storing all hidden states and all activations for the entire sequence before running the backward pass. For a sequence of length $T$, this means $O(T)$ memory for activations and $O(T)$ computation for the backward pass. When $T$ is thousands or tens of thousands of tokens, full BPTT becomes impractical.

Truncated BPTT addresses this by dividing the sequence into fixed-length chunks and running a separate forward/backward pass through each chunk. Within each chunk, gradients flow freely. Across chunks, the hidden state is passed forward but gradients are stopped at the chunk boundary.

The sequence of length $T$ is divided into non-overlapping chunks of length $k$ (the truncation length). For chunk $c$ spanning timesteps $[ck+1, (c+1)k]$:

1. The forward pass starts with the final hidden state from the previous chunk, $\mathbf{h}_{ck}$, carried forward without gradient tracking
2. The forward pass runs through all $k$ steps in the chunk, accumulating hidden states and losses
3. The backward pass propagates gradients back through the $k$ steps in the chunk only
4. The weight gradients are accumulated across all chunks

The stopping of gradients at chunk boundaries is implemented using PyTorch's detach() operation, which creates a new tensor with the same values but no gradient history.

The forward-pass hidden state still carries information from earlier in the sequence. When chunk $c$ begins, the initial hidden state $\mathbf{h}_{ck}$ encodes everything the network saw in chunks 0 through $c-1$. That information is used in the forward pass to make predictions and compute losses. But the backward pass cannot trace gradients through $\mathbf{h}_{ck}$ into earlier chunks, because detach() has severed the computational graph. The result is an asymmetry: the model can use long-range information to make predictions, but cannot learn from gradient signal that would require looking beyond the current chunk boundary.

Truncated BPTT makes a practical trade-off: it gives up the ability to learn dependencies that span more than $k$ timesteps in exchange for bounded memory and predictable computation per update step. In practice, for most RNN training tasks, a truncation length between 20 and 200 captures the dependencies that matter. The exact value is a hyperparameter that depends on the sequence structure, and it interacts with the model's capacity to represent longer-range information in its hidden state.

The following visualization shows how gradient coverage changes with the truncation window relative to the full sequence length. Each bar represents the fraction of the sequence covered by one backward pass:

With $k = 20$, the 200-step sequence requires 10 separate backward passes. With $k = 100$, only 2 passes are needed. The trade-off is that larger chunks require more memory to store activations during the forward pass of each chunk.

The truncation length $k$ controls the trade-off between several competing factors:

• Memory cost: full BPTT requires $O(T)$ memory; truncated BPTT requires $O(k)$ memory per chunk
• Computational cost: each backward pass covers $k$ steps; you run $T/k$ backward passes per sequence
• Gradient quality: longer chunks capture more long-range dependencies; shorter chunks lose information about distant correlations

A common choice is $k = 35$ for language modeling, matching the LSTM training paper by Zaremba et al. (2014). For tasks with explicit long-range dependencies, larger $k$ values or gated architectures such as LSTMs help.

Choosing the truncation length is not purely a memory decision; it should reflect the structure of the data. If the sequence contains strong dependencies that span 100 tokens, a truncation window of 20 will prevent the model from ever learning them through gradient updates. In practice, you can estimate the relevant dependency range by examining what information would be needed to make accurate predictions. For next-word prediction in typical English text, dependencies rarely need to span more than a few dozen words, which explains why modest truncation windows work well for language modeling.

One underappreciated consideration is the interaction between the truncation length and the learning dynamics. Very short truncation windows produce gradients that are effectively local to a small context window. The model learns to make predictions using only nearby context, which is fast but limits what it can ultimately learn. Very long truncation windows give the model more gradient signal but also give it vanishing gradient problems of its own, since the gradient still has to propagate through $k$ matrix multiplications to reach the beginning of the chunk.

The memory comparison between full BPTT and truncated BPTT is significant. For a sequence of length $T$ with hidden size $d$ and batch size $B$:

• Full BPTT activations: $O(T \cdot d \cdot B)$ bytes (all hidden states stored)
• Truncated BPTT activations: $O(k \cdot d \cdot B)$ bytes (only current chunk's states stored)

For long sequences (e.g., $T = 1000$, $k = 35$), truncated BPTT reduces memory by roughly $T / k = 28.6\times$. The computational cost per step is similar, but the memory savings allow much larger batch sizes or hidden state sizes.

Beyond hidden states, the backward pass also requires storing intermediate activations, including the pre-activation values (the inputs to the tanh function) at each step, since the tanh derivative must be computed from these values. The total activation memory scales with the depth of the backward pass, making truncation doubly beneficial: it reduces both the hidden state storage and the intermediate activation storage.

Understanding the computational cost of BPTT matters because it directly constrains what models you can train on what hardware. The cost includes more than FLOPs; memory bandwidth, cache utilization, and the need to store intermediate activations all play important roles.

Let's be concrete. Suppose you have a sequence of length $T$, hidden state dimension $d$, input dimension $k$, and batch size $B$.

Forward pass cost (full sequence):

The dominant operation at each step is the matrix-vector product $\mathbf{W}_{hh}\mathbf{h}_{t-1}$, which costs $O(d^2)$ operations (since $\mathbf{W}_{hh}$ is $d \times d$). Across all $T$ steps and $B$ batches, the forward pass costs $O(T \cdot d^2 \cdot B)$.

Backward pass cost (full BPTT):

The backward pass has roughly the same cost as the forward pass, $O(T \cdot d^2 \cdot B)$. Computing gradient contributions to $\mathbf{W}_{hh}$ involves outer products at each step costing $O(d^2)$, and contributions to $\mathbf{W}_{xh}$ cost $O(d \cdot k)$ per step.

Total cost (full BPTT): $O(T \cdot d^2 \cdot B)$ for the complete forward and backward passes.

Total cost (truncated BPTT with chunk size $k_{\text{chunk}}$): The same asymptotic cost, but practical benefits arise from fitting the computation within GPU cache, allowing larger batches, and parallelism across chunks. Each of the $T / k_{\text{chunk}}$ chunks runs a self-contained forward/backward pass of depth $k_{\text{chunk}}$.

The hidden dimension $d$ dominates the cost in practice. Doubling the hidden size quadruples the computation for the recurrent matrix products. This is why RNNs with large hidden states are expensive to train on long sequences, and why the community moved toward transformer architectures that compute attention in $O(T^2 \cdot d)$ rather than $O(T \cdot d^2)$, a trade-off that favors larger $d$ at the cost of longer sequences.

The sequential nature of the RNN forward pass is a critical bottleneck that the computational cost analysis understates. To compute $\mathbf{h}_t$, you must first have computed $\mathbf{h}_{t-1}$. This strict temporal dependency means that no matter how many parallel processors you have, the forward pass through a sequence of length $T$ requires $T$ sequential steps. Modern GPUs are optimized for massively parallel operations, and sequential dependencies waste this parallelism. A transformer processes all tokens in parallel (though attention still has quadratic cost in sequence length), which is one key reason transformers train dramatically faster than RNNs on modern hardware despite having similar theoretical FLOPs.

For applications that require full BPTT on long sequences, gradient checkpointing offers a practical solution. Instead of storing all intermediate activations during the forward pass, you store only a subset of checkpoint activations and recompute the missing ones during the backward pass when they are needed.

The trade-off is straightforward: checkpointing reduces memory by a factor of roughly $\sqrt{T}$ (with the right checkpointing strategy) at the cost of one additional forward pass through the checkpointed segments. If the sequence is very long and memory is the bottleneck, this trade is usually worth it. PyTorch provides checkpointing support through torch.utils.checkpoint.checkpoint_sequential, which can be applied to RNN segments without modifying the gradient computation.

Gradient checkpointing is especially useful during fine-tuning of pretrained language models on long documents, where full sequences may contain thousands of tokens that would exhaust GPU memory under standard full-BPTT training. The additional compute overhead of recomputation is typically a 30-50% increase in training time, which is often acceptable when the alternative is running out of memory or reducing batch size to a single example.

Derivations and code can make BPTT seem abstract. Working through a tiny numerical example by hand is the fastest way to build intuition about what the equations compute, where the values come from, and why the backward pass has to visit timesteps in reverse order.

Let's trace through BPTT on a minimal example to make the mechanics concrete.

Consider an RNN with hidden dimension $d = 2$, processing a sequence of three steps. The recurrent weight matrix is:

After the forward pass, the hidden states are:

The error signals flowing back from the output layer (local gradients from each timestep's loss) are:

These local gradients come from the output layer: specifically, they are the derivatives of each timestep's loss $\mathcal{L}_t$ with respect to the hidden state $\mathbf{h}_t$ at that step. In practice, this derivative is computed by backpropagating through the output weight matrix $\mathbf{W}_{hy}$ and then through the softmax (if applicable). For this example, we take them as given to focus on the backward-through-time portion.

Step 1: Initialize $\boldsymbol{\delta}_3$ at the last timestep

The backward pass starts from the final timestep. The initial error signal is simply the local gradient:

The tanh derivative (elementwise) at timestep 3 is:

The modified error (tanh derivative applied before backpropagating) is:

Step 2: Propagate backward to $\boldsymbol{\delta}_2$

The recurrence carries the error signal one step further back:

Computing the matrix-vector product:

Note that $\mathbf{W}_{hh}^T$ is $\mathbf{W}_{hh}$ with rows and columns swapped, because we are differentiating $\mathbf{W}_{hh}\mathbf{h}$ with respect to $\mathbf{h}$.

Adding the local gradient:

Step 3: Accumulate the gradient contribution to $\mathbf{W}_{hh}$

At each timestep, the contribution to $\frac{\partial \mathcal{L}}{\partial \mathbf{W}_{hh}}$ is an outer product of the modified error signal and the previous hidden state. For timestep 3:

where:

• $(\boldsymbol{\delta}_3 \odot \mathbf{g}_3)$: the $2 \times 1$ modified error vector at time 3
• $\mathbf{h}_2^T$: the $1 \times 2$ previous hidden state as a row vector
• The outer product gives a $2 \times 2$ gradient matrix matching the shape of $\mathbf{W}_{hh}$

The total gradient for $\mathbf{W}_{hh}$ is the sum of these outer products across all three timesteps:

This example illustrates several important properties of BPTT. First, the backward pass must run in strict reverse order: you cannot compute $\boldsymbol{\delta}_2$ until you have computed $\boldsymbol{\delta}_3$. Second, each timestep's gradient contribution is an outer product that has exactly the shape of the weight matrix, which is why the gradient accumulation is just a simple sum. Third, the magnitude of the gradient contribution from each timestep depends both on the local error signal and on how much that signal has been attenuated by the chain of Jacobians from later steps.

It is instructive to think about what the gradient $\frac{\partial \mathcal{L}}{\partial \mathbf{W}_{hh}}$ says. The entry at position $(i, j)$ in this matrix tells you: if you increase the weight connecting hidden dimension $j$ of the previous state to hidden dimension $i$ of the current state, how does the total loss change? A large positive value means increasing that weight increases the loss. A large negative value means increasing it decreases the loss. A near-zero value means that weight has minimal effect on the loss, perhaps because the gradient signal has been attenuated by vanishing gradients along its path.

When the gradient for $\mathbf{W}_{hh}$ is very small (close to the zero matrix), gradient descent will barely update the recurrent weights. This is the vanishing gradient problem: the model cannot effectively learn from gradient signal that traveled a long path through the unrolled network. When the gradient is very large, a single gradient step will make a drastic change to the weights, destabilizing training. The next chapter examines both pathologies in rigorous detail.

PyTorch's autograd system handles full BPTT automatically when you call .backward() on a loss computed over a sequence. The computational graph is built during the forward pass and traversed during the backward pass. All you need to do is ensure you are not breaking the graph with unnecessary detach() calls.

Let's build a simple character-level language model to demonstrate BPTT:

The model has a standard RNN layer followed by a linear output layer. The number of parameters depends only on the vocabulary size and hidden dimension, not the sequence length, because the weights are shared across all timesteps. This weight sharing is what gives RNNs their efficiency in parameters, but it also creates the gradient challenges we have been examining.

The nn.RNN module in PyTorch implements the vanilla RNN we have been deriving. Internally, it runs the recurrence $\mathbf{h}_t = \tanh(\mathbf{W}_{xh}\mathbf{x}_t + \mathbf{W}_{hh}\mathbf{h}_{t-1} + \mathbf{b})$ at each timestep and accumulates the hidden states. When .backward() is called on the final loss, PyTorch automatically unrolls this computation and computes BPTT gradients through the accumulated graph.

Full BPTT is the default behavior in PyTorch. When you run a forward pass through the entire sequence and then call .backward(), PyTorch automatically propagates gradients through all timesteps:

The gradient norms show how much signal reaches the weight matrices from backpropagating through all seq_length timesteps. For short sequences, these gradients are healthy. For very long sequences (hundreds of steps), they often collapse toward zero, illustrating the vanishing gradient problem.

The key to truncated BPTT is the detach() operation. When you detach the hidden state before passing it to the next chunk, you pass the values forward (preserving memory across chunks) without passing the gradient history (stopping gradient flow at chunk boundaries):

The detach() call is the single most important line in truncated BPTT. Without it, gradients would flow back through the entire sequence history, turning the training into full BPTT but with multiple optimizer steps mid-sequence. With it, each chunk is a self-contained backward pass.

One common mistake when implementing truncated BPTT is calling optimizer.zero_grad() before detaching the hidden state. The order matters: you should detach the hidden state from the previous chunk before zeroing gradients, not after. The code above handles this correctly, but it is easy to get wrong when refactoring.

A useful diagnostic is to track gradient norms over training. Healthy gradients have norms that do not consistently collapse to zero or explode to large values:

The gradient norms reveal a clear pattern: longer sequences produce smaller gradients for the recurrent weight matrix. This is the vanishing gradient problem in action. The gradient signal must travel through more matrix multiplications to reach the earliest timesteps, and each multiplication by $\mathbf{W}_{hh}$ (combined with the tanh derivative) tends to shrink the signal. This degradation motivates the architectural solutions covered in the following chapters on LSTMs and GRUs.

Let's compare the gradient behavior of full BPTT and truncated BPTT on the same sequence:

Truncated BPTT tends to produce larger gradient norms than full BPTT on long sequences, precisely because it avoids the extended chain of matrix products that attenuates the signal in full BPTT. This counter-intuitive result is why truncated BPTT can sometimes train faster: the gradients are noisier (they do not account for long-range dependencies), but they are not vanishing.

The trade-off here is between gradient fidelity and gradient magnitude. Full BPTT computes the exact gradient of the loss with respect to $\mathbf{W}_{hh}$, but on long sequences this exact gradient may be so small that learning is effectively stalled. Truncated BPTT computes an approximate gradient that ignores long-range contributions, but the resulting gradient is often large enough to update the weights. For tasks where long-range dependencies matter, this approximation can harm final model quality; for tasks dominated by local context, it often works as well or better.

The key parameters for BPTT training are:

• chunk_size (truncation window): Controls how many timesteps are unrolled for each backward pass. Typical values are 20-200. Larger values capture longer dependencies but use more memory.
• hidden_size: The dimensionality of the hidden state. Larger hidden states give the model more capacity but increase both forward and backward pass costs as $O(d^2)$.
• gradient clipping threshold: As covered in the Gradient Clipping chapter, clipping by global norm (typical threshold 1.0-5.0) prevents exploding gradients from destabilizing training.

Let's visualize the memory tradeoff concretely:

The chart makes the advantage of truncated BPTT clear. For a sequence of length 1000 with hidden dimension 512 and batch size 32, full BPTT requires roughly 64 MB of activations just for the hidden states. Truncated BPTT with chunk size 50 reduces this to about 3.2 MB per chunk, a $1000/50 = 20\times$ reduction. This is before accounting for the even larger memory needed to store the pre-activation values (inputs to the tanh function), which are also needed for computing the tanh derivatives during the backward pass.

The quality of gradient signal reaching early timesteps determines whether the model can learn long-range dependencies. We can track how gradient norms decay as we move further back in time using PyTorch's retain_grad():

The gradient decay plot reveals why vanilla RNNs struggle with long-range dependencies. Early timesteps receive only a tiny fraction of the gradient signal that the later timesteps receive. This asymmetry means the model updates its weights based primarily on recent history, making it difficult to learn patterns that require looking far back in time.

This plot also reveals something subtle about the nature of learning in RNNs: even if a model could in principle encode long-range dependencies in its hidden state (and a sufficiently large hidden state can represent very complex functions), training with BPTT cannot teach it to do so if the gradients decay before reaching the relevant early timesteps. The representational capacity of the model and the trainability of that capacity are two different things. A vanilla RNN might theoretically be able to remember that a sentence began with "The fluffy orange", but BPTT cannot reliably teach it to do so if those words appeared many steps ago.

The LSTM and GRU architectures, introduced in the following chapters, address this limitation through gating mechanisms that provide more direct gradient paths through time. By using learned gates to selectively write to and read from a cell state, LSTMs can maintain gradients that are closer to 1.0 in magnitude over long distances, allowing gradient signal to travel dozens or hundreds of steps without severe attenuation.

Backpropagation through time was not invented all at once. The algorithm's history reflects the broader history of learning in recurrent networks and reveals why certain design choices were made.

The conceptual framework for training recurrent networks via gradient descent was laid out by David Rumelhart, Geoffrey Hinton, and Ronald Williams in their 1986 paper on backpropagation in feedforward networks. The extension to recurrent networks was recognized soon after: if you unroll a recurrent network across time, it becomes a deep feedforward network, and standard backpropagation applies. Paul Werbos described this approach formally in the late 1980s, and it was further developed and popularized by Williams and Zipser in 1989 and 1995. The term "backpropagation through time" became standard in the early 1990s.

The vanishing gradient problem was identified as the key obstacle by Sepp Hochreiter in his 1991 German-language diploma thesis, a remarkable piece of work that performed a systematic analysis of why standard RNNs failed to learn long-range dependencies. Hochreiter showed mathematically that gradient magnitudes decrease exponentially with distance from the output, making it essentially impossible to train vanilla RNNs on tasks requiring memory beyond about 10-20 timesteps. This analysis directly motivated the design of the Long Short-Term Memory (LSTM) network, published by Hochreiter and Schmidhuber in 1997.

The LSTM's central innovation was the constant error carousel: a cell state that gradients can flow through without passing through the tanh nonlinearity, provided the forget gate remains open. This design creates paths where the product of Jacobians can be kept near 1.0, solving the vanishing gradient problem not by changing the training algorithm but by changing the network architecture.

It is worth appreciating what this history reveals about the relationship between algorithms and architectures in deep learning. BPTT is the correct algorithm for training any differentiable recurrent system. The problem is not the algorithm; the problem is that vanilla RNNs have gradient flow properties that make the algorithm fail to learn long-range dependencies. The solution was not to fix BPTT but to design architectures whose gradient flow properties are compatible with BPTT. This same pattern appears throughout deep learning: batch normalization was designed partly to improve gradient flow in very deep feedforward networks, and residual connections in ResNets create shortcut paths that keep gradients healthy at depth.

Once you understand BPTT for a vanilla RNN, it is informative to consider how the same algorithm behaves for LSTMs and GRUs. These architectures do not require a different training algorithm; they use BPTT exactly as described above. What they change is the structure of the forward computation, which in turn changes the structure of the Jacobians that appear in the backward pass.

An LSTM maintains two state vectors at each timestep: the hidden state $\mathbf{h}_t$ (sometimes called the output state) and the cell state $\mathbf{c}_t$ (the memory). The key gradient path in an LSTM runs through the cell state, not the hidden state. The cell state update equation is:

where $\mathbf{f}_t$ is the forget gate and $\mathbf{i}_t$ is the input gate. The gradient of $\mathbf{c}_t$ with respect to $\mathbf{c}_{t-1}$ is simply $\text{diag}(\mathbf{f}_t)$: a diagonal matrix containing the forget gate values. Because forget gate values are squashed by a sigmoid into $(0, 1)$, this Jacobian is always diagonal and always has entries between 0 and 1. The forget gate is learned and can be set close to 1.0 for dimensions that need long-term memory. When the forget gate is near 1, gradients flow backward through the cell state with almost no attenuation.

This is a fundamentally different gradient regime from the vanilla RNN, where the recurrent Jacobian $\text{diag}(1 - \mathbf{h}_j^2) \cdot \mathbf{W}_{hh}^T$ involves both the tanh derivative (always in $(0, 1)$) and a full matrix $\mathbf{W}_{hh}^T$ that can have large or small spectral properties depending on training. The LSTM essentially replaces that matrix product with a learned diagonal gate, giving the model explicit control over gradient flow.

A Gated Recurrent Unit (GRU) achieves similar benefits with a simpler structure, maintaining only one state vector. The update rule involves a reset gate $\mathbf{r}_t$ and an update gate $\mathbf{z}_t$:

The gradient of $\mathbf{h}_t$ with respect to $\mathbf{h}_{t-1}$ has a component proportional to $\text{diag}(1 - \mathbf{z}_t)$. When the update gate $\mathbf{z}_t$ is close to 0, this component is close to 1, passing gradients backward with minimal attenuation. The GRU achieves a similar constant-error-carousel effect to the LSTM through the interpolation structure of the update gate.

Both LSTMs and GRUs still use BPTT. The benefit they provide is architectural: by creating gradient paths that do not inevitably pass through full matrix products and saturated nonlinearities, they allow BPTT to compute meaningful gradients even at large distances from the output. The training algorithm is unchanged; the model is designed to be more amenable to training.

Training RNNs with BPTT involves several common pitfalls that can cause silently incorrect or very slow training.

When iterating over a dataset multiple times, a common mistake is to carry the hidden state from the end of one epoch into the beginning of the next without detaching it. If the hidden state from epoch $n$'s last batch is passed directly into epoch $n+1$'s first batch without calling .detach(), the computational graph grows across batches. By the time you call .backward() on the epoch $n+1$ loss, PyTorch will try to backpropagate through the hidden state all the way back to wherever the graph started in epoch $n$. This can lead to memory errors (the graph is unexpectedly large) or incorrect gradients (the graph spans multiple batches of unrelated data).

The safe practice is to always initialize the hidden state fresh at the start of each independent sequence in your dataset, or to call .detach() on any hidden state that is being carried across batch boundaries.

When processing a dataset of multiple independent sequences (documents, sentences, examples), you must reset the hidden state between sequences. Each sequence in the dataset is independent; the hidden state from processing document A should not be used as the initial state for processing document B. Failing to reset leads to a subtle bug where the model's predictions at the start of each sequence are influenced by the content of the previous sequence.

This bug is particularly insidious because the training loss may still decrease: the model might learn to use the leaked context from the previous sequence for certain tasks. But the model's behavior at inference time will be wrong if sequences are processed in a different order or in isolation, because the expected hidden state contamination will not be present.

While vanishing gradients are the chronic problem, exploding gradients are the acute one. A single exploding gradient event can corrupt a model's weights catastrophically, requiring restarting training from a checkpoint. Gradient clipping, covered in detail in the Gradient Clipping chapter, is the standard remedy: before applying the optimizer step, you rescale the global gradient norm to a maximum threshold value if it exceeds that threshold.

In PyTorch, gradient clipping is one line:

This should be called after loss.backward() and before optimizer.step(). The typical threshold for RNN training is 1.0 to 5.0, with 1.0 being a conservative and widely-used default.

The learning rate for RNN training typically needs to be smaller than for feedforward networks, because the gradient scale varies significantly with sequence length. When you increase the sequence length (or chunk size), the gradient norms change, and a learning rate that worked well for short sequences may be too large or too small for longer ones. Using adaptive optimizers like Adam helps, because Adam normalizes gradients per-parameter, but even with Adam, careful learning rate tuning is often necessary when changing the sequence length.

A related consideration is learning rate warmup: starting with a very small learning rate and gradually increasing it to the target value over the first few thousand steps. Warmup helps because early in training, the model's hidden states are poorly initialized and the gradient surface is unstable. Large early gradient steps can send weights to poor regions of parameter space from which recovery is slow.

BPTT is the correct training algorithm for RNNs, but it comes with practical constraints that shaped how RNN-based systems were designed for many years.

The most significant limitation is the cost of backpropagating through long sequences. Every additional timestep adds a matrix multiplication to the gradient path, and each multiplication tends to shrink the gradient (if the spectral radius of $\mathbf{W}_{hh}$ is less than 1) or grow it (if greater than 1). This sensitivity to the eigenvalue spectrum of the recurrent weight matrix is mathematically fundamental. BPTT faithfully computes the true gradient; the problem is that the true gradient is often very small or very large for long sequences.

Truncated BPTT is the standard practical solution, but it introduces its own issues. The truncation window acts as a hard limit on the dependencies the gradient can capture. A model trained with chunk size 35 simply cannot learn from gradient signal originating more than 35 steps in the past, no matter how informative that signal would be. This is not a memory issue; it is a deliberate algorithmic choice that trades gradient completeness for computational tractability. For some tasks, such as parsing long documents or understanding multi-paragraph narratives, this limitation is severe.

Memory also presents a real barrier for full BPTT on modern hardware. During the backward pass, PyTorch needs to retain all intermediate activations from the forward pass to compute gradients. For a sequence of length 1000 with a hidden state of dimension 512 and batch size 32, this amounts to hundreds of megabytes for a single sequence, before accounting for optimizer state. Gradient checkpointing offers a middle ground: instead of storing all activations, you recompute them during the backward pass, trading compute for memory. PyTorch supports this via torch.utils.checkpoint.checkpoint_sequential.

The sequential nature of RNN computation is arguably the deeper limitation. Because $\mathbf{h}_t$ depends on $\mathbf{h}_{t-1}$, the forward pass cannot be parallelized across timesteps. For a sequence of length 1000, you must execute 1000 sequential matrix operations, regardless of available hardware parallelism. Modern GPUs can process the matrix operations for a batch of sequences in parallel, but they cannot process different timesteps of the same sequence in parallel. This makes RNNs fundamentally slow to train compared to transformers, which can process all positions in a sequence simultaneously through attention.

Despite these limitations, BPTT remains the conceptual foundation for training all sequence models. The gated architectures (LSTMs, GRUs) covered in the next chapters do not replace BPTT; they improve the gradient flow properties of the network so that BPTT works better. Even modern transformer architectures, which sidestep recurrence entirely through attention, are trained using the same chain-rule-based backpropagation, just through a different computational graph. Understanding BPTT deeply is the foundation for understanding why LSTMs work, why transformers are preferred for long sequences, and what any future sequence model needs to address.

One subtle practical issue deserves mention: the handling of the initial hidden state. Typical practice initializes $\mathbf{h}_0 = \mathbf{0}$. In truncated BPTT, the hidden state at the start of each chunk is the detached final state from the previous chunk. This detachment means the model cannot update weights to make better use of the carried-forward state based on losses from future chunks. In practice this is acceptable, but it does mean the first few steps of each chunk may be slightly miscalibrated as the hidden state re-accumulates context.

An alternative approach sometimes used in practice is to learn the initial hidden state as a parameter. By treating $\mathbf{h}_0$ as a trainable vector, the model can learn a good starting state that minimizes the warmup problem at the beginning of each sequence or chunk. This is a small but non-negligible improvement for tasks where the first few tokens of a sequence are important.

Backpropagation through time trains recurrent neural networks by unrolling the recurrent computation graph into a feedforward graph and applying standard backpropagation through it. The core ideas are:

• Unrolling converts recurrence to depth: An RNN over a sequence of length $T$ becomes a feedforward network with $T$ layers for the purpose of gradient computation, with all layers sharing the same weight matrices.
• Gradient accumulation: Because weights are shared across all timesteps, the gradient for any parameter is the sum of contributions from every timestep where that parameter participated.
• The backward error signal: The $\boldsymbol{\delta}_t$ recurrence propagates gradient information from later timesteps to earlier ones, accumulating contributions at each step.
• Truncated BPTT: Dividing sequences into chunks and using detach() at chunk boundaries bounds memory and computation while preserving the ability to carry forward state information.
• Gradient decay: Vanilla RNNs suffer from vanishing gradients when trained with BPTT on long sequences because the product of Jacobians across many timesteps tends toward zero. This motivates the gated architectures in the following chapters.
• Sequential computation: The strict timestep-by-timestep dependency in RNN forward passes limits hardware utilization and makes RNNs slower to train than architectures that can parallelize across positions, which was a key motivation for the development of the transformer.

The next chapter examines the vanishing gradient problem in detail, analyzing exactly why the product of Jacobians causes gradient signal to decay and what that means for what vanilla RNNs can and cannot learn.

Ready to test your understanding? Take this quick quiz to reinforce what you've learned about backpropagation through time.