Original Substack post

Have you ever noticed that in a deep learning codebase, weights are rarely initialized as plain $\mathcal{N}(0,1)$?

You usually see something like this instead:

$$W_{ij}\sim \mathcal{N}\left(0,\frac{1}{n_{in}}\right)$$

or sometimes:

$$W_{ij}\sim \mathcal{N}\left(0,\frac{2}{n_{in}}\right)$$

In code this often looks like drawing normal random numbers and dividing by the square root of the input dimension.

Why do we do this?

Because if we initialize every weight with variance $1$, then as signals move forward through many layers, and as gradients move backward through those same layers, their scale can multiply again and again. After enough multiplications, the numbers either become huge or almost zero.

That is the exploding and vanishing gradients problem.

This post is the mathematics behind it.

The short version is:

Deep learning is repeated multiplication. Repeated multiplication is controlled by variance, derivatives, Jacobians, and singular values.

Once that sentence becomes obvious, Xavier initialization, He initialization, residual connections, normalization, and gradient clipping all start feeling like the same family of fixes.

---

Start With One Neuron

Consider a single neuron in one layer. Ignore the bias for now.

$$y=w_1{x}_1+w_2{x}_2+\cdots +w_n{x}_n$$

or compactly:

$$y=\sum _{i=1}^n{w}_i{x}_i$$

Here $n$ is the number of input dimensions, usually called $n_{in}$.

Let us make three simple assumptions:

1. The inputs are normalized:

$$\mathbb{E}[x_i]=0,Var(x_i)=1$$

1. The weights are initialized independently:

$$\mathbb{E}[w_i]=0,Var(w_i)={\sigma }_w^2$$

1. The $x_i$ and $w_i$ are independent of each other.

We want to know the variance of $y$.

Since each term is $w_i{x}_i$:

$$Var(w_i{x}_i)=\mathbb{E}[w_i^2{x}_i^2]-\mathbb{E}[w_i{x}_i{]}^2$$

Because the weights and inputs are independent and both have mean zero:

$$\mathbb{E}[w_i{x}_i]=0$$

and:

$$\mathbb{E}[w_i^2{x}_i^2]=\mathbb{E}[w_i^2]\mathbb{E}[x_i^2]$$

Since $\mathbb{E}[w_i^2]=Var(w_i)={\sigma }_w^2$ and $\mathbb{E}[x_i^2]=Var(x_i)=1$:

$$Var(w_i{x}_i)={\sigma }_w^2$$

Now $y$ is the sum of $n$ independent terms, so variances add:

$$Var(y)=\sum _{i=1}^{n}Var(w_i{x}_i)$$ $$Var(y)=n{\sigma }_w^2$$

This is the whole problem in one line.

If we initialize with ${\sigma }_w^2=1$, then:

$$Var(y)=n$$

So if the layer has 768 inputs, the output variance becomes about 768 after just one linear layer. That is already too large.

To keep the output variance around $1$, we want:

$$n{\sigma }_w^2=1$$

so:

$${\sigma }_w^2=\frac{1}{n}$$

and therefore:

$${\sigma }_w=\frac{1}{\sqrt{n}}$$

This is the basic statistical reason behind dividing by $\sqrt{n_{in}}$.

---

Why This Becomes Worse Across Layers

One layer is not the problem. Deep learning means stacking many layers.

Suppose a simple feed-forward network is:

$$h_0=x$$ $$z_l=W_l{h}_{l-1}$$ $$h_l=\varphi (z_l)$$

where $\varphi$ is an activation function like tanh, sigmoid, ReLU, or GELU.

If we temporarily ignore the activation, and each layer multiplies variance by some factor $c$, then after $L$ layers:

$$Var(h_L)\approx c^{L}Var(x)$$

If $c>1$, activations explode.

If $c<1$, activations vanish.

If $c\approx 1$, the network has a chance of preserving signal scale.

This is why tiny errors in scale become massive in deep networks. A factor of $1.2$ seems harmless for one layer, but over 100 layers:

$$1.2^{100}\approx 82,817$$

A factor of $0.8$ seems harmless too, but:

$$0.8^{100}\approx 2.0\times 10^{-10}$$

So a deep network does not merely care whether a layer is mathematically valid. It cares whether repeated composition keeps signals in a usable numerical range.

---

Backpropagation Is The Same Problem In Reverse

Now let us look at gradients.

Suppose the loss is $\mathcal{L}$ and we want gradients with respect to earlier layers.

Backpropagation repeatedly applies the chain rule.

For a network:

$$h_L=f_L(f_{L-1}(\cdots f_1(x)))$$

we get:

$$\frac{\partial \mathcal{L}}{\partial h_{l-1}}={\frac{\partial h_l}{\partial h_{l-1}}}^T\frac{\partial \mathcal{L}}{\partial h_l}$$

The object:

$$J_l=\frac{\partial h_l}{\partial h_{l-1}}$$

is the Jacobian of layer $l$.

So the gradient at an early layer is:

$$\frac{\partial \mathcal{L}}{\partial h_0}=J_1^T{J}_2^T\cdots J_L^T\frac{\partial \mathcal{L}}{\partial h_L}$$

Again, repeated multiplication.

This time we are multiplying Jacobian matrices instead of scalar variances.

If this product grows in norm, gradients explode. If this product shrinks in norm, gradients vanish.

That is the mathematical heart of the issue.

---

The Jacobian Of A Layer

For one layer:

$$z_l=W_l{h}_{l-1}$$ $$h_l=\varphi (z_l)$$

The Jacobian is:

$$J_l=\frac{\partial h_l}{\partial h_{l-1}}$$

Using the chain rule:

$$J_l=D_l{W}_l$$

where:

$$D_l=diag({\varphi }^{\prime }(z_l))$$

So each layer's Jacobian has two parts:

1. $W_l$, the weight matrix.
2. $D_l$, the diagonal matrix of activation derivatives.

This is important because gradients can vanish or explode because of either part.

Bad weight scale can break training.

Bad activation derivatives can also break training.

Usually both interact.

---

Singular Values: The Real Scale Of A Matrix

For scalars, repeated multiplication is easy to understand. If you multiply by $2$ many times, you explode. If you multiply by $0.5$ many times, you vanish.

For matrices, the closest idea is singular values.

A matrix $A$ stretches some directions and shrinks others. Its largest singular value tells us the maximum stretch:

$$\Vert Ax\Vert \le {\sigma }_{max}(A)\Vert x\Vert$$

If the Jacobians have singular values much larger than $1$, gradients can explode.

If the Jacobians have singular values much smaller than $1$, gradients vanish.

For many layers:

​J1T​J2T​⋯JLT​g​

is controlled by the product of the singular values along the directions the gradient travels.

A rough bound is:

​J1T​J2T​⋯JLT​g​≤(l=1∏L​σmax​(Jl​))∥g∥

So even if every layer has maximum singular value only $1.1$, after 100 layers the upper bound contains:

$$1.1^{100}\approx 13,780$$

And if every layer has singular values around $0.9$:

$$0.9^{100}\approx 0.000026$$

This is why people often talk about keeping Jacobian singular values near $1$. The ideal is sometimes called dynamical isometry: signals and gradients pass through depth without being crushed or inflated too much.

---

Eigenvalues: The RNN Version

Eigenvalues become especially intuitive in recurrent neural networks.

A simple RNN has hidden state:

$$h_t=\varphi (W_h{h}_{t-1}+W_x{x}_t+b)$$

The same recurrent matrix $W_h$ is reused at every time step.

During backpropagation through time, the gradient from time $T$ back to time $t$ contains products like:

$$\frac{\partial h_T}{\partial h_t}=\prod _{k=t+1}^T\frac{\partial h_k}{\partial h_{k-1}}$$

For the RNN:

$$\frac{\partial h_k}{\partial h_{k-1}}=D_k{W}_h$$

where:

$$D_k=diag({\varphi }^{\prime }(W_h{h}_{k-1}+W_x{x}_k+b))$$

So:

$$\frac{\partial h_T}{\partial h_t}=(D_T{W}_h)(D_{T-1}{W}_h)\cdots (D_{t+1}{W}_h)$$

If we ignore the changing activation derivative for a moment, the core repeated object is:

$$W_h^{T-t}$$

Now eigenvalues tell the story.

If $W_h$ has an eigenvector $v$ with eigenvalue $\lambda$, then:

$$W_{h}v=\lambda v$$

and after repeated multiplication:

$$W_h^{k}v={\lambda }^{k}v$$

If $|\lambda |>1$, that direction explodes.

If $|\lambda |<1$, that direction vanishes.

If $|\lambda |\approx 1$, that direction can preserve information over time.

This is why vanilla RNNs struggle with long-range dependencies. To remember something from 500 time steps ago, the gradient has to survive hundreds of repeated multiplications.

That is a brutal requirement.

---

Activation Functions Can Kill Gradients Too

Weights are only half the story. Activations matter because their derivatives appear inside every Jacobian.

Sigmoid

The sigmoid is:

$$\sigma (x)=\frac{1}{1+e^{-x}}$$

Its derivative is:

$${\sigma }^{\prime }(x)=\sigma (x)(1-\sigma (x))$$

The maximum value of this derivative is $0.25$.

So every sigmoid layer contributes a derivative factor at most $0.25$.

Across many layers, this is terrible:

$$0.25^{20}\approx 9.1\times 10^{-13}$$

This is why old deep networks with sigmoid activations were difficult to train.

Also, sigmoid saturates. For large positive or negative inputs, the output becomes close to $1$ or $0$, and the derivative becomes almost zero.

So once a unit saturates, learning through it becomes tiny.

Tanh

Tanh is better than sigmoid because it is zero-centered:

$$\tanh (x)\in [-1,1]$$

Its derivative is:

$$\frac{d}{dx}\tanh (x)=1-{\tanh }^2(x)$$

The maximum derivative is $1$, but for large $|x|$, tanh also saturates and the derivative becomes close to zero.

So tanh can still vanish if activations drift into saturation.

ReLU

ReLU is:

$$\varphi (x)=\max (0,x)$$

Its derivative is:

$${\varphi }^{\prime }(x)=\{\begin{array}{ll}1,& x>0\\ 0,& x<0\end{array}$$

ReLU helps because positive units pass gradient with derivative $1$.

But roughly half the units may be inactive at initialization if pre-activations are symmetric around zero. Those inactive units have derivative $0$.

So ReLU does not shrink every active path, but it drops many paths entirely.

This is why ReLU networks usually use He initialization instead of Xavier.

---

Xavier Initialization

Xavier initialization tries to keep variance stable through layers.

For activations that are roughly symmetric and not half-zeroed like ReLU, a common choice is:

$$Var(W_{ij})=\frac{1}{n_{in}}$$

or, balancing both forward and backward flow:

$$Var(W_{ij})=\frac{2}{n_{in}+n_{out}}$$

The second form is often called Glorot or Xavier initialization.

Why include $n_{out}$?

Forward propagation cares about $n_{in}$ because each output sums over input dimensions.

Backward propagation cares about $n_{out}$ because each input receives gradient contributions from output dimensions.

So Xavier balances both directions:

$$Var(h_l)\approx Var(h_{l-1})$$

and:

$$Var\left(\frac{\partial \mathcal{L}}{\partial h_{l-1}}\right)\approx Var\left(\frac{\partial \mathcal{L}}{\partial h_l}\right)$$

It is not magic. It is variance bookkeeping.

---

He Initialization

For ReLU, about half the units are zeroed out at initialization.

So if we used:

$$Var(W_{ij})=\frac{1}{n_{in}}$$

then the variance after ReLU would roughly be cut in half.

To compensate, He initialization uses:

$$Var(W_{ij})=\frac{2}{n_{in}}$$

or:

$$W_{ij}\sim \mathcal{N}\left(0,\frac{2}{n_{in}}\right)$$

This keeps the post-ReLU signal variance closer to stable.

In code, this is why you often see:

$$std=\sqrt{\frac{2}{n_{in}}}$$

for ReLU-style networks.

For transformers using GELU, initialization is more architecture-specific, but the same principle remains: keep signal and gradient scales from drifting too fast with depth.

---

A Tiny Numerical Example

Suppose every layer accidentally multiplies gradient norm by $0.7$.

After 10 layers:

$$0.7^{10}\approx 0.028$$

After 50 layers:

$$0.7^{50}\approx 1.8\times 10^{-8}$$

After 100 layers:

$$0.7^{100}\approx 3.2\times 10^{-16}$$

That gradient is basically gone.

Now suppose every layer multiplies gradient norm by $1.3$.

After 10 layers:

$$1.3^{10}\approx 13.8$$

After 50 layers:

$$1.3^{50}\approx 497,929$$

After 100 layers:

$$1.3^{100}\approx 2.48\times 10^{11}$$

That gradient is unusably large.

So the problem is not that gradients are mysterious. The problem is that exponentials are unforgiving.

---

Why RNNs Were Especially Painful

Feed-forward networks multiply across depth.

RNNs multiply across time.

A sequence length of 1,000 means the model may need useful gradients across 1,000 recurrent steps.

For a vanilla RNN, the gradient from a late output to an early hidden state contains:

$$\prod _{k=t+1}^T{D}_k{W}_h$$

Even if the recurrent matrix is well-scaled at initialization, training changes it. And even if $W_h$ is fine, the activation derivatives can shrink gradients.

This is one reason LSTMs and GRUs were such a big deal. They create gated paths where information can flow additively, not only through repeated nonlinear matrix multiplication.

A simplified LSTM cell state update is:

$$c_t=f_t\odot c_{t-1}+i_t\odot {\tilde{c}}_t$$

The old memory $c_{t-1}$ is multiplied by a forget gate $f_t$, but it also has a relatively direct path forward.

The derivative:

$$\frac{\partial c_t}{\partial c_{t-1}}=f_t$$

If $f_t$ is near $1$, the gradient can pass through many time steps more easily.

This does not make LSTMs perfect, but it explains why they helped long before transformers became dominant.

---

Residual Connections

Residual connections attack the same problem from another angle.

Instead of making a layer learn:

$$h_{l+1}=F(h_l)$$

a residual block does:

$$h_{l+1}=h_l+F(h_l)$$

Now the Jacobian is:

$$\frac{\partial h_{l+1}}{\partial h_l}=I+\frac{\partial F}{\partial h_l}$$

The identity matrix $I$ gives gradients a direct route backward.

So even if $\frac{\partial F}{\partial h_l}$ is small, the gradient is not forced to pass only through $F$.

This is a major reason very deep ResNets and transformers train at all.

Residual connections do not remove the need for good initialization or normalization, but they make the optimization landscape much more forgiving.

---

Normalization

Normalization methods also fight scale drift.

BatchNorm normalizes using batch statistics. LayerNorm normalizes within each example across features.

A simplified LayerNorm is:

$$\mu =\frac{1}{d}\sum _{i=1}^d{x}_i$$ $${\sigma }^2=\frac{1}{d}\sum _{i=1}^d(x_i-\mu {)}^2$$ $$\text{LayerNorm}(x{)}_i={\gamma }_i\frac{x_i-\mu }{\sqrt{{\sigma }^2+ϵ}}+{\beta }_i$$

This keeps activations in a controlled range before the next transformation.

In transformers, LayerNorm is one of the quiet heroes. Without it, residual streams and attention/MLP blocks can drift in scale as depth increases.

There are two common layouts:

Post-norm:

$$h_{l+1}=\text{LayerNorm}(h_l+F(h_l))$$

Pre-norm:

$$h_{l+1}=h_l+F(\text{LayerNorm}(h_l))$$

Pre-norm transformers are usually easier to train at large depth because the residual path remains more direct.

Again, this is gradient-flow engineering.

---

Gradient Clipping

Gradient clipping is mostly a fix for exploding gradients.

If the gradient vector is $g$ and its norm is too large, we rescale it:

$$g\leftarrow g\cdot \frac{\tau }{\Vert g\Vert }$$

when:

$$\Vert g\Vert >\tau$$

Here $\tau$ is the clipping threshold.

So if the gradient norm becomes $1000$ and the threshold is $1$, clipping rescales the whole gradient to norm $1$.

This does not solve vanishing gradients. It also does not fix the underlying reason the gradient exploded. But it prevents a single bad update from destroying the model parameters.

RNN training historically used gradient clipping a lot because exploding gradients were common during backpropagation through time.

Modern large-model training also often clips gradients because occasional spikes are normal at scale.

---

A Concrete Mental Model

Imagine a signal moving through a deep network.

At each layer, three things happen:

1. The weight matrix stretches or shrinks the vector.
2. The activation derivative keeps, shrinks, or kills gradient paths.
3. The next layer repeats the process.

Forward pass:

$$h_0\to h_1\to h_2\to \cdots \to h_L$$

Backward pass:

$$\frac{\partial \mathcal{L}}{\partial h_L}\to \frac{\partial \mathcal{L}}{\partial h_{L-1}}\to \cdots \to \frac{\partial \mathcal{L}}{\partial h_0}$$

The forward pass can lose information if activations collapse or explode.

The backward pass can lose learning signal if Jacobian products shrink or explode.

Training works best when both stay in a reasonable range.

That is why initialization is not a small implementation detail. It determines the starting geometry of learning.

---

Why Plain Normal(0, 1) Is Bad

Let us return to the original question.

If:

$$W_{ij}\sim \mathcal{N}(0,1)$$

and a layer has $n_{in}$ inputs, then:

$$Var(y)=n_{in}$$

For $n_{in}=1024$:

$$Var(y)=1024$$

so the standard deviation is:

$$\sqrt{1024}=32$$

After one layer, activations are already much larger than the normalized input.

If you use tanh or sigmoid, this pushes units into saturation, where derivatives are near zero. So the forward activations become uninformative and the backward gradients vanish.

If you use ReLU, large activations can propagate and make later layers unstable. Gradients can become huge because matrix products keep stretching them.

So Normal$(0,1)$ is usually wrong not because randomness is bad, but because the randomness has the wrong scale for the number of summed inputs.

---

The Main Fixes, In One Place

Initialization controls the starting scale:

$$Var(W_{ij})\approx \frac{1}{n_{in}}$$

or for ReLU:

$$Var(W_{ij})\approx \frac{2}{n_{in}}$$

Activation choice controls derivative behavior. Sigmoid and tanh can saturate. ReLU and GELU usually preserve gradients better, but still need good scale.

Normalization keeps activations from drifting too far during training.

Residual connections create direct gradient paths:

$$\frac{\partial }{\partial h_l}(h_l+F(h_l))=I+F^{\prime }(h_l)$$

Gradient clipping prevents occasional explosions from causing catastrophic updates.

Gated architectures like LSTMs create paths where memory and gradients can survive over time.

All of these are different ways of fighting the same mathematical enemy: bad products of many Jacobians.

---

Final Intuition

Exploding and vanishing gradients are not random bugs in neural networks.

They are what you should expect when you repeatedly multiply by matrices and activation derivatives.

In a deep feed-forward network, the multiplication happens across layers.

In an RNN, it happens across time.

The gradient is basically carried by a product like:

$$J_1^T{J}_2^T\cdots J_L^T$$

or for an RNN:

$$(D_T{W}_h)(D_{T-1}{W}_h)\cdots (D_{t+1}{W}_h)$$

If the product's effective singular values are mostly above $1$, gradients explode.

If they are mostly below $1$, gradients vanish.

If they stay near $1$, learning can flow.

That is why we divide by $\sqrt{n}$.

That is why Xavier and He initialization exist.

That is why residual connections, normalization, gates, and clipping are everywhere.

Deep learning works when the model is expressive enough to learn complicated functions, but numerically disciplined enough that the learning signal can actually travel through it.