Math & Statistics for ML

The big picture: why ML needs math

Machine learning looks like magic from the outside, but inside it is a small number of mathematical ideas used over and over. If you understand these ideas well, every new model or paper becomes "the same pieces in a new arrangement" rather than something to memorize. This page is ordered so that a newcomer can read it top to bottom: we start with describing data, move to uncertainty, then to the machinery of computation and learning, and finish with the math inside large language models.

Where each branch of math shows up

One neural network, all the math at once

Here is a tiny example that touches almost every section of this page. A student studied 7 hours, slept 8 hours and had a previous score of 6 out of 10. We put these three numbers in a row vector x = [7, 8, 6]. A layer with two hidden neurons has a 3×2 weight matrix W. Multiplying x by W produces two numbers, each a dot product of the input with one column of weights.

We then squash each value with the sigmoid σ(z) = 1/(1 + e^−z), giving [0.9998, 0.9608], and a second layer with weights [0.7, 0.3] combines them into 0.9998(0.7) + 0.9608(0.3) ≈ 0.988, which we read as a 98.8% probability of passing. That single forward pass used vectors, matrix multiplication and a non-linear function. Training would then compare 0.988 to the true label with a loss (information theory and probability), compute the gradient of that loss with respect to every weight (calculus), and update the weights (optimization). Stack 1,000 students as rows of a matrix and the same single multiplication processes all of them at once, which is exactly why GPUs are so good at ML.

  data (matrix)        linear algebra          calculus              optimization
 [n x features] ──▶  Z = XW + b  ──▶  A = f(Z) ──▶ loss L ──▶ dL/dW via chain rule ──▶ W = W - lr * dL/dW
                          ▲                          │
                          └──────────── repeat for many mini-batches ◀──────────┘
 probability: what outputs mean      statistics: is the final model really better?

Descriptive statistics: summarizing data

Before training anything, you need to understand your data. Descriptive statistics compress a column of thousands of numbers into a handful of summaries: where the center is, how spread out the values are, what shape the distribution has, and which points are unusual. These summaries drive decisions such as how to scale features, whether to use the mean or median to fill missing values, and which rows might be data-entry errors.

Measures of central tendency: mean, median, mode

We will use one running example: x = [22, 18, 14, 10, 15, 20, 25, 12] (n = 8).

• Median: sort the data, take the middle. Sorted: [10, 12, 14, 15, 18, 20, 22, 25]. With an even count, average the two middle values: (15 + 18)/2 = 16.5. For odd n it is simply the middle element.
• Mode: the most frequent value. A dataset can be unimodal, multimodal (several values tie) or have no mode (all values unique, as here). Mode is the only center measure that works for categorical data ("most common city").
• Weighted mean: Σ w_ix_i / Σ w_i. Used for class-weighted metrics and averaging over batches of unequal size.
• Geometric mean: (Π x_i)^1/n, the right average for growth rates and ratios. BLEU is a geometric mean of n-gram precisions.
• Harmonic mean: n / Σ(1/x_i), the right average for rates. F1 is the harmonic mean of precision and recall.

Robustness. Add one outlier, 250, to the example. The mean jumps from 17 to (136 + 250)/9 ≈ 42.9, far away from every normal value, while the median only moves from 16.5 to 18. This is why the median is preferred for skewed quantities such as income, house prices, request latency and time-on-page, and why "median latency" and "p99 latency" are reported instead of the mean.

Measures of position: quantiles, percentiles, quartiles

Quantiles cut sorted data into equal-sized chunks. The median is the 0.5 quantile (half below, half above). Quartiles cut into four parts (Q1 = 25th percentile, Q2 = median, Q3 = 75th percentile), deciles into ten, percentiles into a hundred, and tertiles into three. They tell you the relative standing of a value: "this response time is at the 99th percentile" means 99% of requests were faster.

When the cut point falls between two data points, most libraries (NumPy's default "linear" method) interpolate. With 0-based indexing, the position of quantile q in sorted data of length n is p = q(n − 1); the value is x_⌊p⌋ + (p − ⌊p⌋)(x_⌊p⌋+1 − x_⌊p⌋).

1. Sort [10, 12, 14, 15, 18, 20, 22, 25], n = 8.
2. Q1 position 0.25 × 7 = 1.75, between index 1 (12) and index 2 (14). Value = 12 + 0.75 × (14 − 12) = 13.5.
3. Q3 position 0.75 × 7 = 5.25, between 20 and 22. Value = 20 + 0.25 × 2 = 20.5.
4. IQR = Q3 − Q1 = 20.5 − 13.5 = 7. This is the spread of the middle 50% of the data.
5. Tertile example For [15, 18, 67, 100, 21, 50], sorted [15, 18, 21, 50, 67, 100], the 1/3 quantile sits at position 5/3 ≈ 1.67, giving 18 + 0.67 × 3 = 20; the 2/3 quantile sits at 3.33, giving 50 + 0.33 × 17 ≈ 55.7.

import numpy as np
x = np.array([22, 18, 14, 10, 15, 20, 25, 12])
q1, med, q3 = np.percentile(x, [25, 50, 75])   # 13.5, 16.5, 20.5
iqr = q3 - q1                                   # 7.0
np.quantile([15, 18, 67, 100, 21, 50], [1/3, 2/3])  # [20.0, 55.67]

Box plots and the IQR outlier rule

A box plot draws the five-number summary: minimum, Q1, median, Q3, maximum. The box spans Q1 to Q3, a line marks the median, and "whiskers" extend to the most extreme points within 1.5 × IQR of the box. Points beyond the whiskers are drawn individually as potential outliers.

          whisker      box (middle 50%)      whisker
   10 |-----------[ 13.5 ===|=== 20.5 ]------------| 25          o 250  (outlier)
                  Q1      median 16.5    Q3
     fences: 3 ............................................ 31

Measures of variability (spread)

Two datasets can share a mean and be completely different: [49, 50, 51] and [0, 50, 100] both average 50. Spread measures capture that difference.

The squared deviations for the example are 25, 1, 9, 49, 4, 9, 64, 25, which sum to 186. Squaring does two things: it stops positive and negative deviations from cancelling out (they always sum to exactly zero around the mean), and it makes a deviation of 8 count four times as much as a deviation of 4.

Why divide by n − 1 (Bessel's correction)?

The sample mean x̄ is computed from the same data, and it is by construction the point that minimizes the sum of squared deviations. So deviations from x̄ are systematically a little smaller than deviations from the true, unknown μ. Dividing by n − 1 instead of n exactly corrects that bias on average. Another way to say it: once you know x̄ and n − 1 of the deviations, the last deviation is determined, so only n − 1 are "free" (degrees of freedom). For n in the thousands the difference is negligible; for n = 5 it is 25%.

x = np.array([22, 18, 14, 10, 15, 20, 25, 12])
x.var()          # 23.25   population (ddof=0) - NumPy default
x.var(ddof=1)    # 26.57   sample - pandas .var() default!
x.std(), x.std(ddof=1)   # 4.82, 5.15

Shape: skewness and kurtosis

Skewness measures asymmetry: E[((X − μ)/σ)³]. Positive (right) skew means a long tail to the right, with mean > median (income, latency, file sizes). Negative (left) skew means a long left tail, mean < median (exam scores on an easy test, age at retirement). Zero skew means symmetric, like the normal distribution. For example, [1, 2, 2, 3, 3, 3, 4, 10] has skewness ≈ 1.8 because of the single 10.

Kurtosis measures tail heaviness: E[((X − μ)/σ)⁴]. A normal distribution has kurtosis 3, so "excess kurtosis" = kurtosis − 3 is usually reported. High excess kurtosis (leptokurtic, e.g. financial returns, Student-t) means extreme values happen far more often than a Gaussian would predict, which matters for risk models and for choosing robust losses.

  right (positive) skew          symmetric              left (negative) skew
      |\                            /\                              /|
      | \__                        /  \                          __/ |
      |    \_____                 /    \                   _____/    |
   mode < median < mean      mean = median = mode      mean < median < mode

Z-scores and standardization

A z-score expresses a value as "how many standard deviations from the mean".

Z-scores let you compare apples and oranges (a height z-score vs. a weight z-score), and the z-score outlier rule flags |z| > 3 as unusual. Standardizing an entire feature (subtract the mean, divide by the standard deviation) gives it mean 0 and standard deviation 1, which is the default preprocessing for linear models, SVMs, k-NN, PCA and neural networks. Min-max normalization, x' = (x − min)/(max − min), instead squeezes values into [0, 1] and is common for pixel intensities.

Outliers: detecting and handling them

• Detection: IQR fences (robust, no normality assumption), z-score > 3 (assumes roughly normal data, and the outlier itself inflates σ), modified z-score using the median and median absolute deviation (robust), Mahalanobis distance for multivariate data, and model-based methods such as Isolation Forest.
• Handling: first ask whether it is an error (a negative age, a 999 sentinel) or a genuine rare event (a real fraud). Fix or drop errors. For genuine extremes, consider capping/winsorizing at the 1st/99th percentile, log-transforming, using robust statistics (median, IQR), robust losses (MAE, Huber), or tree models, which are largely insensitive to monotonic scaling.
• Never drop outliers blindly in fraud or anomaly detection; there the outliers are the signal.

Probability foundations and Bayes' theorem

Probability is the mathematics of uncertainty. Almost every ML output is a probability statement: "92% spam", "the next token is 'Paris' with probability 0.41". To interpret those numbers, combine them, and avoid classic traps, you need a handful of rules.

The vocabulary

• Sample space S: all possible outcomes (a die: {1, 2, 3, 4, 5, 6}).
• Event A: a subset of outcomes you care about ("even" = {2, 4, 6}).
• Probability P(A): a number in [0, 1]. For equally likely outcomes, P(A) = favorable / total = 3/6.
• Complement A^c: "not A". P(A^c) = 1 − P(A).
• Axioms (Kolmogorov): P(A) ≥ 0; P(S) = 1; for mutually exclusive events, P(A ∪ B) = P(A) + P(B). Everything else is derived from these three.

The core rules

Worked example (complement trick). An LLM agent calls a tool that fails 5% of the time, independently. Over 10 calls, P(at least one failure) = 1 − 0.95¹⁰ = 1 − 0.599 = 0.401. A small per-step error rate compounds quickly in multi-step pipelines, which is why agent reliability degrades with chain length.

Joint, marginal and conditional probability

These three are easiest to understand from a table. Over 100 observed days:

• Joint P(Rainy, Heavy) = 30/100 = 0.30: both at once; read from a cell.
• Marginal P(Rainy) = 40/100 = 0.40: one variable alone; read from a total. "Marginalizing out" traffic means summing over it: P(Rainy) = P(Rainy, Heavy) + P(Rainy, Light).
• Conditional P(Heavy | Rainy) = 30/40 = 0.75: restrict to the rainy row.
• Flipped conditional P(Rainy | Heavy) = 30/45 ≈ 0.67. Not the same! P(A | B) ≠ P(B | A) in general, and confusing the two is the "prosecutor's fallacy".

The joint always factorizes both ways: P(A, B) = P(A | B) P(B) = P(B | A) P(A). Rearranging that identity gives Bayes' theorem. In ML terms, a discriminative classifier (logistic regression, most neural nets) models the conditional P(y | x); a generative model (Naive Bayes, GMMs, VAEs, diffusion, LLMs over text) models a joint or the data distribution P(x) and can sample new data.

Independence and conditional independence

A and B are independent if P(A | B) = P(A). In the table, P(Heavy) = 0.45 but P(Heavy | Rainy) = 0.75, so rain and traffic are dependent. Conditional independence means A and B are independent once you know C: P(A, B | C) = P(A | C) P(B | C). Ice-cream sales and drownings are correlated, but conditionally independent given temperature. Naive Bayes assumes all features are conditionally independent given the class, which is almost never true but works surprisingly well for text.

Bayes' theorem

Worked example 1: the medical test (base-rate fallacy)

A disease affects 1% of people. A test catches 95% of sick people (sensitivity) and wrongly flags 5% of healthy people (false-positive rate). You test positive. What is the chance you are sick?

1. Prior P(sick) = 0.01, P(healthy) = 0.99.
2. Likelihoods P(+ | sick) = 0.95, P(+ | healthy) = 0.05.
3. Evidence P(+) = 0.95 × 0.01 + 0.05 × 0.99 = 0.0095 + 0.0495 = 0.059.
4. Posterior P(sick | +) = 0.0095 / 0.059 ≈ 0.161, only about 16%.
5. Second independent positive test Use 0.161 as the new prior: 0.95 × 0.161 / (0.95 × 0.161 + 0.05 × 0.839) ≈ 0.785. Evidence accumulates.

The intuition with counts: out of 10,000 people, 100 are sick and 95 of them test positive; of the 9,900 healthy people, 495 test positive. So 95 of 590 positives are truly sick, about 16%. The rare condition means false alarms from the huge healthy group swamp the true positives. This is exactly why fraud and anomaly detectors with seemingly high accuracy still produce mostly false alerts, and why precision collapses on imbalanced data.

Worked example 2: a spam filter

30% of emails are spam. 80% of spam contains the word "free"; 25% of all emails contain "free". Then P(spam | "free") = 0.8 × 0.3 / 0.25 = 0.96. Naive Bayes extends this to many words by multiplying the per-word likelihoods (in log space, adding log-likelihoods) under the conditional independence assumption.

def bayes(prior, sens, fpr):
    evidence = sens * prior + fpr * (1 - prior)
    return sens * prior / evidence

p1 = bayes(0.01, 0.95, 0.05)   # 0.161
p2 = bayes(p1, 0.95, 0.05)     # 0.785  (second positive test)

Counting: permutations and combinations

Many probability questions reduce to counting equally likely outcomes. Permutations count ordered arrangements: P(n, k) = n!/(n − k)!. Combinations count unordered selections: C(n, k) = n!/(k!(n − k)!). For example, the number of ways to pick 3 features out of 10 is C(10, 3) = 120; the number of 7-heads outcomes in 10 coin flips is also C(10, 7) = 120, so P(exactly 7 heads) = 120/1024 ≈ 0.117.

Random variables, expectation, variance and covariance

A random variable turns random outcomes into numbers so we can do arithmetic on them. Once we have numbers, we can ask for the average outcome (expectation), how much outcomes wobble (variance), and whether two quantities move together (covariance and correlation). These are the workhorses behind every loss function, every initialization scheme and PCA.

Discrete vs. continuous

The cumulative distribution function (CDF) F(x) = P(X ≤ x) works for both. It rises from 0 to 1, and its inverse (the quantile function) answers "what value is exceeded only 10% of the time?". Integration is how we get from a density to a probability, which is the main use of integral calculus in ML: normalizing distributions, computing expectations, and areas such as ROC-AUC.

Expectation

• Linearity: E[aX + bY + c] = aE[X] + bE[Y] + c, always, even if X and Y are dependent. This makes many problems trivial: the expected number of fixed points in a random shuffle of n items is n × (1/n) = 1.
• Function of a variable: E[g(X)] = Σ g(x)p(x). In general E[g(X)] ≠ g(E[X]); Jensen's inequality says for convex g, E[g(X)] ≥ g(E[X]). This is why the ELBO is a lower bound and why log of an average ≥ average of logs.
• Products: E[XY] = E[X]E[Y] only if X and Y are uncorrelated (independence is sufficient).
• In ML: MSE = E[(ŷ − y)²], cross-entropy = E[−log p(y | x)], an RL return is an expected reward. Training minimizes an empirical average that approximates this expectation.

Variance and standard deviation

• Var(aX + b) = a² Var(X): shifting does not change spread; scaling by a scales variance by a².
• Var(X + Y) = Var(X) + Var(Y) + 2 Cov(X, Y). For independent variables the covariance term vanishes.
• Average of n independent copies: Var(x̄) = σ²/n. This one identity explains why larger batches give less noisy gradients, why ensembles help, and why standard errors shrink with √n.
• Sum of n independent unit-variance terms has variance n. This is why attention scores are divided by √d_k and why Xavier/He initialization scales weights by 1/√fan-in.

Covariance and correlation

Worked example. Hours studied X = [1, 2, 3, 4, 5], score Y = [52, 55, 61, 64, 68]. Means: 3 and 60. Deviations of X: [−2, −1, 0, 1, 2]; of Y: [−8, −5, 1, 4, 8]. Products: 16, 5, 0, 4, 16, sum 41. Sample covariance = 41/4 = 10.25. Sample standard deviations: √(10/4) ≈ 1.58 and √(170/4) ≈ 6.52. Correlation = 10.25/(1.58 × 6.52) ≈ 0.99, a near-perfect positive linear relationship.

• ρ = +1 or −1: perfect linear relationship; ρ = 0: no linear relationship.
• Zero correlation does not imply independence. If X is symmetric around 0 and Y = X², then Cov(X, Y) = E[X³] = 0, yet Y is completely determined by X.
• Correlation is not causation. A hidden confounder can drive both variables.
• Spearman correlation is Pearson on ranks; it captures any monotonic relationship and is robust to outliers. Kendall's tau counts concordant vs. discordant pairs.
• Pearson correlation between two mean-centered vectors equals their cosine similarity.

The covariance matrix

For d features, the covariance matrix Σ is d × d with variances on the diagonal and covariances off-diagonal. It is always symmetric and positive semi-definite (all eigenvalues ≥ 0). With a centered data matrix X_c (n × d), Σ = X_c^TX_c/(n − 1). Its eigenvectors are the principal components (see PCA below), its trace is the total variance, and its inverse appears in the multivariate normal density and the Mahalanobis distance.

import numpy as np
X = np.array([1, 2, 3, 4, 5]); Y = np.array([52, 55, 61, 64, 68])
np.cov(X, Y)          # [[2.5, 10.25], [10.25, 42.5]]  (ddof=1 by default!)
np.corrcoef(X, Y)[0, 1]   # 0.994
data = np.random.randn(1000, 3)
Sigma = np.cov(data, rowvar=False)   # 3x3; rows are samples

Probability distributions you must know

A distribution is the "shape" of randomness: which values are possible and how likely each one is. A small family of distributions covers most of ML, and each one comes with a story about how data is generated. Knowing the story tells you which distribution (and therefore which loss function) fits a problem.

Discrete distributions

Worked examples. Binomial: 10 fair coin flips, P(exactly 7 heads) = C(10,7)(0.5)¹⁰ = 120/1024 ≈ 0.117. Poisson: a server averages λ = 3 errors per hour; P(0 errors) = e⁻³ ≈ 0.050 and P(2 errors) = 9e⁻³/2 ≈ 0.224. The Poisson is the limit of a binomial with large n and small p (np = λ); a Poisson whose variance is much larger than its mean signals overdispersion, where a negative binomial fits better.

Continuous distributions

Exponential example. Customers arrive at λ = 2 per hour. Mean wait between arrivals is 1/2 hour, and P(wait > 1 hour) = e⁻² ≈ 0.135. The exponential is memoryless: having already waited 30 minutes does not change the distribution of the remaining wait.

Beta as a conjugate prior. If your prior on a click-through rate is Beta(α, β) and you observe k clicks in n views, the posterior is Beta(α + k, β + n − k). You can think of α and β as "pseudo-counts" of prior successes and failures. The Dirichlet does the same for K categories (it is the conjugate prior of the categorical/multinomial), which is why add-one (Laplace) smoothing in Naive Bayes corresponds to a Dirichlet(1, ..., 1) prior.

The normal distribution in depth

The normal (Gaussian) distribution is the star of statistics because of the central limit theorem (sums and averages of many independent effects become normal) and because it is mathematically convenient: sums of independent normals are normal, it is fully described by two numbers, and it has maximum entropy for a given mean and variance.

• 68-95-99.7 rule: about 68% of values lie within ±1σ, 95% within ±2σ (more precisely ±1.96σ), 99.7% within ±3σ.
• Standard normal Z ~ N(0, 1); any normal converts with z = (x − μ)/σ.
• Worked example: exam scores ~ N(70, 10²). P(60 ≤ X ≤ 80) = P(−1 ≤ Z ≤ 1) ≈ 0.683. The 90th percentile is 70 + 1.28 × 10 = 82.8.
• Anomaly example: legitimate transaction amounts ~ N(50, 20²). A $150 charge has z = 5, and P(Z ≥ 5) ≈ 3 × 10⁻⁷, a strong anomaly signal (if the normal assumption holds, which for money it often does not; heavy tails are common).
• Multivariate normal N(μ, Σ): density ∝ exp(−½(x − μ)^TΣ⁻¹(x − μ)). Contours are ellipses aligned with the eigenvectors of Σ. Used in Gaussian mixture models, Gaussian processes, VAE latents and diffusion.

from scipy import stats
stats.norm.cdf(1) - stats.norm.cdf(-1)     # 0.683
stats.norm.ppf(0.90, loc=70, scale=10)     # 82.8
stats.binom.pmf(7, n=10, p=0.5)            # 0.117
stats.poisson.pmf(2, mu=3)                 # 0.224
stats.expon.sf(1, scale=1/2)               # 0.135  P(T > 1)
rng = np.random.default_rng(0)
rng.normal(0, 1, size=(3, 4)); rng.beta(2, 5, size=10); rng.dirichlet([1, 1, 1])

Sampling, the law of large numbers and the central limit theorem

We almost never see a whole population; we see a sample and use it to estimate population quantities. Two theorems explain why this works. The law of large numbers (LLN) says averages settle down to the truth as the sample grows. The central limit theorem (CLT) says the error of that average is approximately normal, which lets us build confidence intervals and hypothesis tests for almost any kind of data.

Populations, samples, estimators

• Parameter: a fixed but unknown population value (μ, σ, p). Statistic/estimator: a quantity computed from a sample (x̄, s, p̂) that estimates it.
• Bias of an estimator = E[estimate] − true value. x̄ is unbiased for μ; the 1/n variance is biased; the 1/(n−1) variance is unbiased.
• Standard error (SE) = standard deviation of the estimator across repeated samples. For the mean, SE = σ/√n (estimated by s/√n). For a proportion, SE = √(p(1−p)/n).
• Sampling schemes: simple random, stratified (sample within each class to keep proportions, as in train_test_split(stratify=y)), cluster, systematic. Bootstrap resampling (sampling with replacement from your sample) estimates the variability of any statistic without formulas.
• Sampling bias: survivorship bias, selection bias, non-response bias. No amount of data fixes a biased sampling process.

Law of large numbers

LLN justifies empirical risk minimization (the training-set average loss approximates the expected loss), Monte Carlo estimation (average many random samples to estimate an integral), and trusting evaluation metrics computed on large test sets.

Central limit theorem

Worked example. A single die roll is uniform on 1..6 (mean 3.5, σ ≈ 1.71), nothing like a bell. Average 100 rolls, and repeat 10,000 times: the averages form a bell curve centered at 3.5 with standard deviation 1.71/√100 ≈ 0.171. So about 95% of those averages fall in 3.5 ± 0.34.

rng = np.random.default_rng(0)
means = rng.integers(1, 7, size=(10_000, 100)).mean(axis=1)
means.mean(), means.std()      # about 3.5 and 0.171 - and a histogram is bell-shaped

Why this matters in ML

• Mini-batch SGD: the batch-mean gradient is unbiased with variance σ²/B. Quadrupling the batch size only halves the gradient noise, which is part of why learning rate is often scaled with batch size.
• Confidence intervals on metrics: accuracy on 1,000 test examples at 90% has SE = √(0.9 × 0.1/1000) ≈ 0.0095, so a 95% CI of about ±1.9 points. A "0.5 point improvement" on that test set is noise.
• A/B tests: differences in conversion rates are approximately normal by the CLT, which gives the z-test below.
• Initialization: a neuron's pre-activation is a sum of many weighted inputs and is approximately normal; its variance grows with fan-in, motivating 1/fan-in scaling.

Inferential statistics: confidence intervals, hypothesis tests and A/B testing

Inferential statistics answers "is this pattern real, or could it be luck?" It is how you decide whether a new model is actually better, whether a product change moved a metric, or whether a feature is related to the target. ML engineers use it every time they run an online experiment or compare two models on a finite test set.

Confidence intervals

Worked example. Average session length from n = 100 users is 50 seconds, with s = 10. SE = 10/√100 = 1. The 95% CI is 50 ± 1.96, i.e. [48.04, 51.96] seconds.

Correct interpretation: if we repeated the sampling procedure many times, 95% of the intervals built this way would contain the true mean. It is not "there is a 95% probability the true mean is in this particular interval" (that is the Bayesian credible interval's interpretation). Wider intervals come from more variability, smaller samples, or higher confidence levels.

Bootstrap CI: resample the data with replacement B = 1,000+ times, recompute the statistic each time, and take the 2.5th and 97.5th percentiles. This works for medians, F1 scores, AUC and anything else without a neat SE formula.

Bootstrap versus cross-validation

Both reuse the sample, but they answer different questions. The bootstrap estimates the sampling distribution of a statistic (a confidence interval or standard error for accuracy, F1, a median, or the gap between two models). Cross-validation estimates the expected performance of a training procedure by repeatedly holding out folds and retraining. Use the bootstrap on a frozen model and a test set when you want error bars; use CV on the training data when you want to compare algorithms or tune hyperparameters. Do not treat a bootstrap of the training loss as a substitute for a held-out fold: every bootstrap draw still overlaps heavily with the original sample, so it is optimistic for generalization. Out-of-bag scores from bagged trees are the exception that works, because each tree is scored only on rows it did not see.

scores = np.array(per_example_correct)     # 0/1 per test example
boot = [rng.choice(scores, size=len(scores), replace=True).mean() for _ in range(2000)]
lo, hi = np.percentile(boot, [2.5, 97.5])   # 95% bootstrap CI for accuracy

Hypothesis testing step by step

1. State hypotheses Null H₀: no effect (μ_A = μ_B). Alternative H₁: an effect (μ_A ≠ μ_B for two-sided, or > for one-sided).
2. Choose α before looking at data, usually 0.05. This is the Type I error rate you accept.
3. Pick a test statistic that measures the effect in units of noise, e.g. t = (difference) / SE.
4. Compute the p-value P(statistic at least this extreme | H₀ true).
5. Decide If p ≤ α, reject H₀. Otherwise, fail to reject. Always report the effect size and a CI, not just the p-value.

What a p-value is and is not

• It is the probability of seeing data at least as extreme as yours if the null hypothesis were true.
• It is not the probability that the null is true, not the probability the result is due to chance, and not a measure of effect size. With huge n, a trivial 0.01% lift can have p < 0.001; with tiny n, a large real effect can have p = 0.3.

Type I and Type II errors, power

Power (usually targeted at 0.8) increases with sample size, effect size and α, and decreases with variance. There is a trade-off: lowering α to avoid false positives raises β unless you collect more data. In classifier terms, Type I errors are false positives (hurting precision) and Type II errors are false negatives (hurting recall).

The common tests and when to use them

t-test example. A new model's average latency over n = 25 requests is 52 ms with s = 5 ms. Is it different from the 50 ms target? t = (52 − 50)/(5/√25) = 2.0 with 24 degrees of freedom, two-sided p ≈ 0.057. At α = 0.05 we fail to reject, though the evidence is borderline; more data would settle it. The t-distribution has heavier tails than the normal to account for estimating σ from a small sample; with n > 30 or so, t and z are nearly identical.

A/B testing worked end to end

Control (A): 1,000 users, 100 conversions (10.0%). Treatment (B): 1,000 users, 120 conversions (12.0%). Is B better?

1. Pooled rate p̂ = (100 + 120)/2000 = 0.11.
2. Standard error √(0.11 × 0.89 × (1/1000 + 1/1000)) ≈ 0.0140.
3. z statistic (0.12 − 0.10)/0.0140 ≈ 1.43.
4. p-value two-sided ≈ 0.153. Not significant at 0.05, even though the lift is a relative 20%.
5. Same answer via chi-square On the 2×2 table [[100, 900], [120, 880]], expected counts are 110 and 890 per row; χ² = 2(10²/110) + 2(10²/890) ≈ 2.04 = z², same p-value.
6. Sample size needed To detect 10% → 12% with α = 0.05 and 80% power: n ≈ (1.96 + 0.84)² × (0.10 × 0.90 + 0.12 × 0.88)/0.02² ≈ 3,840 users per group. The test above was underpowered.

from statsmodels.stats.proportion import proportions_ztest
z, p = proportions_ztest(count=[120, 100], nobs=[1000, 1000])   # z=1.43, p=0.153
from scipy import stats
chi2, p, dof, expected = stats.chi2_contingency([[100, 900], [120, 880]], correction=False)
t, p = stats.ttest_ind(a, b, equal_var=False)   # Welch t-test
t, p = stats.ttest_rel(model_a_scores, model_b_scores)   # paired

A/B testing pitfalls

• Peeking: checking the p-value daily and stopping when it dips below 0.05 inflates the false-positive rate far above 5%. Fix the sample size in advance or use sequential testing methods.
• Multiple comparisons: test 20 metrics at α = 0.05 and you expect one false "win". Use Bonferroni (α/m) or control the false discovery rate (Benjamini-Hochberg).
• Sample ratio mismatch: if a 50/50 split produced 52/48, randomization or logging is broken; a chi-square test detects it.
• Novelty and primacy effects, weekly seasonality (run full weeks), network effects between users, and Simpson's paradox when segments have different mixes.
• Statistical vs. practical significance: always report the effect size with a confidence interval and translate it into business impact.

Linear algebra: vectors and matrices

Linear algebra is the engine room of ML. Every data point is a vector, every dataset and weight layer is a matrix, and nearly all computation in a neural network is matrix multiplication. GPUs exist in ML because they do matrix multiplications extremely fast in parallel. If you can read shapes and reason about matrix products, you can read most model code.

Scalars, vectors, matrices and your data

• Scalar: a single number (a learning rate, a loss value).
• Vector: an ordered list of numbers, a point or direction in n-dimensional space. A 28×28 grayscale image flattened is a vector in 784 dimensions; a word embedding is a vector in 768 or 4096 dimensions.
• Matrix: a rectangular grid, m rows × n columns; element a_ij is row i, column j. A dataset of 3 houses with 4 columns (area, bedrooms, age, price) is a 3×4 matrix; a grayscale image is an H×W matrix and a color image is three stacked matrices (R, G, B).

Vector operations

The dot product is the single most important operation in ML. A neuron computes w·x + b. A linear classifier's decision is the sign of a dot product. Attention scores are dot products of queries and keys. Recommendation scores are dot products of user and item vectors: with Alice's tastes [action 0.9, romance 0.1, sci-fi 0.8] and a movie's profile [5.0, 0.5, 1.0], the predicted affinity is 4.5 + 0.05 + 0.8 = 5.35. Geometrically, a positive dot product means the vectors point roughly the same way, zero means they are orthogonal (perpendicular), negative means they point in opposite directions.

Norms in ML: L2 norms appear in weight decay (penalizing ||w||₂²), gradient clipping (rescale if ||g|| > threshold) and normalizing embeddings before cosine similarity. L1 norms appear in Lasso (penalizing ||w||₁, which drives weights to exactly zero). The general L_p norm is (Σ|a_i|^p)^1/p; the Frobenius norm of a matrix is the L2 norm of all its entries flattened.

Matrix multiplication

Worked example. [[1, 2], [3, 4]] × [[5, 6], [7, 8]]: top-left = 1·5 + 2·7 = 19, top-right = 1·6 + 2·8 = 22, bottom-left = 3·5 + 4·7 = 43, bottom-right = 3·6 + 4·8 = 50. Result [[19, 22], [43, 50]]. Reversing the order gives [[23, 34], [31, 46]], a different answer.

• Not commutative: AB ≠ BA in general (and BA may not even be defined).
• Associative: (AB)C = A(BC). The order of evaluation can change cost dramatically: for a 1000×1000 A, B and a vector x, computing A(Bx) costs 2 million operations while (AB)x costs a billion.
• Distributive: A(B + C) = AB + AC.
• Element-wise (Hadamard) product A ⊙ B multiplies matching entries; it requires equal shapes and is used in gating (LSTM gates, SwiGLU). In NumPy, A * B is element-wise and A @ B is matrix multiplication.
• Three views of Ax: dot products of each row of A with x; a linear combination of A's columns weighted by x's entries; a transformation of the vector x.

In a neural network, a batch of B inputs with d features is a B×d matrix X. A dense layer with h outputs has a d×h weight matrix W and a bias vector of length h: Z = XW + b is B×h. The bias is broadcast across rows. One matrix multiply processes the whole batch, which is why batching is so efficient on GPUs. Matrix factorization in recommender systems uses the same idea: a sparse users×movies rating matrix R is approximated by P (users×k) times Q (k×movies), and every unseen rating is predicted by one dot product.

import numpy as np
A = np.array([[1, 2], [3, 4]]); B = np.array([[5, 6], [7, 8]])
A @ B        # [[19, 22], [43, 50]]
A * B        # [[ 5, 12], [21, 32]]  element-wise
X = np.random.randn(32, 784)     # batch of 32 flattened images
W = np.random.randn(784, 128) * np.sqrt(2 / 784)   # He init
b = np.zeros(128)
Z = X @ W + b                    # (32, 128); b broadcasts over the batch
H = np.maximum(0, Z)             # ReLU

Transpose, identity, inverse, determinant, trace, rank

Worked inverse. A = [[2, 1], [1, 2]]: det = 2·2 − 1·1 = 3, so A⁻¹ = (1/3)[[2, −1], [−1, 2]]. Check: row 1 of A times column 1 of A⁻¹ = (2·2 + 1·(−1))/3 = 1. For [[1, 2], [2, 4]], det = 4 − 4 = 0: the matrix maps the whole plane onto a line, information is lost, and no inverse exists. Such a matrix is called singular.

Linear independence and span: vectors are linearly independent if none can be written as a combination of the others. The span is every point reachable by combining them. A basis is a minimal set that spans the space; in d dimensions any basis has d vectors. Duplicate or perfectly correlated features make a data matrix rank-deficient, which is the linear-algebra face of multicollinearity.

Special matrices: symmetric (A = A^T, e.g. covariance matrices; they have real eigenvalues and orthogonal eigenvectors), orthogonal (Q^TQ = I; columns are orthonormal; pure rotations or reflections that preserve lengths and angles, used in orthogonal initialization and appearing in SVD), diagonal (scale each axis independently), positive definite (x^TAx > 0 for all non-zero x; a bowl-shaped quadratic; Hessians of strictly convex functions), and sparse (mostly zeros; user-item matrices, one-hot features, TF-IDF matrices).

Solving linear systems and least squares

A system Ax = b has a unique solution x = A⁻¹b when A is square and invertible. In ML we usually have more equations (samples) than unknowns (features), so there is no exact solution; we minimize the squared residual instead. That is linear regression.

Geometrically, Xw* is the orthogonal projection of y onto the column space of X: the residual is perpendicular to every feature column. In practice, you never form the inverse. Libraries solve with QR decomposition or SVD (np.linalg.lstsq), which is faster and far more numerically stable. For millions of rows or features, iterative methods (gradient descent) beat any closed form, since forming X^TX costs O(nd²) and inverting it O(d³).

X = np.c_[np.ones(len(x)), x]          # add bias column
w, *_ = np.linalg.lstsq(X, y, rcond=None)   # stable least squares
# avoid: np.linalg.inv(X.T @ X) @ X.T @ y  (squares the condition number)
np.linalg.solve(A, b)                  # for square systems; never inv(A) @ b

Eigenvalues, SVD and PCA

Matrices can look complicated, but most of them have a few "natural axes" along which their behavior is simple. Eigen-decomposition and the singular value decomposition (SVD) find those axes. Principal component analysis (PCA) uses them to compress data by keeping only the directions where the data varies most. The same ideas power recommendation systems, latent semantic analysis, image compression, and the low-rank tricks used to fine-tune LLMs.

Eigenvalues and eigenvectors

Worked example. A = [[2, 1], [1, 2]]. The characteristic equation is (2 − λ)² − 1 = 0, so λ = 3 or λ = 1. For λ = 3: (A − 3I)v = 0 gives v = [1, 1] (check: A[1, 1] = [3, 3] = 3[1, 1]). For λ = 1: v = [1, −1] (A[1, −1] = [1, −1]). The trace 4 equals 3 + 1 and the determinant 3 equals 3 × 1, which is always true: trace = sum of eigenvalues, determinant = product of eigenvalues.

• Eigendecomposition: A = QΛQ⁻¹ (eigenvectors as columns of Q, eigenvalues on the diagonal of Λ). For symmetric A, Q is orthogonal: A = QΛQ^T (the spectral theorem).
• Powers: A^k = QΛ^kQ⁻¹. Repeated multiplication blows up along eigenvalues with |λ| > 1 and dies along |λ| < 1. This is exactly the vanishing/exploding gradient story in RNNs, where the same weight matrix is applied at every time step.
• Definiteness: all eigenvalues > 0 means positive definite (a bowl); mixed signs means indefinite (a saddle). The Hessian's eigenvalues classify critical points in optimization.
• Other uses: PageRank is the leading eigenvector of a link matrix; spectral clustering uses eigenvectors of a graph Laplacian; the condition number of a symmetric positive definite matrix is λ_max/λ_min and controls how fast gradient descent converges.

Singular value decomposition (SVD)

• Rank = number of non-zero singular values.
• Truncated SVD: keep the top k singular values and vectors, A_k = U_kΣ_kV_k^T. By the Eckart-Young theorem, this is the best possible rank-k approximation in both Frobenius and spectral norm.
• Compression example: a 100×100 grayscale image has 10,000 numbers. Keeping k = 10 needs 100×10 + 10 + 10×100 = 2,010 numbers, about 80% smaller, and often preserves most visible structure because natural images are approximately low rank.
• Uses: pseudoinverse and stable least squares; PCA; latent semantic analysis (SVD of a term-document matrix); collaborative filtering; denoising; and understanding LoRA, which assumes the fine-tuning weight update ΔW is approximately low rank and learns it as B (d×r) times A (r×d) with r as small as 8, cutting trainable parameters from d² to 2dr.

PCA step by step

1. Center (and usually standardize) Subtract each feature's mean. Standardize too if features have different units, otherwise the feature with the largest scale dominates.
2. Covariance matrix Σ = X_c^TX_c/(n − 1), a d×d symmetric matrix.
3. Eigen-decompose The eigenvectors of Σ are the principal components (orthogonal directions); each eigenvalue is the variance of the data along its direction.
4. Sort and choose k Rank components by eigenvalue. Explained variance ratio of component i = λ_i/Σλ. Pick k to reach, say, 95% cumulative variance, or use the "elbow" of a scree plot.
5. Project Z = X_cV_k gives n×k coordinates. Reconstruct approximately with Z V_k^T + mean.

Worked example. Students scored in math, science and English give eigenvalues 5.2, 3.1 and 0.4 (total 8.7). Explained variance: 60%, 36% and 4%. Keeping the first two components retains 96% of the variance while dropping a third of the features. For the 2D covariance matrix [[2, 1], [1, 2]], PC1 is the direction [1, 1]/√2 with variance 3, explaining 3/(3 + 1) = 75% of the total.

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
Xs = StandardScaler().fit_transform(X)
pca = PCA(n_components=0.95).fit(Xs)          # keep 95% of variance
Z = pca.transform(Xs)
pca.explained_variance_ratio_

# PCA from scratch via SVD (numerically preferred over eigh of the covariance)
Xc = X - X.mean(axis=0)
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
explained_var = S**2 / (len(X) - 1)
Z = Xc @ Vt[:k].T

t-SNE and UMAP are non-linear alternatives for visualization. t-SNE matches Gaussian neighborhoods in high dimensions to Student-t neighborhoods in 2D; it produces clean clusters, but cluster sizes and distances between clusters are not meaningful. UMAP preserves more global structure and scales better. Neither should be used as input features for a downstream model without care, since the embedding is not stable or invertible.

Tensors, shapes and broadcasting

In deep learning code, "tensor" simply means a multi-dimensional array. Most bugs in model code are shape bugs, so reading and reasoning about shapes is a core practical skill. Broadcasting is the set of rules that lets you combine arrays of different shapes without writing loops, and it is the source of both elegant code and silent errors.

Common tensor shapes

Broadcasting rules

Compare shapes from the rightmost dimension leftward. Two dimensions are compatible if they are equal, or one of them is 1, or one is missing (treated as 1). The size-1 dimension is virtually stretched (no memory copy) to match.

• (32, 784) + (784,) → (32, 784): the bias is added to each of 32 rows.
• (B, 1, H, W) * (B, C, H, W) → (B, C, H, W): one mask applied to all channels.
• (3, 1) + (1, 4) → (3, 4): an outer-sum table, handy for pairwise distances.
• (32, 784) + (50,) → error: 784 vs. 50.

x = np.random.randn(5, 3)
mu, sd = x.mean(axis=0), x.std(axis=0)     # shapes (3,)
z = (x - mu) / sd                          # (5,3) - (3,) broadcasts: column-wise standardize
# pairwise squared distances between rows of A (n,d) and B (m,d) without loops
D2 = (A**2).sum(1)[:, None] + (B**2).sum(1)[None, :] - 2 * A @ B.T   # (n, m)

# the classic silent bug
y_true = np.array([1., 2., 3.])            # (3,)
y_pred = np.array([[1.], [2.], [3.]])      # (3,1)
(y_pred - y_true).shape                    # (3,3) !  MSE is now wrong, no error raised

Reshaping operations

• reshape / view: same data, new shape; the total size must stay equal. (B, 28, 28) → (B, 784). Use -1 to infer one dimension.
• transpose / permute: reorder axes. (B, H, W, C) → (B, C, H, W). In attention, (B, T, h, d_k) → (B, h, T, d_k) so each head is a separate batch of matrix multiplies.
• squeeze / unsqueeze (NumPy expand_dims, or indexing with None): remove or add size-1 axes.
• concatenate vs. stack: concatenate joins along an existing axis; stack creates a new axis. Stacking 32 vectors of shape (F,) gives (32, F).
• einsum: a compact notation for tensor contractions. np.einsum('bhqd,bhkd->bhqk', Q, K) computes attention scores for all heads at once.
• Reductions with axis: x.sum(axis=0) collapses rows (per-feature totals); keepdims=True keeps the axis so the result still broadcasts correctly, which is essential in softmax.

Distance and similarity metrics

Much of ML is "find things that are similar": nearest-neighbor classification, clustering, retrieval for RAG, recommendation, deduplication and anomaly detection. The choice of distance or similarity metric defines what "similar" means, so the wrong metric silently ruins an otherwise good pipeline.

The main metrics (A = [1, 3], B = [4, 7])

Cosine vs. Euclidean. A 100-word and a 1,000-word document on the same topic have count vectors pointing in almost the same direction but with very different lengths: cosine says "very similar", Euclidean says "far apart". For L2-normalized vectors the two agree in ranking, because ||a − b||² = 2 − 2 cos θ. That identity is why vector databases often normalize embeddings and then use fast inner-product search.

Mahalanobis intuition. It measures distance in units of standard deviation along each principal direction of the data. A point 3 units away along a tightly clustered direction is far more anomalous than a point 3 units away along a loose direction, and Mahalanobis captures that while Euclidean does not. With Σ = I it reduces to Euclidean distance.

Distributional "distances": KL divergence (asymmetric), Jensen-Shannon divergence (symmetric, bounded), Wasserstein or earth mover's distance (used in WGANs and drift detection), and the Kolmogorov-Smirnov statistic. These compare whole distributions rather than two points.

from numpy.linalg import norm
a, b = np.array([1, 3]), np.array([4, 7])
norm(a - b), norm(a - b, 1), norm(a - b, np.inf)   # 5.0, 7.0, 4.0
cos = a @ b / (norm(a) * norm(b))                  # 0.981
# batch cosine similarity: normalize once, then one matrix multiply
E = emb / norm(emb, axis=1, keepdims=True)
sims = E @ E.T                                     # (n, n)

The curse of dimensionality

In high dimensions, distances between random points concentrate: the nearest and farthest neighbors end up at nearly the same distance, so "nearest" loses meaning. Volume grows exponentially, so data becomes sparse, and a fixed-size neighborhood contains almost no points. Remedies include dimensionality reduction (PCA, learned embeddings), feature selection, using cosine similarity on meaningful embeddings, and models that do not rely on raw distances.

Calculus: derivatives, gradients and the chain rule

Linear algebra lets a network compute; calculus lets it learn. A derivative tells you how much an output changes when you nudge an input. With millions of weights, the gradient collects one such sensitivity per weight into a vector that points uphill on the loss surface. Step the other way and the loss goes down. The chain rule lets us compute all those sensitivities efficiently through many layers, which is backpropagation.

Derivatives: the slope detector

Optimization by hand. A coffee shop's daily profit at price p is −20p² + 200p − 300. Setting the derivative −40p + 200 to zero gives p = $5 and a profit of $200; at $4 or $6 profit is $180. The second derivative −40 < 0 confirms a maximum. Training a model is the same idea with millions of "prices", except we cannot solve the equation in closed form, so we walk towards it iteratively.

ML's must-know derivatives

• Sigmoid: σ(x) = 1/(1 + e^−x), σ'(x) = σ(x)(1 − σ(x)). Derivation: σ = (1 + e^−x)⁻¹, chain rule gives e^−x/(1 + e^−x)², which factors into σ(1 − σ). Its maximum is 0.25 at x = 0; at x = 5 it is only 0.0066. That flatness causes vanishing gradients.
• Tanh: tanh'(x) = 1 − tanh²(x), maximum 1 at 0.
• ReLU: derivative 1 for x > 0, 0 for x < 0; undefined at exactly 0, where frameworks use a subgradient (0). Hitting exactly 0 is rare, so this does not matter in practice.
• Softmax + cross-entropy: for logits z, softmax p and one-hot target y, ∂L/∂z = p − y. Beautifully simple: the gradient is "predicted minus true". The same is true for sigmoid + binary cross-entropy and for linear output + MSE (up to a constant).
• MSE: d/dŷ (ŷ − y)² = 2(ŷ − y).
• L2 penalty: d/dw (λw²) = 2λw, which shrinks each weight towards 0 proportionally (weight decay). L1 penalty: d/dw (λ|w|) = λ sign(w), a constant push that can drive weights exactly to 0.

Partial derivatives and the gradient

With several inputs, a partial derivative asks "how does the output change if I nudge just this one input and hold the others fixed?". The gradient stacks all partial derivatives into a vector.

The directional derivative in unit direction u is ∇f · u, which is maximized when u points along ∇f (by the Cauchy-Schwarz inequality). The gradient is perpendicular to the contour lines (level sets) of the function. In a network with a billion weights, the gradient is a vector with a billion entries, one per weight.

Worked training step. A one-weight house-price model predicts ŷ = w · area + b. For a 1,500 sq ft house that sold for $300,000, with w = 100, b = 50,000, the prediction is $200,000 and the squared-error loss is (−100,000)² = 10¹⁰. Then ∂L/∂w = 2(ŷ − y) · area = −3 × 10⁸ and ∂L/∂b = 2(ŷ − y) = −2 × 10⁵. Both are negative, so increasing w and b reduces loss. Notice how the gradient for w is 1,500 times larger than for b purely because of the feature's scale, which is why unscaled features make gradient descent zig-zag and need a tiny learning rate.

The chain rule: the heart of backpropagation

For a composition L = f(g(h(x))), dL/dx = f' · g' · h'. A neural network is a long composition, so the gradient of the loss with respect to an early weight is a product of many local derivatives. Computing these products from the output backwards, and reusing shared intermediate results, is backpropagation (next section).

Jacobian and Hessian

• Jacobian J: for a vector function f: Rⁿ → R^m, the m×n matrix of all first partial derivatives, J_ij = ∂f_i/∂x_j. Example: f(x, y) = [x²y, 5x + sin y] has J = [[2xy, x²], [5, cos y]]. The chain rule for vector functions is a product of Jacobians; backprop computes vector-Jacobian products without ever materializing the full Jacobian. The softmax Jacobian is diag(p) − pp^T.
• Hessian H: for a scalar function, the n×n matrix of second partial derivatives, H_ij = ∂²f/∂x_i∂x_j. It describes curvature. For f = w₁² + 3w₁w₂ + w₂², H = [[2, 3], [3, 2]], with eigenvalues 5 and −1. Mixed signs mean the origin (where the gradient is zero) is a saddle point, not a minimum.
• Second-derivative test: at a point with zero gradient, H positive definite (all eigenvalues > 0) means local minimum, negative definite means local maximum, mixed signs mean saddle.
• Cost: for n = 10⁹ parameters, the Hessian has 10¹⁸ entries, so it is never formed. Methods use Hessian-vector products, diagonal approximations or quasi-Newton approximations (L-BFGS) instead.

Taylor series: why small gradient steps work

Example: e^x = 1 + x + x²/2 + x³/6 + ... At x = 1: one term gives 1, two give 2, three give 2.5, five give 2.708, versus the true 2.718.

• Gradient descent comes from the first-order approximation: with Δ = −η∇f, f decreases by about η||∇f||², as long as η is small enough that the linear approximation holds.
• Newton's method minimizes the second-order approximation exactly: Δ = −H⁻¹∇f. It converges in very few steps near a minimum but is too expensive for large networks and is attracted to saddle points.
• Learning-rate limit: along a direction with curvature λ, gradient descent is stable only if η < 2/λ. The sharpest direction (largest Hessian eigenvalue) caps the learning rate for the whole model.

Integration in one paragraph

Integration is the reverse of differentiation and computes areas: ∫₀³ x² dx = [x³/3]₀³ = 9. In ML it appears mostly in probability: a density integrates to 1, P(a ≤ X ≤ b) is an area, E[X] = ∫ x f(x) dx, and the Bayesian evidence P(data) = ∫ P(data | θ)P(θ) dθ is usually intractable, which is why we use Monte Carlo sampling and variational approximations. ROC-AUC is the integral of the ROC curve, computed numerically with the trapezoidal rule.

import torch
w = torch.tensor([1.0, 2.0], requires_grad=True)
f = w[0]**2 + 3*w[0]*w[1] + w[1]**2
f.backward()
w.grad                       # tensor([8., 7.])

# numerical gradient check (central difference) - use to test hand-written gradients
def num_grad(f, x, h=1e-5):
    g = np.zeros_like(x)
    for i in range(len(x)):
        e = np.zeros_like(x); e[i] = h
        g[i] = (f(x + e) - f(x - e)) / (2 * h)
    return g

Backpropagation and vanishing/exploding gradients

Backpropagation is the algorithm that computes the gradient of the loss with respect to every weight in a network. It is nothing more than the chain rule applied systematically from the output layer backwards, caching intermediate values from the forward pass so that nothing is computed twice. It costs roughly the same as two or three forward passes, no matter how many weights there are, which is what makes training billion-parameter models feasible.

Worked example: a tiny two-layer network

One input x = 0.5, a hidden weight w₁ = 0.8, an output weight w₂ = 1.2, sigmoid activations, target y = 1 and loss L = ½(ŷ − y)².

1. Forward: hidden z₁ = x · w₁ = 0.4, h = σ(0.4) = 0.599.
2. Forward: output z₂ = h · w₂ = 0.718, ŷ = σ(0.718) = 0.672.
3. Loss L = ½(0.672 − 1)² = 0.0537.
4. Output error ∂L/∂ŷ = ŷ − y = −0.328; ∂ŷ/∂z₂ = ŷ(1 − ŷ) = 0.220; so δ₂ = ∂L/∂z₂ = −0.0722.
5. Gradient for w₂ ∂L/∂w₂ = δ₂ · h = −0.0722 × 0.599 = −0.0432.
6. Pass error backwards ∂L/∂h = δ₂ · w₂ = −0.0867; δ₁ = ∂L/∂z₁ = −0.0867 × h(1 − h) = −0.0867 × 0.240 = −0.0208.
7. Gradient for w₁ ∂L/∂w₁ = δ₁ · x = −0.0104.
8. Update (learning rate 0.5) w₁ = 0.8 + 0.0052 = 0.805; w₂ = 1.2 + 0.0216 = 1.222. Both increase, pushing ŷ towards 1.

Notice that the early weight's gradient (−0.0104) is smaller than the later one's (−0.0432), because it was multiplied by an extra sigmoid derivative (at most 0.25) and weight. Repeat that over 30 layers and early layers receive almost no signal.

Backprop in matrix form

For a dense layer Z = XW + b with upstream gradient G = ∂L/∂Z (shape B×h):

# one training step of a 2-layer MLP with softmax cross-entropy, in NumPy
Z1 = X @ W1 + b1;  A1 = np.maximum(0, Z1)
Z2 = A1 @ W2 + b2
P = np.exp(Z2 - Z2.max(1, keepdims=True)); P /= P.sum(1, keepdims=True)
loss = -np.log(P[np.arange(B), y]).mean()

G2 = P.copy(); G2[np.arange(B), y] -= 1; G2 /= B      # dL/dZ2 = (p - onehot)/B
dW2 = A1.T @ G2;  db2 = G2.sum(0)
G1 = (G2 @ W2.T) * (Z1 > 0)                          # ReLU derivative
dW1 = X.T @ G1;   db1 = G1.sum(0)
for p, g in [(W1, dW1), (b1, db1), (W2, dW2), (b2, db2)]:
    p -= lr * g

Automatic differentiation (PyTorch autograd, JAX, TensorFlow) builds a computational graph during the forward pass and applies exactly these rules in reverse. Reverse mode is efficient when there are many inputs (weights) and one scalar output (the loss); forward mode is efficient in the opposite case. The trade-off is memory: activations from the forward pass must be stored for the backward pass, which is why activation checkpointing (recompute instead of store) exists for large models.

Vanishing and exploding gradients

Backprop multiplies one local factor per layer. If those factors are typically below 1 the product shrinks exponentially (vanishing); if above 1 it grows exponentially (exploding).

Symptoms: vanishing shows as early layers' weights barely changing, a loss that plateaus early, and tiny gradient norms in lower layers. Exploding shows as loss spikes, NaN or inf values, and huge gradient norms.

Optimization: gradient descent, momentum, Adam and schedules

Optimization is the procedure that turns gradients into better weights. The basic recipe is always the same, take a step against the gradient, but the details (how much data per step, how to smooth noisy gradients, how to adapt step sizes per parameter, and how to change the learning rate over time) decide whether training is fast, stable and ends in a solution that generalizes.

The update rule and the learning rate

• Too small (η = 0.0001): about 10,000 steps to make the same progress; may stall on plateaus.
• Just right (η = 0.1): smooth convergence.
• Too large (η = 1.5): each step multiplies w by 1 − 2(1.5) = −2, so 5 → −10 → 20 → −40: divergence. For this loss the stability limit is η < 1, which is 2/curvature.

Batch, stochastic and mini-batch gradient descent

In practice "SGD" almost always means mini-batch SGD. The mini-batch gradient is an unbiased estimate of the full gradient with variance proportional to 1/B. Some gradient noise is helpful: it tends to steer training towards flatter minima that generalize better. A common heuristic when increasing the batch size by k is to increase the learning rate by about k (linear scaling), with warmup.

Momentum and Nesterov

In a long narrow valley, plain GD bounces between the walls and creeps along the floor. Momentum damps the bouncing and accelerates along the floor. Nesterov momentum evaluates the gradient at the "look-ahead" point w − ηβv, which corrects overshooting a little earlier.

Adaptive methods: AdaGrad, RMSProp, Adam, AdamW

• AdaGrad: divides the step by the square root of the sum of all past squared gradients. Great for sparse features (rare words get larger steps), but the denominator only grows, so learning eventually stalls.
• RMSProp: uses an exponential moving average of squared gradients instead, s = 0.9s + 0.1g², and steps by ηg/(√s + ε). Parameters with large, noisy gradients get smaller steps.
• Adam combines momentum (first moment) and RMSProp (second moment) with bias correction:

Why bias correction? m and v start at 0, so early averages are biased towards 0. At t = 1, m₁ = 0.1g but m̂₁ = g; v₁ = 0.001g² but v̂₁ = g². So the first update is about η · sign(g): every parameter moves by roughly η regardless of its gradient's scale, which shows Adam's scale invariance.

AdamW decouples weight decay from the adaptive update: w = w − η(m̂/(√v̂ + ε) + λw). With plain Adam, an L2 term added to the loss gets divided by √v̂, so heavily updated weights are barely regularized. AdamW is the default for transformers. Other optimizers you may hear about: Adafactor (memory-light second moments), LAMB/LARS (layer-wise scaling for huge batches), Lion (sign-based), and second-order-inspired methods such as Shampoo.

Learning-rate schedules

Why warmup? At the start, weights are random, Adam's second-moment estimates are unreliable, and gradients can be large. A full learning rate at step 0 often causes an early divergence. Ramping up over a few hundred or thousand steps lets statistics settle. Why decay? Late in training, a smaller step reduces the noise floor so the model can settle into a minimum instead of bouncing around it.

import torch
opt = torch.optim.AdamW(model.parameters(), lr=3e-4, betas=(0.9, 0.95), weight_decay=0.1)
sched = torch.optim.lr_scheduler.OneCycleLR(opt, max_lr=3e-4, total_steps=10_000, pct_start=0.03)
for x, y in loader:
    loss = loss_fn(model(x), y)
    opt.zero_grad(); loss.backward()
    torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)
    opt.step(); sched.step()

Convexity and the loss landscape

A function is convex if the line segment between any two points on its graph lies above the graph: f(θx + (1−θ)y) ≤ θf(x) + (1−θ)f(y). Equivalently (for twice-differentiable functions) its Hessian is positive semi-definite everywhere. For convex functions every local minimum is a global minimum, and gradient descent with a suitable learning rate is guaranteed to find it. Linear regression with MSE, logistic regression with cross-entropy, ridge/lasso and SVMs are convex. Neural networks are not.

• Local minima: in high dimensions, most local minima of large networks have loss values close to the global minimum; bad local minima are rarely the practical problem.
• Saddle points and plateaus: far more common than local minima in high dimensions; momentum and gradient noise help escape them.
• Sharp vs. flat minima: flat minima (low curvature) tend to generalize better; small batches, weight decay and techniques such as sharpness-aware minimization favor them.
• Condition number: the ratio of largest to smallest curvature. High condition numbers produce zig-zagging; feature scaling, normalization layers and adaptive optimizers all effectively reduce it.

Regularization, seen as optimization

Regularization adds a complexity penalty to the loss: L_total = L_data + λR(w). With L2 (ridge, weight decay), R = ||w||₂², which shrinks all weights smoothly. With L1 (lasso), R = ||w||₁, which pushes small weights exactly to zero, producing sparse models and built-in feature selection. Example with λ = 0.01 and weights [3.0, 0.1, 0.5, 2.0]: the L1 penalty is 0.01 × 5.6 = 0.056, the L2 penalty is 0.01 × 13.26 ≈ 0.133, and L2 punishes the large weight 3.0 disproportionately. Geometrically, the L1 constraint region is a diamond whose corners lie on the axes, so the loss contours usually touch it at a corner where some weights are zero; the L2 region is a circle with no corners. Other regularizers include dropout, early stopping, data augmentation and label smoothing. The MLE/MAP section below shows that L2 and L1 are Gaussian and Laplace priors in disguise.

Information theory: entropy, cross-entropy and KL divergence

Information theory measures surprise and uncertainty. It gives us the most widely used loss in deep learning (cross-entropy), the standard way to compare two distributions (KL divergence), the evaluation metric for language models (perplexity), and the split criterion in decision trees (information gain). The key idea is simple: rare events carry more information than common ones.

Information content and entropy

• Fair coin: H = −(0.5 log₂ 0.5 + 0.5 log₂ 0.5) = 1 bit, the maximum for two outcomes.
• Biased coin (90% heads): H = −(0.9 log₂ 0.9 + 0.1 log₂ 0.1) ≈ 0.47 bits: more predictable.
• Certain outcome: H = 0. Uniform over K outcomes: H = log₂ K (a fair die: 2.58 bits). Entropy is maximized by the uniform distribution.
• Entropy is the minimum average number of bits needed to encode samples from the distribution (Shannon's source coding theorem), which is why it connects to compression.
• In ML: the entropy of a model's predicted distribution measures its uncertainty (useful for active learning and uncertainty-based routing); entropy regularization encourages exploration in RL; decision trees choose splits that reduce entropy.

Cross-entropy: the classification loss

Worked example. The label is "cat". If the model predicts [cat 0.9, dog 0.1], the loss is −ln 0.9 = 0.105. If it predicts [cat 0.2, dog 0.8], the loss is −ln 0.2 = 1.609. If it predicts cat with 0.01, the loss is 4.6; as the probability goes to 0 the loss goes to infinity. Cross-entropy punishes confident mistakes very heavily, which pushes models towards calibrated probabilities.

Binary cross-entropy: L = −[y log p + (1 − y) log(1 − p)], used with a sigmoid output. Categorical cross-entropy: L = −Σ y_k log p_k, used with softmax. Class labels as strings are first converted to integer indices or one-hot vectors. The minus sign exists because log-probabilities are negative; negating turns "maximize the log-probability of the correct class" into "minimize a positive loss".

KL divergence

Worked example (nats). P = [0.5, 0.5], Q = [0.9, 0.1]. D_KL(P||Q) = 0.5 ln(0.5/0.9) + 0.5 ln(0.5/0.1) = −0.294 + 0.805 = 0.511. D_KL(Q||P) = 0.9 ln 1.8 + 0.1 ln 0.2 = 0.529 − 0.161 = 0.368. Different values, so the order matters.

• Cross-entropy = entropy + KL. Since H(P) of the labels is fixed, minimizing cross-entropy is exactly minimizing KL from the data distribution to the model.
• Forward KL D(P||Q) (used by MLE) is "mean-seeking" or "mass-covering": Q is heavily penalized for assigning near-zero probability where P has mass, so it spreads out to cover all modes. Reverse KL D(Q||P) (used in variational inference) is "mode-seeking": Q prefers to fit one mode well and ignore others.
• Where KL appears: VAE loss (KL between the encoder's Gaussian and a standard normal prior, closed form ½Σ(μ² + σ² − log σ² − 1)); knowledge distillation (student matches the teacher's temperature-softened distribution); RLHF and DPO (a KL penalty keeps the tuned policy close to the reference model to prevent reward hacking); data-drift monitoring.
• Jensen-Shannon divergence: JS(P, Q) = ½KL(P||M) + ½KL(Q||M) with M = (P + Q)/2. Symmetric and bounded by ln 2; the original GAN objective implicitly minimizes it.

Mutual information and information gain

Worked example. In the weather/traffic table (rain 0.4, heavy traffic 0.45, joint rainy-heavy 0.30), the mutual information is about 0.18 bits: knowing the weather removes about 0.18 of the 0.99 bits of uncertainty about traffic. Information gain in decision trees is exactly mutual information between a split and the label: Gain = H(parent) − Σ(n_i/n)H(child_i). The Gini impurity 1 − Σp_k² is a cheaper alternative with similar behavior. MI is also used for feature selection and in contrastive learning (the InfoNCE loss is a lower bound on mutual information).

Perplexity

Worked example. A language model assigns probabilities 0.5, 0.25 and 0.1 to three true next tokens. Mean negative log-likelihood = (0.693 + 1.386 + 2.303)/3 = 1.461 nats, so PPL = e^1.461 ≈ 4.31: the model was as unsure as choosing uniformly among about 4.3 tokens. A model with a training loss of 2.0 nats per token has perplexity e² ≈ 7.4. Perplexities are only comparable between models that use the same tokenizer and test set.

def entropy(p):  p = np.asarray(p); return -(p * np.log2(p)).sum()
entropy([0.5, 0.5]), entropy([0.9, 0.1])        # 1.0, 0.469 bits
P, Q = np.array([0.5, 0.5]), np.array([0.9, 0.1])
kl = (P * np.log(P / Q)).sum()                  # 0.511 nats
import torch.nn.functional as F
loss = F.cross_entropy(logits, targets)         # takes raw logits; applies log-softmax internally
kd = F.kl_div(F.log_softmax(s/T, -1), F.softmax(t/T, -1), reduction='batchmean') * T*T

Maximum likelihood (MLE) and maximum a posteriori (MAP)

MLE and MAP are the two philosophies behind almost every loss function. MLE picks the parameters that make the observed data most probable. MAP does the same but also multiplies in a prior belief about the parameters. Once you see this, many "recipes" stop being arbitrary: cross-entropy is MLE for a categorical model, MSE is MLE for Gaussian noise, and weight decay is a Gaussian prior.

MLE

Coin example. 70 heads in 100 flips. Log-likelihood = 70 log p + 30 log(1 − p). Setting the derivative 70/p − 30/(1 − p) to zero gives p̂ = 0.7, the observed fraction. Gaussian example. For data from N(μ, σ²), the MLE of μ is the sample mean and the MLE of σ² is the 1/n (biased) variance. MLE estimators are consistent and asymptotically efficient, but they can overfit with small data: a coin that landed heads 3 times out of 3 gets p̂ = 1.

Losses are negative log-likelihoods

The choice of loss is therefore an implicit assumption about how the targets are generated. If your regression residuals are heavy-tailed, MSE's Gaussian assumption is wrong and a Laplace (MAE) or Huber loss is more appropriate.

MAP

• Gaussian prior P(w) ∝ exp(−||w||²/(2τ²)) gives −log P(w) = ||w||²/(2τ²) + const: this is L2 regularization, with λ = 1/(2τ²) (up to the noise-variance scaling). A tight prior (small τ) means strong regularization.
• Laplace prior P(w) ∝ exp(−|w|/b) gives L1 regularization; the Laplace density's sharp peak at zero is why L1 produces exact zeros.
• Beta prior on a coin: with prior Beta(2, 2) (a gentle belief that the coin is roughly fair) and 7 heads in 10 flips, the posterior is Beta(9, 5). MAP = (7 + 1)/(10 + 2) = 0.667 versus MLE = 0.7; the posterior mean is 9/14 ≈ 0.643. With 3 heads out of 3, MAP gives 4/5 = 0.8 instead of MLE's overconfident 1.0.
• As data grows, the likelihood dominates the prior, and MAP converges to MLE.

Numerical stability: log-sum-exp, stable softmax and precision

Math on paper uses real numbers with infinite precision; computers use floating-point numbers with a limited range and precision. Naive implementations of softmax, log-likelihoods and variance can overflow to infinity, underflow to zero or produce NaN. A few standard tricks prevent this, and knowing them is a strong signal in interviews.

Floating-point limits

In FP32, e⁸⁹ already overflows to infinity, and e⁻¹⁰⁴ underflows to 0. In FP16, e¹² ≈ 162,755 already overflows. Logits of 100 or more are common, so naive exponentiation is dangerous.

The log-sum-exp trick

Example. x = [1000, 999]. Naively, e¹⁰⁰⁰ overflows and the answer is inf. With the trick: 1000 + log(e⁰ + e⁻¹) = 1000 + log(1.368) = 1000.313.

Stable softmax and log-softmax

def softmax(x, axis=-1):
    z = x - x.max(axis=axis, keepdims=True)       # shift for stability
    e = np.exp(z)
    return e / e.sum(axis=axis, keepdims=True)

def log_softmax(x, axis=-1):
    m = x.max(axis=axis, keepdims=True)
    return x - m - np.log(np.exp(x - m).sum(axis=axis, keepdims=True))

softmax(np.array([1000., 999.]))     # [0.731, 0.269]   naive version returns [nan, nan]
from scipy.special import logsumexp, expit   # expit = numerically safe sigmoid

This is why frameworks fuse softmax and cross-entropy into one operation that takes raw logits: computing log(softmax(x)) in two steps can produce log(0) = −inf when a probability underflows, while x_y − logsumexp(x) never does. The same applies to sigmoid + BCE (BCEWithLogitsLoss), which uses the identity log(1 + e^z) = max(z, 0) + log(1 + e^−|z|).

Work in log space

The probability of a 1,000-token sentence is a product of 1,000 probabilities around 0.1, i.e. about 10⁻¹⁰⁰⁰, far below FP32's smallest positive number. Sums of logs (−2,303 nats here) are perfectly representable. Every language model, HMM, Naive Bayes implementation and beam search works with log-probabilities for this reason.

Other stability tricks

• Epsilon in denominators and logs: normalization layers use √(σ² + ε) with ε ≈ 10⁻⁵; Adam uses √v̂ + ε; clip probabilities to [ε, 1 − ε] before taking logs.
• Loss scaling in FP16 mixed precision: multiply the loss by a large factor (for example 1,024) so small gradients do not underflow, then unscale before the optimizer step; dynamic scalers back off when inf/NaN appears.
• Stable variance: Welford's online algorithm instead of E[X²] − E[X]², which suffers catastrophic cancellation when the mean is large.
• log1p and expm1: log(1 + x) and e^x − 1 for tiny x lose all precision if computed directly.
• Attention scaling and masking: dividing by √d_k keeps logits moderate; masks use a large negative number (or −inf with care) so masked positions get zero weight; fused attention kernels (FlashAttention-style) use an online softmax with running max for stability and memory efficiency.
• Condition numbers: prefer solvers (QR, Cholesky, SVD) to explicit inverses; add a small ridge λI to near-singular matrices.

The activation function and loss function zoos

Activation functions add the non-linearity that lets neural networks represent curved decision boundaries; loss functions define what "good" means. Each comes with a mathematical shape, a derivative that determines how gradients flow, and a set of typical failure modes. Picking them correctly is a mix of math and convention, and interviewers love to probe the reasoning.

Why non-linearity is essential

Without activations, a stack of layers collapses: W₃(W₂(W₁x)) = (W₃W₂W₁)x, a single linear map that cannot even solve XOR. With a non-linear function between layers, a network with enough hidden units can approximate any continuous function on a bounded domain (the universal approximation theorem); depth makes such approximations far more parameter-efficient.

Activation functions

Choosing: hidden layers of CNNs and MLPs: ReLU (or variants). Transformers: GELU or SwiGLU. RNN/LSTM: tanh for states, sigmoid for gates. Output layer: sigmoid for binary or multi-label, softmax for multi-class, identity for regression, softplus or exp for positive targets, tanh for bounded targets in [−1, 1].

Gated linear units. SwiGLU computes (xW₁ ⊙ SiLU(xW₂))W₃: one branch decides how much of the other to let through. It uses three weight matrices instead of two, so the hidden size is usually reduced to about 8d/3 to keep parameters constant.

Regression losses

Example. Predictions [$320K, $190K] vs. actual [$300K, $200K]: errors +20K and −10K. MSE = (400M + 100M)/2 = 250M (RMSE ≈ 15.8K); MAE = (20K + 10K)/2 = 15K. Add one wild error of 200K and the MSE explodes while the MAE grows modestly.

Classification and ranking losses

import torch.nn as nn
nn.MSELoss(); nn.L1Loss(); nn.HuberLoss(delta=1.0)
nn.BCEWithLogitsLoss(pos_weight=torch.tensor([9.0]))   # logits in; weight positives for imbalance
nn.CrossEntropyLoss(label_smoothing=0.1)                # logits in; integer class targets

The math of generative AI: softmax temperature, attention and positional encoding

Large language models are built from the pieces above: embeddings (vectors), matrix multiplications, softmax, dot-product similarity, normalization, residual connections and cross-entropy training. This section walks through the specific formulas that make a transformer work, at the level needed to explain them in an interview. The architecture itself is covered in depth on the Transformers & LLMs page.

Tokens, embeddings and logits

Text is split into tokens, each mapped to an integer ID, and each ID selects a row of an embedding matrix E of shape V × d (vocabulary size by model dimension). A 100,000-token vocabulary with d = 4,096 is already 410 million parameters. After the transformer layers, the final hidden vector h is multiplied by an output matrix (often tied to E) to produce one logit per vocabulary entry: a raw, unnormalized score. Softmax turns the logits into a probability distribution over the next token, and training minimizes cross-entropy against the actual next token.

Softmax with temperature

Temperature changes the probabilities, but not their ranking, so greedy decoding gives the same token at any positive temperature. T = 0 is implemented as argmax (dividing by zero is avoided). Negative temperatures would invert the ranking and are not used. Typical settings: 0 to 0.3 for code, extraction and factual answers; about 0.7 for general chat; 0.9 to 1.2 for creative writing. Temperature is also used in knowledge distillation (a teacher's softmax at T = 2 to 5 reveals "dark knowledge" about which wrong classes are similar) and in contrastive losses (a small τ such as 0.07 sharpens similarity distributions).

Sampling strategies as probability operations

• Greedy: pick argmax p. Deterministic; can loop or be bland.
• Top-k: keep the k highest-probability tokens, renormalize them to sum to 1, sample. A fixed k is too many when the model is confident and too few when it is uncertain.
• Top-p (nucleus): keep the smallest set of tokens whose cumulative probability reaches p (for example 0.9), renormalize, sample. The set size adapts to the model's confidence, so top-p behaves like a "dynamic k".
• Combined: filters are applied in sequence, so the surviving set is the intersection. With probabilities C = 0.22, F = 0.18, D = 0.16, A = 0.13, ... top-k = 2 keeps {C, F}; top-p = 0.95 would keep eight tokens; the intersection is {C, F}, so top-k is the binding constraint.
• Beam search: keep the b best partial sequences by total log-probability. Good for translation; tends to be bland for open-ended text. Length normalization divides the log-probability by length to avoid favoring short outputs.
• Repetition penalties: reduce the logits of tokens already generated. Min-p sampling keeps tokens whose probability is at least a fraction of the top token's probability.

Scaled dot-product attention

1. Project Each token produces a query ("what am I looking for?"), a key ("what do I contain?") and a value ("what information do I pass on?").
2. Score Dot product of every query with every key: high when they point in similar directions.
3. Scale Divide by √d_k to keep scores in a range where softmax still has useful gradients.
4. Mask (optional) Add −∞ (in practice a large negative number) to disallowed positions, e.g. future tokens in a decoder, so their weights become exactly 0 after softmax.
5. Normalize Row-wise softmax: each token's weights over all tokens are positive and sum to 1.
6. Mix Multiply by V: each output is a weighted average of value vectors.

Worked example (3 tokens, d_k = 2). Q = K = [[1, 0], [0, 1], [1, 1]], V = [[2, 0], [0, 3], [1, 1]]. For the third token, the raw scores against the three keys are [1, 1, 2]. Divided by √2 they become [0.707, 0.707, 1.414]. Softmax gives [0.248, 0.248, 0.503]. The output is 0.248[2, 0] + 0.248[0, 3] + 0.503[1, 1] = [1.00, 1.25]: mostly its own value, blended with the others.

Why √d_k? If the components of q and k are independent with mean 0 and variance 1, then q·k = Σ q_ik_i is a sum of d_k terms each with variance 1, so its variance is d_k and its standard deviation √d_k. With d_k = 64, raw scores have a typical magnitude around 8; softmax of such values is nearly one-hot, its Jacobian diag(p) − pp^T is nearly zero, and gradients vanish. Dividing by √d_k = 8 restores unit variance. The factor is not a guess; it is exactly the standard deviation of the dot product.

Why softmax before multiplying by V? The scores decide where to look; V decides what to retrieve. Normalizing the weights first makes the output a convex combination of value vectors, keeping its scale stable regardless of sequence length. Cost: computing QK^T takes O(T²d) time and O(T²) memory per head, the reason long contexts are expensive and why efficient kernels and approximations exist.

def attention(Q, K, V, mask=None):
    d_k = Q.shape[-1]
    scores = Q @ K.swapaxes(-1, -2) / np.sqrt(d_k)         # (..., T, T)
    if mask is not None:
        scores = np.where(mask, scores, -1e9)             # causal / padding mask
    w = softmax(scores, axis=-1)                          # rows sum to 1
    return w @ V, w

T = 4
causal = np.tril(np.ones((T, T), dtype=bool))             # token i sees tokens 0..i

Multi-head attention

A single softmax can only produce one weighting pattern per token; multiple heads let the model attend to different relationships at once (neighboring tokens, subject-verb links, coreference, copying earlier tokens). The four projection matrices W_Q, W_K, W_V, W_O contribute 4d² parameters per layer; the feed-forward block (d → 4d → d) contributes about 8d², which is why FFNs hold roughly two-thirds of a transformer's parameters. Heads must be initialized randomly and independently, otherwise they would receive identical gradients and remain copies. Grouped-query attention (GQA) and multi-query attention (MQA) share key/value projections across groups of heads, shrinking the KV cache used during generation, whose size is 2 × layers × KV heads × d_k × sequence length × bytes per value (about 0.5 MB per token for a 7B model with 32 layers, 32 heads of 128 dimensions, in FP16).

Positional encoding

Self-attention is permutation-equivariant: shuffle the input tokens and the outputs shuffle the same way. It has no built-in notion of order, yet "dog bites man" and "man bites dog" mean different things. Positional information must be injected.

Example: for the first pair of dimensions (i = 0), position 1 gets [sin 1, cos 1] = [0.841, 0.540] and position 2 gets [sin 2, cos 2] = [0.909, −0.416]. Using sine and cosine together means each pair is a point on a circle; moving k positions forward is a fixed rotation of that point, so PE(pos + k) is a linear function of PE(pos). That lets attention learn relative offsets. The multiple frequencies act like the hands of a clock (seconds, minutes, hours) so every position gets a unique combination.

Rotary position embedding (RoPE)

RoPE splits each query and key vector into pairs of dimensions and rotates pair i by the angle mθ_i, where m is the token's position and θ_i = 10000^−2i/d.

Because rotations preserve length, RoPE does not change vector norms; it only changes angles. It is applied to Q and K (not V) inside every attention layer rather than added once at the input. Low-frequency pairs encode long-range position, high-frequency pairs encode local order. Context-extension methods (position interpolation, NTK-aware scaling, YaRN) rescale the angles so a model trained on, say, 4K tokens can handle much longer inputs with a little fine-tuning.

def rope(x, pos, base=10000.0):
    # x: (T, d) with d even; pos: (T,)
    d = x.shape[-1]
    theta = base ** (-np.arange(0, d, 2) / d)            # (d/2,)
    ang = pos[:, None] * theta[None, :]                  # (T, d/2)
    cos, sin = np.cos(ang), np.sin(ang)
    x1, x2 = x[:, 0::2], x[:, 1::2]
    out = np.empty_like(x)
    out[:, 0::2] = x1 * cos - x2 * sin
    out[:, 1::2] = x1 * sin + x2 * cos
    return out

Normalization and the rest of the block

A transformer block is: x = x + Attention(Norm(x)); x = x + FFN(Norm(x)) (the pre-norm arrangement). The residual additions keep a direct gradient path through dozens of layers, and the normalizations keep activations at a stable scale. For generation, the model runs this stack for each new token, reusing cached keys and values from earlier tokens.

Quick revision

• Mean uses every value and is pulled by outliers; the median is robust and preferred for skewed data such as income or latency.
• Quartiles split sorted data into four parts; IQR = Q3 − Q1; the box-plot outlier rule flags points beyond Q1 − 1.5·IQR or Q3 + 1.5·IQR.
• Variance is the mean squared deviation; standard deviation is its square root, in the data's units.
• Sample variance divides by n − 1 (Bessel's correction) because deviations from the sample mean are systematically too small.
• Right skew: long right tail, mean > median. Excess kurtosis > 0: heavier tails than a normal distribution.
• z = (x − μ)/σ; fit scalers on training data only to avoid leakage.
• P(A ∪ B) = P(A) + P(B) − P(A ∩ B); P(A ∩ B) = P(A)P(B | A); independence means P(A ∩ B) = P(A)P(B).
• P(A | B) ≠ P(B | A); Bayes: posterior = likelihood × prior / evidence.
• 1% prevalence, 95% sensitivity, 5% false-positive rate: a positive test means only about 16% chance of disease (base-rate fallacy).
• Linearity of expectation holds even for dependent variables; Var(aX + b) = a²Var(X).
• Var(X + Y) = Var X + Var Y + 2Cov(X, Y); the variance of a mean of n i.i.d. values is σ²/n.
• Correlation is covariance scaled to [−1, 1]; zero correlation does not imply independence (X and X²).
• Bernoulli (p, p(1−p)), binomial (np, np(1−p)), Poisson (λ, λ), exponential (1/λ, 1/λ²), uniform ((a+b)/2, (b−a)²/12).
• Normal: 68-95-99.7 within 1, 2, 3 standard deviations; 1.96 for a two-sided 95% interval.
• Beta is the conjugate prior for a Bernoulli/binomial rate; Dirichlet for categorical; posterior = prior pseudo-counts + observed counts.
• LLN: sample means converge to the true mean. CLT: sample means are approximately normal with standard error σ/√n, whatever the data's shape.
• A 95% CI means 95% of intervals built this way contain the true value, not a 95% probability for this interval.
• A p-value is P(data at least this extreme | H₀), not P(H₀ | data).
• Type I = false positive (α); Type II = false negative (β); power = 1 − β, increased by larger n, larger effects, lower variance.
• Use Welch's t-test for two means, paired tests for the same units, chi-square for categorical association, z-test for proportions.
• A/B pitfalls: peeking, multiple comparisons, underpowered tests, sample-ratio mismatch, novelty effects.
• Bootstrap puts error bars on a statistic (resample with replacement); CV estimates how well a training procedure generalises (retrain on folds). They are not interchangeable.
• Student-t variance is ν/(ν−2) only for ν > 2; the mean exists only for ν > 1.
• Label smoothing: true class gets (1 − ε) + ε/K, every other class gets ε/K.
• Matrix product (m×n)(n×p) = (m×p); not commutative but associative; (AB)^T = B^TA^T.
• Only square, full-rank matrices (det ≠ 0) are invertible; rank = number of independent rows or columns.
• Normal equation w = (X^TX)⁻¹X^Ty; in practice use lstsq/QR/SVD, never an explicit inverse.
• Av = λv; trace = sum of eigenvalues, determinant = product; symmetric matrices have real eigenvalues and orthogonal eigenvectors.
• SVD A = UΣV^T works for any matrix; truncated SVD is the best low-rank approximation.
• PCA: center/standardize, covariance, eigenvectors = directions of max variance, keep top k by explained variance.
• Broadcasting compares shapes right to left; dimensions must match or be 1; (N,) minus (N, 1) silently becomes (N, N).
• Cosine similarity ignores magnitude; for unit vectors, squared Euclidean distance = 2 − 2cos θ.
• The gradient points in the direction of steepest ascent; gradient descent steps along −∇L.
• Chain rule: multiply local derivatives along a path, sum over paths; backprop applies it from the loss backwards.
• σ'(x) = σ(1 − σ) ≤ 0.25; the softmax + cross-entropy gradient with respect to logits is p − y.
• Hessian eigenvalues all positive = local minimum; mixed signs = saddle point; learning rate must satisfy η < 2/λ_max.
• Vanishing gradients: fix with ReLU-family activations, He/Xavier init, residuals, normalization; exploding: gradient clipping.
• Momentum averages gradients; RMSProp scales by recent gradient magnitude; Adam combines both with bias correction; AdamW decouples weight decay.
• Warmup avoids early divergence; cosine decay lets the model settle; tune learning rates on a log scale.
• Convex losses have one global minimum; neural network losses are non-convex but mostly have good minima and many saddles.
• Entropy H = −Σp log p; cross-entropy = entropy + KL; KL ≥ 0 and is not symmetric.
• Perplexity = e^{cross-entropy per token}; a loss of 2.0 nats means perplexity about 7.4.
• MSE = Gaussian MLE, MAE = Laplace MLE, cross-entropy = categorical MLE; L2 = Gaussian prior, L1 = Laplace prior (MAP).
• Stable softmax subtracts the max; log-sum-exp = m + log Σe^{x − m}; work with log-probabilities; pass logits to fused losses.
• BF16 has FP32's range with less precision, so LLMs train in BF16 without loss scaling; FP16 needs loss scaling.
• Temperature divides logits before softmax: T < 1 sharpens, T > 1 flattens; it never changes the argmax.
• Attention = softmax(QK^T/√d_k)V; dividing by √d_k keeps score variance at 1 so softmax does not saturate.
• Attention is order-blind; sinusoidal, learned, ALiBi or RoPE encodings add position; RoPE rotates Q and K so scores depend on relative distance.

Glossary

A/B test

A randomized experiment comparing a control and a treatment on a pre-chosen metric using a hypothesis test.

Adam

An optimizer that combines momentum (a running mean of gradients) with per-parameter step sizes from a running mean of squared gradients, plus bias correction.

AdamW

Adam with weight decay applied directly to the weights rather than through the loss; the default for transformers.

Backpropagation

Computing gradients of the loss for all parameters by applying the chain rule backwards through the computation graph.

Bayes' theorem

P(A | B) = P(B | A)P(A)/P(B): how to update a prior belief into a posterior after seeing evidence.

Bessel's correction

Dividing by n − 1 instead of n to get an unbiased sample variance.

Bias (of an estimator)

The difference between an estimator's expected value and the true parameter.

Bootstrap

Estimating the variability of a statistic by recomputing it on many resamples drawn with replacement.

Broadcasting

Rules that let array operations combine different shapes by virtually stretching size-1 dimensions.

Central limit theorem

The distribution of a sample mean approaches a normal distribution with standard deviation σ/√n as n grows.

Chain rule

The derivative of a composition is the product of the derivatives of its parts.

Condition number

Ratio of the largest to smallest singular value (or curvature); high values mean numerical instability and slow gradient descent.

Conditional probability

P(A | B) = P(A ∩ B)/P(B): the probability of A within the world where B happened.

Confidence interval

A range built by a procedure that captures the true parameter in a stated fraction (e.g. 95%) of repeated samples.

Conjugate prior

A prior that yields a posterior in the same family, e.g. Beta for a binomial rate.

Convex function

A function whose chords lie above its graph; every local minimum is global.

Correlation

Covariance divided by the product of standard deviations; measures linear association on a [−1, 1] scale.

Cosine similarity

The dot product of two vectors divided by the product of their lengths; the cosine of the angle between them.

Covariance

E[(X − μ_X)(Y − μ_Y)]: whether two variables move together, in their product units.

Covariance matrix

A symmetric positive semi-definite matrix of pairwise feature covariances with variances on the diagonal.

Cross-entropy

−Σp log q: the average surprise of true outcomes under predicted distribution q; the standard classification loss.

Determinant

A scalar giving the volume scaling of a linear map; zero means the matrix is singular.

Dot product

The sum of element-wise products of two vectors; measures alignment and underlies neurons and attention.

Eigenvector / eigenvalue

A direction v that a matrix only scales, Av = λv; λ is the scale factor.

Entropy

−Σp log p: the average uncertainty or information content of a distribution.

Expectation

The probability-weighted average of a random variable.

Gradient

The vector of partial derivatives; points in the direction of steepest increase.

Gradient clipping

Rescaling the gradient when its norm exceeds a threshold, to prevent exploding updates.

Hessian

The matrix of second partial derivatives, describing the curvature of a function.

Hypothesis test

A procedure that decides whether data provide enough evidence to reject a null hypothesis at level α.

i.i.d.

Independent and identically distributed: samples drawn independently from the same distribution.

Independence

Two events or variables where knowing one gives no information about the other: P(A ∩ B) = P(A)P(B).

Interquartile range (IQR)

Q3 − Q1, the spread of the middle half of the data.

Jacobian

The matrix of first partial derivatives of a vector-valued function.

KL divergence

Σp log(p/q): the extra surprise from using q instead of the true p; non-negative and asymmetric.

Kurtosis

A measure of tail heaviness; the normal distribution has kurtosis 3 (excess kurtosis 0).

Law of large numbers

Sample averages converge to the expected value as the sample size grows.

Learning rate

The step-size multiplier in gradient descent.

Likelihood

P(data | parameters) viewed as a function of the parameters.

Log-sum-exp trick

Computing log Σe^x_i as m + log Σe^x_i−m with m = max x to avoid overflow.

Logit

A raw, unnormalized model score that softmax or sigmoid turns into a probability; also log(p/(1 − p)).

MAP estimation

Choosing parameters that maximize likelihood times prior; equivalent to regularized MLE.

Marginal probability

The probability of one variable obtained by summing (or integrating) a joint distribution over the others.

Maximum likelihood (MLE)

Choosing parameters that make the observed data most probable.

Median

The middle value of sorted data; the 0.5 quantile.

Momentum

An optimizer technique that accumulates an exponentially weighted sum of past gradients to smooth and accelerate updates.

Mutual information

How much knowing one variable reduces uncertainty about another; zero only under independence.

Norm

A measure of vector length, e.g. L1 (sum of absolute values) or L2 (Euclidean length).

Null hypothesis

The default claim of no effect or no difference that a test tries to reject.

Orthogonal matrix

A square matrix with orthonormal columns (Q^TQ = I); preserves lengths and angles.

p-value

The probability, assuming the null hypothesis, of a result at least as extreme as the one observed.

Perplexity

The exponential of the average per-token cross-entropy; the effective number of choices a language model is unsure between.

Positional encoding

Information added to token representations so attention can use word order.

Power (statistical)

The probability that a test detects a real effect of a given size; 1 − β.

Principal component analysis (PCA)

Linear dimensionality reduction that projects data onto the top eigenvectors of its covariance matrix.

Quantile

A cut point below which a given fraction of sorted data falls (percentiles, quartiles, deciles).

Rank

The number of linearly independent rows or columns of a matrix.

RoPE

Rotary position embedding: rotating query and key vector pairs by position-dependent angles so attention depends on relative position.

Saddle point

A point with zero gradient that is a minimum in some directions and a maximum in others.

Scaled dot-product attention

softmax(QK^T/√d_k)V: weighting value vectors by normalized query-key similarity.

Singular value decomposition (SVD)

Factorizing any matrix as UΣV^T with orthogonal U, V and non-negative singular values.

Skewness

A measure of asymmetry of a distribution; positive means a longer right tail.

Softmax

A function that turns a vector of logits into positive probabilities summing to 1: e^z_i/Σe^z_j.

Standard deviation

The square root of variance, expressed in the data's units.

Standard error

The standard deviation of an estimator across repeated samples, e.g. σ/√n for the mean.

Temperature

A divisor applied to logits before softmax that controls how sharp or flat the output distribution is.

Tensor

A multi-dimensional array, such as a (batch, tokens, features) activation.

Type I / Type II error

Rejecting a true null (false positive) / failing to reject a false null (false negative).

Variance

The expected squared deviation from the mean.

z-score

The number of standard deviations a value lies from the mean.

Interview questions