The Vanishing and Exploding Gradient Problem
Imagine training a 100-layer network and watching the loss refuse to move. The weights in the early layers barely change, epoch after epoch. Meanwhile, the later layers overfit wildly. This is not a bug in your code. It is the vanishing gradient problem, one of the oldest and most consequential obstacles in deep learning. Its mirror image, the exploding gradient problem, produces the opposite symptom. Losses suddenly turn to `NaN` after looking perfectly healthy for a hundred steps. Both problems stem from the same arithmetic fact. Backpropagation computes gradients for early layers by multiplying together many partial derivatives. There is one per layer, via the chain rule. When those per-layer factors are consistently less than one, the product shrinks geometrically with depth. When they are consistently greater than one, the product grows geometrically with depth. Depth, the very property that makes neural networks expressive, is also what makes this multiplication dangerous. This is not a niche concern. It is the central reason plain deep networks were considered untrainable through the 1990s and early 2000s. It directly motivated the design of Activation Functions: Sigmoid, Tanh, ReLU alternatives. It also shaped careful Weight Initialization Strategies, normalization layers, and residual connections. Recurrent networks unrolled over long sequences suffer an especially acute version of the problem. This happens since the "depth" is the sequence length, which can be in the thousands. This article gives a precise mathematical account of why gradients vanish or explode. It works through the chain-rule computation layer by layer. It quantifies the effect with an explicit worked example. It also surveys the fixes that make modern deep networks trainable. By the end you should be able to look at an activation function or a weight matrix. You should be able to predict, before ever running the code, whether it will cause trouble at depth 100.
- Backpropagation multiplies per-layer Jacobians via the chain rule; a deep network's gradient is a long product of matrices.
- If the typical singular value of each layer's Jacobian is `< 1`, the product shrinks exponentially with depth — vanishing gradients.
- If it is `> 1`, the product grows exponentially — exploding gradients, often surfacing as `NaN` losses.
- Saturating activations (sigmoid, tanh) have derivatives `≤ 0.25` or `≤ 1`, making vanishing the default failure mode for deep sigmoid/tanh networks.
- Fixes include ReLU-family activations, careful weight initialization (Xavier, He), normalization layers, residual/skip connections, and gradient clipping.
- Recurrent networks face the same problem across time steps, which motivated gated architectures like LSTM and GRU.
A deep network's gradient is a product of many per-layer Jacobians, so its magnitude grows or shrinks exponentially with depth unless each factor stays close to unit scale.
The idea
Consider a feedforward network with `L` layers, where layer `ℓ` computes `hℓ = f(Wℓ hℓ-1 + bℓ)` for some activation function `f`. Training uses Backpropagation to compute the gradient of the loss `E` with respect to an early-layer weight `W1`. The chain rule expresses this as a product over every intervening layer: ``` ∂E/∂W1 = ∂E/∂hL · ∂hL/∂hL-1 · ∂hL-1/∂hL-2 · ... · ∂h2/∂h1 · ∂h1/∂W1 ``` Each factor `∂hℓ/∂hℓ-1` is a Jacobian matrix, equal to `diag(f′(zℓ)) · Wℓ`, where `zℓ = Wℓ hℓ-1 + bℓ` is the pre-activation. The gradient reaching layer 1 is the product of `L-1` such Jacobians. This is the crux of the entire phenomenon: a deep network's early-layer gradient is not a single derivative but a matrix product with as many factors as there are layers between the loss and that weight. The vanishing gradient problem occurs when the typical operator norm (loosely, the largest singular value) of each Jacobian factor is smaller than 1. The product of `L-1` numbers each less than 1 decays geometrically, so `‖∂E/∂W1‖ ≈ ρ^(L-1) ‖∂E/∂WL‖` for some `ρ < 1`. With `ρ = 0.25` and `L = 30`, the ratio `ρ^29` is on the order of `10^-18` — the early-layer gradient is, for all numerical purposes, zero. The layer's weights stop updating and the network effectively stops learning in its early layers, even though later layers may continue to adapt. The exploding gradient problem is the same mechanism with `ρ > 1`. The product grows as `ρ^(L-1)`, and gradients can overflow floating-point range, producing `NaN` or `Inf` values that corrupt every subsequent parameter update. Both failure modes are two sides of a single multiplicative process: repeated composition of near-linear maps whose "gain" is not exactly 1.
How it works
To make this concrete, decompose each Jacobian factor into two pieces: the derivative of the activation function, and the weight matrix itself.
Activation derivative saturation. For sigmoid `σ(z) = 1/(1+e^-z)`, the derivative is `σ′(z) = σ(z)(1-σ(z))`, which is maximized at `z=0` with value `0.25` and decays to `0` as `|z|` grows. For tanh, `tanh′(z) = 1 - tanh²(z)`, maximized at `1` at `z=0` but also decaying to `0` at the tails. In both cases, once pre-activations drift away from zero (which happens naturally as signals propagate through layers), the derivative shrinks toward zero — this is "saturation." Multiplying `diag(f′(z))` into the Jacobian directly caps its scale, and for sigmoid the cap is a hard `0.25` even at the best possible operating point.
Weight matrix scale. The other factor, `Wℓ`, contributes its own singular values. If weights are initialized with too small a variance, `Wℓ` shrinks the signal (compounding the activation-derivative shrinkage); if initialized too large, `Wℓ` can amplify the signal enough to counteract or overpower the activation derivative, potentially producing explosion instead.
Step-by-step mechanism during backpropagation: ``` ALGORITHM: Backpropagation gradient flow (illustrating vanishing/exploding) Input: L layers, weights W_1..W_L, activations f, loss E, forward pre-activations z_1..z_L and outputs h_1..h_L already computed. 1. Compute the output-layer gradient: delta_L = dE/dh_L * f'(z_L) // elementwise product 2. For ell = L-1 down to 1:
# backpropagate through the affine map of layer ell+1 g = W_{ell+1}^T * delta_{ell+1} # this is dE/dh_ell
# backpropagate through the activation of layer ell delta_ell = g * f'(z_ell) # elementwise product
# accumulate norm ratio to monitor vanishing/exploding ratio = norm(delta_ell) / norm(delta_{ell+1}) record(ell, ratio) 3. For each layer ell, the weight gradient is: dE/dW_ell = delta_ell * h_{ell-1}^T Output: dE/dW_1 ... dE/dW_L
# Vanishing: ratio << 1 consistently => delta_1 exponentially small
# Exploding: ratio >> 1 consistently => delta_1 exponentially large / NaN ``` The key diagnostic in this pseudocode is the `ratio` computed at each layer: it is exactly the local "gain" of the backward pass through one layer. If you log `ratio` across layers during training and see it consistently below 1 (say, averaging `0.2`–`0.3` for sigmoid networks), the gradient signal is vanishing. If it consistently exceeds 1 (common with poorly scaled weights or `tanh` at small pre-activation with large weight variance), the gradient signal is exploding. The diagram below shows this as a computational graph, with the multiplicative shrinkage annotated at each backward edge.
Complexity & correctness
Time and space complexity. Backpropagation itself costs `O(L · d²)` time per example for `L` fully connected layers of width `d` (matrix-vector products dominate), and `O(L · d)` space to store activations for the backward pass — this is unchanged by whether gradients vanish or explode. The vanishing/exploding phenomenon is not a complexity issue in the algorithmic sense; it is a numerical conditioning issue. The cost of computing the gradient is identical whether that gradient is informative or numerically useless. The practical cost it imposes is on convergence: an ill-conditioned gradient signal means gradient descent needs vastly more iterations (or never converges within floating-point precision) to update early layers meaningfully.
Correctness argument — why the product decays/grows geometrically. Write the backward recursion as `δℓ = diag(f′(zℓ)) · Wℓ+1ᵀ · δℓ+1`. Unrolling this recursion from layer `L` down to layer `1` gives ``` δ1 = [ Π_{ℓ=1}^{L-1} diag(f'(z_ℓ)) · W_{ℓ+1}ᵀ ] · δL ``` Take operator norms and use submultiplicativity, `‖AB‖ ≤ ‖A‖‖B‖`: ``` ‖δ1‖ ≤ ‖δL‖ · Π_{ℓ=1}^{L-1} ‖diag(f'(z_ℓ))‖ · ‖W_{ℓ+1}‖ ``` If every factor `‖diag(f′(zℓ))‖ · ‖Wℓ+1‖` is bounded above by some constant `ρ < 1` (true, for instance, whenever activations saturate and weights are not compensating), then ``` ‖δ1‖ ≤ ρ^{L-1} · ‖δL‖ → 0 as L → ∞ ``` This is a rigorous upper bound showing vanishing is guaranteed whenever the per-layer gain is uniformly below 1 — it is not merely empirical folklore. Symmetrically, if each factor is bounded below by some `ρ > 1` (true when weight matrices are large relative to the activation derivative, e.g., initialized with too-large variance), a matching lower bound on the smallest singular value gives `‖δ1‖ ≥ ρ^{L-1} ‖δL‖ → ∞`, proving explosion. The two regimes are separated by the razor's-edge condition `ρ ≈ 1` at every layer — this is precisely the design target of variance-preserving initialization schemes (foreshadowing Xavier and He Initialization).
A useful sanity computation. For sigmoid activations, `‖diag(f′(z))‖ ≤ 0.25` always, regardless of the weights. So no matter how weights are initialized, the activation-derivative factor alone caps the per-layer gain unless `‖Wℓ‖ > 4` on average — and such large weights themselves risk instability and poor optimization elsewhere (e.g., saturating the next layer's sigmoid even harder, since larger pre-activations push `|z|` further from 0). This is why sigmoid/tanh networks are fundamentally harder to train at depth than ReLU networks: ReLU's derivative is exactly 1 for all `z > 0`, so the activation-derivative factor does not itself contribute shrinkage on the active half of its domain — the entire per-layer gain is then controlled by `Wℓ` alone, which is exactly what principled initialization schemes are designed to keep near 1.
Honest tradeoffs. ReLU removes activation-derivative shrinkage but introduces a new failure mode — "dying ReLUs," where a unit's pre-activation is permanently negative and its gradient is exactly zero forever, a distinct (non-multiplicative) pathology. Gradient clipping controls explosion but does nothing for vanishing, and can bias the direction of very large gradients if applied naively (clip by norm, not by value, to preserve direction). Residual connections mitigate vanishing by providing an additive identity path with derivative exactly 1, so the product-of-Jacobians argument no longer applies in the same way — but they do not eliminate the underlying tension between depth and conditioning, only reroute around it.
Worked example
Take a toy 5-layer sigmoid network, scalar for clarity (`d=1`), so all "matrices" are numbers. Suppose weights are all `w = 1.5` and pre-activations happen to sit near `z = 0.5` at every layer (a realistic post-initialization regime), where `σ′(0.5) ≈ 0.235`. The per-layer backward gain is `g = σ′(z) · w ≈ 0.235 × 1.5 = 0.353`.
| Layer (from output) | Backward gradient magnitude | Cumulative gain `g^k` |
|---|---|---|
| δ5 (output) | 1.000 | `g^0 = 1.000` |
| δ4 | 0.353 | `g^1 = 0.353` |
| δ3 | 0.124 | `g^2 = 0.125` |
| δ2 | 0.044 | `g^3 = 0.044` |
| δ1 (input layer) | 0.0155 | `g^4 = 0.0156` |
By the time the gradient reaches layer 1, it has shrunk to about 1.5% of its output-layer magnitude after only 4 multiplications. Extrapolate this same `g ≈ 0.353` to a 50-layer network: `0.353^49 ≈ 3.4 × 10^-22`, a number that underflows to exact zero in single-precision floating point. The early layers receive literally no gradient signal — training freezes there regardless of how long you run gradient descent. Now flip the scenario: keep `σ′(z) ≈ 0.235` but set `w = 8` (an aggressively large initialization). The gain becomes `g = 0.235 × 8 = 1.88`. Over 20 layers, `1.88^19 ≈ 1.4 × 10^5` — a 100,000-fold amplification of the output gradient by the time it reaches early layers, easily enough to push weight updates into instability or numerical overflow within a handful of training steps. This single scalar example, with only the constant `w` and operating point `z` changed, reproduces both named failure modes from one shared formula `g^{L-1}`.
Common pitfalls & traps
| Pitfall | Why it happens | Fix |
|---|---|---|
| Blaming the optimizer for a dead early layer | Gradients underflow to ~0 before reaching that layer; no learning rate fixes a zero gradient | Check per-layer gradient norms; switch activation/init, or add skip connections |
| Clipping gradients by raw value, not norm | Element-wise clipping distorts the gradient's direction, biasing updates | Clip by global norm (`clip = g · min(1, threshold/‖g‖)`) |
| Using sigmoid/tanh in very deep stacks "because it worked shallow" | Saturation caps derivative at ≤0.25 (sigmoid) or ≤1 (tanh); compounds multiplicatively with depth | Use ReLU/GELU family for hidden layers at depth |
| Assuming ReLU "solves" vanishing gradients | ReLU removes derivative-shrinkage but weight-matrix scale can still shrink or explode the signal; dying ReLUs are a separate zero-gradient trap | Pair ReLU with variance-aware init and monitor for dead units |
| Confusing exploding gradients with a bad learning rate | Both cause huge, unstable updates and diverging loss, but the root cause differs | Log raw gradient norm before scaling by learning rate to distinguish the two |
| Treating RNN "depth" as small because the architecture looks shallow | An RNN unrolled over `T` time steps has the same multiplicative structure with `T` factors, often `T ≫ L` | Use gated units (LSTM/GRU) or truncated backprop through time |
Pioneers & history
The vanishing gradient problem was first rigorously identified and analyzed in the context of recurrent networks, where the effective "depth" from unrolling over time made the problem especially severe and immediately visible in experiments through the early 1990s.
Origins · 1 minSepp Hochreiter and how this was discovered
Where this matters
Every framework that trains networks beyond a handful of layers must contend with this problem, whether the developer thinks about it explicitly or not. Modern image classifiers such as ResNet and its descendants use residual/skip connections specifically because they provide a gradient path with derivative exactly 1, bypassing the geometric decay argument regardless of how deep the network becomes. Transformer architectures, now the backbone of large language models, use residual connections combined with layer normalization for the same reason — without them, stacking dozens of attention/feedforward blocks would be untrainable. Recurrent architectures used in speech recognition and older sequence-to-sequence translation systems relied on LSTM and GRU gating specifically to give gradients a near-linear, non-shrinking path across long time horizons. Any deep learning framework's default weight initializers (Xavier/Glorot, He) are direct engineering responses to the mathematics in this article — they choose initial weight variance precisely to keep the per-layer gain close to 1 at the start of training. Practical tools like gradient clipping in speech and NLP training pipelines, and batch/layer normalization nearly everywhere, exist to keep the multiplicative chain from drifting away from unit gain as training proceeds.
Ask the tutor: How does He initialization derive its specific variance formula from this gain-product argument, and why does it differ from Xavier initialization for ReLU networks?
Is the vanishing gradient problem only about sigmoid and tanh activations?
No — the underlying cause is the product of many per-layer Jacobians shrinking with depth, and activation saturation is only one contributor to that shrinkage. Even with ReLU, whose derivative is exactly 1 on the active region, weight matrices with small singular values can still cause the overall product to shrink. Conversely, ReLU networks with poorly scaled large weights can still explode. The activation function determines one factor in the per-layer gain; the weight matrix scale determines the other. Both must be controlled — this is precisely why initialization schemes like Xavier and He are paired with specific activation choices rather than treated independently.
How do I detect vanishing or exploding gradients in practice?
Log the norm of the gradient at each layer (e.g., `‖∂E/∂Wℓ‖`) during training and plot it against layer index. A healthy network shows roughly comparable gradient norms across layers, perhaps within an order of magnitude. Vanishing shows as gradient norms decaying by orders of magnitude from output layers to input layers, often reaching numbers so small they round to zero in floating point. Exploding shows as gradient norms growing sharply toward the input layers or suddenly spiking to very large values or `NaN` partway through training. Many frameworks expose per-parameter gradient statistics directly, making this a cheap diagnostic to run routinely.
Why are recurrent neural networks especially prone to this problem?
An RNN applies the same weight matrix and activation repeatedly across time steps, and backpropagation through time unrolls this into a computational graph as deep as the sequence is long. A sequence of length 200 is, from the gradient's perspective, a 200-layer network sharing the same weights at every layer. Because natural language and time-series data routinely involve dependencies spanning hundreds of steps, the multiplicative shrinkage argument applies with a much larger effective depth than typical feedforward networks. This is exactly why gated architectures like LSTM and GRU, which include an additive memory-cell path with a near-constant gradient gain, were developed specifically for sequential data.
Does batch normalization solve the vanishing/exploding gradient problem?
Batch normalization helps substantially but does not eliminate the underlying multiplicative structure. By normalizing pre-activations to have controlled mean and variance at each layer, it keeps activation derivatives away from the extreme-saturation regions and stabilizes the scale of signals flowing forward and gradients flowing backward. This makes the per-layer gain more consistently close to 1 throughout training, even as weights change. However, it is typically combined with good initialization and, in very deep networks, with residual connections, because normalization alone does not guarantee the product of per-layer gains stays bounded across arbitrarily many layers.
Can vanishing gradients happen even with ReLU and good initialization?
Yes, though it is far less common and typically less severe. If a large fraction of ReLU units become inactive (dying ReLU), their local gradient is exactly zero, which can produce vanishing gradients for those specific pathways even though the aggregate per-layer Jacobian gain might look acceptable elsewhere. In very deep plain (non-residual) ReLU networks, small systematic biases in weight scale can still compound across dozens of layers, since even a factor of `0.95` per layer becomes `0.95^100 ≈ 0.006` at depth 100. This is why extremely deep architectures still favor residual connections over relying on initialization and activation choice alone.
What is the practical difference between fixing vanishing gradients versus fixing exploding gradients?
Exploding gradients are usually addressed reactively during training: gradient norm clipping caps the update magnitude whenever gradients spike, without changing the architecture. Vanishing gradients generally require proactive architectural or initialization changes, because once a gradient underflows to numerical zero there is no way to "clip it back up" — the information is already lost. This asymmetry is why the deep learning toolkit for exploding gradients (clipping, careful learning-rate schedules) is largely orthogonal to the toolkit for vanishing gradients (activation choice, initialization variance, normalization, residual connections), even though both stem from the same multiplicative chain-rule mechanism.
How does this connect to weight initialization strategies?
The per-layer backward gain is the product of the activation derivative and the weight matrix's singular values. Since the activation derivative is fixed once you choose an activation function, the remaining lever available at the start of training is the initial variance of the weights. Xavier initialization chooses weight variance to keep signal and gradient variance approximately constant across layers for tanh/sigmoid-like activations, while He initialization adjusts this analysis for ReLU's asymmetric, half-zero derivative. Neither guarantees stability throughout all of training, since weights evolve, but both directly target the `ρ ≈ 1` condition derived in this article's correctness argument.