Backpropagation through convolutional and recurrent networks

Backpropagation is one idea, the chain rule applied in reverse over a computation graph, but it looks different depending on how a layer reuses its parameters. A convolutional layer applies the same small kernel at every spatial position. A recurrent layer applies the same weight matrix at every time step. In both cases a single parameter influences the loss through many paths, and its gradient is the sum over all of them. Most of what is special about CNN and RNN backpropagation follows from that one observation.

This guide derives the gradients for convolution, pooling, and a vanilla recurrent network by hand, checks them with small NumPy implementations and finite differences, and then uses the recurrent derivation to explain why gradients vanish or explode over long sequences and how LSTMs change that picture. It assumes you are comfortable with backpropagation in a multilayer perceptron.

A short recap: deltas flow backward

The backpropagation guide covers the fully connected case in detail. The one idea we need is the delta recursion. For a dense layer $z=Wa+b$ with upstream gradient $\delta=\partial L/\partial z$, the local rules are:

The pattern is general. Every layer receives the gradient of the loss with respect to its output, multiplies by its local Jacobian to produce the gradient with respect to its input, and accumulates gradients for its own parameters. When a variable feeds several downstream computations, the multivariable chain rule adds their contributions:

Convolutions and recurrences are nothing more than graphs in which certain variables, the kernel weights or the recurrent matrix, have many children. Keep that sum in mind; every formula below is an instance of it.

Convolution, cross-correlation, and output sizes

In signal processing, the discrete convolution of an image $X$ with a kernel $K$ flips the kernel before sliding it: $(X*K)_{i,j}=\sum_{u,v}X_{i-u,\,j-v}K_{u,v}$. Deep learning libraries almost always compute cross-correlation instead, with no flip:

Because the kernel is learned, the distinction does not change what a network can represent: a network trained with one convention would simply learn the flipped kernel under the other. It does matter when deriving gradients, because the backward pass for the input turns out to be a true convolution, which is a cross-correlation with a kernel rotated by 180 degrees. Throughout this article "conv" means the cross-correlation that frameworks implement.

For an input of height $H$, kernel height $k$, zero padding $p$ on each side, and stride $s$, the output height is:

The same formula applies to width. A $32\times32$ input with a $5\times5$ kernel, no padding, and stride one gives $28\times28$. With a $3\times3$ kernel, padding one, and stride two it gives $\lfloor 31/2\rfloor+1=16$. When $H+2p-k$ is not divisible by $s$, the floor silently discards the last rows of the input; those pixels receive zero gradient. Dumoulin and Visin's guide to convolution arithmetic works through all of these cases, including transposed convolution, with diagrams.

The multi-channel forward pass

A real convolutional layer takes $C$ input channels and produces $F$ output channels, or feature maps. Its kernel tensor has shape $(F,C,k_H,k_W)$ and its bias has one entry per output channel. With stride one and an already padded input, the forward pass for one example is:

With stride $s$ the input index becomes $X_{c,\,si+u,\,sj+v}$. Each output channel sums over all input channels, so a $3\times3$ kernel is really a $C\times3\times3$ block of weights. A typical network stacks these layers with nonlinearities and pooling, then hands a flattened representation to dense layers, as in LeCun et al. (1998).

flowchart LR
    X["Input image<br/>3 x 32 x 32"] --> C1["Conv 3x3<br/>16 filters"]
    C1 --> R1["ReLU"]
    R1 --> P1["Max pool 2x2"]
    P1 --> C2["Conv 3x3<br/>32 filters"]
    C2 --> R2["ReLU"]
    R2 --> P2["Max pool 2x2"]
    P2 --> FL["Flatten"]
    FL --> D["Dense layer"]
    D --> L["Softmax and<br/>cross-entropy loss"]

Parameter sharing

The same kernel weights are reused at every spatial position. A $3\times3$ convolution from 3 to 64 channels has $64\cdot3\cdot9+64=1792$ parameters, regardless of image size. A dense layer connecting a $3\times32\times32$ input to a $64\times32\times32$ output would need $3072\cdot65536=201{,}326{,}592$ weights. Sharing encodes an assumption, that a useful local pattern is useful wherever it appears, and it makes the layer translation equivariant: shifting the input shifts the feature map. That assumption is also why the kernel gradient will be a sum over positions.

Receptive field growth

Each output of a $k\times k$ layer sees a $k\times k$ window of its input. Stacking layers widens the window. For layer $\ell$ with kernel size $k_\ell$ and stride $s_\ell$, the receptive field in input pixels grows as:

Three stride-one $3\times3$ layers therefore see $7\times7$ pixels, the same as one $7\times7$ layer, but with $27C^2$ instead of $49C^2$ weights for $C$ channels in and out, and two extra nonlinearities. Strides and pooling multiply the growth, which is how deep networks reach whole-image context.

Deriving the convolution gradients

Let $G_{f,i,j}=\partial L/\partial Y_{f,i,j}$ be the upstream gradient. We need three things: the gradient for the kernel, the bias, and the input. Work with stride one first.

Kernel. The weight $K_{f,c,u,v}$ appears in every output $Y_{f,i,j}$ of channel $f$, multiplied by $X_{c,i+u,j+v}$. Summing over its children:

Compare this with the forward formula. It is a cross-correlation of the input channel $X_c$ with the upstream gradient map $G_f$, which plays the role of the kernel. For a batch, add the same quantity over examples.

Bias. $b_f$ enters every output of channel $f$ with coefficient one, so $\partial L/\partial b_f=\sum_{i,j}G_{f,i,j}$, summed over the batch as well.

Input. An input pixel $X_{c,m,n}$ contributes to output $Y_{f,i,j}$ whenever $i+u=m$ and $j+v=n$ for some valid kernel offset. Writing $i=m-u$ and $j=n-v$:

where $G$ is treated as zero outside its valid range. The minus signs make this a true convolution. To see it as a cross-correlation, define the rotated kernel $\tilde K_{f,c,u,v}=K_{f,c,\,k_H-1-u,\,k_W-1-v}$ and pad $G$ with $k-1$ zeros on every side to obtain $\tilde G$, where $\tilde G_{f,a,b}=G_{f,\,a-(k_H-1),\,b-(k_W-1)}$. Substituting $u'=k_H-1-u$ and $v'=k_W-1-v$:

So the input gradient is a "full" cross-correlation of the padded upstream gradient with the 180-degree rotated kernel, summed over output channels. The full padding makes the result the same size as the input: an $H-k+1$ map padded by $k-1$ on each side has size $H+k-1$, and correlating with a size-$k$ kernel returns $H$.

flowchart TD
    X["Input X"] --> FWD["Forward: correlate X with K, add b"]
    K["Kernel K"] --> FWD
    FWD --> Y["Output Y"]
    Y --> G["Upstream gradient G = dL/dY"]
    G -.-> DK["dL/dK: correlate X with G"]
    X -.-> DK
    G -.-> DB["dL/db: sum G over positions"]
    G -.-> DX["dL/dX: full correlation of padded G<br/>with 180 degree rotated K"]
    K -.-> DX

Stride and padding. With stride $s$, the kernel gradient uses $X_{c,\,si+u,\,sj+v}$ in place of $X_{c,i+u,j+v}$. For the input gradient, the cleanest view is to dilate $G$ by inserting $s-1$ zeros between its entries, then apply the stride-one formula; this is exactly a transposed convolution. Padding is easiest to handle by computing the gradient with respect to the padded input and cropping the border away, since the padded zeros are constants. The implementation below does exactly that.

A worked numeric example

Take a single channel, a $3\times3$ input, a $2\times2$ kernel, and bias $0.5$:

The output is $2\times2$. For the top-left entry, $Y_{0,0}=1\cdot1+2\cdot(-1)+0\cdot2+1\cdot0+0.5=-0.5$. The full output is $Y=\begin{bmatrix}-0.5&4.5\\3.5&0.5\end{bmatrix}$. Suppose the upstream gradient is:

The kernel gradient correlates $X$ with $G$. Its top-left entry overlays $G$ on the top-left $2\times2$ block of $X$: $1\cdot1+2\cdot0+0\cdot(-1)+1\cdot2=3$. Repeating at the other three offsets gives:

For the input gradient, the center pixel $X_{1,1}$ is touched by all four outputs, each through a different kernel entry: $G_{0,0}K_{1,1}+G_{0,1}K_{1,0}+G_{1,0}K_{0,1}+G_{1,1}K_{0,0}=0+0+1+2=3$. The corner pixel $X_{0,0}$ is touched only by $Y_{0,0}$ through $K_{0,0}$, so its gradient is $1$. The complete result is:

The code below reproduces each number and checks the input gradient against central finite differences of $L=\sum_{i,j}G_{i,j}Y_{i,j}$, a loss chosen so that $\partial L/\partial Y=G$ exactly.

import numpy as np

X = np.array([[1., 2., 0.],
              [0., 1., 3.],
              [2., 1., 1.]])
K = np.array([[1., -1.],
              [2.,  0.]])
b = 0.5

def corr2d(A, B):
    """Valid cross-correlation of A with B (no kernel flip)."""
    H = A.shape[0] - B.shape[0] + 1
    W = A.shape[1] - B.shape[1] + 1
    return np.array([[np.sum(A[i:i+B.shape[0], j:j+B.shape[1]] * B)
                      for j in range(W)] for i in range(H)])

Y = corr2d(X, K) + b
G = np.array([[1., 0.],          # upstream gradient dL/dY
              [-1., 2.]])

dK = corr2d(X, G)                          # correlate input with upstream grad
db = G.sum()
dX = corr2d(np.pad(G, 1), np.rot90(K, 2))  # full conv with rotated kernel

print("Y =\n", Y)
print("dK =\n", dK)
print("db =", db)
print("dX =\n", dX)

# Brute-force check: L = sum(G * Y) has dL/dY = G exactly.
eps = 1e-6
num_dX = np.zeros_like(X)
for i in range(3):
    for j in range(3):
        Xp = X.copy(); Xp[i, j] += eps
        Xm = X.copy(); Xm[i, j] -= eps
        num_dX[i, j] = (np.sum(G * corr2d(Xp, K)) - np.sum(G * corr2d(Xm, K))) / (2 * eps)
print("max |dX - numeric| =", np.abs(dX - num_dX).max())

On my run the finite-difference discrepancy was about $3\times10^{-10}$, which is round-off, not a derivation error.

A multi-channel implementation and the im2col view

The naive implementation below follows the forward formula with explicit loops. Its backward pass does not use the rotated-kernel formula directly. Instead, it walks over outputs and asks which patch and which weights produced each one, then scatters the gradient back. This "scatter" form and the "full correlation" form compute the same sums in a different order; the scatter version is easier to get right and handles padding by cropping.

import numpy as np

def conv2d_forward(X, W, b, pad=0):
    """X: (N, C, H, W_in), W: (F, C, kH, kW), b: (F,). Stride 1."""
    N, C, H, Wd = X.shape
    F, _, kH, kW = W.shape
    Xp = np.pad(X, ((0, 0), (0, 0), (pad, pad), (pad, pad)))
    Ho, Wo = H + 2 * pad - kH + 1, Wd + 2 * pad - kW + 1
    Y = np.zeros((N, F, Ho, Wo))
    for n in range(N):
        for f in range(F):
            for i in range(Ho):
                for j in range(Wo):
                    Y[n, f, i, j] = np.sum(Xp[n, :, i:i+kH, j:j+kW] * W[f]) + b[f]
    return Y, (X, W, pad)

def conv2d_backward(dY, cache):
    X, W, pad = cache
    N, C, H, Wd = X.shape
    F, _, kH, kW = W.shape
    _, _, Ho, Wo = dY.shape
    Xp = np.pad(X, ((0, 0), (0, 0), (pad, pad), (pad, pad)))
    dXp = np.zeros_like(Xp)
    dW = np.zeros_like(W)
    db = dY.sum(axis=(0, 2, 3))
    for n in range(N):
        for f in range(F):
            for i in range(Ho):
                for j in range(Wo):
                    g = dY[n, f, i, j]
                    # Y[n,f,i,j] touched this patch with these weights:
                    dW[f] += g * Xp[n, :, i:i+kH, j:j+kW]
                    dXp[n, :, i:i+kH, j:j+kW] += g * W[f]
    dX = dXp[:, :, pad:pad+H, pad:pad+Wd] if pad else dXp
    return dX, dW, db

def num_grad(f, A, eps=1e-6):
    g = np.zeros_like(A)
    it = np.nditer(A, flags=["multi_index"])
    for _ in it:
        idx = it.multi_index
        old = A[idx]
        A[idx] = old + eps; fp = f()
        A[idx] = old - eps; fm = f()
        A[idx] = old
        g[idx] = (fp - fm) / (2 * eps)
    return g

def rel_err(a, b):
    return np.abs(a - b).max() / max(1e-12, np.abs(a).max() + np.abs(b).max())

rng = np.random.default_rng(0)
X = rng.standard_normal((2, 3, 5, 5))
W = rng.standard_normal((4, 3, 3, 3))
b = rng.standard_normal(4)
R = rng.standard_normal((2, 4, 5, 5))   # random projection makes L a scalar

def loss():
    Y, _ = conv2d_forward(X, W, b, pad=1)
    return np.sum(Y * R)                 # so dL/dY = R

Y, cache = conv2d_forward(X, W, b, pad=1)
dX, dW, db = conv2d_backward(R, cache)
print("dX rel err:", rel_err(dX, num_grad(loss, X)))
print("dW rel err:", rel_err(dW, num_grad(loss, W)))
print("db rel err:", rel_err(db, num_grad(loss, b)))

All three relative errors were below $10^{-9}$ in my run. The random projection $R$ is a useful trick: a loss of the form $\sum Y\odot R$ exercises every output with a different weight, so an index error cannot hide behind symmetric upstream gradients.

Convolution as a matrix multiply

The im2col transformation copies every receptive-field patch into a column of a matrix $A$ with $Ck_Hk_W$ rows and one column per output position. Reshape the kernel into a matrix $W_m$ of shape $(F,\,Ck_Hk_W)$. The entire convolution becomes $Y_m=W_mA$, a dense layer applied to patches. The dense-layer rules from the recap carry over immediately:

The only new step is col2im, which adds each column of $\partial L/\partial A$ back into its patch. Because neighboring patches overlap, a pixel's gradient is the sum of contributions from every patch that contained it, which is the same sum as before.

import numpy as np

def im2col(X, kH, kW):
    """X: (C, H, W) -> columns of shape (C*kH*kW, Ho*Wo), stride 1, no padding."""
    C, H, W = X.shape
    Ho, Wo = H - kH + 1, W - kW + 1
    cols = np.empty((C * kH * kW, Ho * Wo))
    for i in range(Ho):
        for j in range(Wo):
            cols[:, i * Wo + j] = X[:, i:i+kH, j:j+kW].ravel()
    return cols, Ho, Wo

def col2im(dcols, shape, kH, kW):
    """Adjoint of im2col: scatter-add every column back into its patch."""
    C, H, W = shape
    Ho, Wo = H - kH + 1, W - kW + 1
    dX = np.zeros(shape)
    for i in range(Ho):
        for j in range(Wo):
            dX[:, i:i+kH, j:j+kW] += dcols[:, i * Wo + j].reshape(C, kH, kW)
    return dX

rng = np.random.default_rng(1)
C, H, Wd, F, k = 2, 4, 4, 3, 3
X = rng.standard_normal((C, H, Wd))
Wt = rng.standard_normal((F, C, k, k))

cols, Ho, Wo = im2col(X, k, k)
Wm = Wt.reshape(F, -1)                  # (F, C*k*k)
Y = (Wm @ cols).reshape(F, Ho, Wo)      # convolution as one matmul

# Reference: direct loops.
Yref = np.array([[[np.sum(X[:, i:i+k, j:j+k] * Wt[f]) for j in range(Wo)]
                  for i in range(Ho)] for f in range(F)])
print("forward match:", np.allclose(Y, Yref))

# Backward is the dense-layer rule Y = W A  =>  dW = dY A^T,  dA = W^T dY.
dY = rng.standard_normal((F, Ho, Wo))
dYm = dY.reshape(F, -1)
dW = (dYm @ cols.T).reshape(Wt.shape)
dX = col2im(Wm.T @ dYm, X.shape, k, k)

# Same answers as the loop-based rule.
dW_ref = np.zeros_like(Wt); dX_ref = np.zeros_like(X)
for f in range(F):
    for i in range(Ho):
        for j in range(Wo):
            dW_ref[f] += dY[f, i, j] * X[:, i:i+k, j:j+k]
            dX_ref[:, i:i+k, j:j+k] += dY[f, i, j] * Wt[f]
print("dW match:", np.allclose(dW, dW_ref), " dX match:", np.allclose(dX, dX_ref))

Many production libraries have used im2col-style lowering because large matrix multiplies are extremely well optimized, though modern kernels also use direct, FFT, and Winograd algorithms. The memory cost is visible in the shapes: each input pixel is copied up to $k_Hk_W$ times.

Pooling layers route the gradient

Pooling has no parameters, so its backward pass only has to answer one question: how much did each input change the output? For max pooling over a window, the output equals the largest input. Infinitesimally perturbing a non-maximal input does not change the output, so its derivative is zero. Perturbing the maximal input changes the output one for one:

The upstream gradient is routed to the winner and the other inputs receive nothing. The forward pass must therefore remember where each maximum was. Average pooling over an $s\times s$ window has $Y=\frac{1}{s^2}\sum X$, so each input receives an equal share $G_{i,j}/s^2$.

import numpy as np

def maxpool_forward(X, s=2):
    """X: (H, W) with H, W divisible by s. Non-overlapping s x s windows."""
    H, W = X.shape
    blocks = X.reshape(H // s, s, W // s, s)
    Y = blocks.max(axis=(1, 3))
    mask = blocks == Y[:, None, :, None]   # where each max came from
    return Y, mask

def maxpool_backward(dY, mask):
    # Route each upstream gradient to its argmax (ties would all receive it;
    # frameworks usually pick one).
    s = mask.shape[1]
    d = mask * dY[:, None, :, None]
    return d.reshape(mask.shape[0] * s, mask.shape[2] * s)

def avgpool_backward(dY, s=2):
    # Every input in the window gets an equal 1/s^2 share.
    return np.kron(dY, np.ones((s, s))) / (s * s)

X = np.array([[1., 3., 2., 0.],
              [4., 2., 1., 5.],
              [0., 1., 2., 2.],
              [3., 0., 1., 7.]])
Y, mask = maxpool_forward(X)
dY = np.array([[10., 20.],
               [30., 40.]])
print("Y =\n", Y)
print("max-pool dX =\n", maxpool_backward(dY, mask))
print("avg-pool dX =\n", avgpool_backward(dY))

The max-pooled output is $\begin{bmatrix}4&5\\3&7\end{bmatrix}$, and the gradients $10,20,30,40$ land exactly at the positions of $4,5,3,7$. Max is not differentiable at ties, where any convex combination of the tied subgradients is valid; with real-valued activations exact ties are uncommon, but they do occur after ReLU produces many zeros. The comment in the code marks that edge case.

Recurrent networks and backpropagation through time

A vanilla recurrent network reads a sequence $x_1,\ldots,x_T$ and keeps a hidden state:

The loss is a sum of per-step terms, $L=\sum_{t=1}^{T}L_t(y_t)$, with $h_0$ given. Unrolling the loop over time produces an ordinary deep feedforward graph with one layer per time step, in which every layer shares the same weights. Backpropagation on that graph is called backpropagation through time, a name popularized by Werbos (1990).

flowchart LR
    H0["h0"] --> H1["h1"]
    H1 --> H2["h2"]
    H2 --> H3["h3"]
    X1["x1"] --> H1
    X2["x2"] --> H2
    X3["x3"] --> H3
    H1 --> L1["Loss L1"]
    H2 --> L2["Loss L2"]
    H3 --> L3["Loss L3"]
    L3 -.-> H3
    L2 -.-> H2
    L1 -.-> H1
    H3 -.-> H2
    H2 -.-> H1

The recursion for the hidden-state gradient

Define $\delta_t=\partial L/\partial h_t$, the gradient of the whole loss with respect to the hidden state. The state $h_t$ has two children: the output $y_t$, and the next pre-activation $a_{t+1}$. The chain rule adds both paths:

The Jacobian of one step is $\partial h_{t+1}/\partial h_t=\operatorname{diag}(1-h_{t+1}\odot h_{t+1})\,W_{hh}$, because $\tanh'(a)=1-\tanh^2(a)$. It is convenient to name the gradient with respect to the pre-activation:

So the backward pass runs from $t=T$ down to $1$, carrying $W_{hh}^\top g_{t+1}$ from the future and adding the local output gradient.

Summing over time

Each weight matrix is used at every step, so its gradient is the sum of the dense-layer rule applied at each step:

and similarly $\partial L/\partial W_{hy}=\sum_t(\partial L_t/\partial y_t)h_t^\top$. This is the recurrent analogue of the kernel gradient summing over spatial positions.

Implementation and gradient check

The code uses a squared-error loss at every step to keep the output layer simple. With softmax and cross-entropy, only the line computing dy would change.

import numpy as np

def rnn_forward(xs, h0, P):
    """xs: (T, D). Returns per-step hidden states, outputs, and loss."""
    Wxh, Whh, bh, Why, by = P["Wxh"], P["Whh"], P["bh"], P["Why"], P["by"]
    hs, ys = {-1: h0}, {}
    for t in range(len(xs)):
        hs[t] = np.tanh(Wxh @ xs[t] + Whh @ hs[t - 1] + bh)
        ys[t] = Why @ hs[t] + by
    return hs, ys

def rnn_loss(xs, targets, h0, P):
    _, ys = rnn_forward(xs, h0, P)
    return sum(0.5 * np.sum((ys[t] - targets[t]) ** 2) for t in range(len(xs)))

def rnn_bptt(xs, targets, h0, P):
    hs, ys = rnn_forward(xs, h0, P)
    grads = {k: np.zeros_like(v) for k, v in P.items()}
    dh_next = np.zeros_like(h0)              # dL/dh_t arriving from step t+1
    for t in reversed(range(len(xs))):
        dy = ys[t] - targets[t]              # dL_t/dy_t for squared error
        grads["Why"] += np.outer(dy, hs[t])
        grads["by"] += dy
        dh = P["Why"].T @ dy + dh_next       # two paths: output at t, and h_{t+1}
        da = dh * (1.0 - hs[t] ** 2)         # through tanh
        grads["Wxh"] += np.outer(da, xs[t])  # shared weights: sum over time
        grads["Whh"] += np.outer(da, hs[t - 1])
        grads["bh"] += da
        dh_next = P["Whh"].T @ da            # send gradient back to h_{t-1}
    return grads

rng = np.random.default_rng(0)
D, Hd, O, T = 3, 5, 2, 6
P = {"Wxh": 0.5 * rng.standard_normal((Hd, D)),
     "Whh": 0.5 * rng.standard_normal((Hd, Hd)),
     "bh": 0.1 * rng.standard_normal(Hd),
     "Why": 0.5 * rng.standard_normal((O, Hd)),
     "by": 0.1 * rng.standard_normal(O)}
xs = rng.standard_normal((T, D))
targets = rng.standard_normal((T, O))
h0 = np.zeros(Hd)

grads = rnn_bptt(xs, targets, h0, P)
eps = 1e-6
for name, A in P.items():
    num = np.zeros_like(A)
    it = np.nditer(A, flags=["multi_index"])
    for _ in it:
        idx = it.multi_index
        old = A[idx]
        A[idx] = old + eps; lp = rnn_loss(xs, targets, h0, P)
        A[idx] = old - eps; lm = rnn_loss(xs, targets, h0, P)
        A[idx] = old
        num[idx] = (lp - lm) / (2 * eps)
    err = np.abs(num - grads[name]).max() / (np.abs(num).max() + np.abs(grads[name]).max())
    print(f"{name:>3} rel err: {err:.2e}")

Every parameter's relative error was around $2\times10^{-10}$ to $4\times10^{-10}$ in my run. The dh_next variable is the whole of BPTT: delete it and the code still runs, still produces plausible gradients, and silently trains a model that ignores the past. A gradient check catches that bug immediately.

Why gradients vanish or explode

Suppose the loss depends only on the final state. Unrolling the recursion, the gradient that reaches an early state is a product of step Jacobians:

Take spectral norms. Since $0\lt\tanh'(a)\le1$, let $\gamma\le1$ bound the diagonal factors and $\sigma_{\max}$ be the largest singular value of $W_{hh}$:

If $\gamma\sigma_{\max}\lt1$, the contribution of step $k$ shrinks at least exponentially with the distance $T-k$. That is a sufficient condition for vanishing. For the linear case, decompose $W_{hh}$ in its eigenbasis: the component of the gradient along an eigenvector with eigenvalue $\lambda$ is scaled by $\lambda^{T-k}$, so any $|\lambda|\gt1$ produces exponential growth. With the nonlinearity, $\sigma_{\max}\gt1/\gamma$ is necessary for explosion but not sufficient, because saturated units shrink the diagonal factors. Bengio, Simard and Frasconi (1994) analyzed why storing information robustly in such systems tends to make long-range gradients vanish, and Pascanu, Mikolov and Bengio (2013) gave the singular-value argument above along with a geometric picture of exploding gradients.

The consequence for learning is easy to state. The total gradient for $W_{hh}$ is a sum over time of terms, and the terms linking distant steps are exponentially smaller than the terms linking nearby steps. Short-range dependencies dominate the update, so the model has little signal for learning long-range ones.

An experiment

The script below builds $W_{hh}$ as a scaled orthogonal matrix, so every singular value equals the scale. It runs 50 steps, places a unit-norm gradient on $h_T$, and prints the norm of $\partial L/\partial h_t$ as it travels back.

import numpy as np

def grad_norms(scale, T=50, H=32, seed=0, activation="tanh"):
    """Loss depends only on h_T. Return ||dL/dh_t|| for t = T, T-1, ..., 0."""
    rng = np.random.default_rng(seed)
    Q, _ = np.linalg.qr(rng.standard_normal((H, H)))
    Whh = scale * Q                      # every singular value equals `scale`
    Wxh = 0.25 * rng.standard_normal((H, 4))
    xs = rng.standard_normal((T, 4))
    hs = [np.zeros(H)]
    for t in range(T):
        a = Whh @ hs[-1] + Wxh @ xs[t]
        hs.append(np.tanh(a) if activation == "tanh" else a)
    g = np.ones(H) / np.sqrt(H)          # dL/dh_T, unit norm
    norms = [np.linalg.norm(g)]
    for t in range(T, 0, -1):
        d = g * (1 - hs[t] ** 2) if activation == "tanh" else g
        g = Whh.T @ d                    # dL/dh_{t-1}
        norms.append(np.linalg.norm(g))
    return norms                         # norms[k] = ||dL/dh_{T-k}||

for act, scales in [("tanh", [0.5, 1.0, 3.0]), ("linear", [0.9, 1.1])]:
    for s in scales:
        n = grad_norms(s, activation=act)
        print(f"{act:>6} scale={s:<4} steps back: 0 -> {n[0]:.2e}  "
              f"10 -> {n[10]:.2e}  25 -> {n[25]:.2e}  50 -> {n[50]:.2e}")

With this seed, the output was:

The linear rows match the theory exactly: an orthogonal matrix preserves norms, so the gradient scales as $0.9^{50}\approx5.15\times10^{-3}$ and $1.1^{50}\approx117$. The tanh row at scale one still vanishes, because every non-zero activation makes $\tanh'\lt1$, and those factors compound. At scale three the gradient grows despite partial saturation. The exact values depend on the seed and inputs; the exponential trend is the robust part.

Truncated BPTT and gradient clipping

Full BPTT over a long stream needs memory for every hidden state. Truncated BPTT splits the sequence into windows of length $k$: it carries the hidden state forward across windows but stops the backward pass at the window boundary. Gradients from dependencies longer than $k$ steps are simply dropped, which biases the estimate. That is often acceptable because such terms were already tiny, but it also means the model is never directly trained to use information older than the window.

Exploding gradients have a blunt but effective fix. Clip by global norm: if the norm of the concatenated gradient exceeds a threshold $\tau$, rescale it.

Rescaling keeps the direction and limits the step size, which helps when the loss surface has the steep "walls" Pascanu et al. describe.

import numpy as np

def clip_by_global_norm(grads, max_norm):
    total = np.sqrt(sum(np.sum(g * g) for g in grads.values()))
    scale = min(1.0, max_norm / (total + 1e-12))
    return {k: g * scale for k, g in grads.items()}, total

grads = {"Whh": np.full((2, 2), 30.0), "Wxh": np.full((2, 1), 40.0)}
clipped, before = clip_by_global_norm(grads, max_norm=5.0)
after = np.sqrt(sum(np.sum(g * g) for g in clipped.values()))
print(f"norm before: {before:.2f}, after: {after:.2f}")

This prints a norm of $82.46$ before clipping and $5.00$ after. Clipping does nothing for vanishing gradients. For that, the architecture has to change.

LSTM and GRU: an additive path for the gradient

The Long Short-Term Memory network of Hochreiter and Schmidhuber adds a cell state $c_t$ whose update is mostly additive. In the widely used form with a forget gate, introduced by Gers, Schmidhuber and Cummins in 2000, write $z_t=[h_{t-1};x_t]$ and $\sigma$ for the logistic sigmoid:

flowchart LR
    CP["Previous cell c t-1"] --> MULF["Multiply by forget gate f"]
    MULF --> ADD["Add"]
    HX["Previous h and input x"] --> GF["Forget gate f sigmoid"]
    HX --> GI["Input gate i sigmoid"]
    HX --> GC["Candidate tanh"]
    HX --> GO["Output gate o sigmoid"]
    GF --> MULF
    GI --> MULI["Multiply i and candidate"]
    GC --> MULI
    MULI --> ADD
    ADD --> CN["New cell c t"]
    CN --> TH["tanh"]
    TH --> MULO["Multiply by output gate o"]
    GO --> MULO
    MULO --> HN["New hidden h t"]

The key term is the direct derivative of one cell state with respect to the previous one, ignoring the indirect paths through the gates:

Along the cell path, the gradient from step $T$ back to step $k$ is multiplied by $\prod_{t=k+1}^{T}\operatorname{diag}(f_t)$. There is no repeated multiplication by a weight matrix and no repeated $\tanh'$ factor. If the network learns forget gates near one for a particular unit, that unit's gradient can travel many steps with little decay. The gates are themselves learned, so the network decides when to keep and when to erase. Initializing the forget-gate bias to a positive value such as one is a common heuristic for starting with long memory. The full gradient also includes paths through $h_{t-1}$ into every gate, so LSTMs reduce the vanishing problem rather than eliminate it.

The gated recurrent unit of Cho et al. (2014) merges the cell and hidden state with an update gate $z_t$ and a reset gate $r_t$:

The direct path contributes $\operatorname{diag}(1-z_t)$ to the Jacobian, which plays the same role as the forget gate. GRUs have fewer parameters per unit; neither architecture dominates the other across all tasks.

Common bugs and how to catch them

Flipping the kernel twice, or not at all. If your forward pass is cross-correlation, the input gradient needs the rotated kernel and the kernel gradient does not. Mixing conventions gives gradients that are correct for symmetric kernels and wrong otherwise. Test with an asymmetric kernel like the one in the worked example.

Overwriting instead of accumulating. Shared parameters need += across positions, examples, and time steps. Assigning with = keeps only the last contribution. The same applies to input gradients where patches overlap.

Forgetting the padding crop or stride dilation. Gradients for padded zeros must be discarded. With stride greater than one, scattering into the wrong positions often still produces the right shape.

Dropping the recurrent term. Omitting dh_next, or adding it after the tanh derivative instead of before, yields gradients that look reasonable. Only a numerical check reveals the error.

Reusing a stale hidden state across sequences. Reset or deliberately carry $h_0$ between batches. When carrying it in truncated BPTT, stop gradients at the boundary, or memory grows with every batch in autodiff frameworks.

Testing with a symmetric loss. A loss like $\sum Y$ gives the same upstream gradient everywhere, which can hide index transpositions. Use a random projection, as the checks above do, in double precision, with central differences and a relative error.

Practice exercises

1. Extend conv2d_forward and conv2d_backward to stride two. Derive the kernel gradient with the index $si+u$, implement it, and confirm it with the same finite-difference check. Then verify that the input gradient equals the stride-one formula applied to a zero-dilated $G$.

2. Compute the receptive field of a network with layers: $3\times3$ conv stride one, $2\times2$ max pool stride two, $3\times3$ conv stride one, $2\times2$ max pool stride two, $3\times3$ conv stride one. Check your answer by placing a nonzero gradient at a single output position and counting which input pixels receive a nonzero gradient.

3. Replace the squared-error loss in the RNN with softmax cross-entropy over a small vocabulary. Only the output gradient should change. Rerun the gradient check.

4. Modify the vanishing-gradient experiment to track the norm of each per-step term $g_t h_{t-1}^\top$ in the $W_{hh}$ gradient when the loss is placed on the last step. Plot or print how much of the total gradient comes from the last five steps.

5. Implement one LSTM step forward and backward in NumPy, gradient-check it, then repeat the experiment above with forget-gate biases of $-2$, $0$, and $2$. Describe how the norm of $\partial L/\partial c_t$ changes with distance, and explain why in terms of $\prod_t\operatorname{diag}(f_t)$.

Each exercise ends with a measurement. That habit, deriving a formula and then checking it numerically on a case small enough to reason about, is what makes backpropagation through any new layer routine rather than mysterious.