Vanishing and Exploding Gradients

Vanishing and exploding gradients are two failure modes of training deep networks with backpropagation. The error signal that reaches an early layer is the output error multiplied by one factor per layer in between. When those factors tend to shrink vectors, the gradient decays geometrically with depth, and the first convolutional layers of an image classifier barely change; when they tend to stretch vectors, it grows geometrically, and a single update can throw the weights far away. The problem was analyzed for recurrent networks in the early 1990s and for deep feedforward networks in 2010; saturating activation functions, the scale of the initial weights, and later the absence of normalization and shortcuts were each identified as causes, and each has a remedy.

Definition

A network has a vanishing gradient when the norm of the gradient of the loss with respect to a layer’s activations (or weights) is many orders of magnitude smaller in early layers than in late ones, typically shrinking by a roughly constant factor per layer. It has an exploding gradient when that norm grows by a roughly constant factor per layer, so that early layers, or the whole gradient, become very large. In recurrent networks unrolled over time the same happens across time steps instead of layers. Both are properties of a network at particular weights, not of the architecture alone: the same network can have a healthy gradient at one scale of initialization and a vanishing or exploding one at another.

Intuition

Think of a deep stack of $3 \times 3$ convolutions that classifies images of horizontal and vertical stripes. To learn better edge detectors in the first layer, that layer needs to know how the loss changes when its feature maps change. The answer is assembled backward: the classifier’s error is passed through the last layer, then the one before, and so on. Each layer rescales the signal by roughly the same amount, so twenty layers multiply it by that amount twenty times, the way compound interest works. A factor of $0.5$ per layer leaves about $10^{-6}$ of the signal after twenty layers; a factor of $1.5$ multiplies it by about 3,000. Only factors close to 1 keep the gradient on the same scale at every depth.

Two things set the factor. The weights decide how much a layer stretches or shrinks a vector. The activation’s slope decides how much survives the nonlinearity: a sigmoid unit passes at most a quarter of the signal (its slope never exceeds $1/4$), and a saturated unit, far out on the flat part of the curve, passes almost nothing.

Formal Definition

In the notation of the backpropagation article, layer $\ell$ computes $\mathbf{z}^{(\ell)} = \mathbf{W}^{(\ell)}\mathbf{a}^{(\ell-1)} + \mathbf{b}^{(\ell)}$ and $\mathbf{a}^{(\ell)} = \phi(\mathbf{z}^{(\ell)})$, and the error signal $\boldsymbol{\delta}^{(\ell)} = \partial L / \partial \mathbf{z}^{(\ell)}$ obeys $\boldsymbol{\delta}^{(\ell-1)} = \mathbf{D}^{(\ell-1)} (\mathbf{W}^{(\ell)})^\top \boldsymbol{\delta}^{(\ell)}$ with $\mathbf{D}^{(\ell)} = \operatorname{diag}\big(\phi'(\mathbf{z}^{(\ell)})\big)$. Unrolling from the output layer $L$ down to layer $k$,

$$\boldsymbol{\delta}^{(k)} = \mathbf{D}^{(k)} (\mathbf{W}^{(k+1)})^\top \mathbf{D}^{(k+1)} (\mathbf{W}^{(k+2)})^\top \cdots \mathbf{D}^{(L-1)} (\mathbf{W}^{(L)})^\top \boldsymbol{\delta}^{(L)} = \left(\frac{\partial \mathbf{z}^{(L)}}{\partial \mathbf{z}^{(k)}}\right)^{\!\top} \boldsymbol{\delta}^{(L)},$$

a product of $L - k$ layer Jacobians. Since the norm of a product is at most the product of the norms,

$$\big\lVert \boldsymbol{\delta}^{(k)} \big\rVert \le \prod_{\ell=k+1}^{L} \big\lVert \mathbf{W}^{(\ell)} \big\rVert_2 \, \big\lVert \mathbf{D}^{(\ell-1)} \big\rVert_2 \; \big\lVert \boldsymbol{\delta}^{(L)} \big\rVert \le (s\gamma)^{L-k} \, \big\lVert \boldsymbol{\delta}^{(L)} \big\rVert,$$

where $s$ bounds the largest singular value of every weight matrix and $\gamma$ bounds $|\phi'|$. If $s\gamma < 1$, the error signal vanishes at least geometrically with depth. Pascanu, Mikolov and Bengio (2013) proved this sufficient condition for a recurrent network, where the same matrix appears at every step, and noted that inverting it gives a necessary condition for explosion: the largest singular value must exceed $1/\gamma$, which is 1 for tanh and 4 for the sigmoid. The weight gradients follow from the error signals, $\partial L / \partial \mathbf{W}^{(\ell)} = \boldsymbol{\delta}^{(\ell)} (\mathbf{a}^{(\ell-1)})^\top$.

The per-layer factor on average. The bound is a worst case. Glorot and Bengio (2010) computed what happens on average for random weights in the linear regime ($\phi' \approx 1$): with independent weights of variance $\operatorname{Var}[W^{(\ell)}]$ and layers of width $n_\ell$, the variance of the back-propagated gradient at layer $i$ of a $d$-layer network is that at the top multiplied by $\prod_{\ell=i}^{d} n_{\ell+1} \operatorname{Var}[W^{(\ell)}]$. Each layer contributes a factor $n_\text{out} \operatorname{Var}[w]$, where $n_\text{out}$ is the number of units a unit feeds; for a convolution it is $k^2 C_\text{out}$. Including the activation’s slope under the same independence assumption gives the factor below. He et al. (2015) derived it for ReLU, where it reads $\tfrac{1}{2} \hat{n} \operatorname{Var}[w]$ with $\hat{n} = k^2 C_\text{out}$ for a convolution; the general form is our extension, which the example below checks:

$$\chi = n_\text{out}\, \operatorname{Var}[w]\; \mathbb{E}\big[\phi'(z)^2\big], \qquad \frac{\lVert \boldsymbol{\delta}^{(\ell-1)} \rVert}{\lVert \boldsymbol{\delta}^{(\ell)} \rVert} \approx \sqrt{\chi},$$

and the gradient vanishes when $\chi < 1$ and explodes when $\chi > 1$, by $\chi^{(L-k)/2}$ over $L-k$ layers. For the sigmoid, $\mathbb{E}[\sigma'(z)^2] \le 1/16$. For ReLU, $\phi'(z)^2$ is 1 for active units and 0 otherwise, so $\mathbb{E}[\phi'(z)^2]$ is the fraction of active units, about $1/2$ for zero-mean pre-activations, and $\chi \approx 1$ requires $\operatorname{Var}[w] \approx 2 / n$. The derivations that choose $\operatorname{Var}[w]$ from this condition, Glorot and Bengio’s and He et al.’s, belong to weight initialization.

Properties

Weight scale matters as much as the activation. Glorot and Bengio’s commonly used heuristic of the time, uniform weights in $[-1/\sqrt{n}, 1/\sqrt{n}]$, gives $n \operatorname{Var}[W] = 1/3$ per layer, so even in the linear regime the back-propagated gradient shrinks with depth. ReLU removes the activation’s share of the shrinkage on active paths but not the weights’ share, as the example below shows. Scaling up a sigmoid network’s weights does not rescue it either: larger pre-activations push units into the flat tails, where the slope is far below $1/4$.

Saturation, not just small slopes. Glorot and Bengio (2010) trained fully connected networks with up to five hidden layers of 1,000 units on Shapeset-3×2, synthetic 32×32 images of triangles, parallelograms and ellipses. With sigmoid units the top hidden layer was pushed to its lower saturation value 0 early in training; a four-layer network escaped this regime after about 100 epochs, while the five-layer one never did. They traced this to the sigmoid’s non-zero mean and recommended against sigmoid hidden units when starting from small random weights. Tanh networks with the standard initialization saturated layer by layer, starting with the first. On Shapeset-3×2 their Table 1 reports test errors of 82.61% for sigmoid networks, 27.15% for tanh with the standard initialization and 15.60% for tanh with their normalized initialization; the table’s caption says five hidden layers, although the text gives four as the best depth for the sigmoid on this dataset.

Vanishing back-propagated gradients do not imply vanishing weight gradients. Since $\partial L / \partial \mathbf{W}^{(\ell)} = \boldsymbol{\delta}^{(\ell)} (\mathbf{a}^{(\ell-1)})^\top$, activations that shrink going forward can compensate error signals that shrink going backward. Glorot and Bengio showed that in the linear regime, with equal widths and the same initialization in every layer, the variance of the weight gradient is the same for all layers even while the back-propagated gradient vanishes or explodes, and observed exactly that at initialization: the back-propagated gradients got smaller toward the input, the weight gradients did not. The example below shows the same in a ReLU network: the weight gradients are uniform across layers but uniformly tiny.

The gradient changes during training. Glorot and Bengio found that the downward shrinkage of back-propagated gradients after the standard initialization reversed very quickly once training started, and that weight gradients, initially of similar size across layers, drifted apart as training went on, with larger gradients in the lower layers. Gradient health is a property of the current weights, and they call monitoring activations and gradients across layers and training iterations a powerful investigative tool.

Exploding gradients and steep walls. In a recurrent network the Jacobian factors are all built from the same recurrent matrix, so the product behaves like a matrix power. For the linear case Pascanu et al. (2013) show that a spectral radius below 1 is sufficient for long-term gradient components to vanish and one above 1 necessary for them to explode. They also hypothesize that where the gradient explodes along some direction, the curvature does too, so the loss surface has a steep wall; an ordinary gradient step that hits the wall is thrown far across the valley. In our reading, feedforward image networks show the same arithmetic over depth when $\chi > 1$, with the difference that each layer has its own matrix.

History

Sepp (Josef) Hochreiter’s 1991 diploma thesis at the Technical University of Munich, written in German and supervised by Jürgen Schmidhuber, analyzed the error flow of gradient-based learning in recurrent networks. Its summary of that analysis (Section 2.3) observes that the output’s dependence on a weight contains a product of many factors, each an activation derivative times a weight. Because this product falls or rises exponentially with the weights, a weight enters the gradient either too strongly, causing instability, or too weakly, so that it hardly changes; and initializing weights near zero makes the product very small and learning very slow. Chapter 4 of the thesis builds units with constant error backflow to carry information over long time spans. Bengio, Simard and Frasconi (1994) showed, for recurrent networks, a trade-off between storing information robustly against noise and being trainable by gradient descent: in the robust regime the derivative of the current state with respect to a distant past state decreases exponentially with the time between them. Pascanu et al. (2013) describe the LSTM of Hochreiter and Schmidhuber (1997) as a fix for the vanishing case, built on a linear unit with a self-connection of weight 1 and learned gates, and note that it does not explicitly address exploding gradients.

For feedforward networks, Glorot and Bengio (2010) asked why standard gradient descent from random initialization did so poorly on deep networks, which until then had been trained successfully only with new initialization or training mechanisms such as unsupervised pre-training. Their answers, saturation and per-layer Jacobians far from 1, led to the normalized initialization.

Degradation is a different problem. By 2015, He et al. (2016) considered vanishing and exploding gradients largely addressed by normalized initialization and batch normalization, which let networks with tens of layers start converging. What remained was the degradation problem: deeper plain networks reached higher training error than shallower ones. The authors argued that this optimization difficulty is unlikely to be caused by vanishing gradients: their plain networks used batch normalization, which keeps forward signals at non-zero variance, and they state that they checked the backward gradients and found healthy norms (the CVPR paper reports this check in words, without a plot). They conjectured that deep plain networks may have exponentially low convergence rates; in a footnote they report that three times as many training iterations did not remove the degradation. Residual connections were their answer to it.

Remedies

• Non-saturating activations. ReLU and its variants, covered in the activation function article, remove the activation’s share of the shrinkage on active paths.
• Variance-preserving initialization. Choosing $\operatorname{Var}[w]$ so that $\chi \approx 1$, Glorot and Bengio’s normalized initialization for tanh-like units and He et al.’s for ReLU, keeps the gradient’s scale constant at initialization; the weight initialization article derives both.
• Normalization layers. Batch normalization, one of the normalization layers, re-standardizes every channel during training; see batch normalization.
• Residual connections. A block’s Jacobian is $I + \partial \mathcal{F} / \partial \mathbf{x}$, so the identity term carries the gradient through any number of blocks (residual connection).
• Gradient clipping for exploding gradients. Pascanu et al. (2013) rescale the gradient $\mathbf{g}$ whenever its norm exceeds a threshold $\tau$, $\mathbf{g} \leftarrow (\tau / \lVert \mathbf{g} \rVert)\,\mathbf{g}$, which keeps the direction and bounds the step. They describe it as similar to Tomas Mikolov’s earlier element-wise clipping, and suggest setting $\tau$ from the average gradient norm over many updates; in their experience, values from half to ten times that average still converged.

Examples

The PyTorch code below (run with PyTorch 2.14.1) builds a plain 20-layer convolutional network with 16 channels of $3 \times 3$ filters, no normalization and no shortcuts, followed by global average pooling and a linear classifier. The inputs are 64 synthetic $16 \times 16$ images of horizontal or vertical stripes with a random period of 3 to 6 pixels, a random phase and noise. Weights are drawn from $\mathcal{N}(0, g^2 / n_\text{in})$ with fan-in $n_\text{in} = 9\,C_\text{in}$, so the gain $g = \sqrt{2}$ is He et al.’s choice for ReLU. Here $n_\text{out} = n_\text{in}$, so the predicted per-layer factor is $\sqrt{\chi} = g \sqrt{\overline{\phi'(z)^2}}$, with the mean squared slope measured on the network’s own pre-activations. The script backpropagates one cross-entropy loss and prints the norm of $\partial L / \partial \mathbf{z}^{(\ell)}$ at several layers, its geometric-mean ratio per layer over layers 2 to 19, the prediction, and the weight-gradient norms.

import math
import torch
import torch.nn.functional as F

torch.set_num_threads(1)
torch.set_default_dtype(torch.float64)

def stripes(n, size=16, gen=None):
    """Class 0: horizontal stripes, class 1: vertical stripes (random period and phase, plus noise)."""
    y = torch.randint(0, 2, (n,), generator=gen)
    period = torch.randint(3, 7, (n, 1, 1), generator=gen).double()
    phase = torch.rand(n, 1, 1, generator=gen) * 2 * math.pi
    r = torch.arange(size).double()
    rows = torch.sin(2 * math.pi * r[None, :, None] / period + phase).expand(n, size, size)
    cols = torch.sin(2 * math.pi * r[None, None, :] / period + phase).expand(n, size, size)
    x = torch.where(y[:, None, None] == 0, rows, cols) + 0.3 * torch.randn(n, size, size, generator=gen)
    x = (x - x.mean()) / x.std()                      # standardize pixels
    return x[:, None], y                              # (N, 1, H, W) images, labels

def layer_gradients(act, gain, depth=20, C=16, seed=9):
    """Plain CNN of `depth` 3x3 conv layers (no shortcuts, no normalization) at initialization.
    Weights ~ N(0, gain^2 / fan_in), fan_in = 9 * C_in, zero biases. Returns, per layer,
    ||dL/dz|| (gradient w.r.t. the pre-activation feature map), ||dL/dW||, and mean phi'(z)^2."""
    gen = torch.Generator().manual_seed(seed)
    x, y = stripes(64, gen=gen)
    phi = {"sigmoid": torch.sigmoid, "relu": torch.relu}[act]
    h, zs, Ws = x, [], []
    for l in range(depth):
        c_in = 1 if l == 0 else C
        W = (torch.randn(C, c_in, 3, 3, generator=gen) * gain / math.sqrt(9 * c_in)).requires_grad_()
        z = F.conv2d(h, W, padding=1)
        z.retain_grad()
        zs.append(z); Ws.append(W)
        h = phi(z)
    head = torch.randn(2, C, generator=gen) / math.sqrt(C)
    logits = h.mean(dim=(2, 3)) @ head.T              # global average pooling + linear classifier
    F.cross_entropy(logits, y).backward()
    dz = [z.grad.norm().item() for z in zs]
    dW = [W.grad.norm().item() for W in Ws]
    if act == "sigmoid":
        d2 = [(torch.sigmoid(z) * (1 - torch.sigmoid(z))).pow(2).mean().item() for z in zs]
    else:
        d2 = [(z > 0).double().mean().item() for z in zs]
    return dz, dW, d2

he = math.sqrt(2)                                     # He et al.'s gain for ReLU
configs = [("sigmoid", 1.0), ("sigmoid", 4.0), ("relu", 0.5 * he), ("relu", he), ("relu", 1.5 * he)]
print(f"torch {torch.__version__}")
print("config          ||dL/dz|| at layer 1      10      19      20 | ratio/layer  predicted | ||dL/dW|| at layer 2      10      20")
for act, gain in configs:
    dz, dW, d2 = layer_gradients(act, gain)
    ratio = (dz[1] / dz[18]) ** (1 / 17)              # geometric mean of ||dL/dz_l|| / ||dL/dz_l+1||, layers 2..19
    pred = math.sqrt(gain ** 2 * sum(d2[1:18]) / 17)  # sqrt(fan_in * Var[w] * mean phi'(z)^2)
    print(f"{act:7s} gain {gain:4.2f} {dz[0]:18.1e} {dz[9]:7.1e} {dz[18]:7.1e} {dz[19]:7.1e} | {ratio:11.2f} {pred:10.2f} |"
          f" {dW[1]:18.1e} {dW[9]:7.1e} {dW[19]:7.1e}")

print("\nratio/layer over seeds 0-9 (min-max):")
for act, gain in configs:
    rs = []
    for s in range(10):
        dz, _, _ = layer_gradients(act, gain, seed=s)
        rs.append((dz[1] / dz[18]) ** (1 / 17))
    print(f"{act:7s} gain {gain:4.2f}: {min(rs):.2f}-{max(rs):.2f}")

Output:

torch 2.14.1
config          ||dL/dz|| at layer 1      10      19      20 | ratio/layer  predicted | ||dL/dW|| at layer 2      10      20
sigmoid gain 1.00            3.3e-16 3.4e-10 3.3e-04 1.4e-03 |        0.22       0.23 |            1.8e-13 5.0e-08 2.7e-01
sigmoid gain 4.00            3.5e-08 5.8e-06 5.6e-04 8.6e-04 |        0.59       0.65 |            5.1e-06 7.6e-04 2.5e-01
relu    gain 0.71            7.8e-09 6.1e-06 2.2e-03 3.8e-03 |        0.50       0.49 |            5.7e-07 7.9e-07 2.1e-07
relu    gain 1.41            4.2e-03 6.6e-03 4.5e-03 3.9e-03 |        1.00       0.98 |            3.2e-01 5.8e-01 3.7e-01
relu    gain 2.12            1.3e+01 5.5e-01 9.9e-03 5.7e-03 |        1.49       1.48 |            1.1e+03 2.9e+03 2.8e+03

ratio/layer over seeds 0-9 (min-max):
sigmoid gain 1.00: 0.22-0.24
sigmoid gain 4.00: 0.57-0.64
relu    gain 0.71: 0.48-0.55
relu    gain 1.41: 0.95-1.11
relu    gain 2.12: 1.43-1.66

The sigmoid network’s per-layer factor sits just below $g/4 = 0.25$, the value if every unit had the sigmoid’s maximum slope, and unlike the ReLU networks its weight gradients vanish along with the error signal, because sigmoid activations do not shrink going forward and cannot compensate. Quadrupling the sigmoid weights would allow a factor of 1 if every unit sat at $z = 0$, but the predicted factor implies a mean squared slope of $(0.65/4)^2 \approx 0.026$, well below the maximum $1/16$, because larger pre-activations reach the flat tails. With ReLU the per-layer factor follows the gain alone. The half-scale ReLU network is Glorot and Bengio’s observation in miniature: its weight gradients are of the same order in every layer while $\partial L / \partial \mathbf{z}$ vanishes, yet all of them are about a millionth of those at He et al.’s scale. The variance prediction is within 2% of the measured factor for ReLU; for the sigmoid it overestimates the factor, by about 10% at gain 4, where saturation is strongest. The table shows seed 9; over seeds 0–9 the factor varies by a few hundredths to about 0.2, and the five printed ranges do not overlap.

Common Misconceptions

• “Small back-propagated gradients mean the early weights get small updates relative to the late ones.” The weight gradient multiplies the error signal by the layer’s input activations; when activations shrink forward as fast as error signals shrink backward, all layers get gradients of the same size, though possibly all far too small, as in the half-scale ReLU network above.
• “Gradient clipping also fixes vanishing gradients.” Clipping only shrinks gradients that are too large; it cannot enlarge a gradient that has vanished. Pascanu et al. (2013) proposed clipping for explosion and a separate regularization term for the vanishing case.
• “It is a property of the architecture.” The same 20-layer ReLU network vanishes, stays balanced or explodes depending only on the weight scale.

Where It Is Used

The analysis explains several standard design choices in vision networks: ReLU-family activations instead of sigmoids in hidden layers, fan-in-scaled initialization (see training a neural network for its place in the workflow), batch normalization after convolutions, and residual blocks in very deep networks such as ResNet. Monitoring per-layer gradient norms is a debugging check when a deep network fails to train, and gradient clipping, proposed in this form for recurrent networks, bounds the step when a gradient spikes.