57 - Backpropagation Through Time: Vanishing and Exploding Gradients

An RNN is trained like any network: you compute the error, carry it backwards with the chain rule and update the weights. The trick is that the network, unrolled in time, is as deep as the sequence is long, so the gradient has to cross T identical steps. The gradient is a product of T matrices: if these matrices are in some sense "smaller than 1" the product tends to zero (vanishing gradient: the network does not learn distant dependencies); if they are "larger than 1" it can grow without control (exploding gradient: training becomes unstable). This article derives the formula, checks it numerically, measures both phenomena and shows the remedy for the explosion, gradient clipping, together with what it does not cure. It is the sequel to Recurrent neural networks (RNN).

The words you need

The loss L is a single number that measures how wrong the network was, and training means making it go down. The derivative tells you how much a number changes if you nudge another one by a tiny amount: if ∂L/∂w = 3, raising w by 0.01 raises the loss by about 0.03. The gradient is the set of all the derivatives of the loss with respect to all the weights, and the update is w ← w − η · ∂L/∂w, where η is the learning rate, the length of the step.

The chain rule says that if A influences B and B influences C, the effect of A on C is the product of the two effects: in a chain of gears, if the first turns once and the second makes two turns, and the third makes three for every turn of the second, the third makes six turns; the factors multiply. When input and output are vectors the "derivative" is a matrix, the Jacobian: row i, column j tells you how much component i of the output changes as component j of the input varies; here we need one per step, ∂h_t/∂h_{t-1}. The norm ‖v‖ is the "length" of a vector (or the "size" of a matrix) and lets you say with a single number whether a gradient is large or small. Backpropagation is the algorithm that applies the chain rule backwards, from the loss to the weights, reusing the intermediate computations.

The analogy: the telephone game

Ten people in a line pass a message by whispering it. If each one attenuates it a little, silence reaches the last one: that is the vanishing gradient. If each one amplifies it a little, a distorted shout reaches the last one: that is the exploding gradient. In the network the "message" is the error signal that has to reach a step far back in time to tell it "correct your contribution", and every step multiplies it by a matrix.

The maths: deriving BPTT

Take the simplest case, many-to-one, with a single loss at the end. With a_t = W_xh x_t + W_hh h_{t-1} and h_t = tanh(a_t) the loss is:

L = ½ (w · h_T − y)²

The weight W_hh appears at every step. By the chain rule its effect on the loss is the sum of the contributions of every step t:

∂L/∂W_hh = Σ_{t=1..T}  (∂L/∂h_T) · (∂h_T/∂h_t) · (∂h_t/∂W_hh)|direct

Here "direct" means the derivative of h_t with respect to W_hh holding h_{t-1} fixed. The critical piece is ∂h_T/∂h_t, that is, how much a small variation of the state at time t reaches all the way to the end. Since h_j depends on h_{j-1}, it is a product of Jacobians:

∂h_T/∂h_t = Π_{j=t+1..T}  ∂h_j/∂h_{j-1}
∂h_j/∂h_{j-1} = diag(1 − h_j²) · W_hh        (derivative of the tanh: 1 − tanh²)

Let us take norms. The derivative of the tanh is at most 1, and is almost 0 when the tanh is saturated, that is, when |h| is close to 1. Call γ the maximum of the derivative (here 1) and σ the largest singular value of W_hh; each factor has norm at most γσ and therefore:

‖∂h_T/∂h_t‖ ≤ (γ · σ)^(T − t)

From here two regimes follow, analysed by Bengio, Simard and Frasconi (1994) and in a clean form by Pascanu, Mikolov and Bengio (2013). If γσ < 1 the product decays exponentially with the distance T − t, and the gradient of the distant steps is practically zero: it is a guarantee, it always happens. If σ is large the product can grow exponentially, but that is a necessary and not a sufficient condition, because it depends on the directions and on the saturation of the tanh.

The algorithm comes from reading the formula from right to left. Starting from dh = ∂L/∂h_T, at every step backwards you go through the tanh (da = dh · (1 − h²)), you accumulate the contribution to W_hh (da · h_{t-1}ᵀ) and you send the gradient to the previous step (dh = W_hhᵀ · da). This is backpropagation through time (BPTT).

When there is a loss at every step

So far the loss was a single one, at the end. In many-to-many tasks (labelling every word, predicting every next value) every step has its own loss L_t and the total is L = L_1 + … + L_T. The derivative with respect to W_hh becomes a double sum:

∂L/∂W_hh = Σ_{t=1..T}  Σ_{k=1..t}  (∂L_t/∂h_t) · (∂h_t/∂h_k) · (∂h_k/∂W_hh)|direct

The error made at time t sends a signal back towards all the steps k ≤ t, and every step k contributes to the gradient of the shared weight; the term ∂h_t/∂h_k is again the product of Jacobians, with t − k factors. Two practical consequences follow. The "near" contributions (small t − k) arrive strong and the distant ones weak, so the network learns short dependencies first and long ones last, or never. The cost of full BPTT grows with T, and so does memory, because every intermediate state has to be kept until the backward pass.

Why the bound on the norm is true

The inequality ‖∂h_T/∂h_t‖ ≤ (γσ)^(T−t) is proved in two lines. For any matrix A and vector v we have ‖A v‖ ≤ ‖A‖ · ‖v‖, where ‖A‖ is the spectral norm, the maximum stretch that A can produce; it follows that the norm of a product is at most the product of the norms, ‖A B‖ ≤ ‖A‖ · ‖B‖. Each Jacobian is diag(1 − h²) · W_hh: the diagonal has norm at most γ (the largest value of 1 − h², that is, 1) and W_hh has norm σ.

So each factor has norm at most γσ, and a product of T − t factors has norm at most (γσ)^(T−t). If γσ < 1 this tends to zero exponentially, and it is a mathematical guarantee, not a tendency. If instead γσ > 1 the inequality says nothing, because it is only an upper bound and the true product can stay small: that is why the explosion is a necessary and not a sufficient condition.

In the drawing the red arrows go from right to left: the error signal is born in the loss and has to reach h1 through four multiplications. If each factor is 0.5, on arrival the signal is multiplied by 0.5⁴ = 0.0625; with a hundred steps it would become 0.5¹⁰⁰ ≈ 8·10⁻³¹, a number the computer represents without problems but that for the optimiser is equivalent to zero.

Code by hand: BPTT with a numerical check

The safest way to know whether a derivation is right is to compare it with finite differences, (L(w+ε) − L(w−ε)) / 2ε. The listing implements BPTT for the many-to-one network with three neurons and six steps, then nudges every element of W_hh by ±10⁻⁶ and compares.

import numpy as np

rng = np.random.default_rng(1)
n, d, T = 3, 2, 6
Wxh, Whh = rng.normal(0, .6, (n, d)), rng.normal(0, .6, (n, n))
w_out = rng.normal(0, .6, n)
xs, y = rng.normal(size=(T, d)), 0.7

def forward(Whh):
    hs = [np.zeros(n)]
    for x in xs:
        hs.append(np.tanh(Wxh @ x + Whh @ hs[-1]))
    loss = 0.5 * (w_out @ hs[-1] - y) ** 2
    return loss, hs

def bptt(Whh):
    loss, hs = forward(Whh)
    dh = (w_out @ hs[-1] - y) * w_out           # dL/dh_T
    dWhh = np.zeros_like(Whh)
    for t in range(T, 0, -1):                   # walk BACKWARDS in time
        da = dh * (1 - hs[t] ** 2)              # through the tanh
        dWhh += np.outer(da, hs[t - 1])         # this step's contribution to dL/dWhh
        dh = Whh.T @ da                         # pass the gradient to h_{t-1}
    return dWhh

# Check against finite differences
analytic = bptt(Whh)
numeric = np.zeros_like(Whh)
eps = 1e-6
for i in range(n):
    for j in range(n):
        P, M = Whh.copy(), Whh.copy()
        P[i, j] += eps; M[i, j] -= eps
        numeric[i, j] = (forward(P)[0] - forward(M)[0]) / (2 * eps)
print(np.abs(analytic - numeric).max())    # ~1e-10: the manual BPTT is correct

# Output:
# 2.942041055220557e-11

forward is the forward pass we have already seen, with the final loss ½(w·h_T − y)² added, and it returns all the states because the backward pass needs them. bptt starts from the derivative of the loss with respect to the last state (dh) and walks from t = T down to t = 1. Inside the loop are the three lines of the derivation: da = dh * (1 - hs[t]**2) goes through the tanh, dWhh += np.outer(da, hs[t-1]) accumulates this step's contribution to the shared weight (it sums because the same weight is used at every step; with = in place of += the gradient would be that of the first step backwards only) and dh = Whh.T @ da sends the gradient to the previous step.

The maximum error is about 3·10⁻¹¹: the hand-written BPTT coincides with the numerical derivative, so the derivation is right. The loop that runs t backwards also makes the cost clear: memory grows with T because all the states hs have to be kept.

The smallest case: a single neuron

Before matrices, a single neuron with no input, h_t = tanh(w · h_{t−1}). The derivative of h_T with respect to h_0 is the product of the factors w · (1 − h_t²), one per step: the chain rule in numerical form. The function scalar_gradient multiplies one factor per iteration and the final loop tries it with four weights at 1, 5, 10 and 20 steps.

import numpy as np

# One-neuron RNN without input:  h_t = tanh(w * h_{t-1}).
# d h_T / d h_0 is the product of the per-step factors  w * (1 - h_t^2).
def scalar_gradient(w, steps, h0=0.5):
    h, grad = h0, 1.0
    for _ in range(steps):
        h = np.tanh(w * h)
        grad *= w * (1 - h ** 2)         # chain rule: multiply one factor per step
    return grad

for w in (0.5, 0.9, 1.1, 2.0):
    print(f"w={w}:", [f"{scalar_gradient(w, s):.2e}" for s in (1, 5, 10, 20)])

# Output:
# w=0.5: ['4.70e-01', '2.88e-02', '9.00e-04', '8.79e-07']
# w=0.9: ['7.40e-01', '3.30e-01', '1.66e-01', '5.42e-02']
# w=1.1: ['8.24e-01', '3.79e-01', '1.43e-01', '2.01e-02']
# w=2.0: ['8.40e-01', '1.68e-03', '2.14e-07', '3.48e-15']

A surprising result: even with w = 2.0, which "looks" like it amplifies, the gradient vanishes (from 0.84 to 3.5·10⁻¹⁵). The reason is saturation: with a large w the state is pushed towards ±1, where 1 − h² is almost zero, and the factor w(1 − h²) falls below 1. In this single-neuron experiment with tanh the explosion does not appear for any of the four weights; to explode you need matrices with favourable directions (W_hh large in some direction while the state is not saturated). It is one reason why in practice vanishing is always seen, while exploding appears now and then, as sudden spikes.

Seeing the decay: the norm of the product of Jacobians

Let us measure ‖∂h_T/∂h_0‖ directly for a recurrent network with a random matrix rescaled so that the largest eigenvalue in modulus is 0.9 or 1.1, both linear and with tanh. jacobian_norm builds W and divides it by its largest eigenvalue in modulus (the "spectral radius"), so that it equals exactly radius: 0.9 contracts slightly, 1.1 expands slightly. Then it simulates T steps and progressively multiplies the Jacobian of the step (diag(slope) @ W) into the product J; with activation="linear" we remove the tanh to isolate the effect of the matrix alone. At the end it returns np.linalg.norm(J, 2), the 2-norm, that is, the maximum factor by which J amplifies a vector: the right measure for "how much signal gets through".

import numpy as np

def jacobian_norm(T, n=16, radius=1.0, activation="tanh", seed=0):
    """Spectral norm of d h_T / d h_0 for a recurrent net h_t = f(W h_{t-1})."""
    rng = np.random.default_rng(seed)
    W = rng.normal(size=(n, n))
    W = radius * W / max(abs(np.linalg.eigvals(W)))     # rescale: largest |eigenvalue| = radius
    h = rng.normal(size=n) * 0.5
    J = np.eye(n)                                        # running product of Jacobians
    for _ in range(T):
        h = np.tanh(W @ h) if activation == "tanh" else W @ h
        slope = 1 - h**2 if activation == "tanh" else np.ones(n)   # f'(a)
        J = np.diag(slope) @ W @ J                       # chain rule, one step
    return np.linalg.norm(J, 2)

for act in ("linear", "tanh"):
    for radius in (0.9, 1.1):
        vals = [f"{jacobian_norm(T, radius=radius, activation=act):.1e}" for T in (1, 10, 50, 100)]
        print(f"{act:6s} radius={radius}", vals)

# Output:
# linear radius=0.9 ['1.4e+00', '4.4e-01', '6.6e-03', '3.4e-05']
# linear radius=1.1 ['1.7e+00', '3.3e+00', '1.5e+02', '1.8e+04']
# tanh   radius=0.9 ['1.2e+00', '3.6e-01', '5.3e-03', '2.7e-05']
# tanh   radius=1.1 ['1.5e+00', '2.0e+00', '1.4e-02', '4.2e-07']

In the linear case a spectral radius just below 1 makes the gradient disappear in 100 steps (from 1.4 to 3.4·10⁻⁵), and just above 1 makes it explode (1.8·10⁴). With the tanh something instructive happens: even with radius 1.1 the gradient first grows (2.0 at T=10) and then vanishes (4.2·10⁻⁷ at T=100), because the states saturate towards ±1 and the derivative 1 − h² approaches zero. The tanh brakes the explosion, but in the intervals where the state is not saturated the gradient can still rise. It is a single seed, so the order of magnitude counts and not the exact figure. In practice both things are observed: gradient spikes alternating with long flat stretches.

Eigenvalues, and why orthogonal initialisation helps

In the linear case the computation is exact. If h_t = W h_{t−1}, then h_T = W^T h_0. A "normal" matrix can be written with its eigenvalues λ_i (the factors by which it stretches some special directions): after T steps the component along direction i is multiplied by λ_i^T. With |λ| = 0.9 we get 0.9^100 = 2.7·10⁻⁵, with |λ| = 1.1 we get 1.1^100 = 1.4·10⁴: two numbers only apparently close to 1 produce nine orders of magnitude of difference, because T is an exponent.

The suggested remedy is to choose W_hh so that all its "stretch factors" are exactly 1. An orthogonal matrix (QᵀQ = I) preserves the length of every vector. The listing checks it: np.linalg.qr produces an orthogonal Q from a random matrix (every square matrix decomposes into an orthogonal times a triangular one), we check that ‖Qx‖ = ‖x‖, we multiply by Q a thousand times and compare with a random matrix of spectral radius 1.

import numpy as np

rng = np.random.default_rng(0)
n = 16

# 1) the eigenvalue picture, in the LINEAR case: h_T = W^T h_0
W = rng.normal(size=(n, n))
W = 0.9 * W / max(abs(np.linalg.eigvals(W)))
print("0.9^100 =", 0.9 ** 100, "   1.1^100 =", 1.1 ** 100)

# 2) an orthogonal matrix keeps the length of every vector: |Q x| = |x|
Q, _ = np.linalg.qr(rng.normal(size=(n, n)))          # QR gives an orthogonal Q
x = rng.normal(size=n)
print("|x| =", np.linalg.norm(x).round(6), "  |Qx| =", np.linalg.norm(Q @ x).round(6))

# 3) the product of Jacobians of a LINEAR recurrent net with orthogonal W: norm stays 1 forever
J = np.eye(n)
for _ in range(1000):
    J = Q @ J
print("norm of Q^1000:", np.linalg.norm(J, 2).round(6))

# versus a random matrix rescaled to spectral radius 1.0 (not normal -> norms can still drift)
R = rng.normal(size=(n, n)); R = R / max(abs(np.linalg.eigvals(R)))
J = np.eye(n)
for T in (10, 100, 1000):
    J = np.linalg.matrix_power(R, T)
    print(f"random, radius 1: norm of R^{T} = {np.linalg.norm(J, 2):.3e}")

# Output:
# 0.9^100 = 2.6561398887587544e-05    1.1^100 = 13780.61233982238
# |x| = 4.320329   |Qx| = 4.320329
# norm of Q^1000: 1.0
# random, radius 1: norm of R^10 = 1.803e+00
# random, radius 1: norm of R^100 = 1.467e+00
# random, radius 1: norm of R^1000 = 1.651e+00

The first line shows the effect of the exponent. After a thousand multiplications by Q the norm stays 1.0: no decay and no explosion, however many steps. The random matrix with spectral radius 1 is not as clean, with norms around 1.5-1.8: non-normal matrices can amplify temporarily even with spectral radius 1. That is why a common practice is to initialise the recurrent matrix as orthogonal (in PyTorch there is torch.nn.init.orthogonal_; here we did not use it).

An honest caveat: orthogonality holds for the product of matrices. The true Jacobian also contains diag(1 − h²), which is less than 1, so with the tanh the decay comes back (slowly, but it comes back); and during training the weights drift away from orthogonality. It is a good starting point, not a guarantee.

The remedy for the explosion: gradient clipping

If the gradient has an enormous norm, a single update can throw the weights very far away and ruin training (the loss becomes nan). Clipping by norm, proposed by Pascanu et al. (2013), keeps the direction but limits the length: g ← g · min(1, τ / ‖g‖), where τ is a threshold to be chosen by observing the norms that appear during training.

import numpy as np

def clip_by_global_norm(grads, max_norm):
    """Rescale all gradients together if their joint norm exceeds max_norm."""
    total = np.sqrt(sum((g ** 2).sum() for g in grads))
    factor = min(1.0, max_norm / (total + 1e-6))
    return [g * factor for g in grads], total

grads = [np.array([30.0, -40.0]), np.array([0.0, 0.0])]     # exploding: norm = 50
clipped, before = clip_by_global_norm(grads, max_norm=5.0)
print(before, clipped[0])                                    # 50.0, about [3, -4]: same direction, norm 5

# Output:
# 50.0 [ 2.99999994 -3.99999992]

The function computes the global norm (the square root of the sum of squares of all the gradients, as if they were one long vector) and rescales everything by the same factor min(1, threshold/norm). Two choices are worth noting. The factor never exceeds 1, so small gradients are not enlarged. A global norm is used instead of clipping each tensor separately, so the proportions among the gradients of the different weights stay intact and the descent direction does not change. With norm 50 and threshold 5 everything is multiplied by 0.1: the output is (3, −4) with an uncertainty of 10⁻⁸ due to the 1e-6 term added to avoid division by zero.

In PyTorch it is one line between backward() and step(): torch.nn.utils.clip_grad_norm_(params, max_norm), which returns the norm before clipping (useful to log to see whether the threshold triggers often, as you can read in the source). There is also clip_grad_value_, which clips each component separately: it does something different and does not preserve the direction of the gradient. We ran them (PyTorch 2.12, CPU, 2 threads, twice with the same output) on the same two gradients as the NumPy example, global norm 50.

import torch

torch.set_num_threads(2)

def gradients():
    """The same two gradients as the NumPy example: global norm 50."""
    p = [torch.zeros(2, requires_grad=True), torch.zeros(2, requires_grad=True)]
    p[0].grad = torch.tensor([30.0, -40.0])
    p[1].grad = torch.tensor([0.0, 0.0])
    return p

p = gradients()
before = torch.nn.utils.clip_grad_norm_(p, max_norm=5.0)    # returns the norm BEFORE clipping
g = p[0].grad
print("clip_grad_norm_: norm before =", before.item(), " after =", round(g.norm().item(), 6))
print("  gradient:", [round(v, 6) for v in g.tolist()], " direction:", [round(v, 4) for v in (g / g.norm()).tolist()])

p = gradients()
ret = torch.nn.utils.clip_grad_value_(p, clip_value=5.0)   # clips every component to [-5, 5]
g = p[0].grad
print("clip_grad_value_: returns", ret, " norm after =", round(g.norm().item(), 6))
print("  gradient:", g.tolist(), " direction:", [round(v, 4) for v in (g / g.norm()).tolist()])

# Output:
# clip_grad_norm_: norm before = 50.0  after = 5.0
#   gradient: [3.0, -4.0]  direction: [0.6, -0.8]
# clip_grad_value_: returns None  norm after = 7.071068
#   gradient: [5.0, -5.0]  direction: [0.7071, -0.7071]

What the output shows. clip_grad_norm_ returns 50.0, the norm before clipping, and brings the gradient to (3, −4): norm 5 and the direction (0.6, −0.8) unchanged, like our function (up to 32-bit rounding). clip_grad_value_ with threshold 5 brings the same gradient to (5, −5): the norm becomes 7.07, that is, above the threshold, and the direction (0.7071, −0.7071) is no longer the original one; it also returns None, so it does not give you a norm to log. The same 30 and −40 produce two different vectors because the two clippings answer different questions: one limits the length, the other limits each component.

The vanishing gradient is not cured by clipping

Clipping only cuts norms that are too large: if the gradient is already small, scaling it does not bring it back to life. Let us measure it. The network must remember the first element of a random sequence and return it at the end, after T − 1 steps of noise (y = x[:, 0]). Always predicting 0 gives a mean squared error of about 1, the variance of a standard normal: "1" means "I have learned nothing". Every step generates a fresh batch, so the network cannot memorise examples.

The listing writes BPTT by hand and the Adam optimiser in full (it is the same structure as the previous article, with 32 units). The score is not the loss of the last batch, which is noisy, but the error on 2000 new sequences; and the experiment is repeated over 5 seeds, because with a single seed you could tell any story. Clipping is switched off (clip=None) in the first three lines, to isolate the vanishing. The last line redoes the hardest case with clip=1.0.

import numpy as np

def train(T, steps=400, clip=None, seed=0, hidden=32):
    """RNN(1 -> hidden) + linear head, target = FIRST element of the sequence, BPTT by hand, Adam."""
    rng = np.random.default_rng(seed)
    k = 1 / np.sqrt(hidden)
    p = {"Wx": rng.uniform(-k, k, hidden), "Wh": rng.uniform(-k, k, (hidden, hidden)),
         "b": np.zeros(hidden), "w": rng.uniform(-k, k, hidden), "c": np.zeros(())}
    m = {n: np.zeros_like(v) for n, v in p.items()}
    s = {n: np.zeros_like(v) for n, v in p.items()}
    lr, b1, b2, eps = 1e-2, 0.9, 0.999, 1e-8
    for step in range(1, steps + 1):
        x = rng.standard_normal((64, T))
        y = x[:, 0]                                     # target = FIRST element: needs T-1 steps of memory
        hs = [np.zeros((64, hidden))]
        for t in range(T):
            hs.append(np.tanh(x[:, t:t + 1] * p["Wx"] + hs[-1] @ p["Wh"].T + p["b"]))
        err = hs[-1] @ p["w"] + p["c"] - y
        loss = (err ** 2).mean()
        g = {n: np.zeros_like(v) for n, v in p.items()}
        dout = 2 * err / len(err)
        g["w"], g["c"] = hs[-1].T @ dout, dout.sum()
        dh = dout[:, None] * p["w"]
        for t in range(T, 0, -1):                       # BPTT: T steps backwards
            da = dh * (1 - hs[t] ** 2)
            g["Wx"] += (da * x[:, t - 1:t]).sum(axis=0)
            g["Wh"] += da.T @ hs[t - 1]
            g["b"] += da.sum(axis=0)
            dh = da @ p["Wh"]
        if clip:                                        # clip by global norm, as in the block above
            total = np.sqrt(sum((v ** 2).sum() for v in g.values()))
            f = min(1.0, clip / (total + 1e-6))
            g = {n: v * f for n, v in g.items()}
        for n in p:
            m[n] = b1 * m[n] + (1 - b1) * g[n]
            s[n] = b2 * s[n] + (1 - b2) * g[n] ** 2
            p[n] = p[n] - lr * (m[n] / (1 - b1 ** step)) / (np.sqrt(s[n] / (1 - b2 ** step)) + eps)
    # honest score: fresh sequences, not the last training batch
    x = rng.standard_normal((2000, T))
    h = np.zeros((2000, hidden))
    for t in range(T):
        h = np.tanh(x[:, t:t + 1] * p["Wx"] + h @ p["Wh"].T + p["b"])
    return ((h @ p["w"] + p["c"] - x[:, 0]) ** 2).mean()

# Predicting 0 always gives MSE ~ 1.0 (the variance of a standard normal): the "I learned nothing" baseline.
for T in (5, 20, 60):
    scores = [train(T, seed=s) for s in range(5)]
    print(f"T={T:3d}  test MSE over 5 seeds: min {min(scores):.3f}  median {sorted(scores)[2]:.3f}  max {max(scores):.3f}")

# Clipping on the hardest case: does it rescue the long dependency?
scores = [train(60, clip=1.0, seed=s) for s in range(5)]
print(f"T= 60  with clip=1.0:  min {min(scores):.3f}  median {sorted(scores)[2]:.3f}  max {max(scores):.3f}")

# Output:
# T=  5  test MSE over 5 seeds: min 0.002  median 0.003  max 0.005
# T= 20  test MSE over 5 seeds: min 0.026  median 0.280  max 0.973
# T= 60  test MSE over 5 seeds: min 0.072  median 0.230  max 1.014
# T= 60  with clip=1.0:  min 0.186  median 0.805  max 1.027

At distance 5 the network learns almost perfectly (median error 0.003). At distance 20 the result depends heavily on the seed: from 0.026 to 0.973, that is, from "it learned" to "it learned nothing". At distance 60 with 400 steps the median error is 0.23, with one seed reaching 0.072 and one staying at 1.014. So the distant dependency is not impossible to learn, but it is much harder and less reliable: the signal reaching the first step is weak and the outcome depends on the initialisation. The last line is the lesson on clipping: with threshold 1.0 the median error at T=60 rises to 0.805, worse than without. It did not cure the vanishing. It did worse, because with a low threshold it triggers often and slows every update, as if the learning rate were smaller. There are 5 seeds, 400 steps and a single task: they indicate an order of magnitude, not a law.

The other remedies, and their limits

• Truncated BPTT. You unroll only for k steps (for example 35) and cut the gradient before that; the state is carried forward but without gradient (h.detach() in PyTorch). It costs less memory, but by construction the network cannot learn dependencies longer than k steps.
• Orthogonal or identity initialisation of W_hh (Le, Jaitly and Hinton, 2015, for the identity with ReLU): you start with a product of Jacobians that does not decay. It helps at the start, it guarantees nothing during training.
• Gated cells (LSTM, GRU). They change the structure so that the gradient has an almost linear path through time. It is the structural remedy and the subject of the next articles.

The three remedies act at different levels: the first limits how far the gradient can travel, the second improves the starting point, the third changes the architecture. Only the last solves the problem at its root, and in fact it is the one used in practice. Clipping, unlike these three, does not touch the vanishing: it limits the length of the gradient that travels, not how far it gets.

How to choose the clipping threshold

The practical procedure is: train once without clipping for a few hundred steps, record the gradient norm before cutting at every step, look at its distribution and choose τ a little above the 95th-99th percentile. A hypothetical example, with numbers invented for teaching purposes: if over 10,000 steps the norm lies in 95% of cases between 0.2 and 1.5, with very rare spikes at 80, 300 and 2000, a sensible threshold is τ = 2 or τ = 5. With τ = 5 a spike of 300 is reduced by a factor of 60 and the normal steps are not touched; with τ = 0.1 all of them would be reduced, and you would be training with a smaller learning rate in effect, which is exactly what the experiment with clip=1.0 showed on a small scale.

The traps you really fall into

The most frequent is a loss that becomes nan after many quiet epochs: it is almost always a gradient spike, and recording the pre-clipping norm at every step reveals it (a jump from 1 to 10⁴ is the signature). The second is a threshold that is too low, as above. The third is thinking that "the loss goes down" means "the network has long memory": it can reduce the error by exploiting only the last elements, and you need a test that really requires the distant dependency, like the first-element one. The last concerns truncated BPTT: if you reuse the state between batches without detaching it from the graph (detach()), PyTorch keeps the whole history in memory and after a few batches you run out of RAM or get an error on the second backward().

Checked with four calculations

A linear RNN with maximum eigenvalue 0.95 gives after 200 steps a factor 0.95^200 ≈ 3.5·10⁻⁵: the gradient has vanished. Clipping does not solve the vanishing because it acts only if the norm exceeds the threshold, and a gradient close to zero already has a norm below the threshold. Truncated BPTT limits how far the gradient travels backwards, clipping limits how large the travelling gradient is: different remedies for different problems. In the BPTT code += is used for dWhh because the same weight appears at every step and the total gradient is the sum of the contributions of all the steps.

What you can do now

You can derive BPTT, write it in NumPy and test it with finite differences. You can explain why the gradient of an RNN is a product of Jacobians diag(1 − h²) · W_hh, measure the decay or growth of its norm and tell vanishing, which clipping does not cure, from the explosion, which clipping contains. The next step is to apply it all to real data: RNN on a time series; then the structure that changes the problem at its root, the LSTM cell. For the prerequisites, backpropagation step by step and gradient descent.

Other lessons in the series that link to this one: LSTM Put to the Test: The Adding Problem, Stacked Layers and Padding, Attention, Derived in Full: Query, Key, Value and the Transformer Block Written From Scratch.