Maths Catch-Up · Variances Add, Standard Deviations Don't: Where the √n Comes From

The previous part used σ_z = √n · σ_W · σ_x to decide what W1 should be multiplied by, but the formula itself was borrowed. This part derives it from the definition of expectation, and answers a more valuable question along the way: why does diagnosing a network mean watching the variance, rather than the mean or the maximum?

Expectation, where it all starts

Expectation is the average, written E[X] or μ.

X = [2, 4, 4, 4, 5, 5, 7, 9]
E[X] = (2+4+4+4+5+5+7+9) / 8 = 5

E[X²] means “square each number first, then average” — the order matters, E[X²] ≠ (E[X])². The difference between those two is exactly the variance, as we are about to see.

What A and B stand for

The A, B and X in the formulas below are all random variables — quantities whose value is uncertain but whose distribution is known. They are placeholders, like x and y in algebra. In this setting they stand for:

zⱼ = x₁w₁ + x₂w₂ + … + x₃₀w₃₀

So each of the three rules answers one concrete question:

Using letters rather than the concrete symbols is the point: any pair works, as long as they are independent with zero mean. Writing x₁w₁ would suggest it only holds for that particular term.

“Random variable” does not mean regenerated repeatedly

This one is easy to muddle. A random variable is like a die: before the roll you don’t know the face, but you know each has probability 1/6. That does not mean W is regenerated during training.

W1 = torch.randn((30, 200))     # runs once, 6000 numbers
W1.std()                        # ≈ 1.0, which is σ_W itself

W1 is generated once, then used and updated by gradients — that is just the ordinary training loop. Probability applies because that single generation already contains 6000 samples from one distribution. After the forward pass, z holds 32×200 = 6400 values, which are naturally a random mixture across samples and neurons. What the formula predicts is how spread out those 6400 numbers are:

z.std()      # 4.916  ← this is what the formula computes

One generation, one forward pass, and the formula already matches.

σ_W is a parameter of the distribution, not one weight

A single fixed number has no variance — −0.722 is just −0.722. σ_W describes “how large numbers of this kind typically are”. That is why σ_W carries no subscript: all 30 weights come from one distribution and share one σ_W. It is also why the terms inside the summation below carry no subscript.

Two properties get used repeatedly:

E[A + B] = E[A] + E[B]        always
E[AB]    = E[A] · E[B]        only if A and B are independent

The condition on the second one matters — every derivation below rests on it.

Variance: how spread out

Two sets of numbers, both averaging 100:

set A: 99, 100, 101, 100, 100      ← packed together
set B: 40, 160, 100,  20, 180      ← spread wide

Same average, obviously not the same thing. Variance is what quantifies the difference.

Var(X) = E[(X − μ)²]

Three steps: subtract the mean → square → average. On set A (μ = 100):

Var = (1+0+1+0+0) / 5 = 0.4        σ = √0.4 = 0.632

Why square

Because adding the deviations directly makes them cancel:

E[X − μ] = E[X] − μ = 0        always zero

Squaring makes every deviation positive. Absolute value would work too, but it isn’t differentiable, and it lacks the properties below.

A useful equivalent form

Var(X) = E[X²] − (E[X])²

Derivation: E[(X−μ)²] = E[X² − 2μX + μ²] = E[X²] − 2μ² + μ² = E[X²] − μ².

What makes it work is that μ is a fixed number (the average), not a random variable, so it moves freely in and out of the expectation: E[2μX] = 2μ·E[X] = 2μ², and E[μ²] = μ².

When μ = 0, Var(X) = E[X²]. Numbers from randn have mean 0, which simplifies everything below substantially — not a coincidence, but one of the reasons initialisation deliberately centres at zero.

Checking on set A:

E[X²]   = (99² + 100² + 101² + 100² + 100²) / 5 = 10000.4
(E[X])² = 100² = 10000
Var     = 10000.4 − 10000 = 0.4      ✅ matches the definition

“Square then average” is always ≥ “average then square”, so variance is never negative; they are equal only when every number is identical.

Standard deviation: back to the original units

σ = √Var(X)

Variance is in squared units and cannot be compared against the data directly. Standard deviation shares the data’s units and can be read straight off.

Derive in variance, report in standard deviation. Because only variance adds (section four below); standard deviation does not. That is the whole reason for the division of labour.

Three rules of arithmetic

Rule 1: multiplying by a constant

Var(cX) = c² · Var(X)        σ_cX = |c| · σ_X

Derivation (with μ = 0): E[(cX)²] = c²E[X²].

This is why W1 = randn(...) * 0.2 puts the standard deviation at exactly 0.2 — the constant multiplies straight onto it, no trial and error needed. The “hand-tuned 0.2” of the previous part was not found by trying; it was computed.

Rule 2: addition (requires independence)

Var(A + B) = Var(A) + Var(B)

Derivation (both means zero):

E[(A+B)²] = E[A²] + 2·E[AB] + E[B²]
                     ↑
        under independence = E[A]E[B] = 0 × 0 = 0

The cross term vanishing is the whole thing. Counter-example, when A is perfectly correlated with itself:

Var(A + A) = Var(2A) = 4·Var(A) ≠ 2·Var(A)

Measured, Var(a+a) = 3.999, not 2.

Why the cross term vanishes

E[AB] = 0 needs two conditions at once, each doing one step:

Intuitively, E[AB] is the covariance — whether A and B tend to rise and fall together. Under independence the four sign combinations are equally likely and cancel:

E[A] = 0 but E[A²] ≠ 0, without contradiction

An easy point to muddle: if E[A] = 0, why isn’t the E[A²] in the expansion also zero?

Because A takes both signs, while A² is never negative.

First powers cancel, squares do not. And that is precisely why variance exists — if E[A²] were zero too, the variance would be zero and the whole derivation pointless. Only the middle term disappears from the expansion; the two squared terms on either side are the variances themselves.

Rule 3: multiplication (requires independence + zero mean)

Var(AB) = σ_A² · σ_B²

Derivation:

Var(AB) = E[A²B²] − (E[AB])² = E[A²]E[B²] − (E[A]E[B])²

With zero means the second term drops, and E[A²] = σ_A², giving σ_A²σ_B².

With non-zero means the formula is much messier: Var(AB) = σ_A²σ_B² + σ_A²μ_B² + σ_B²μ_A². Taking μ_A = 2, μ_B = 3, σ = 1: theory gives 1 + 9 + 4 = 14, measurement gives 13.991 ✅

All three, checked at once

a, b = torch.randn(2_000_000), torch.randn(2_000_000)

Var(a)     = 0.9997     theory 1
Var(3a)    = 8.9969     theory 9        ← rule 1
Var(a+b)   = 1.9985     theory 2        ← rule 2
Var(ab)    = 1.0007     theory 1        ← rule 3
Var(a+a)   = 3.9986     ← not independent, rule 2 does not apply

Variances add, standard deviations don’t

This is the hinge of the whole derivation.

The signs cancel, so adding 30 random numbers gives a total of 5.5× rather than 30×. A random walk is the same fact: flip a coin 100 times and you end up about √100 = 10 steps out.

So the derivation has to be done in variance space, taking the square root only at the end. The previous part’s line “the √n comes from variances adding and the standard deviation taking a root” is this section in full.

Back to the network

One neuron: zⱼ = Σᵢ xᵢwᵢⱼ, with i from 1 to n.

Step 1 — a single term (rule 3)

Var(xᵢwᵢⱼ) = σ_x² σ_W²

Step 2 — n terms added (rule 2)

Var(zⱼ) = n · σ_x² σ_W²

There are two equal signs in that step, and they are entirely different in kind:

① Var(Σ…) = ΣVar(…) — Var passing through the summation, which is rule 2 applied 29 times and only holds under independence ② ΣVar(xᵢwᵢ) = Σσ_x²σ_W² — pure algebraic substitution of step 1’s result. The right-hand side has no subscript i, because all 30 terms share a distribution and a variance

Then “the same constant added n times = times n”. Equal sign ① is the only step in the whole derivation with any content.

Checking ①:

x = torch.randn(100000, 30) * 0.911
W = torch.randn(30)
terms = x * W                                  # (100000, 30), every term

print("left  Var(∑):", terms.sum(1).var().item())
print("right ∑Var  :", terms.var(0).sum().item())    # nearly identical ✅

# the dependent counter-example
a  = torch.randn(100000)
st = torch.stack([a, a], dim=1)                # two identical columns
print("left :", st.sum(1).var().item())         # ≈ 4
print("right:", st.var(0).sum().item())         # ≈ 2   ❌ not equal

Step 3 — take the root

σ_z = √n · σ_W · σ_x

Checked: σ_x = 0.909, theory √30 × 1 × 0.911 = 4.99, measured 4.958 ✅

Two assumptions

Both assumptions decay once training starts: gradients push W until it is no longer independent of the data, and the mean is no longer 0. So this formula governs step 0, not step 1000 — which is exactly why BatchNorm exists.

The normal distribution, and forecasting saturation

Many small independent factors added together tend towards a normal distribution (the central limit theorem). zⱼ is 30 terms added, so it is necessarily close to normal — whatever the distributions of x and W were.

That property buys one concrete ability: knowing σ_z tells you how many neurons will die, without running the model.

Predicted 54.8%, measured 62.3% in the previous part — the same order. The gap has two sources: |z| > 3 and |h| > 0.99 are not the same threshold, and a single sample fluctuates.

Running it backwards is more interesting: asking for saturation under 5% is asking that 3 be at least 2 σ out, that is σ_z ≲ 1.5. The hand-tuned 0.2 from the previous part gives σ_z ≈ 1.0 — inside that range.

Why diagnosis watches the variance

Three properties nothing else has:

① Predictable — there is a closed form, so the next layer can be computed without running the model. The figure above is exactly that.

② Compounds layer by layer — multiply the per-layer ratios and you have the total effect. 1.5 per layer is 57× after ten (explosion); 0.5 per layer is 0.001 (vanishing). This is the entire basis for the previous part’s “per-layer σ ratio ≈ 1”.

③ Works backwards too — gradients propagating backwards obey the same rules, because ∂L/∂x = ∂L/∂z · Wᵀ is also a matrix product. One metric monitors both directions.

Variance is the only quantity that connects “what to initialise with” to “how deep this can stack”.

Looking back at the previous part’s metrics, they are all faces of the same thing:

Two practical details

n or n−1

t.std()                  # correction=1 by default, divides by n−1
t.std(correction=0)      # divides by n

With thousands of elements the difference is negligible — not worth agonising over.

What “the variance of a matrix” means

W1.shape          # (30, 200)
W1.std()          # all 6000 elements as one pile of numbers, ≈ 1.0
W1[:, 0].std()    # column 0, also ≈ 1.0
W1[0, :].std()    # row 0, also ≈ 1.0

Because randn generates each element independently, it slices the same way in every direction. “The standard deviation of a matrix” is a lazy phrase; what it really means is “the standard deviation of the distribution each element is drawn from”.

Check it yourself

import torch
a, b = torch.randn(1_000_000), torch.randn(1_000_000)

# the three rules
print((3*a).var(), '   theory 9')
print((a+b).var(), '   theory 2')
print((a*b).var(), '   theory 1')
print((a+a).var(), '   ← not independent: 4, not 2')

# √n growth
for n in [1, 4, 30, 200]:
    s = torch.randn(200000, n).sum(1)
    print(f'n={n:<4} std={s.std():.3f}  √n={n**0.5:.3f}')

# the network setting
x = torch.randn(5000, 30) * 0.911
W = torch.randn(30, 200)
print(f'σ_z measured {(x@W).std():.3f}   theory {(30**0.5)*0.911:.3f}')

# 68/95/99.7
u = torch.randn(1_000_000)
for k in [1, 2, 3]:
    print(f'within ±{k}σ: {100*(u.abs()<k).float().mean():.2f}%')

Nothing new was introduced here. This part only pays back the formula the previous one borrowed: the √n comes from “variances add, standard deviations don’t”, and variances add because the cross term vanishes under independence. Followed all the way down, the whole thing rests on one condition, E[AB] = E[A]E[B].

And precisely because it rests on that condition, it starts failing the moment training begins — gradients push W until it is no longer independent of the data. The next section is BatchNorm, whose idea is rather than getting it right, force it: subtract the mean and divide by the standard deviation at every layer’s output, regardless of what happened upstream.