Backpropagation Through Time › Deriving BPTT Gradients

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:

$$\frac{\partial \mathcal{L}}{\partial \mathbf{W}_{hh}} = \sum_{t=1}^{T} \frac{\partial \mathcal{L}_t}{\partial \mathbf{W}_{hh}}$$

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:

$$\frac{\partial \mathcal{L}_t}{\partial \mathbf{W}_{hh}} = \sum_{k=1}^{t} \frac{\partial \mathcal{L}_t}{\partial \mathbf{h}_t} \cdot \frac{\partial \mathbf{h}_t}{\partial \mathbf{h}_k} \cdot \frac{\partial \mathbf{h}_k}{\partial \mathbf{W}_{hh}}$$

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$:

$$\frac{\partial \mathbf{h}_t}{\partial \mathbf{h}_k} = \prod_{j=k+1}^{t} \frac{\partial \mathbf{h}_j}{\partial \mathbf{h}_{j-1}}$$

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:

$$\frac{\partial \mathbf{h}_j}{\partial \mathbf{h}_{j-1}} = \text{diag}(1 - \mathbf{h}_j^2) \cdot \mathbf{W}_{hh}^T$$

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:

$$\frac{\partial \mathcal{L}}{\partial \mathbf{W}_{hh}} = \sum_{t=1}^{T} \sum_{k=1}^{t} \frac{\partial \mathcal{L}_t}{\partial \mathbf{h}_t} \cdot \left(\prod_{j=k+1}^{t} \frac{\partial \mathbf{h}_j}{\partial \mathbf{h}_{j-1}}\right) \cdot \frac{\partial \mathbf{h}_k}{\partial \mathbf{W}_{hh}}$$

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$:

$$\boldsymbol{\delta}_t = \frac{\partial \mathcal{L}}{\partial \mathbf{h}_t} = \frac{\partial \mathcal{L}_t}{\partial \mathbf{h}_t} + \left(\frac{\partial \mathbf{h}_{t+1}}{\partial \mathbf{h}_t}\right)^T \boldsymbol{\delta}_{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:

$$\frac{\partial \mathcal{L}}{\partial \mathbf{W}_{hh}} = \sum_{t=1}^{T} (\boldsymbol{\delta}_t \odot \mathbf{g}_t) \cdot \mathbf{h}_{t-1}^T$$ $$\frac{\partial \mathcal{L}}{\partial \mathbf{W}_{xh}} = \sum_{t=1}^{T} (\boldsymbol{\delta}_t \odot \mathbf{g}_t) \cdot \mathbf{x}_t^T$$

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.