BackPropagation Through Time (BPTT)

Backpropagation Through Time (BPTT) is the application of the canonical backpropagation algorithm to a Recurrent Neural Network (RNN) after the recurrent computation has been explicitly unrolled over a finite number of time steps.

Prerequisite: backpropagation as automatic differentiation

For a deeper understanding of how backpropagation works at the level of computational graphs, local derivatives, gradient accumulation, and topological ordering, see the Micrograd series. BPTT uses the same mechanism, with the additional structure introduced by the temporal unrolling of a recurrent computation.

In essence, the recurrent layer is unfolded along the temporal axis and treated as a deep feedforward computational graph whose temporal copies share the same parameters. Backpropagation is then applied to this unfolded graph. The number of time steps over which gradients are propagated is therefore a relevant modeling choice: in full BPTT it coincides with the sequence length, whereas in truncated BPTT it becomes an explicit truncation horizon.

Question

How should the hidden state before the beginning of the sequence be initialized, namely the initial state ${\mathbf{h}}_0$?

Note

At the first time step ($t=1$), no previous hidden state is available. It is therefore customary to set

$${\mathbf{h}}_0=\vec{0}$$

that is, a zero vector. This choice is analogous to zero-padding in convolutional neural networks: in CNNs, zeros are added at the spatial boundary; here, mutatis mutandis, the same idea is applied at the temporal boundary of the sequence.

Conceptually, BPTT is equivalent to backpropagating the loss through a deep feedforward graph made of $T$ temporal copies of the same recurrent cell, from $t=1$ to $t=T$. Unlike a standard MLP with independent layers, however, these temporal copies all share the same parameters.

---

In depth derivation

For an explicit derivation of Backpropagation Through Time (BPTT), the RNN is represented in its unfolded form along the temporal axis, as shown in the figure above. The diagram highlights both the forward computation, from inputs to local losses, and the backward propagation of gradients through time.

The compact RNN block shown on the left of the figure should be read as a schematic view of a single recurrent cell. The hidden-to-output matrix ${\mathbf{W}}_{yh}$ maps the hidden state to the logit vector ${\mathbf{o}}_t$; the prediction ${\hat{\mathbf{y}}}_t$ is obtained only after applying the output activation, typically a softmax in multiclass classification.

Component / Symbol	Description
$C$	Number of classes in the classification task.
${\mathbf{W}}_{xh}\in {\mathbb{R}}^{n_{\text{neurons}}\times n_{\text{inputs}}}$, ${\mathbf{W}}_{hh}\in {\mathbb{R}}^{n_{\text{neurons}}\times n_{\text{neurons}}}$, ${\mathbf{W}}_{yh}\in {\mathbb{R}}^{C\times n_{\text{neurons}}}$	Weight matrices. Although they are drawn at each time step, they are shared across time. This parameter sharing reduces the total number of parameters and expresses the recurrent nature of the model.
${\mathbf{x}}_t\in {\mathbb{R}}^{n_{\text{inputs}}\times 1}$	Input column vector at time step $t$.
${\mathbf{h}}_t\in {\mathbb{R}}^{n_{\text{neurons}}\times 1}$	Hidden state at time step $t$.
${\mathbf{o}}_t\in {\mathbb{R}}^{C\times 1}$	Raw output vector, or logit vector, before the softmax.
${\hat{\mathbf{y}}}_t\in {\mathbb{R}}^{C\times 1}$	Predicted probability vector. It is obtained by applying the softmax function to the logits ${\mathbf{o}}_t$. For each component $i$, ${\hat{\mathbf{y}}}_t^{(i)}=\frac{e^{{\mathbf{o}}_t^{(i)}}}{{\sum }_{j=1}^C{e}^{{\mathbf{o}}_t^{(j)}}}$ The softmax ensures that ${\hat{\mathbf{y}}}_t$ is a probability vector, which makes it suitable for multiclass classification and compatible with cross-entropy loss.
$L_t$	Local loss at time step $t$.
Propagation Phases	The computation has two phases. In the forward pass, shown by gray arrows, the input ${\mathbf{x}}_t$ and the previous hidden state ${\mathbf{h}}_{t-1}$ are combined to produce the new hidden state ${\mathbf{h}}_t$ and the output. In the backward pass, shown by red arrows, the gradient of the total loss is propagated backward through the unfolded time steps and accumulated for the shared parameters.

RNN as a Dynamical System

A standard RNN can be viewed as a dynamical system defined by a state transition function, which updates the internal state of the network, and an output function, which maps the hidden state to the model output.

Single-Sequence Formulation

For a single sequence, the system can be written in explicit form as

$$\begin{array}{rl}{\mathbf{h}}_t& =f_h({\mathbf{x}}_t,{\mathbf{h}}_{t-1})\\ {\hat{\mathbf{y}}}_t& =f_o({\mathbf{h}}_t).\end{array}$$

For a conventional RNN, these functions are usually instantiated as

$$\boxed{\begin{array}{rl}{\mathbf{a}}_t& ={\mathbf{W}}_{xh}{\mathbf{x}}_t+{\mathbf{W}}_{hh}{\mathbf{h}}_{t-1}+{\mathbf{b}}_h\\ {\mathbf{h}}_t& =f_h({\mathbf{x}}_t,{\mathbf{h}}_{t-1})={\varphi }_h({\mathbf{a}}_t)\\ {\mathbf{o}}_t& ={\mathbf{W}}_{yh}{\mathbf{h}}_t+{\mathbf{b}}_y\\ {\hat{\mathbf{y}}}_t& =f_o({\mathbf{h}}_t)={\varphi }_o({\mathbf{o}}_t).\end{array}}$$

The vector ${\mathbf{a}}_t={\mathbf{W}}_{xh}{\mathbf{x}}_t+{\mathbf{W}}_{hh}{\mathbf{h}}_{t-1}+{\mathbf{b}}_h$ is the hidden pre-activation: the recurrent counterpart of the pre-activation $z$ in an MLP, with the hidden state obtained as its element-wise nonlinearity, ${\mathbf{h}}_t={\varphi }_h({\mathbf{a}}_t)$. It is named explicitly because the backward pass differentiates the loss with respect to it.

The functions ${\varphi }_h(\cdot )$ and ${\varphi }_o(\cdot )$ are nonlinear functions. In particular:

• ${\varphi }_h(\cdot )$ is an element-wise hidden activation function, commonly $\tanh (\cdot )$.
• ${\varphi }_o(\cdot )$ is the output activation. In multiclass classification, it is typically the softmax, because its output can be interpreted as a probability distribution over classes.

Symbol	Dimension	Description
${\mathbf{x}}_t$	${\mathbb{R}}^{n_{\text{inputs}}\times 1}$	Input vector at time $t$.
${\mathbf{h}}_{t-1}$, ${\mathbf{h}}_t$	${\mathbb{R}}^{n_{\text{neurons}}\times 1}$	Previous and current hidden states.
${\mathbf{a}}_t$	${\mathbb{R}}^{n_{\text{neurons}}\times 1}$	Hidden pre-activation, ${\mathbf{a}}_t={\mathbf{W}}_{xh}{\mathbf{x}}_t+{\mathbf{W}}_{hh}{\mathbf{h}}_{t-1}+{\mathbf{b}}_h$, with ${\mathbf{h}}_t={\varphi }_h({\mathbf{a}}_t)$.
${\mathbf{o}}_t$	${\mathbb{R}}^{C\times 1}$	Logit vector before the softmax.
${\hat{\mathbf{y}}}_t$	${\mathbb{R}}^{C\times 1}$	Predicted probability vector after the softmax.
${\mathbf{W}}_{xh}$	${\mathbb{R}}^{n_{\text{neurons}}\times n_{\text{inputs}}}$	Input-to-hidden weight matrix.
${\mathbf{W}}_{hh}$	${\mathbb{R}}^{n_{\text{neurons}}\times n_{\text{neurons}}}$	Hidden-to-hidden recurrent weight matrix.
${\mathbf{W}}_{yh}$	${\mathbb{R}}^{C\times n_{\text{neurons}}}$	Hidden-to-output/logit weight matrix.
${\mathbf{b}}_h$	${\mathbb{R}}^{n_{\text{neurons}}\times 1}$	Hidden-state bias vector.
${\mathbf{b}}_y$	${\mathbb{R}}^{C\times 1}$	Output/logit bias vector.
${\varphi }_h(\cdot )$	—	Hidden activation function, for instance $\tanh$.
${\varphi }_o(\cdot )$	—	Output activation function, typically softmax in classification.
$L_t$	$\mathbb{R}$	Local scalar loss at time $t$.

Mini-Batch Formulation

Important

For computational efficiency, training is usually performed on mini-batches of sequences. In this setting, the input vector of a single time step, ${\mathbf{x}}_t\in {\mathbb{R}}^{n_{\text{inputs}}}$, is generalized to an input matrix ${\mathbf{X}}_t\in {\mathbb{R}}^{B\times n_{\text{inputs}}}$, where $B$ is the mini-batch size.

Similarly, the previous hidden states are collected into a matrix ${\mathbf{H}}_{t-1}\in {\mathbb{R}}^{B\times n_{\text{neurons}}}$, where each row contains the hidden state associated with one sequence in the batch.

The recurrent layer can then be evaluated for the whole mini-batch as

$${\mathbf{H}}_t={\varphi }_h({\mathbf{X}}_t{\mathbf{W}}_{xh}^{\top }+{\mathbf{H}}_{t-1}{\mathbf{W}}_{hh}^{\top }+{\mathbf{b}}_h),{\mathbf{O}}_t={\mathbf{H}}_t{\mathbf{W}}_{yh}^{\top }+{\mathbf{b}}_y,{\hat{\mathbf{Y}}}_t=\mathrm{softmax}({\mathbf{O}}_t).$$

Here:

• ${\mathbf{X}}_t\in {\mathbb{R}}^{B\times n_{\text{inputs}}}$ is the mini-batch input matrix at time $t$.
• ${\mathbf{H}}_{t-1}\in {\mathbb{R}}^{B\times n_{\text{neurons}}}$ is the matrix of previous hidden states.
• ${\mathbf{H}}_t\in {\mathbb{R}}^{B\times n_{\text{neurons}}}$ contains the hidden states computed at time $t$.
• ${\mathbf{O}}_t\in {\mathbb{R}}^{B\times C}$ contains the logits for each example in the batch.
• ${\hat{\mathbf{Y}}}_t\in {\mathbb{R}}^{B\times C}$ contains the predicted probabilities after the softmax.

The main tensor dimensions are summarized below.

Symbol	Dimension	Description
${\mathbf{X}}_t$	${\mathbb{R}}^{B\times n_{\text{inputs}}}$	Mini-batch input at time $t$.
${\mathbf{H}}_{t-1}$	${\mathbb{R}}^{B\times n_{\text{neurons}}}$	Mini-batch hidden state at time $t-1$.
${\mathbf{W}}_{xh}$	${\mathbb{R}}^{n_{\text{neurons}}\times n_{\text{inputs}}}$	Input-to-hidden weight matrix.
${\mathbf{W}}_{hh}$	${\mathbb{R}}^{n_{\text{neurons}}\times n_{\text{neurons}}}$	Hidden-to-hidden recurrent weight matrix.
${\mathbf{W}}_{yh}$	${\mathbb{R}}^{C\times n_{\text{neurons}}}$	Hidden-to-output/logit weight matrix.
${\mathbf{b}}_h$	${\mathbb{R}}^{n_{\text{neurons}}\times 1}$, broadcast to ${\mathbb{R}}^{B\times n_{\text{neurons}}}$	Hidden-state bias.
${\mathbf{b}}_y$	${\mathbb{R}}^{C\times 1}$, broadcast to ${\mathbb{R}}^{B\times C}$	Output/logit bias.
${\hat{\mathbf{Y}}}_t$	${\mathbb{R}}^{B\times C}$	Mini-batch predicted probabilities at time $t$.

Why the Transposed Weight Matrices?

The transposes come from using two compatible conventions at the same time: column vectors for the mathematical derivation, and row-wise mini-batches for implementation.

• Single-sequence view. The input ${\mathbf{x}}_t$ is a column vector of shape $[n_{\text{inputs}}\times 1]$. The matrix ${\mathbf{W}}_{xh}$ has shape $[n_{\text{neurons}}\times n_{\text{inputs}}]$, so the product ${\mathbf{W}}_{xh}{\mathbf{x}}_t$ is dimensionally valid.

• Mini-batch view. The input ${\mathbf{X}}_t$ stores one example per row, so it has shape $[B\times n_{\text{inputs}}]$. To apply the same linear map to every row, the weight matrix must appear on the right with shape $[n_{\text{inputs}}\times n_{\text{neurons}}]$.

Since ${\mathbf{W}}_{xh}$ has been defined in the column-vector convention as $[n_{\text{neurons}}\times n_{\text{inputs}}]$, the mini-batch computation uses its transpose:

$${\mathbf{X}}_t{\mathbf{W}}_{xh}^{\top }\in {\mathbb{R}}^{B\times n_{\text{neurons}}}.$$

The same reasoning applies to the recurrent term:

$${\mathbf{H}}_{t-1}{\mathbf{W}}_{hh}^{\top }\in {\mathbb{R}}^{B\times n_{\text{neurons}}}.$$

This is also the convention used by PyTorch. A layer such as nn.Linear(in_features, out_features) stores weight with shape [out_features, in_features], and the forward computation is performed as input @ weight.T + bias. The same pattern appears in nn.RNN, where the recurrent update is written in batched form using terms such as x_t @ W_ih.T and h_{t-1} @ W_hh.T.

Compact Notation for Input and Previous Hidden State

In optimized implementations, it is common to concatenate the input at time $t$ and the previous hidden state into a single matrix:

$$[{\mathbf{X}}_t{\mathbf{H}}_{t-1}]\in {\mathbb{R}}^{B\times (n_{\text{inputs}}+n_{\text{neurons}})}.$$

This matrix is paired with an extended weight matrix

$${\mathbf{W}}_{\text{cat}}=\left[\begin{array}{c}{\mathbf{W}}_{xh}^{\top }\\ {\mathbf{W}}_{hh}^{\top }\end{array}\right]\in {\mathbb{R}}^{(n_{\text{inputs}}+n_{\text{neurons}})\times n_{\text{neurons}}}.$$

The recurrent transition can then be written compactly as

$${\mathbf{H}}_t={\varphi }_h\left([{\mathbf{X}}_t{\mathbf{H}}_{t-1}]{\mathbf{W}}_{\text{cat}}+{\mathbf{b}}_h\right).$$

This notation simplifies the implementation and reduces the recurrent transition to a single matrix multiplication.

Note

At the first time step ($t=1$), the previous hidden state is initialized as

$${\mathbf{H}}_0=0,$$

where $0\in {\mathbb{R}}^{B\times n_{\text{neurons}}}$ is a matrix of zeros.

Important

This structure makes it possible to derive gradients systematically with respect to the shared parameters ${\mathbf{W}}_{xh}$, ${\mathbf{W}}_{hh}$, ${\mathbf{W}}_{yh}$, ${\mathbf{b}}_h$, and ${\mathbf{b}}_y$. This is the core mechanism of BPTT.

---

Gradient Derivation: Backpropagation Through Time

BPTT Derivation: Returning to a Single Sequence

The forward computation has been described both for a single sequence and for a mini-batch. In the mini-batch formulation, matrices such as ${\mathbf{X}}_t\in {\mathbb{R}}^{B\times n_{\text{inputs}}}$ and ${\mathbf{H}}_t\in {\mathbb{R}}^{B\times n_{\text{neurons}}}$ are used, with one sequence per row.

For the derivation of Backpropagation Through Time, the notation returns to a single sequence. Thus, ${\hat{\mathbf{y}}}_t$, ${\mathbf{y}}_t$, ${\mathbf{o}}_t$, and ${\mathbf{h}}_t$ are vectors rather than mini-batch matrices. This is purely a notational choice; the mini-batch case follows from the single-sequence case by the linearity of differentiation.

If the mini-batch loss is defined as the average of the per-sequence losses,

$${\mathcal{L}}_{\text{batch}}=\frac{1}{B}\sum _{i=1}^B{\mathcal{L}}^{(i)},$$

with ${\mathcal{L}}^{(i)}$ the loss computed on the $i$-th sequence of the mini-batch, then for any parameter $\theta$ the linearity of $\partial /\partial \theta$ gives

$$\frac{\partial {\mathcal{L}}_{\text{batch}}}{\partial \theta }=\frac{1}{B}\sum _{i=1}^B\frac{\partial {\mathcal{L}}^{(i)}}{\partial \theta }.$$

The single-sequence derivation that follows therefore applies unchanged at the mini-batch level: each $\partial {\mathcal{L}}^{(i)}/\partial \theta$ is computed by the formulas derived below, applied to sequence $i$ in isolation, and the resulting per-sequence gradients are averaged. The same identity holds with sum reduction (i.e., without the $1/B$ factor), which is the convention used by some implementations.

The total loss over the unfolded sequence is defined as

$$\begin{array}{rl}\mathcal{L}(\hat{\mathbf{y}},\mathbf{y})& =\sum _{t=1}^T{L}_t({\hat{\mathbf{y}}}_t,{\mathbf{y}}_t)\\ & =-\sum _{t=1}^T{\mathbf{y}}_t^{\top }\log {\hat{\mathbf{y}}}_t\\ & =-\sum _{t=1}^T{\mathbf{y}}_t^{\top }\log \left[\mathrm{softmax}({\mathbf{o}}_t)\right].\end{array}$$

Here, ${\mathbf{y}}_t$ is the one-hot target vector at time $t$, and ${\hat{\mathbf{y}}}_t$ is the predicted probability distribution.

The sequence loss sums the local losses over the $T$ time steps. Averaging instead, $\frac{1}{T}{\sum }_t{L}_t$, is an equally common convention: it rescales every gradient derived below by $1/T$ and changes nothing structural, exactly as the mini-batch reduction can be taken as a sum or a mean. The summed form is kept here to match the standard BPTT presentation.

Cross-Entropy and One-Hot Targets

For a single time step, the multiclass cross-entropy can be written as

$$L_t=-{\mathbf{y}}_t^{\top }\log {\hat{\mathbf{y}}}_t,$$

where ${\mathbf{y}}_t\in \{0,1{\}}^C$ is a one-hot target vector and ${\hat{\mathbf{y}}}_t\in {\mathbb{R}}^C$ is the predicted probability vector after softmax. The dot product hides the sum over classes:

$${\mathbf{y}}_t^{\top }\log {\hat{\mathbf{y}}}_t=\sum _{c=1}^C{y}_t^{(c)}\log {\hat{y}}_t^{(c)}.$$

Since ${\mathbf{y}}_t$ is one-hot, only the true class contributes. Therefore,

$$L_t=-\log {\hat{y}}_t^{(c^{\star })}.$$

---

Error Signals and the Backward Recursion

This section derives the gradients of $\mathcal{L}$ with respect to the five shared parameters. The derivation follows the same pattern as backpropagation for MLPs: introduce error signals at intermediate points of the computational graph, derive a backward recursion that computes them in a single sweep, and express each parameter gradient as an outer product of an error signal and a forward input. The whole procedure runs in $\mathcal{O}(T)$ time on a sequence of length $T$.

What distinguishes BPTT from MLP backpropagation is the temporal axis. The same parameters are reused at every time step, so each parameter contributes to every local loss $L_1,\dots ,L_T$. The recursion that propagates the error runs from time $T$ down to time $1$, and the parameter gradients are sums over the $T$ temporal slices. Apart from this extra summation, every step has a direct counterpart in the MLP case.

Relation to the MLP error signal

The error signal ${\delta }_j^l=\partial \mathcal{L}/\partial z_j^l$ introduced in MLP backpropagation is the gradient of the loss with respect to the pre-activation of a unit. The same role is played here by ${\mathbf{\delta }}_t^a$.

The letters used for activations and pre-activations are swapped between the two sections:

• the MLP notes use $z$ for the pre-activation and $a$ for the activation
• the RNN notes use ${\mathbf{a}}_t$ for the pre-activation and ${\mathbf{h}}_t$ for the activation.

The structural role of the error signal is identical in both: a gradient against the pre-activation, computed by a backward recursion and combined with forward inputs to produce parameter gradients.

Three error signals

For a single sequence, define the error signal at three points of the computational graph:

$${\mathbf{\delta }}_t^o:=\frac{\partial \mathcal{L}}{\partial {\mathbf{o}}_t},{\mathbf{\delta }}_t^h:=\frac{\partial \mathcal{L}}{\partial {\mathbf{h}}_t},{\mathbf{\delta }}_t^a:=\frac{\partial \mathcal{L}}{\partial {\mathbf{a}}_t}.$$

These are the gradients of the total loss with respect to the logits, the hidden state, and the pre-activation at time $t$. All three are vectors of the appropriate dimension.

The first one is immediate. Under the standard softmax-plus-cross-entropy combination, the derivative simplifies to

$${\mathbf{\delta }}_t^o={\hat{\mathbf{y}}}_t-{\mathbf{y}}_t.$$

This is the result of a clean cancellation between the softmax derivative and the cross-entropy derivative. The takeaway is that the gradient at the output is the prediction error, signed and component-wise: a vector of mismatches between the predicted probabilities ${\hat{\mathbf{y}}}_t$ and the one-hot targets ${\mathbf{y}}_t$, with no residual Jacobian to track. A broader treatment of the same identity is in Softmax and Cross-Entropy.

Optional: derivation of ${\mathbf{\delta }}_t^o={\hat{\mathbf{y}}}_t-{\mathbf{y}}_t$

Since the total loss $\mathcal{L}={\sum }_s{L}_s$ depends on ${\mathbf{o}}_t$ only through the local loss $L_t$ at the same time step,

$${\mathbf{\delta }}_t^o=\frac{\partial \mathcal{L}}{\partial {\mathbf{o}}_t}=\frac{\partial L_t}{\partial {\mathbf{o}}_t}.$$

The local loss is $L_t=-{\sum }_c{y}_t^{(c)}\log {\hat{y}}_t^{(c)}$, with ${\hat{y}}_t^{(c)}=e^{o_t^{(c)}}/{\sum }_j{e}^{o_t^{(j)}}$.

Softmax derivative. Differentiating the softmax with respect to $o_t^{(i)}$,

$$\frac{\partial {\hat{y}}_t^{(c)}}{\partial o_t^{(i)}}={\hat{y}}_t^{(c)}({\delta }_{ci}-{\hat{y}}_t^{(i)}),$$

where ${\delta }_{ci}$ is the Kronecker delta (which equals $1$ when $c=i$ and $0$ otherwise; unrelated to the error signal $\mathbf{\delta }$).

Cross-entropy derivative. Differentiating $L_t$ with respect to ${\hat{y}}_t^{(c)}$,

$$\frac{\partial L_t}{\partial {\hat{y}}_t^{(c)}}=-\frac{y_t^{(c)}}{{\hat{y}}_t^{(c)}}.$$

Chain rule. Combining the two via the multivariate chain rule across the $C$ logits,

$$\begin{array}{rl}\frac{\partial L_t}{\partial o_t^{(i)}}& =\sum _c\frac{\partial L_t}{\partial {\hat{y}}_t^{(c)}}\frac{\partial {\hat{y}}_t^{(c)}}{\partial o_t^{(i)}}\\ & =\sum _c\left(-\frac{y_t^{(c)}}{{\hat{y}}_t^{(c)}}\right){\hat{y}}_t^{(c)}({\delta }_{ci}-{\hat{y}}_t^{(i)})\\ & =-\sum _c{y}_t^{(c)}({\delta }_{ci}-{\hat{y}}_t^{(i)})\\ & =-y_t^{(i)}+{\hat{y}}_t^{(i)}\sum _c{y}_t^{(c)}\\ & ={\hat{y}}_t^{(i)}-y_t^{(i)},\end{array}$$

where the last step uses ${\sum }_c{y}_t^{(c)}=1$ because ${\mathbf{y}}_t$ is one-hot. Stacking the scalar derivatives into a vector gives the boxed identity ${\mathbf{\delta }}_t^o={\hat{\mathbf{y}}}_t-{\mathbf{y}}_t$.

Where the simplification comes from. The factor ${\hat{y}}_t^{(c)}$ at the front of the softmax derivative cancels the $1/{\hat{y}}_t^{(c)}$ that the cross-entropy derivative leaves behind. Without this cancellation, the gradient at the output would still depend on the softmax outputs in a non-trivial way; with it, the gradient is the prediction error in its rawest form, which is exactly what makes the downstream backward recursion clean.

Backward recursion for ${\mathbf{\delta }}_t^h$

The hidden state ${\mathbf{h}}_t$ affects the total loss through two paths:

1. Locally, through the output at time $t$: ${\mathbf{h}}_t\to {\mathbf{o}}_t\to {\hat{\mathbf{y}}}_t\to L_t$.
2. Recursively, through the future hidden states: ${\mathbf{h}}_t\to {\mathbf{a}}_{t+1}\to {\mathbf{h}}_{t+1}\to \cdots$, which in turn affect $L_{t+1},L_{t+2},\dots ,L_T$.

Each path contributes to ${\mathbf{\delta }}_t^h$ via the multivariate chain rule. The general identity, for any vector-valued forward step $\mathbf{u}=\mathbf{f}(\mathbf{v})$ and any scalar loss $\mathcal{L}$, is

$$\frac{\partial \mathcal{L}}{\partial \mathbf{v}}={\left(\frac{\partial \mathbf{u}}{\partial \mathbf{v}}\right)}^{\top }\frac{\partial \mathcal{L}}{\partial \mathbf{u}}.$$

Here $\partial \mathbf{u}/\partial \mathbf{v}$ is the Jacobian matrix of $\mathbf{f}$, and the transpose makes the dimensions match: a column gradient with respect to $\mathbf{u}$ on the right side produces, after multiplication by the transposed Jacobian, a column gradient with respect to $\mathbf{v}$ on the left side.

Applied to the two paths through which ${\mathbf{h}}_t$ affects $\mathcal{L}$:

• Local path $({\mathbf{h}}_t\to {\mathbf{o}}_t)$. From the output equation ${\mathbf{o}}_t={\mathbf{W}}_{yh}{\mathbf{h}}_t+{\mathbf{b}}_y$, the Jacobian is $\partial {\mathbf{o}}_t/\partial {\mathbf{h}}_t={\mathbf{W}}_{yh}$. Its contribution to ${\mathbf{\delta }}_t^h$ is therefore ${\mathbf{W}}_{yh}^{\top }{\mathbf{\delta }}_t^o$.
• Recursive path $({\mathbf{h}}_t\to {\mathbf{a}}_{t+1})$. From the recurrent transition ${\mathbf{a}}_{t+1}={\mathbf{W}}_{xh}{\mathbf{x}}_{t+1}+{\mathbf{W}}_{hh}{\mathbf{h}}_t+{\mathbf{b}}_h$, the Jacobian is $\partial {\mathbf{a}}_{t+1}/\partial {\mathbf{h}}_t={\mathbf{W}}_{hh}$. Its contribution is ${\mathbf{W}}_{hh}^{\top }{\mathbf{\delta }}_{t+1}^a$, where ${\mathbf{\delta }}_{t+1}^a$ is the future error signal that has already been computed by the backward sweep.

Summing the two contributions yields the recursion for ${\mathbf{\delta }}_t^h$:

$$\boxed{{\mathbf{\delta }}_t^h={\mathbf{W}}_{yh}^{\top }{\mathbf{\delta }}_t^o+{\mathbf{W}}_{hh}^{\top }{\mathbf{\delta }}_{t+1}^a.}$$

The recursion is initialized at the right end of the sequence. No future hidden state exists beyond $t=T$, so the recursive term vanishes by setting ${\mathbf{\delta }}_{T+1}^a:=0$, and

$${\mathbf{\delta }}_T^h={\mathbf{W}}_{yh}^{\top }{\mathbf{\delta }}_T^o.$$

From ${\mathbf{\delta }}_t^h$ to ${\mathbf{\delta }}_t^a$

The pre-activation error is one element-wise nonlinearity away from the hidden-state error, and the corresponding Jacobian has a very particular structure that is worth deriving in detail.

The map ${\mathbf{h}}_t={\varphi }_h({\mathbf{a}}_t)$ acts component by component: the $i$-th coordinate of ${\mathbf{h}}_t$ depends only on the $i$-th coordinate of ${\mathbf{a}}_t$,

$$h_t^{(i)}={\varphi }_h\left(a_t^{(i)}\right),i=1,\dots ,n_{\text{neurons}}.$$

Changing $a_t^{(j)}$ with $j\ne i$ has no effect on $h_t^{(i)}$, because the activation function never mixes coordinates. The partial derivatives of ${\mathbf{h}}_t$ with respect to ${\mathbf{a}}_t$ are therefore

$$\frac{\partial h_t^{(i)}}{\partial a_t^{(j)}}=\{\begin{array}{ll}{\varphi }_h^{\prime }\left(a_t^{(i)}\right)& \text{if }i=j,\\ 0& \text{if }i\ne j.\end{array}$$

Assembling these scalar derivatives into the Jacobian matrix gives a square matrix of shape $n_{\text{neurons}}\times n_{\text{neurons}}$ whose only non-zero entries lie on the diagonal:

$$\frac{\partial {\mathbf{h}}_t}{\partial {\mathbf{a}}_t}=\mathrm{diag}({\varphi }_h^{\prime }({\mathbf{a}}_t))=\left[\begin{array}{cccc}{\varphi }_h^{\prime }(a_t^{(1)})& 0& \cdots & 0\\ 0& {\varphi }_h^{\prime }(a_t^{(2)})& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\varphi }_h^{\prime }(a_t^{(n_{\text{neurons}})})\end{array}\right].$$

The diagonal structure is a direct consequence of the activation acting independently on each coordinate; there is no coupling between coordinates to fill in the off-diagonal entries. The same construction, an element-wise nonlinearity contributing a diagonal $\mathrm{diag}(f^{\prime })$ factor to a Jacobian product, recurs throughout deep learning: it governs the exploding gradient across time and the gradient highway through feedforward residual blocks.

A useful general identity follows: multiplying any vector $\mathbf{v}$ by a diagonal matrix $\mathrm{diag}(\mathbf{d})$ is the same as taking the element-wise product,

$$\mathrm{diag}(\mathbf{d})\mathbf{v}=\mathbf{d}\odot \mathbf{v}.$$

Applying the chain rule with the Jacobian above, and using that the transpose of a diagonal matrix is the diagonal matrix itself, the pre-activation error becomes

$${\mathbf{\delta }}_t^a={\left(\frac{\partial {\mathbf{h}}_t}{\partial {\mathbf{a}}_t}\right)}^{\top }{\mathbf{\delta }}_t^h=\mathrm{diag}({\varphi }_h^{\prime }({\mathbf{a}}_t)){\mathbf{\delta }}_t^h={\varphi }_h^{\prime }({\mathbf{a}}_t)\odot {\mathbf{\delta }}_t^h.$$

Boxing the final form:

$$\boxed{{\mathbf{\delta }}_t^a={\varphi }_h^{\prime }({\mathbf{a}}_t)\odot {\mathbf{\delta }}_t^h,}$$

where $\odot$ denotes the element-wise (Hadamard) product. For the canonical choice ${\varphi }_h=\tanh$, the derivative is ${\varphi }_h^{\prime }({\mathbf{a}}_t)=1-{\mathbf{h}}_t\odot {\mathbf{h}}_t$, which can be evaluated for free from the hidden state ${\mathbf{h}}_t$ saved during the forward pass.

The backward sweep in algorithmic form

Combining the two updates yields a single loop that walks the unrolled graph from right to left:

$${\mathbf{\delta }}_{T+1}^a:=0,\text{for }t=T,T-1,\dots ,1:\{\begin{array}{l}{\mathbf{\delta }}_t^h={\mathbf{W}}_{yh}^{\top }{\mathbf{\delta }}_t^o+{\mathbf{W}}_{hh}^{\top }{\mathbf{\delta }}_{t+1}^a,\\ {\mathbf{\delta }}_t^a={\varphi }_h^{\prime }({\mathbf{a}}_t)\odot {\mathbf{\delta }}_t^h.\end{array}$$

Each iteration performs a constant number of matrix-vector products. The full backward pass therefore visits each of the $T$ time steps exactly once, for a total cost of $\mathcal{O}(T\cdot n_{\text{neurons}}^2)$, linear in the sequence length.

Dynamic programming along the temporal axis

Standard backpropagation in an MLP is already a dynamic-programming algorithm: each layer’s error signal ${\mathbf{\delta }}^l$ is computed once from ${\mathbf{\delta }}^{l+1}$ and reused in every parameter gradient at layer $l$. Backprop is linear in the depth precisely because of this reuse. The novelty of BPTT is not dynamic programming itself, but the observation that the same DP structure survives the temporal axis: across time steps, just as across depth, error signals are computed once and reused.

This deserves emphasis because the naive chain-rule expansion of $\partial \mathcal{L}/\partial {\mathbf{W}}_{hh}$ produces a double sum over temporal paths that looks quadratic in $T$. Sharing ${\mathbf{W}}_{hh}$ across all time steps could in principle have spoiled the DP structure, forcing each path to be computed independently. The recursion on ${\mathbf{\delta }}_t^h$ and ${\mathbf{\delta }}_t^a$ above is the proof that this does not happen: the temporal dependency stays linear in $T$, exactly as the depth-wise dependency stays linear in $L$ for an MLP. Same principle, extended to a new axis.

This is also exactly what PyTorch’s autograd computes when loss.backward() is called on the unrolled cell of the minimal PyTorch implementation: the recursion above runs inside the autograd engine, one matrix-vector product per time step. The hand-written derivation in this note and the production-grade autograd path implement the same algorithm.

---

Parameter Gradients as Outer-Product Sums

Once ${\mathbf{\delta }}_t^o$ and ${\mathbf{\delta }}_t^a$ are available for every $t$, the gradients with respect to the five shared parameters are immediate, because the same template applies to every weight matrix in the cell. It is worth deriving the template once in general before listing the five formulas.

Why the gradient with respect to a weight matrix is an outer product. Consider a generic linear step $\mathbf{u}=\mathbf{Wv}+\mathbf{c}$ inside the forward pass. The gradient of $\mathcal{L}$ with respect to the scalar entry $W_{ij}$ uses the chain rule:

$$\frac{\partial \mathcal{L}}{\partial W_{ij}}=\frac{\partial \mathcal{L}}{\partial u_i}\frac{\partial u_i}{\partial W_{ij}}.$$

The first factor is the $i$-th component of the error signal at the output of the linear step, ${\delta }_u^{(i)}:=\partial \mathcal{L}/\partial u_i$. The second factor follows from $u_i={\sum }_k{W}_{ik}{v}_k+c_i$: only the term $W_{ij}{v}_j$ involves the specific entry $W_{ij}$, so

$$\frac{\partial u_i}{\partial W_{ij}}=v_j.$$

Multiplying the two factors gives the scalar derivative

$$\frac{\partial \mathcal{L}}{\partial W_{ij}}={\delta }_u^{(i)}{v}_j,$$

and assembling these scalars over the index pair $(i,j)$ recovers the outer product

$$\frac{\partial \mathcal{L}}{\partial \mathbf{W}}={\mathbf{\delta }}_u{\mathbf{v}}^{\top }.$$

The same template, applied to a bias term $\mathbf{c}$, gives $\partial \mathcal{L}/\partial \mathbf{c}={\mathbf{\delta }}_u$, since $\partial u_i/\partial c_i=1$.

Applying the template to the RNN cell. Each of the five parameters of the cell appears in exactly one linear step of the forward pass, multiplying exactly one input vector at each time step $t$:

• ${\mathbf{W}}_{yh}$ multiplies ${\mathbf{h}}_t$ in the output equation, and its error signal at the output is ${\mathbf{\delta }}_t^o$;
• ${\mathbf{W}}_{xh}$ multiplies ${\mathbf{x}}_t$ in the recurrent transition, and its error signal at the output is ${\mathbf{\delta }}_t^a$;
• ${\mathbf{W}}_{hh}$ multiplies ${\mathbf{h}}_{t-1}$ in the recurrent transition, and its error signal at the output is ${\mathbf{\delta }}_t^a$;
• ${\mathbf{b}}_y$ and ${\mathbf{b}}_h$ inherit the corresponding output error signals directly.

Summing the per-step outer products over $t$ gives the five gradients:

$$\boxed{\begin{array}{rlrl}\frac{\partial \mathcal{L}}{\partial {\mathbf{W}}_{yh}}& =\sum _{t=1}^T{\mathbf{\delta }}_t^o{\mathbf{h}}_t^{\top },& \frac{\partial \mathcal{L}}{\partial {\mathbf{b}}_y}& =\sum _{t=1}^T{\mathbf{\delta }}_t^o,\\ \frac{\partial \mathcal{L}}{\partial {\mathbf{W}}_{xh}}& =\sum _{t=1}^T{\mathbf{\delta }}_t^a{\mathbf{x}}_t^{\top },& \frac{\partial \mathcal{L}}{\partial {\mathbf{W}}_{hh}}& =\sum _{t=1}^T{\mathbf{\delta }}_t^a{\mathbf{h}}_{t-1}^{\top },\\ \frac{\partial \mathcal{L}}{\partial {\mathbf{b}}_h}& =\sum _{t=1}^T{\mathbf{\delta }}_t^a.& & \end{array}}$$

The structure is the same as ordinary backpropagation in an MLP; the only addition is the temporal sum. The two input-side weight matrices have a parallel role: ${\mathbf{W}}_{xh}$ multiplies ${\mathbf{x}}_t$ in the forward pass, ${\mathbf{W}}_{hh}$ multiplies ${\mathbf{h}}_{t-1}$, and the corresponding gradient is the outer product of ${\mathbf{\delta }}_t^a$ with whichever input the matrix saw at time $t$. No separate derivation is needed for the two cases: they are the same operation applied to two different inputs.

Equivalence with the long-form chain rule

The compact formulas above are not an alternative to the chain rule; they are its natural form once the backward recursion has been applied. Expanding ${\mathbf{\delta }}_t^a$ via the recursion reproduces the long-form expression

$$\frac{\partial \mathcal{L}}{\partial {\mathbf{W}}_{hh}}=\sum _{t=1}^T\sum _{k=1}^t\frac{\partial L_t}{\partial {\hat{\mathbf{y}}}_t}\frac{\partial {\hat{\mathbf{y}}}_t}{\partial {\mathbf{o}}_t}\frac{\partial {\mathbf{o}}_t}{\partial {\mathbf{h}}_t}\left(\prod _{j=k}^{t-1}\frac{\partial {\mathbf{h}}_{j+1}}{\partial {\mathbf{h}}_j}\right)\frac{\partial {\mathbf{h}}_k}{\partial {\mathbf{W}}_{hh}},$$

which sums over all causal paths from ${\mathbf{W}}_{hh}$ to each local loss. The inner index $k$ corresponds to the time step at which ${\mathbf{W}}_{hh}$ is “tagged” along the path: the parameter participates in the computation of every hidden state ${\mathbf{h}}_k$, and each such participation defines one path that ends at the loss $L_t$. The same sum applies to ${\mathbf{W}}_{xh}$ with $\partial {\mathbf{h}}_k/\partial {\mathbf{W}}_{xh}$ in place of $\partial {\mathbf{h}}_k/\partial {\mathbf{W}}_{hh}$.

This expansion is correct but should not be implemented as written. Computing each path independently leads to $\mathcal{O}(T^2)$ operations. The backward recursion in the previous section is the same sum evaluated in $\mathcal{O}(T)$ by reusing each ${\mathbf{\delta }}_t^a$ across all the paths that pass through time $t$.

---

The Product of Recurrent Jacobians

The unrolled form contains one term that is worth isolating, because it governs the stability of the entire algorithm:

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

This is the chain of recurrent Jacobians multiplied together along the temporal gap between time $k$ and time $t$. When the upper index is smaller than the lower one, the product is the identity by convention.

The same product appears, in disguise, inside the backward recursion: every step of the backward sweep applies ${\mathbf{W}}_{hh}^{\top }$ to ${\mathbf{\delta }}^a$ once, so propagating an error from time $t$ back to time $k$ involves $t-k$ applications of a Jacobian of ${\mathbf{W}}_{hh}$. For example, for $t=3$ and $k=1$,

$$\frac{\partial {\mathbf{h}}_3}{\partial {\mathbf{h}}_1}=\frac{\partial {\mathbf{h}}_3}{\partial {\mathbf{h}}_2}\cdot \frac{\partial {\mathbf{h}}_2}{\partial {\mathbf{h}}_1},$$

a product of two recurrent Jacobians, each of shape $n_{\text{neurons}}\times n_{\text{neurons}}$.

Stability is dictated by this product

The chain ${\prod }_{j=k}^{t-1}\partial {\mathbf{h}}_{j+1}/\partial {\mathbf{h}}_j$ is the term through which BPTT carries gradient information across time. When its norm decays geometrically along the temporal gap, the gradient vanishes; when it grows geometrically, the gradient explodes. The full stability analysis, including the vanishing and exploding gradient regimes in vanilla RNNs, is the subject of BPTT Problems, and the optimizer-level remedy for the explosion is the subject of Gradient Clipping.

---

Conclusion

Every gradient computed by BPTT has the same shape:

$$\frac{\partial \mathcal{L}}{\partial \theta }=\sum _{t=1}^T\left(\text{local error signal at time }t\right)\otimes \left(\text{local forward input at time }t\right),$$

where the error signal is provided by the backward recursion on ${\mathbf{\delta }}_t^h$ and ${\mathbf{\delta }}_t^a$, and the forward input is whatever quantity the parameter $\theta$ multiplied during the forward pass (${\mathbf{x}}_t$, ${\mathbf{h}}_{t-1}$, ${\mathbf{h}}_t$, or the constant vector $1$ for biases). A single backward sweep through the unrolled graph in $\mathcal{O}(T)$ time computes every parameter gradient.

The same recurrent-Jacobian product that makes this computation possible also exposes vanilla RNNs to the stability issues analyzed in BPTT Problems: the algorithm’s efficiency and its instability share the same mathematical source.

For sequences where even $\mathcal{O}(T)$ becomes prohibitive in memory (because every ${\mathbf{h}}_t$ and ${\mathbf{a}}_t$ must be retained for the backward pass), the practical training algorithm is Truncated BPTT, treated in BPTT Variants.