Recurrent Neural Network

A recurrent neural network (RNN) processes a sequence one element at a time and carries a state vector from each step to the next. The state at step $t$ is computed from the state at $t-1$ and the current input by the same function, with the same weights, at every step. A recurrent network can therefore read sequences of any length: the columns of a word image, the frames of a video, the words of a caption. Unrolled over $T$ steps it is a $T$-layer neural network whose layers share their weights, and it is trained by applying backpropagation to that unrolled graph, a procedure called backpropagation through time. Because the gradient that reaches early steps is a product of one factor per step, it tends to shrink or grow geometrically with the time span. Gated units, the long short-term memory (LSTM) and the gated recurrent unit (GRU), were built to keep that product near 1. In vision, recurrent layers, usually LSTMs, were combined with convolutional networks for captioning, video recognition, text recognition and video prediction. Captioning and text recognition later moved to attention and transformer decoders, and recurrent update operators appear in iterative refinement such as RAFT’s optical flow.

The defining LSTM paper (Hochreiter and Schmidhuber, 1997) could not be read for this article. Its host for the authors’ copy failed a certificate check and the publisher’s site is not used. The LSTM is therefore described here from readable papers that restate it: Pascanu, Mikolov and Bengio (2013), Jozefowicz, Zaremba and Sutskever (2015), Cho et al. (2014), the LRCN paper (Donahue et al., 2015) and Chapter 10 of the Deep Learning book (Goodfellow, Bengio and Courville, 2016).

Core Idea

A multilayer perceptron applied to a sequence of fixed length would need separate weights for every position, so it could not handle other lengths or use what it learned at one position at another. A recurrent network shares one transition across time. Goodfellow et al. (§10.1) compare this with a convolutional network, which shares one filter across image positions; a one-dimensional convolution over time also shares parameters, but each output then depends only on a few neighbouring inputs, whereas each recurrent state depends on everything before it through the chain of earlier states.

The state $\mathbf{h}_t$ is a learned summary of $\mathbf{x}_1, \dots, \mathbf{x}_t$. It must be of fixed size, so the network has to decide what to keep. Long-range information survives only if the transition preserves it over many steps, which is where training becomes hard.

Building Blocks

• The transition. The simple (“vanilla”, or Elman-style) recurrent layer computes $\mathbf{h}_t = \tanh(W \mathbf{h}_{t-1} + U \mathbf{x}_t + \mathbf{b})$, with a hidden-to-hidden matrix $W$, an input matrix $U$ and a bias (Goodfellow et al., Eqs. 10.8–10.9). Pascanu et al. (2013) attribute the recurrent network to work of the 1980s, citing Rumelhart et al. (1986), Elman (1990) and Werbos (1988).
• Outputs. A readout $\mathbf{o}_t = V \mathbf{h}_t + \mathbf{c}$, often followed by a softmax, can be emitted at every step (one label per frame or per character) or only at the end (one label per sequence).
• The initial state. $\mathbf{h}_0$ is set to zero, given by the user or learned (Pascanu et al., §1.1); in captioning and other encoder–decoder models it carries the encoding of the input.
• Gated cells. LSTM and GRU cells replace the single $\tanh$ unit with a small circuit of sigmoid gates that control what is written to, kept in and read from the state.

How Data Flows

Goodfellow et al. describe these patterns (§10.2, §10.2.4 and §10.4), all of which appear in vision:

• Sequence to sequence of the same length. An output at every step: per-column character scores in a text recognizer, per-frame labels in a video.
• Sequence to one output. The network reads the whole sequence and a readout uses the final state (or an average of per-step outputs, as LRCN does for video clips).
• One input to a sequence. A caption decoder receives an image encoding and emits words one at a time, each conditioned on the words before it.
• Sequence to sequence of a different length. An encoder RNN summarizes the input and a decoder RNN generates the output, the encoder–decoder design of Cho et al. (2014) and of Sutskever, Vinyals and Le (2014).

Computation is sequential: step $t$ cannot start before step $t-1$ has finished. The forward pass over $T$ steps costs $O(T)$ time that parallelization cannot reduce, and the states must be stored for the backward pass, so memory is also $O(T)$ (Goodfellow et al., §10.2).

Teacher forcing. A network whose next input is its own previous output (a caption decoder that is fed the word it just produced) can be trained by feeding the correct previous output instead, so that every step is trained from the ground truth. In a network whose only recurrence runs from outputs to the next step, the steps then decouple and can be trained in parallel; once the hidden units also depend on earlier steps, as in an LSTM caption decoder, backpropagation through time is still needed. At test time the model’s own outputs are fed back, a mismatch between training and use (Goodfellow et al., §10.2.1).

Mathematical Formulation

Backpropagation through time. For the simple recurrence with a loss summed over steps, $L = \sum_t L_t$, the gradient with respect to a shared weight is a sum over time. Pascanu et al. (2013, Eqs. 3–5) write it as

$$\frac{\partial L}{\partial \theta} = \sum_{t} \sum_{k \le t} \frac{\partial L_t}{\partial \mathbf{h}_t} \, \frac{\partial \mathbf{h}_t}{\partial \mathbf{h}_k} \, \frac{\partial^+ \mathbf{h}_k}{\partial \theta}, \qquad \frac{\partial \mathbf{h}_t}{\partial \mathbf{h}_k} = \prod_{i=k+1}^{t} \frac{\partial \mathbf{h}_i}{\partial \mathbf{h}_{i-1}},$$

where $\partial^+$ is the “immediate” derivative that treats $\mathbf{h}_{k-1}$ as a constant. Each term measures how the parameters at step $k$ affect the loss at step $t$; Pascanu et al. call the terms with $t - k$ large the long-term contributions. In practice the sum is computed backward through time with one error vector. For $\mathbf{h}_t = \tanh(\mathbf{a}_t)$, $\mathbf{a}_t = W \mathbf{h}_{t-1} + U \mathbf{x}_t$ and a per-step loss with error $\mathbf{e}_t = \partial L_t / \partial \mathbf{o}_t$,

$$\boldsymbol{\delta}_t = \frac{\partial L}{\partial \mathbf{h}_t} = V^\top \mathbf{e}_t + W^\top \Big(\boldsymbol{\delta}_{t+1} \odot \big(1 - \mathbf{h}_{t+1}^2\big)\Big), \qquad \frac{\partial L}{\partial W} = \sum_t \big(\boldsymbol{\delta}_t \odot (1 - \mathbf{h}_t^2)\big)\, \mathbf{h}_{t-1}^\top,$$

and similarly for $U$ (with $\mathbf{x}_t$ in place of $\mathbf{h}_{t-1}$) and $V$ (with $\mathbf{e}_t \mathbf{h}_t^\top$). This is ordinary backpropagation on the unrolled graph with the weight gradients summed over the copies; no special algorithm is needed (Goodfellow et al., §10.2.2). The example below implements it and checks it against reverse-mode automatic differentiation as implemented in PyTorch.

Why gradients vanish or explode over time. For the simple recurrence each factor is $\partial \mathbf{h}_i / \partial \mathbf{h}_{i-1} = D_i W$ with $D_i = \operatorname{diag}(1 - \mathbf{h}_i^2)$, so $\partial \mathbf{h}_t / \partial \mathbf{h}_k$ is a product of $t-k$ Jacobians built from the same matrix. Pascanu et al. (§2.1) show for the linear case that a spectral radius $\rho(W) < 1$ is sufficient for the long-term components to vanish and $\rho(W) > 1$ necessary for them to explode. For a nonlinearity with $|\phi'| \le \gamma$ they show that a largest singular value $\lambda_1 < 1/\gamma$ is sufficient for vanishing, since then $\lVert \partial \mathbf{h}_t / \partial \mathbf{h}_k \rVert \le (\lambda_1 \gamma)^{t-k}$; inverting the argument, $\lambda_1 > 1/\gamma$ is necessary for explosion, with $\gamma = 1$ for tanh and $1/4$ for the sigmoid. The vanishing gradient problem article treats the same product over the depth of feedforward networks.

Bengio, Simard and Frasconi (1994) showed a deeper obstacle, as Goodfellow et al. (§10.7) summarize it: to store information robustly against small perturbations, a recurrent network must operate where gradients vanish, so the gradient of a long-term interaction is exponentially smaller than that of a short-term one. Goodfellow et al. add that in those experiments the probability of training a simple RNN successfully by SGD fell rapidly to zero for dependencies spanning only 10 or 20 steps. Jozefowicz et al. (2015, §2) put the practical consequence this way: vanishing does not make the whole gradient small, only its long-term components, so the network learns short-term dependencies and misses the long ones.

The LSTM cell. In the equations used by PyTorch (torch.nn.LSTM, release 2.14.1) and, up to notation, by the LRCN paper,

$$\begin{aligned} \mathbf{i}_t &= \sigma(W_{ii}\mathbf{x}_t + W_{hi}\mathbf{h}_{t-1} + \mathbf{b}_i), & \mathbf{f}_t &= \sigma(W_{if}\mathbf{x}_t + W_{hf}\mathbf{h}_{t-1} + \mathbf{b}_f), \\ \mathbf{g}_t &= \tanh(W_{ig}\mathbf{x}_t + W_{hg}\mathbf{h}_{t-1} + \mathbf{b}_g), & \mathbf{o}_t &= \sigma(W_{io}\mathbf{x}_t + W_{ho}\mathbf{h}_{t-1} + \mathbf{b}_o), \\ \mathbf{c}_t &= \mathbf{f}_t \odot \mathbf{c}_{t-1} + \mathbf{i}_t \odot \mathbf{g}_t, & \mathbf{h}_t &= \mathbf{o}_t \odot \tanh(\mathbf{c}_t), \end{aligned}$$

with input gate $\mathbf{i}$, forget gate $\mathbf{f}$, candidate (LRCN’s “input modulation gate”) $\mathbf{g}$ and output gate $\mathbf{o}$. The cell state $\mathbf{c}_t$ is updated additively. Along it, $\partial \mathbf{c}_t / \partial \mathbf{c}_{t-1}$ contains the term $\operatorname{diag}(\mathbf{f}_t)$, so the direct path from $\mathbf{c}_k$ to $\mathbf{c}_t$ multiplies the gradient by $\prod_{i=k+1}^{t} \mathbf{f}_i$ elementwise, and a forget gate near 1 passes it almost unchanged. Pascanu et al. (§3.1) describe the LSTM as relying on a linear unit with a self-connection of value 1, guarded by learned input and output gates, and note several variations of this structure. Making the weight on this self-loop a gate computed from the context, rather than a constant, is credited by Goodfellow et al. (§10.10.1) to Gers, Schmidhuber and Cummins (2000). Jozefowicz et al. explain the benefit by the additive update: a network that computes an increment $\Delta S_t$ and adds it to the state has the same expressive power as one that computes $S_t$ directly, but its long-range gradient cannot vanish. They also note that an LSTM with $n$ cells has a $2n$-dimensional state, $(\mathbf{h}_t, \mathbf{c}_t)$: a slow state $\mathbf{c}$ and a fast state $\mathbf{h}$. LSTMs come in slightly different forms. The version Jozefowicz et al. specify follows Graves (2013) without peephole connections, and as printed in their §2 it computes the output gate with $\tanh$ rather than a sigmoid. Goodfellow et al. (Eqs. 10.40–10.44) apply a sigmoid to the candidate and mention a variant in which the cell state is an extra input to the gates.

The forget-gate bias. Jozefowicz et al. (§2.2) point out that the common initialization with small random weights sets the forget gate to about 0.5, which introduces a vanishing factor of 0.5 per step. They recommend a forget-gate bias of 1 or 2 so the gate starts near 1, an idea they trace to Gers et al. (2000). The example below measures what this does at initialization.

The GRU. Cho et al. (2014, §2.3) proposed a simpler unit “motivated by the LSTM”, with a reset gate $\mathbf{r}$ and an update gate $\mathbf{z}$ and no separate cell:

$$\mathbf{r}_t = \sigma(W_r \mathbf{x}_t + U_r \mathbf{h}_{t-1}), \quad \mathbf{z}_t = \sigma(W_z \mathbf{x}_t + U_z \mathbf{h}_{t-1}), \quad \tilde{\mathbf{h}}_t = \tanh\big(W \mathbf{x}_t + U(\mathbf{r}_t \odot \mathbf{h}_{t-1})\big), \quad \mathbf{h}_t = \mathbf{z}_t \odot \mathbf{h}_{t-1} + (1 - \mathbf{z}_t) \odot \tilde{\mathbf{h}}_t.$$

A reset gate near 0 lets the unit drop the previous state; the update gate decides how much of it carries over, which Cho et al. compare to the LSTM’s memory cell and to a leaky-integration unit. The paper itself does not use the name “gated recurrent unit”; Jozefowicz et al. and Goodfellow et al. call this unit the GRU. Sources differ in two conventions. In Cho et al., Jozefowicz et al. and PyTorch, $\mathbf{z}$ is the weight on the old state; in RAFT’s convolutional GRU (Teed and Deng, 2020, Eq. 6) the update reads $\mathbf{h}_t = (1 - \mathbf{z}_t) \odot \mathbf{h}_{t-1} + \mathbf{z}_t \odot \tilde{\mathbf{h}}_t$, so $\mathbf{z}$ is the weight on the new candidate. And PyTorch’s torch.nn.GRU applies the reset gate after the recurrent matrix, $\tanh(W_{in}\mathbf{x}_t + \mathbf{b}_{in} + \mathbf{r}_t \odot (W_{hn}\mathbf{h}_{t-1} + \mathbf{b}_{hn}))$, rather than to $\mathbf{h}_{t-1}$ before it as in Cho et al.’s Eq. 8. The two placements are different functions, so weights cannot simply be moved between implementations that follow different papers.

Why It Works

Weight sharing across time gives a recurrent network the same advantage weight sharing across space gives a convolutional network: what is learned at one position applies at every other, and the parameter count does not grow with the sequence (Goodfellow et al., §10.1). LRCN’s authors make the same point for video: the last predictions of a $T$-step recurrent network are computed by a $T$-layer nonlinear function, yet the weights are reused at every step, so the model learns generic step-to-step dynamics instead of dynamics tied to the frame index (Donahue et al., 2015, §3).

The state lets the network keep information for as long as it needs. Gates make that retention depend on the input, so the time scale of each unit can change with the sequence rather than being fixed by its weights (Goodfellow et al., §10.10.1).

Training Considerations

• Gradient clipping. Pascanu et al. (2013, §2.3) hypothesize that when the gradient explodes along a direction, the curvature does too, so the loss surface has a steep wall, and an ordinary gradient step that hits the wall throws the parameters far across the valley. Their Algorithm 1 rescales the gradient $\hat{\mathbf{g}}$ whenever its norm exceeds a threshold, $\hat{\mathbf{g}} \leftarrow (\text{threshold} / \lVert \hat{\mathbf{g}} \rVert)\, \hat{\mathbf{g}}$, which keeps its direction. They describe it as similar to Mikolov’s earlier element-wise clipping and suggest choosing the threshold from the average gradient norm over many updates; in their experience values from half to ten times that average still converged, though at different speeds (§3.2). Goodfellow et al. (§10.11.1) note that norm clipping per minibatch biases the averaged gradient, a heuristic bias that is empirically useful, and that element-wise clipping does not keep the gradient direction but still descends. Clipping only bounds large gradients; it does nothing for vanishing ones. In PyTorch the norm form is torch.nn.utils.clip_grad_norm_.
• A regularizer for vanishing gradients. Pascanu et al. (§3.3) also add a term that penalizes changes in the norm of the error signal as it travels back one step, $\Omega = \sum_k \big(\lVert \tfrac{\partial L}{\partial \mathbf{h}_{k+1}} \tfrac{\partial \mathbf{h}_{k+1}}{\partial \mathbf{h}_k} \rVert / \lVert \tfrac{\partial L}{\partial \mathbf{h}_{k+1}} \rVert - 1\big)^2$. On the temporal order problem, a synthetic task that requires remembering two symbols across a long sequence, clipping mattered more as sequences grew longer (Fig. 7), and a 50-unit tanh network trained on lengths from 50 to 200 succeeded in all 5 runs and made no errors on lengths up to 5,000 (§4.1). Without clipping, training with the regularizer alone usually failed because of exploding gradients. On character-level Penn Treebank, predicting the next character, their Table 1 reports test bits per character of 1.50 without clipping, 1.42 with clipping and 1.41 with clipping and the regularizer.
• Initialization. Pascanu et al.’s “smart tanh” initialization used a sparse recurrent matrix with spectral radius 0.95, and Goodfellow et al. (§10.8) report that an initial spectral radius of 1.2 with sparse initialization performs well when echo-state-network techniques are used to initialize a fully trained network. For LSTMs, the forget-gate bias above. The weight initialization article covers the feedforward case.
• Unrolling a fixed number of steps. Long sequences are often cut into chunks. Jozefowicz et al. (§3.5) unrolled their networks for 35 steps and carried the hidden state from one minibatch to the next, so gradients flow back at most 35 steps while the state sees the whole sequence; this is commonly called truncated backpropagation through time (our term, not theirs).
• Gates matter unequally. In Jozefowicz et al.’s ablations the forget gate was the most important LSTM gate (removing it gave drastically worse results on their arithmetic and XML tasks, though it mattered little for language modelling), the input gate second, and the output gate the least important, except on the music data with dropout, where LSTMs without an input or without an output gate did best (§4). They evaluated 10,000 architectures (§3.4); the three best, all similar to the GRU, outperformed the LSTM on some problems, but none consistently beat the LSTM and the GRU in all conditions (§5). Their summary (§4) states that the GRU outperformed the LSTM on all tasks except language modelling and that the LSTM with a forget-gate bias of 1 outperformed both on almost all tasks; their tables are more mixed. The GRU was ahead of the LSTM on the arithmetic, XML and Penn Treebank accuracies of Table 1 and on Nottingham without dropout, but behind it on Nottingham with dropout and on Piano-midi (Table 2). The biased LSTM beat both on arithmetic, on Nottingham with dropout, on Piano-midi and in the Penn Treebank models with dropout (Table 3), but trailed the GRU on XML, on Table 1’s Penn Treebank accuracy, on Nottingham without dropout and in the 5M-parameter Penn Treebank model without dropout. Their results came from random hyperparameter search with clipping thresholds of 1 to 20 among the hyperparameters (§3.7).

Variants

• Bidirectional RNNs. Two recurrent networks, one reading forward and one backward, with outputs combined at each step, so each output can depend on the whole input. Goodfellow et al. (§10.3) attribute them to Schuster and Paliwal (1997) and list handwriting and speech recognition among their successes. Text recognizers such as CRNN and the BLSTM systems of handwriting recognition use bidirectional LSTMs because a character’s identity depends on strokes on both sides.
• Stacked (deep) RNNs. The state at each step can be split into layers, each feeding the next; Goodfellow et al. (§10.5) cite Graves et al. (2013) as the first to show a significant benefit of this decomposition, and a later paper by Pascanu and colleagues (2014) for adding depth inside the input-to-hidden, hidden-to-hidden and hidden-to-output maps, which lengthens the path between time steps unless skip connections are added. LRCN’s caption model stacks four LSTMs.
• Leaky units and echo state networks. A unit can mix its old value with the new one at a fixed rate, $\mathbf{h}_t = \alpha \mathbf{h}_{t-1} + (1 - \alpha)\phi(\cdot)$ (Pascanu et al., §3.1). Echo state networks do not learn the recurrent and input weights at all: they fix them and train only a linear readout (Goodfellow et al., §10.8). The two sources describe the fixed weights differently. Pascanu et al. (§3.1) write that the spectral radius of the recurrent matrix is smaller than 1 by construction, so information dies out exponentially; Goodfellow et al. (§10.8) write that echo state networks fix the spectral radius at a value such as 3, where saturating nonlinearities like tanh keep the forward dynamics from exploding, and that recent work on them advocated radii much larger than 1.
• Convolutional recurrent units. For images and video the state is a feature map and the matrix products become convolutions, which keeps the spatial layout and shares weights across positions. The convolutional LSTM of Shi et al. (2015) could not be read for this article. PixelRNN’s Row LSTM shows the construction: its four gates are computed by a convolution over the previous row’s hidden state plus a precomputed convolution over the input, $[\mathbf{o}_i, \mathbf{f}_i, \mathbf{i}_i, \mathbf{g}_i] = \sigma(K^{ss} \circledast \mathbf{h}_{i-1} + K^{is} \circledast \mathbf{x}_i)$, where $\sigma$ is the logistic sigmoid for the three gates and tanh for $\mathbf{g}_i$, with the LSTM’s cell and output updates applied elementwise (van den Oord, Kalchbrenner and Kavukcuoglu, 2016, Eq. 3). Video Pixel Networks combine their encoder outputs over time with a convolutional LSTM of 256 channels that preserves resolution, citing Hochreiter and Schmidhuber and Shi et al. for it (Kalchbrenner et al., 2017, §1, §5.1). RAFT’s update operator is a convolutional GRU.
• Two-dimensional recurrence. PixelRNN’s Row LSTM and Diagonal BiLSTM scan an image row by row or diagonal by diagonal; the autoregressive image model article compares them with masked convolutions.

Use in Computer Vision

Image captioning. Show and Tell (Vinyals et al., 2015) feeds a CNN’s image encoding to an LSTM language model only once, at the first step; the authors report that feeding the image at every step gave inferior results because the network could exploit noise in the image and overfit more easily (§3.1). LRCN’s captioner instead gives a stack of four LSTMs with 1,000 units each the image features at every step, together with the previous word (Donahue et al., 2015, §5). The image captioning article follows these decoders to attention and transformer captioners.

Video recognition and description. LRCN (Donahue et al., 2015, §4) runs a CNN on each of 16 frames and feeds the features to a single-layer LSTM with 256 hidden units that predicts the class at every step; the predictions are averaged over the clip, and at test time over clips taken with a stride of 8 frames. On UCF-101, averaged over the three splits, the LSTM on the first fully connected layer (LRCN-fc6) reached 68.19% with RGB input and 77.46% with optical-flow images, against 67.70% and 72.19% for the same CNN classifying single frames, and 82.66% against 78.84% for the 1:2 weighted average of RGB and flow (Table 1). The temporal modeling article places recurrence among the other ways of combining frames, and video captioning covers recurrent encoder–decoders for video description.

Text recognition. A word or line image can be read as a sequence of columns. CRNN reads the columns of a CNN feature map with a deep bidirectional LSTM and transcribes them with connectionist temporal classification, so no character segmentation is needed; handwriting recognizers built on BLSTMs and multidimensional LSTMs do the same for lines of handwriting.

Video prediction and generation. Recurrent cores, often convolutional LSTMs, summarize past frames in video prediction models, from Srivastava et al.’s LSTM encoder–decoders to Video Pixel Networks and stochastic predictors.

Iterative refinement. RAFT estimates optical flow by repeatedly updating a flow field with a convolutional GRU whose weights are tied across iterations (Teed and Deng, 2020, §3.3). The authors designed the update operator to mimic the steps of an optimization algorithm, with tied weights and bounded activations to encourage convergence to a fixed point. Here the recurrence runs over refinement iterations, not over time in the input.

Limitations

• Sequential computation. Each step waits for the previous one, so training and inference over long sequences cannot be parallelized across time (Goodfellow et al., §10.2). Video and captioning models built on recurrent layers therefore run step by step even on hardware that could process all frames at once.
• Long paths between positions. Information from step $k$ reaches step $t$ only through $t - k$ transitions, and the gradient through a product of as many Jacobians. Gating reduces the decay but does not remove it, as the example below shows at initialization.
• A fixed-size state. Everything the network remembers must fit in $\mathbf{h}_t$ (and $\mathbf{c}_t$). For an encoder–decoder this is a bottleneck, which motivated attention over all encoder states in the Bahdanau attention paper.
• Replacement by attention. The Attention Is All You Need paper, as that article summarizes its Table 1, compares a recurrent layer, with $O(n d^2)$ computation, $O(n)$ sequential operations and a maximum path length of $O(n)$ between positions, against self-attention, with $O(n^2 d)$ computation but $O(1)$ sequential operations and path length. The transformer dropped recurrence altogether. Captioners replaced their LSTM decoders with transformers, as the image captioning article describes, and text recognizers moved from CTC to attention and transformer decoders, as the CRNN article describes.

Example: Gradients Through Time at Initialization

The code below measures how the gradient travels through time in untrained recurrent networks reading an image. The input is a strip of 16 handwritten digits from scikit-learn’s bundled $8 \times 8$ digit images, placed side by side into an $8 \times 128$ image, standardized, and read one column per step, as a text recognizer reads a word ($\mathbf{x}_t \in \mathbb{R}^8$, $T = 128$). No network is trained.

• Simple RNNs. Linear ($\mathbf{h}_t = W\mathbf{h}_{t-1} + U\mathbf{x}_t$) and tanh networks with 64 units, $\mathbf{h}_0 = 0$ and no bias. $W$ is a Gaussian matrix rescaled to spectral radius exactly $\rho \in \{0.5, 0.9, 1.0, 1.2, 3.0\}$; 1.2 is the initial radius Goodfellow et al. (§10.8) report performing well, and 3 the echo-state-network value they cite. $U$ has entries of variance $1/8$, and the linear and tanh networks of a given seed share $W$ and $U$.
• LSTM. 64 cells with PyTorch’s equations and its default initialization, uniform in $\pm 1/\sqrt{64}$, all biases zero except the forget-gate bias $b_f \in \{0, 1\}$ (Jozefowicz et al., §2.2).
• Measures. The largest singular value $s$ of $\partial \mathbf{h}_T / \partial \mathbf{h}_0$, accumulated step by step as $D_T W \cdots D_1 W$, at several $T$; its per-step rate $(s_{128} / s_{32})^{1/96}$ for seed 0 and its range over seeds 0–9; and $\lVert W \rVert_2$, which enters Pascanu et al.’s bound $s \le \lVert W \rVert_2^T$ for tanh. For the LSTM: the full $\partial \mathbf{c}_T / \partial \mathbf{c}_0$, through the gates and $\mathbf{h}$ (derivative with respect to $\mathbf{c}_0$ with $\mathbf{h}_0$ held fixed), the direct cell path $\max_i \prod_t f_{t,i}$, and the geometric-mean forget-gate value.
• Checks. Backpropagation through time written in NumPy, for the loss $\tfrac12 \sum_{t=1}^{127} \lVert V\mathbf{h}_t - \mathbf{x}_{t+1} \rVert^2$ (predict the next column) of the tanh network with $\rho = 1$, is compared with PyTorch’s autograd; the accumulated Jacobians are compared with torch.autograd.functional.jacobian over 8 steps; the NumPy LSTM is compared with torch.nn.LSTMCell; and the LSTM’s Jacobian with gates computed from a detached $\mathbf{h}$ is compared with $\operatorname{diag}(\prod_t \mathbf{f}_t)$.

The script runs in float64 with PyTorch 2.14.1 and NumPy 2.5.4 (scikit-learn 1.9.1 provides the digits) in about 1.5 seconds (approximate, measured while other agents shared the machine); the output was identical over two runs.

import numpy as np
import torch
from sklearn.datasets import load_digits

torch.set_num_threads(1)
torch.set_default_dtype(torch.float64)
print(f"torch {torch.__version__}, numpy {np.__version__}")

# ---- data: a strip of 16 handwritten digits, read column by column ----
digits = load_digits().images[:16]                     # 16 bundled 8x8 images, no download
strip = np.concatenate(list(digits), axis=1)          # 8 x 128 image
X = ((strip - strip.mean()) / strip.std()).T          # X[t] = column t, x_t in R^8, T = 128
T, d = X.shape
n = 64
print(f"data: strip {strip.shape} from {len(digits)} digits; mean {X.mean():.1e}, std {X.std():.12f}")

def make_rnn(rho, seed):
    rng = np.random.default_rng(seed)
    G = rng.standard_normal((n, n)) / np.sqrt(n)
    W = rho * G / np.max(np.abs(np.linalg.eigvals(G)))  # spectral radius exactly rho
    U = rng.standard_normal((n, d)) / np.sqrt(d)
    return W, U

def rnn_jacobian_norms(W, U, linear, Ts=(1, 8, 32, 64, 128)):
    """s(dh_T/dh_0) (largest singular value) for h_t = phi(W h_{t-1} + U x_t), h_0 = 0."""
    h, J, out = np.zeros(n), np.eye(n), {}
    for t in range(1, T + 1):
        a = W @ h + U @ X[t - 1]
        h = a if linear else np.tanh(a)
        D = np.ones(n) if linear else 1 - h ** 2
        J = (D[:, None] * W) @ J                          # dh_t/dh_0 = D_t W dh_{t-1}/dh_0
        if t in Ts:
            out[t] = np.linalg.norm(J, 2)
    return out

sig = lambda z: 1 / (1 + np.exp(-z))

def make_lstm(bf, seed):
    rng = np.random.default_rng(1000 + seed)
    k = 1 / np.sqrt(n)                                     # PyTorch's default U(-1/sqrt(n), 1/sqrt(n))
    Wih = rng.uniform(-k, k, (4 * n, d)); Whh = rng.uniform(-k, k, (4 * n, n))
    b = np.zeros(4 * n); b[n:2 * n] = bf                   # gate order i, f, g, o as in PyTorch
    return Wih, Whh, b

def lstm_run(Wih, Whh, b, Ts=(1, 8, 32, 64, 128), steps=T):
    """Returns s(dc_T/dc_0) (full Jacobian of the (h, c) recurrence), the direct cell path
    max_i prod_t f_{t,i}, the geometric-mean forget gate, and the final (h, c)."""
    h, c = np.zeros(n), np.zeros(n)
    P = np.eye(2 * n); fprod = np.ones(n); logf = 0.0; out, direct = {}, {}
    Wi, Wf, Wg, Wo = np.split(Whh, 4)
    for t in range(1, steps + 1):
        a = Wih @ X[t - 1] + Whh @ h + b
        i, f, g, o = sig(a[:n]), sig(a[n:2*n]), np.tanh(a[2*n:3*n]), sig(a[3*n:])
        c_new = f * c + i * g
        tc = np.tanh(c_new)
        dc_dh = (c * f * (1 - f))[:, None] * Wf + (g * i * (1 - i))[:, None] * Wi + (i * (1 - g**2))[:, None] * Wg
        dc_dc = np.diag(f)
        dh_dh = (tc * o * (1 - o))[:, None] * Wo + (o * (1 - tc**2))[:, None] * dc_dh
        dh_dc = (o * (1 - tc**2) * f)[:, None] * np.eye(n)
        Jt = np.block([[dh_dh, dh_dc], [dc_dh, dc_dc]])
        P = Jt @ P
        h, c = o * tc, c_new
        fprod *= f; logf += np.log(f).mean()
        if t in Ts:
            out[t] = np.linalg.norm(P[n:, n:], 2)          # dc_t/dc_0
            direct[t] = fprod.max()
    return out, direct, np.exp(logf / steps), h, c, P, fprod

Ts = (1, 8, 32, 64, 128)
rate = lambda o: (o[128] / o[32]) ** (1 / 96)
print("\nmodel            ||W||_2 | s(dh_T/dh_0) at T =    1        8       32       64      128 | rate/step  seeds 0-9")
rhos = [0.5, 0.9, 1.0, 1.2, 3.0]
bound_ok, radius_err, lin_rate_err = True, 0.0, 0.0
for kind in ["linear", "tanh"]:
    for rho in rhos:
        W, U = make_rnn(rho, 0)
        o = rnn_jacobian_norms(W, U, kind == "linear")
        s2 = np.linalg.norm(W, 2)
        rs = []
        for seed in range(10):
            Ws, Us = make_rnn(rho, seed)
            radius_err = max(radius_err, abs(np.max(np.abs(np.linalg.eigvals(Ws))) - rho))
            os_ = rnn_jacobian_norms(Ws, Us, kind == "linear")
            rs.append(rate(os_))
            if kind == "tanh":
                s2s = np.linalg.norm(Ws, 2)
                bound_ok &= all(os_[t] <= s2s ** t * (1 + 1e-12) for t in Ts)
            else:
                lin_rate_err = max(lin_rate_err, abs(rate(os_) / rho - 1))
        vals = " ".join(f"{o[t]:8.1e}" for t in Ts)
        print(f"{kind:6s} rho {rho:3.1f} {s2:9.2f} | {' ' * 22}{vals} | {rate(o):9.3f}  {min(rs):.3f}-{max(rs):.3f}")

print("\nLSTM (64 cells)  | s(dc_T/dc_0) at T =    1        8       32       64      128 | rate/step  seeds 0-9 | direct-path rate  geo-mean f")
for bf in [0.0, 1.0]:
    o, dr, gm, *_ = lstm_run(*make_lstm(bf, 0))
    rs, drs = [], []
    for seed in range(10):
        os_, ds_, *_ = lstm_run(*make_lstm(bf, seed))
        rs.append(rate(os_)); drs.append(rate(ds_))
    vals = " ".join(f"{o[t]:8.1e}" for t in Ts)
    print(f"forget bias {bf:3.1f}  | {' ' * 19}{vals} | {rate(o):9.3f}  {min(rs):.3f}-{max(rs):.3f} | {rate(dr):16.3f} {gm:11.3f}")

# ---- checks against PyTorch autograd ----
print("\nchecks (relative error unless stated)")
Xt = torch.tensor(X)
# (1) BPTT in NumPy vs autograd: tanh RNN, rho = 1, loss = 1/2 sum_t ||V h_t - x_{t+1}||^2
W, U = make_rnn(1.0, 0)
V = np.random.default_rng(7).standard_normal((d, n)) / np.sqrt(n)
h0 = np.zeros(n)
H = [h0]
for t in range(T):
    H.append(np.tanh(W @ H[-1] + U @ X[t]))
gW, gU, gV, dh_next = np.zeros_like(W), np.zeros_like(U), np.zeros_like(V), np.zeros(n)
for t in range(T, 0, -1):                                   # backpropagation through time
    dh = dh_next.copy()
    if t <= T - 1:
        e = V @ H[t] - X[t]                                 # X[t] is x_{t+1}
        gV += np.outer(e, H[t]); dh += V.T @ e
    da = dh * (1 - H[t] ** 2)
    gW += np.outer(da, H[t - 1]); gU += np.outer(da, X[t - 1]); dh_next = W.T @ da
gh0 = dh_next
Wt, Ut, Vt, h0t = (torch.tensor(a, requires_grad=True) for a in (W, U, V, h0))
h, L = h0t, 0.0
for t in range(T):
    h = torch.tanh(Wt @ h + Ut @ Xt[t])
    if t + 1 <= T - 1:
        L = L + 0.5 * ((Vt @ h - Xt[t + 1]) ** 2).sum()
L.backward()
for name, g_np, g_t in [("W", gW, Wt.grad), ("U", gU, Ut.grad), ("V", gV, Vt.grad), ("h_0", gh0, h0t.grad)]:
    print(f"  BPTT vs autograd, dL/d{name:4s}: {np.abs(g_np - g_t.numpy()).max() / np.abs(g_t.numpy()).max():.1e}")

# (2) Jacobian accumulation vs autograd, 8 steps
Xs8 = Xt[:8]
def tanh_steps(h):
    for t in range(8):
        h = torch.tanh(torch.tensor(W) @ h + torch.tensor(U) @ Xs8[t])
    return h
Jt = torch.autograd.functional.jacobian(tanh_steps, torch.zeros(n)).numpy()
hh, Jn = np.zeros(n), np.eye(n)
for t in range(8):
    hh = np.tanh(W @ hh + U @ X[t]); Jn = ((1 - hh ** 2)[:, None] * W) @ Jn
print(f"  tanh dh_8/dh_0 NumPy vs autograd: {np.abs(Jn - Jt).max() / np.abs(Jt).max():.1e}")

Wih, Whh, b = make_lstm(1.0, 0)
cell = torch.nn.LSTMCell(d, n)
with torch.no_grad():
    cell.weight_ih.copy_(torch.tensor(Wih)); cell.weight_hh.copy_(torch.tensor(Whh))
    cell.bias_ih.copy_(torch.tensor(b)); cell.bias_hh.zero_()
def lstm_c(c0, detach_gates=False):
    h, c = torch.zeros(n), c0
    for t in range(8):
        hg = h.detach() if detach_gates else h
        a = torch.tensor(Wih) @ Xs8[t] + torch.tensor(Whh) @ hg + torch.tensor(b)
        i, f, g, o = torch.sigmoid(a[:n]), torch.sigmoid(a[n:2*n]), torch.tanh(a[2*n:3*n]), torch.sigmoid(a[3*n:])
        c = f * c + i * g
        h = o * torch.tanh(c)
    return c
_, _, _, h8, c8, P8, fprod8 = lstm_run(Wih, Whh, b, Ts=(), steps=8)
Jc = torch.autograd.functional.jacobian(lstm_c, torch.zeros(n)).numpy()
print(f"  LSTM dc_8/dc_0 NumPy vs autograd: {np.abs(P8[n:, n:] - Jc).max() / np.abs(Jc).max():.1e}")
Jd = torch.autograd.functional.jacobian(lambda c0: lstm_c(c0, True), torch.zeros(n)).numpy()
print(f"  LSTM direct path (gates detached) vs diag(prod f): {np.abs(Jd - np.diag(fprod8)).max() / fprod8.max():.1e}")
with torch.no_grad():
    ht, ct = torch.zeros(1, n), torch.zeros(1, n)
    for t in range(8):
        ht, ct = cell(Xs8[t][None], (ht, ct))
print(f"  NumPy LSTM vs torch.nn.LSTMCell, max |diff| in (h_8, c_8): {max(np.abs(h8 - ht[0].numpy()).max(), np.abs(c8 - ct[0].numpy()).max()):.1e}")
print(f"  spectral radius of W vs rho, max |diff| (all rows, seeds 0-9): {radius_err:.1e}")
print(f"  linear rate vs rho, max relative difference (seeds 0-9): {lin_rate_err:.1e}")
print(f"  tanh: s(dh_T/dh_0) <= ||W||_2^T for every rho, T and seed: {bound_ok}")

Output:

torch 2.14.1, numpy 2.5.4
data: strip (8, 128) from 16 digits; mean 8.7e-18, std 1.000000000000

model            ||W||_2 | s(dh_T/dh_0) at T =    1        8       32       64      128 | rate/step  seeds 0-9
linear rho 0.5      0.98 |                        9.8e-01  1.1e-02  5.9e-10  1.5e-19  7.8e-39 |     0.500  0.500-0.501
linear rho 0.9      1.77 |                        1.8e+00  1.2e+00  8.6e-02  3.2e-03  3.7e-06 |     0.900  0.899-0.901
linear rho 1.0      1.97 |                        2.0e+00  2.7e+00  2.5e+00  2.7e+00  2.7e+00 |     1.001  0.999-1.001
linear rho 1.2      2.36 |                        2.4e+00  1.2e+01  8.6e+02  3.1e+05  3.6e+10 |     1.201  1.199-1.201
linear rho 3.0      5.90 |                        5.9e+00  1.8e+04  4.7e+15  9.2e+30  3.1e+61 |     3.002  2.998-3.003
tanh   rho 0.5      0.98 |                        7.8e-01  1.0e-03  6.1e-15  2.0e-30  2.6e-61 |     0.329  0.302-0.329
tanh   rho 0.9      1.77 |                        1.4e+00  8.2e-02  3.1e-07  2.4e-15  1.8e-31 |     0.559  0.524-0.567
tanh   rho 1.0      1.97 |                        1.6e+00  1.7e-01  6.4e-06  6.8e-13  1.4e-26 |     0.609  0.584-0.623
tanh   rho 1.2      2.36 |                        1.9e+00  5.1e-01  6.4e-04  3.0e-09  1.5e-18 |     0.704  0.672-0.717
tanh   rho 3.0      5.90 |                        4.7e+00  6.7e+01  6.8e+04  2.2e+08  6.3e+14 |     1.270  1.203-1.290

LSTM (64 cells)  | s(dc_T/dc_0) at T =    1        8       32       64      128 | rate/step  seeds 0-9 | direct-path rate  geo-mean f
forget bias 0.0  |                     6.1e-01  3.6e-02  1.6e-06  1.9e-12  1.2e-24 |     0.648  0.620-0.649 |            0.507       0.496
forget bias 1.0  |                     8.1e-01  4.0e-01  2.4e-02  4.2e-04  8.2e-08 |     0.877  0.852-0.879 |            0.737       0.727

checks (relative error unless stated)
  BPTT vs autograd, dL/dW   : 4.2e-16
  BPTT vs autograd, dL/dU   : 1.6e-16
  BPTT vs autograd, dL/dV   : 2.1e-16
  BPTT vs autograd, dL/dh_0 : 4.0e-16
  tanh dh_8/dh_0 NumPy vs autograd: 9.2e-16
  LSTM dc_8/dc_0 NumPy vs autograd: 4.4e-16
  LSTM direct path (gates detached) vs diag(prod f): 4.3e-16
  NumPy LSTM vs torch.nn.LSTMCell, max |diff| in (h_8, c_8): 5.9e-17
  spectral radius of W vs rho, max |diff| (all rows, seeds 0-9): 2.7e-14
  linear rate vs rho, max relative difference (seeds 0-9): 1.1e-03
  tanh: s(dh_T/dh_0) <= ||W||_2^T for every rho, T and seed: True

Linear networks follow Pascanu et al.’s linear analysis: the per-step rate equals $\rho$ to within about 0.1% for every seed, so the gradient vanishes for $\rho < 1$ and explodes for $\rho > 1$. At $\rho = 1$ the norm does not stay at 1 but settles near 2.7: a random $W$ is not symmetric, its largest singular value is about twice its spectral radius (the $\lVert W \rVert_2$ column), and the product can amplify before the spectral radius takes over. The spectral radius sets the long-run rate; the singular values set the short-run size.

Tanh networks decay faster than linear ones at every radius. In our reading of the construction, each factor $D_t$ scales rows by $1 - h^2 \le 1$, and the input drives some units toward saturation, where that factor is small. At $\rho = 1.2$ the gradient across 128 columns still vanishes, at about 0.7 per step on this input. At $\rho = 3$ it grows by about 1.27 per step over the 96 steps measured (1.20–1.29 over seeds), reaching $6 \times 10^{14}$ after 128 steps, although every state stays inside $[-1, 1]$. Goodfellow et al.’s §10.8 statement concerns the forward dynamics of echo state networks, which saturation keeps bounded; the same section notes that backpropagation can retain unbounded dynamics while the forward pass is bounded, and here it does. The run is consistent with Pascanu et al.’s necessary condition: explosion requires $\lambda_1 > 1$ for tanh, and here $\lambda_1 = 5.9$. Their sufficient bound $s \le \lVert W \rVert_2^T$ holds everywhere but is far from tight; at $\rho = 1$ it allows growth by $1.97^{128}$ where the measured value is $10^{-26}$.

LSTM. The direct cell path decays at the forget gate’s own rate: 0.51 per step with $b_f = 0$, where the gate averages 0.50, Jozefowicz et al.’s factor of 0.5; and 0.74 with $b_f = 1$, close to $\sigma(1) \approx 0.73$. The full Jacobian $\partial \mathbf{c}_T / \partial \mathbf{c}_0$ decays more slowly than the direct path (0.65 and 0.88 per step), because $\mathbf{c}_0$ also reaches $\mathbf{c}_T$ through $\mathbf{h}$ and the gates. With the bias, the cell gradient after 128 steps is $8 \times 10^{-8}$ instead of $10^{-24}$, and it decays more slowly than in any of the tanh networks with $\rho \le 1.2$ (against $10^{-26}$ for tanh at $\rho = 1$). It still vanishes. In our reading, a randomly initialized LSTM does not yet hold information over long spans; that requires forget gates near 1, which training must produce, and the bias only moves the starting point closer.

The checks agree to round-off ($10^{-16}$): the NumPy implementation of backpropagation through time gives the same gradients as autograd for every parameter and for $\mathbf{h}_0$, the accumulated Jacobians match autograd’s, and the hand-written LSTM matches PyTorch’s cell. The detached-gate Jacobian equals $\operatorname{diag}(\prod_t \mathbf{f}_t)$ to round-off, which is the direct cell path isolated by stopping the gradient through the gates.

The measurements are of untrained networks on one input strip; trained networks change their Jacobians as their weights move, as the vanishing gradient problem article notes for feedforward networks.

Important Models

Model	Year	Recurrent unit	Use
LSTM (Hochreiter and Schmidhuber; forget gate Gers et al.)	1997, 2000	Gated cell with additive state	Sequences with long-range dependencies
RNN Encoder–Decoder (Cho et al.)	2014	Reset and update gates (GRU)	Phrase scoring for machine translation
LRCN (Donahue et al.)	2015	LSTM on CNN features	Activity recognition, captioning, video description
Show and Tell (Vinyals et al.)	2015	LSTM language model	Image captioning
CRNN (Shi, Bai and Yao)	2015	Deep bidirectional LSTM	Scene text recognition with CTC
ConvLSTM (Shi et al.)	2015	LSTM with convolutional transitions	Spatiotemporal prediction
PixelRNN (van den Oord et al.)	2016	Row LSTM, Diagonal BiLSTM	Autoregressive image modelling
Video Pixel Networks (Kalchbrenner et al.)	2017	Resolution-preserving ConvLSTM	Video likelihood modelling
RAFT (Teed and Deng)	2020	Convolutional GRU over iterations	Optical flow refinement