The maths of deep learning

import numpy as np

a = np.array([1., 2., 3.])
b = np.array([4., 5., 6.])

print(a @ b)                       # 32.0 = 1*4 + 2*5 + 3*6
print(np.sum(a * b))               # the same thing, spelled out

# Geometric reading: |a||b|cos(theta)
cos = (a @ b) / (np.linalg.norm(a) * np.linalg.norm(b))
print(np.degrees(np.arccos(cos)))  # the angle between them

print(np.array([1., 0.]) @ np.array([0., 1.]))   # 0.0 — orthogonal

# A linear weighted combination IS a dot product, and it is what every layer
# of a network computes: weights dotted with inputs, plus a bias.
weights = np.array([0.5, 0.3, 0.2])
inputs = np.array([10., 20., 30.])
print(weights @ inputs)            # 17.0

# Matrix multiplication applies many dot products at once: element (i, j) of
# A @ B is row i of A dotted with column j of B.
A = np.random.randn(3, 4)
B = np.random.randn(4, 5)
print((A @ B).shape)               # (3, 5) — inner dimensions must match

# The attention score matrix, in one line:
Q, K = np.random.randn(6, 64), np.random.randn(6, 64)
scores = Q @ K.T / np.sqrt(64)     # (6, 6): every position against every other

# A @ B is NOT B @ A. Order matters, and reversing it is a common bug that
# sometimes still runs when the shapes happen to be compatible.

import numpy as np

print(np.log(np.e), np.log(1), np.log2(1024))    # 1.0  0.0  10.0
print(np.log(0.001 * 0.002), np.log(0.001) + np.log(0.002))   # equal

# Why log-probabilities: the product underflows, the sum does not.
probs = np.full(2000, 0.5)
print(np.prod(probs))              # 0.0 — underflowed
print(np.sum(np.log(probs)))       # -1386.29 — fine

def softmax(z):
    z = np.asarray(z, dtype=np.float64)
    z = z - z.max()                # shift for stability; result unchanged
    e = np.exp(z)
    return e / e.sum()

logits = np.array([2.0, 1.0, 0.1])
p = softmax(logits)
print(p, p.sum())                  # sums to 1

# Softmax is shift-invariant but NOT scale-invariant: adding a constant to
# every logit changes nothing, which is why subtracting the max is free.
print(np.allclose(softmax(logits), softmax(logits + 100)))     # True

# max vs argmax — the value vs the position. argmax is what turns a
# distribution into a prediction.
print(logits.max(), logits.argmax())        # 2.0  0
print(np.argsort(logits)[::-1][:2])         # top-2 indices
print(np.argmax([[1, 5], [7, 2]], axis=1))  # per row: [1 0]

import numpy as np

def entropy(p):
    p = np.asarray(p, dtype=np.float64)
    p = p[p > 0]                          # 0 log 0 is defined as 0
    return float(-(p * np.log(p)).sum())

print(entropy([1.0, 0.0, 0.0]))           # 0.0    — certain
print(entropy([1/3, 1/3, 1/3]))           # 1.0986 — maximal for 3 outcomes
print(np.log(3))                          # same
print(np.log(50257))                      # 10.82  — an untrained LM

def cross_entropy(p, q):
    p, q = np.asarray(p), np.clip(np.asarray(q), 1e-12, 1)
    return float(-(p * np.log(q)).sum())

truth = np.array([0., 1., 0.])            # one-hot: the true token is index 1
print(cross_entropy(truth, [0.1, 0.8, 0.1]))   # 0.223 — good prediction
print(cross_entropy(truth, [0.4, 0.2, 0.4]))   # 1.609 — poor
print(cross_entropy(truth, [0.5, 0.001, 0.499]))  # 6.91 — confidently wrong

# With one-hot truth it reduces to -log q[correct]:
print(-np.log(0.8))

# Perplexity is its exponential — the effective number of choices.
losses = [0.223, 1.609, 6.908]
for L in losses:
    print(f"loss {L:.3f}  perplexity {np.exp(L):8.2f}")

import numpy as np

x = np.array([2., 4., 4., 4., 5., 5., 7., 9.])

print(x.mean(), np.median(x))
print(x.var(), x.var(ddof=1))        # population vs sample (n-1)
print(x.std(), x.std(ddof=1))

# Standard error: the precision of the MEAN, not the spread of the data.
sem = x.std(ddof=1) / np.sqrt(len(x))
print(f"mean {x.mean():.2f} +/- {1.96 * sem:.2f}  (95% CI)")

# Sampling variability, demonstrated: the same population, different samples.
rng = np.random.default_rng(0)
population = rng.normal(100, 15, 100_000)
for n in [10, 50, 200, 1000]:
    means = [rng.choice(population, n).mean() for _ in range(500)]
    print(f"n={n:5}  sd of means {np.std(means):5.2f}  "
          f"predicted {15 / np.sqrt(n):5.2f}")
# The standard error formula predicts it exactly. This is why 'I found this
# neuron on 20 stimuli' is not a result.

# Bootstrap: a confidence interval without distributional assumptions.
def bootstrap_ci(data, statistic=np.mean, n_boot=10_000, alpha=0.05):
    rng = np.random.default_rng(0)
    boots = [statistic(rng.choice(data, len(data), replace=True))
             for _ in range(n_boot)]
    return np.percentile(boots, [100 * alpha / 2, 100 * (1 - alpha / 2)])

print(bootstrap_ci(x))

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
group_a = rng.normal(0.5, 1.0, 40)
group_b = rng.normal(0.0, 1.0, 40)

# Welch's t-test — does not assume equal variances. Prefer it as the default.
t, p = stats.ttest_ind(group_a, group_b, equal_var=False)
print(f"t = {t:.3f}, p = {p:.4f}")

# Effect size: what the difference actually is, independent of n.
pooled = np.sqrt((group_a.var(ddof=1) + group_b.var(ddof=1)) / 2)
d = (group_a.mean() - group_b.mean()) / pooled
print(f"Cohen's d = {d:.3f}")     # ~0.2 small, ~0.5 medium, ~0.8 large

# Non-parametric alternative when the data are clearly not normal — which
# post-GELU activations usually are not.
print(stats.mannwhitneyu(group_a, group_b))

# Multiple comparisons: the correction that makes neuron scans valid.
from statsmodels.stats.multitest import multipletests
p_values = rng.uniform(0, 1, 5000)          # pure noise
print(f"uncorrected 'significant': {(p_values < 0.05).sum()}")   # ~250
reject, p_adj, *_ = multipletests(p_values, alpha=0.05, method="fdr_bh")
print(f"after FDR correction: {reject.sum()}")                   # ~0

import numpy as np, torch

# f(x) = x^3 - 3x, so f'(x) = 3x^2 - 3, zero at x = +/-1.
f = lambda x: x ** 3 - 3 * x
fp = lambda x: 3 * x ** 2 - 3

# Numerical check by central difference — always verify a hand-derived gradient.
def numerical(f, x, h=1e-5):
    return (f(x + h) - f(x - h)) / (2 * h)

for x in [-2.0, -1.0, 0.0, 1.0, 2.0]:
    print(f"x={x:+.1f}  analytic {fp(x):+7.3f}  numerical {numerical(f, x):+7.3f}")

# The chain rule, by hand: y = (3x + 1)^2
#   du/dx = 3,  dy/du = 2u  =>  dy/dx = 6(3x + 1)
chain = lambda x: 6 * (3 * x + 1)

# And by autograd, which is the same computation done for you.
x = torch.tensor(2.0, requires_grad=True)
y = (3 * x + 1) ** 2
y.backward()
print(x.grad.item(), chain(2.0))          # 42.0  42.0

# Product rule: d/dx [x^2 * sin(x)] = 2x sin(x) + x^2 cos(x)
x = torch.tensor(1.0, requires_grad=True)
(x ** 2 * torch.sin(x)).backward()
print(x.grad.item(),
      2 * np.sin(1) + np.cos(1))

# Vanishing gradients, in one line: a product of many small derivatives.
print(np.prod([0.25] * 50))        # ~8e-31 — this is why deep sigmoid nets
                                   # would not train, and why residual
                                   # connections were the fix.