Every loss function is a probability statement. Every metric is a random variable. Learn the rules once and most of ML stops being magic.
This page is the theory under the models. It covers the rules of probability, the inequalities that bound error, and the limit theorems behind every A/B test. Then it moves to information theory. Entropy, cross-entropy and KL explain why we train with the losses we use. The last sections cover sketches and sampling, the tools you reach for at scale.
Each topic has a plain definition, the math, the intuition, and why it matters on the job. Each ends with short interview checks. For worked puzzles like Bayes on a medical test, see Statistics and Probability. This page does not repeat them.
Probability space. A sample space Ω of outcomes, a set of events, and a measure P on those events. P follows three axioms. Everything else is derived from them.
The three axioms
Non-negative. P(A) ≥ 0 for every event A.
Normalized. P(Ω) = 1.
Countably additive. For disjoint A1, A2, … we have P(∪i Ai) = ∑i P(Ai).
From these you get the working rules. P(Ac) = 1 − P(A). P(A ∪ B) = P(A) + P(B) − P(A ∩ B). The union bound P(∪ Ai) ≤ ∑ P(Ai) holds with no assumptions. It powers many proofs in learning theory.
Conditional probability
Conditioning shrinks the world to the event B and renormalizes.
P(A | B) = P(A ∩ B) / P(B), P(B) > 0
Chain rule: P(A1, …, An) = P(A1) P(A2 | A1) … P(An | A1, …, An−1)
Total prob.: P(A) = ∑i P(A | Bi) P(Bi) for a partition {Bi}
Bayes: P(B | A) = P(A | B) P(B) / P(A)
The chain rule is exact for any order. An autoregressive language model is just the chain rule over tokens. It models P(xt | x<t) and multiplies.
Independence and conditional independence
Independent. A ⊥ B when P(A ∩ B) = P(A) P(B). Knowing B tells you nothing about A.
Conditionally independent. A ⊥ B | C when P(A ∩ B | C) = P(A | C) P(B | C). Once you know C, B adds nothing about A.
Mutual independence. Every subset factors, not just every pair. Pairwise independence is weaker.
Neither kind of independence implies the other. Two examples make this clear.
Conditional but not marginal. Two features both depend on the class label. Given the label, they are independent. Without it, they are correlated through the label. Naive Bayes assumes exactly this structure.
Marginal but not conditional. Flip two fair coins X and Y. Let Z = X XOR Y. X and Y are independent. Given Z, knowing X fixes Y. This is “explaining away”. Conditioning on a common effect creates dependence.
The XOR example also shows pairwise without mutual independence. X, Y and Z are each pairwise independent. But any two of them fix the third.
Conditioning on a collider creates bias. Suppose you only study users who clicked. Clicks depend on both relevance and position. So inside your data, relevance and position look dependent even if they are not. This is selection bias, and it is the same math as the XOR coins.
Why it matters in practice
Language models are the chain rule. Perplexity and log likelihood come straight from it.
Naive Bayes and graphical models encode conditional independence to cut parameters from exponential to linear.
Training on logged data conditions on what the old system showed. That is collider bias. Inverse propensity weights fix it.
Union bounds give family-wise error control. Bonferroni is the union bound applied to many tests.
Interview check
Does pairwise independence imply mutual independence? No. Take two fair bits and their XOR. Each pair is independent. The triple is not, since any two fix the third.
Give a case where conditioning makes two independent variables dependent. Talent and looks are independent in the population. Among famous actors, who need one or the other, they look negatively correlated. Conditioning on the common effect did that.
What does naive Bayes assume? Features are independent given the class. It does not assume they are independent overall.
Random variables, expectation and variance
Random variable. A function X from outcomes to numbers. Its distribution is fully described by the CDF F(x) = P(X ≤ x). A discrete X has a PMF p(x). A continuous X has a density f(x) with F(x) = ∫−∞x f(t) dt.
Expectation and variance
E[X] = ∑x x p(x) or ∫ x f(x) dx
E[g(X)] = ∑x g(x) p(x) (law of the unconscious statistician)
Var(X) = E[(X − E X)2] = E[X2] − (E X)2
Cov(X, Y) = E[(X − E X)(Y − E Y)] = E[XY] − E[X] E[Y]
ρ(X, Y) = Cov(X, Y) / (σX σY) ∈ [−1, 1]
Expectation is linear with no conditions. E[aX + bY] = a E[X] + b E[Y], even if X and Y are dependent. Variance is not linear.
Var(aX + bY) = a2 Var(X) + b2 Var(Y) + 2ab Cov(X, Y)
Var(wTX) = wT Σ w for a random vector X with covariance Σ
The second line explains why every covariance matrix is positive semi-definite. A variance can never be negative, so wTΣw ≥ 0 for every w.
Zero correlation is not independence
Let X be uniform on [−1, 1] and Y = X2. Then Cov(X, Y) = E[X3] = 0. But Y is a function of X. Correlation only sees linear dependence. Independence implies zero covariance, never the reverse. The one big exception is a jointly Gaussian pair.
Averages of correlated variables
Take n variables, each with variance σ2 and pairwise correlation ρ. Their mean has this variance.
Var(X̄) = ρ σ2 + (1 − ρ) σ2 / n
As n grows, the second term goes to zero. The first does not. This one line explains bagging and random forests. Averaging trees kills the (1 − ρ) part. Random feature subsets lower ρ itself.
Short derivation of the second. Write Var(X) = E[X2] − (E X)2. Apply the tower rule to E[X2] = E[ Var(X|Y) + E[X|Y]2 ]. Then subtract (E[E[X|Y]])2. The last two pieces form Var(E[X|Y]).
The intuition is simple. Total spread splits into spread inside each group plus spread between group means.
Why it matters in practice
Uncertainty in deep ensembles. Total predictive variance splits into mean member variance (aleatoric) plus variance of member means (epistemic).
Variance reduction in A/B tests. Stratification and CUPED remove the “between” term. Only within-stratum noise is left.
Portfolio of models. Var(wTX) tells you when blending two scores helps. Low correlation is the win.
Feature checks. A near-zero Pearson correlation does not mean a feature is useless. Use mutual information for non-linear links.
Interview check
Why does a random forest beat one deep tree? Averaging cuts variance to ρσ2 + (1 − ρ)σ2/n. Bootstrap rows and random features lower ρ. Bias stays about the same.
State the law of total variance and give one use. Var(X) = E[Var(X|Y)] + Var(E[X|Y]). Stratified sampling removes the second term. That is why it beats simple random sampling.
Can two variables have zero covariance and still be dependent? Yes. X uniform on [−1, 1] and Y = X2.
Key inequalities
Tail bound. An upper limit on how likely a random variable is to land far from where you expect. Stronger assumptions buy tighter bounds.
Markov
For X ≥ 0 and any a > 0, we have P(X ≥ a) ≤ E[X] / a.
Proof in one line. a · 1{X ≥ a} ≤ X pointwise. Take expectations. It needs only a mean and non-negativity. So it is weak but always available.
Chebyshev
Apply Markov to (X − μ)2. For any k > 0 we get P(|X − μ| ≥ kσ) ≤ 1 / k2.
Applied to a sample mean, P(|X̄ − μ| ≥ t) ≤ σ2 / (n t2). That is already a proof of the weak law of large numbers.
Jensen
For convex f, f(E[X]) ≤ E[f(X)]. For concave f the sign flips. Equality holds when X is constant or f is linear where X lives.
X takes x1 or x2 with equal odds. The chord midpoint sits above the curve. The vertical gap is E[f(X)] − f(E[X]) ≥ 0.
The ELBO. log p(x) = log Eq[p(x,z)/q(z)] ≥ Eq[log p(x,z) − log q(z)]. That is the bound VAEs and EM maximize.
EM never decreases likelihood. The E-step makes the Jensen bound tight. The M-step raises the bound.
Metric bias. E[log X] ≤ log E[X]. Averaging log-revenue per user is not the log of mean revenue.
Hoeffding
Let X1, …, Xn be independent with Xi ∈ [ai, bi]. For any t > 0:
P( |X̄ − E X̄| ≥ t ) ≤ 2 exp( −2 n2 t2 / ∑i (bi − ai)2 )
For Xi ∈ [0, 1]: P( |X̄ − μ| ≥ t ) ≤ 2 e−2 n t2
Solve for n: n ≥ ln(2/δ) / (2 t2)
Chebyshev decays like 1/n. Hoeffding decays exponentially in n. The price is bounded variables. The proof bounds the moment generating function, then applies Markov to eλX̄. That trick is the Chernoff method.
A worked number. You want a click rate within ±0.01 with 95% confidence. Hoeffding says n ≥ ln(40) / (2 × 0.0001) ≈ 18,445. That holds for any distribution on [0, 1]. The CLT gives a smaller n, about 9,604 at p = 0.5, but only as an approximation.
From Hoeffding to generalization
Take a finite hypothesis class H and a 0-1 loss. Apply Hoeffding to each h. Then take a union bound over all |H| of them. With probability at least 1 − δ, every h satisfies:
|R(h) − R̂n(h)| ≤ √( (ln|H| + ln(2/δ)) / (2n) )
This is the template for all of learning theory. VC dimension and Rademacher complexity replace ln|H| for infinite classes.
Why it matters in practice
Bandits. UCB1 adds √(2 ln t / na) to each arm. That bonus is a Hoeffding confidence radius.
Sizing eval sets. Hoeffding gives a distribution-free sample size for a target accuracy margin.
Variational inference and EM. Both rest on the Jensen bound.
Sketches. Count-min error bounds come from Markov plus independent rows. See the sketching section.
Alert thresholds. Chebyshev gives a safe bound when you cannot assume normality.
Interview check
Prove KL ≥ 0. Apply Jensen to the concave log. Ep[log(q/p)] ≤ log Ep[q/p] = log ∑ q = 0.
Chebyshev or Hoeffding for a bounded metric? Hoeffding. Its tail falls as e−2nt2 instead of 1/n.
Where does the UCB bonus come from? Set the Hoeffding tail e−2nt2 equal to t−4 and solve for the radius. You get √(2 ln t / n).
Laws of large numbers and the CLT
Limit theorems. The law of large numbers says the sample mean converges to the true mean. The central limit theorem says how the error is shaped on the way there.
Modes of convergence
Almost surely. P(lim Xn = X) = 1. The strongest kind.
In probability. For every ε > 0, P(|Xn − X| > ε) → 0.
In distribution. The CDFs converge at every point where the limit CDF is continuous.
Almost sure implies in probability. In probability implies in distribution. None of the reverse holds in general.
Precise statements
Let X1, X2, … be i.i.d. with mean μ. Let X̄n be the mean of the first n.
Weak law. If E|X| < ∞, then X̄n → μ in probability.
Strong law. If E|X| < ∞, then X̄n → μ almost surely. A finite variance is not needed.
Central limit theorem. If also Var(X) = σ2 < ∞, then √n (X̄n − μ) / σ → N(0, 1) in distribution.
Multivariate CLT. For i.i.d. vectors with covariance Σ, √n (X̄n − μ) → N(0, Σ).
Lindeberg-Feller. The CLT also holds for independent, non-identical terms. No single term may dominate the total variance.
How fast: Berry-Esseen
The CLT is a limit. Berry-Esseen gives a rate. Let ρ = E|X − μ|3. Then:
supx | P( √n (X̄n − μ)/σ ≤ x ) − Φ(x) | ≤ C ρ / (σ3 √n), C < 0.48
Skewed data has a large ρ / σ3. So it needs far more samples before normal confidence intervals are honest. Revenue per user is the classic case.
The delta method
Often you need a smooth function of a mean, not the mean itself. If √n (X̄n − μ) → N(0, σ2) and g is differentiable with g′(μ) ≠ 0:
√n ( g(X̄n) − g(μ) ) → N( 0, g′(μ)2 σ2 )
Ratio metric R = Ȳ / X̄ (clicks per impression, per user):
Var(R) ≈ (1/n) [ Var(Y)/μX2 − 2 μY Cov(X,Y)/μX3 + μY2 Var(X)/μX4 ]
When the CLT fails. A Cauchy variable has no mean. Its sample mean is Cauchy again for every n, so it never settles. Power-law data with tail index below 2 has infinite variance. The mean converges, but slower than 1/√n and not to a Gaussian. Cap, log-transform, or use medians and the bootstrap.
Why it matters in practice
A/B tests use the CLT for z-tests on means. Berry-Esseen tells you when a small or skewed test is unsafe.
Ratio metrics like CTR per user need the delta method. Treating each impression as independent gives intervals that are far too narrow.
SGD noise. A mini-batch gradient is a sample mean. Its noise shrinks like 1/√B, which guides batch size and learning rate scaling.
Monte Carlo estimates of any expectation have error O(1/√n), no matter the dimension.
Interview check
Does the strong law need finite variance? No. A finite mean is enough. The CLT is what needs finite variance.
Your CTR is clicks over impressions, randomized by user. How do you get its variance? Aggregate to user-level sums. Apply the delta method to the ratio of means. Or bootstrap over users.
What controls how fast the CLT kicks in? The skew. Berry-Esseen bounds the error by Cρ/(σ3√n).
The multivariate Gaussian
Multivariate Gaussian. A distribution on vectors in ℝd, fixed by a mean μ and a positive definite covariance Σ. Every linear combination of its coordinates is a 1D Gaussian.
The quadratic form in the exponent is the squared Mahalanobis distance. Contours of equal density are ellipses. Their axes are the eigenvectors of Σ. Their radii scale with the square roots of its eigenvalues.
The closure properties
Split x into blocks x1 and x2. Split μ and Σ the same way, with blocks Σ11, Σ12, Σ21, Σ22.
Affine: A x + b ~ N( Aμ + b, A Σ AT )
Marginal: x1 ~ N( μ1, Σ11 ) (just read off the block)
Conditional: x1 | x2 ~ N( μ1|2, Σ1|2 )
μ1|2 = μ1 + Σ12 Σ22−1 (x2 − μ2)
Σ1|2 = Σ11 − Σ12 Σ22−1 Σ21 (Schur complement)
Product: N(x; a, A) N(x; b, B) ∝ N(x; c, C), C = (A−1 + B−1)−1, c = C(A−1a + B−1b)
Read the conditional slowly. The new mean is a linear regression of x1 on x2. The coefficient matrix is Σ12Σ22−1. The new covariance does not depend on the observed value x2. Observing data always shrinks the covariance, since the subtracted term is PSD.
In 2D with correlation ρ, this becomes E[x1 | x2] = μ1 + ρ (σ1/σ2)(x2 − μ2). The variance becomes σ12(1 − ρ2). That is regression to the mean in one line.
The precision matrix
Let Λ = Σ−1. A zero in Σ means two coordinates are marginally independent. A zero in Λ means they are independent given all the others. Sparse precision matrices are Gaussian graphical models. The graphical lasso estimates them with an L1 penalty.
For a jointly Gaussian vector, uncorrelated implies independent. This fails if only the marginals are Gaussian. Take X ~ N(0,1) and Y = SX with S a random sign. Y is N(0,1) and uncorrelated with X. But |Y| = |X|.
Why it is everywhere
The CLT. Sums of many small effects become Gaussian. Noise and averaged gradients both qualify.
Maximum entropy. Among all densities with a given mean and covariance, the Gaussian has the most entropy. So it assumes the least beyond those two moments.
Closed algebra. Affine maps, marginals, conditionals and products all stay Gaussian. Inference is linear algebra.
Conjugacy. A Gaussian prior with a Gaussian likelihood gives a Gaussian posterior.
Squared loss. Minimizing squared error is MLE under Gaussian noise. So every least-squares model is a Gaussian model.
Sampling and the reparameterization trick
Factor Σ = L LT with Cholesky. Draw ε ~ N(0, I). Then x = μ + Lε has the right distribution, by the affine rule. A VAE uses this to backprop through sampling. The randomness sits in ε, so gradients flow to μ and L.
Why it matters in practice
Gaussian processes are the conditional formula applied to a kernel matrix. Bayesian optimization of hyperparameters runs on them.
Kalman filters track state with the product and conditional rules at every step.
Anomaly detection. Squared Mahalanobis distance follows a χ2d law under the model. That gives a principled threshold.
Thompson sampling with a Bayesian linear model keeps a Gaussian posterior over weights. Each update is closed form.
VAEs and diffusion models are built on Gaussian noise and reparameterized samples.
Interview check
Give the conditional mean and covariance of x1 given x2. μ1 + Σ12Σ22−1(x2 − μ2), and Σ11 − Σ12Σ22−1Σ21. The covariance ignores the observed value.
What does a zero in the precision matrix mean? Those two variables are conditionally independent given the rest.
Two Gaussian variables with zero correlation. Are they independent? Only if they are jointly Gaussian. The random-sign example breaks it.
Exponential family and sufficient statistics
Exponential family. Distributions whose log density is linear in a fixed set of statistics T(x). Bernoulli, categorical, Gaussian, Poisson, exponential, gamma, beta and Dirichlet all belong.
Differentiate A under the integral. You get two facts that drive everything else.
∇A(η) = E[ T(X) ]
∇2A(η) = Cov[ T(X) ] ≽ 0
So A is convex. The log likelihood ηT∑T(xi) − nA(η) is concave in η. MLE has no bad local optima. Set the gradient to zero and you get moment matching.
∇A(η̂) = (1/n) ∑i T(xi) model moments = data moments
Three members worked out
Bernoulli. px(1−p)1−x = exp( x log(p/(1−p)) − log(1 + eη) ). So η is the log-odds, and ∇A = 1/(1+e−η) is the sigmoid.
Categorical. ηk are logits. A(η) = log ∑k eηk is logsumexp. Its gradient is softmax.
Poisson. λxe−λ/x! gives η = log λ, T(x) = x, A = eη, h(x) = 1/x!. Hence the log link for counts.
Sigmoid and softmax are not design choices. They are the mean maps of the Bernoulli and the categorical. A generalized linear model sets η = wTx. The canonical link is the inverse of ∇A. The gradient of the negative log likelihood is then always (prediction − target) · x.
Sufficient statistics
A statistic T(X) is sufficient for θ when the data hold no extra information about θ once you know T. Formally, P(X | T(X), θ) does not depend on θ.
Fisher-Neyman factorization. T is sufficient if and only if p(x | θ) = g(T(x), θ) h(x).
For n Bernoulli draws, the count of ones is sufficient.
For a Gaussian with unknown mean and variance, (∑xi, ∑xi2) is sufficient.
The Pitman-Koopman-Darmois theorem says the converse. Among families with fixed support, only exponential families have sufficient statistics of fixed size as n grows.
The practical meaning is compression. You can fold a stream of a billion events into a handful of running sums. You can still fit the model exactly.
Conjugate priors and maximum entropy
Every exponential family has a conjugate prior of the form p(η) ∝ exp(ηTχ − νA(η)). The posterior just adds data statistics to χ and counts to ν. Beta-Bernoulli and Dirichlet-categorical are the famous cases.
There is also a dual view. Fix the expected values of T(x). The maximum entropy distribution that meets them is exactly the exponential family with those statistics.
Why it matters in practice
Choosing the loss. Counts get Poisson loss with a log link. Binary labels get log loss. Positive skewed targets get gamma or Tweedie loss.
Streaming and federated stats. Sufficient statistics are sums. They merge across shards and devices with no loss.
Bandits. Beta-Bernoulli Thompson sampling updates two counters per arm.
Convex training. Logistic and Poisson regression are convex because A is convex.
Interview check
Why is logistic regression convex? The NLL is A(η) − yη with η = wTx. A is convex, and convexity survives a linear map.
Derive softmax from first principles. Write the categorical in exponential form with logits η. Then A is logsumexp, and the mean ∇A is softmax.
What is sufficient for a Poisson rate from n draws? The total ∑xi. The MLE is that total divided by n.
Markov chains and PageRank
Markov chain. A random process where the next state depends only on the current state. The past matters only through the present.
Transition matrix and stationary distribution
P(Xt+1 = j | Xt = i, Xt−1, …) = P(Xt+1 = j | Xt = i) = Pij
Each row of P sums to 1.
Distribution after t steps: πt = π0 Pt
Stationary distribution: π P = π, ∑i πi = 1
So π is a left eigenvector of P with eigenvalue 1. A row-stochastic matrix always has one.
A three-state chain. Solving πP = π gives π = (2, 4, 3) / 9. B gets the most mass because both A and C feed it.
When the stationary distribution is unique and reached
Irreducible. Every state can reach every other state. Then a finite chain has exactly one stationary π.
Aperiodic. The chain does not cycle on a fixed clock. A self-loop anywhere in an irreducible chain is enough.
Ergodic theorem. If a finite chain is irreducible and aperiodic, then π0Pt → π from any start. Time averages also converge to π-averages.
Mixing speed. The error shrinks like |λ2|t, where λ2 is the second largest eigenvalue in modulus. A big spectral gap 1 − |λ2| means fast mixing.
Detailed balance and MCMC
If πiPij = πjPji for all i and j, then π is stationary. Sum both sides over i to see it. Such chains are called reversible.
Metropolis-Hastings builds a chain with detailed balance for any target p known up to a constant. Propose x′ from q(x′ | x). Accept with probability min(1, p(x′)q(x | x′) / (p(x)q(x′ | x))). The unknown normalizer cancels in the ratio. That is why MCMC works for Bayesian posteriors.
PageRank is a stationary distribution
Model a surfer who follows a random out-link with probability d. Otherwise, the surfer jumps to a random page. The rank of a page is the long-run share of time spent there.
G = d P + (1 − d) (1/n) 1 1T d ≈ 0.85
PageRank = the π with π G = π
The teleport term makes every entry of G positive. So the chain is irreducible and aperiodic, and π is unique. The second eigenvalue of G is at most d. So power iteration converges like dt. About 50 steps give 4 digits at d = 0.85. Pages with no out-links are set to jump uniformly.
import numpy as np
def pagerank(adj: np.ndarray, d: float = 0.85, tol: float = 1e-10) -> np.ndarray:
"""Power iteration for PageRank. adj[i, j] = 1 if page i links to page j."""
n = adj.shape[0]
out = adj.sum(axis=1, keepdims=True)
# Dangling pages (no out-links) jump uniformly.
P = np.where(out > 0, adj / np.maximum(out, 1), 1.0 / n)
G = d * P + (1 - d) / n # row-stochastic "Google matrix"
pi = np.full(n, 1.0 / n)
while True:
new = pi @ G # one step of the chain
if np.abs(new - pi).sum() < tol:
return new
pi = new
adj = np.array([[0, 1, 1, 0],
[0, 0, 1, 0],
[1, 0, 0, 0],
[0, 0, 1, 0]], dtype=float)
pi = pagerank(adj)
print(np.round(pi, 4), round(pi.sum(), 6))
# [0.3725 0.1958 0.3941 0.0375] 1.0
Page 3 has no in-links, so it keeps only the teleport mass (1 − d)/n = 0.0375. Page 2 has three in-links and wins.
Why it matters in practice
Graph ranking. PageRank and personalized PageRank score nodes in social and item graphs. Personalized PageRank teleports to the user’s own nodes.
Random-walk embeddings. DeepWalk and node2vec train on walks drawn from a Markov chain.
Bayesian inference. MCMC samples posteriors when no closed form exists.
User journeys. Session funnels and churn are often modeled as chains with absorbing states.
Language models are high-order chains. An n-gram model is a Markov chain on the last n − 1 tokens.
Interview check
Why does PageRank need the teleport term? Without it the web graph has sinks and cycles. The chain would be reducible or periodic, and π could be non-unique or never reached.
How fast does power iteration converge? Like |λ2|t. For the Google matrix |λ2| ≤ d.
Why does Metropolis-Hastings not need the normalizing constant? The acceptance rule uses only the ratio p(x′)/p(x). The constant cancels.
Entropy, cross-entropy and KL divergence
Entropy. The average surprise of a random outcome. It is also the fewest bits per symbol any code can use, on average, to send draws from that distribution.
The last step holds because H(p̂) does not depend on θ. So three views agree. Maximize likelihood. Minimize cross-entropy. Minimize forward KL from data to model.
For classification the target is one-hot y and the model gives softmax probabilities s = softmax(z).
L(z, y) = − ∑k yk log sk = − zc + log ∑k ezk
∂L/∂zk = sk − yk
The gradient is prediction minus target, the exponential family pattern again. Compute it from logits with logsumexp. Never take log of a softmax output, since it underflows.
import numpy as np
def kl(p: np.ndarray, q: np.ndarray) -> float:
"""KL(p || q) in nats. Assumes q > 0 wherever p > 0."""
m = p > 0
return float(np.sum(p[m] * np.log(p[m] / q[m])))
def entropy(p: np.ndarray) -> float:
m = p > 0
return float(-np.sum(p[m] * np.log(p[m])))
p = np.array([0.70, 0.20, 0.10])
q = np.array([0.40, 0.40, 0.20])
cross = float(-np.sum(p * np.log(q)))
print(f"H(p)={entropy(p):.4f} H(p,q)={cross:.4f} KL(p||q)={kl(p, q):.4f}")
print(f"H(p,q) - H(p) = {cross - entropy(p):.4f}")
print(f"KL(q||p)={kl(q, p):.4f} (not equal: KL is asymmetric)")
# H(p)=0.8018 H(p,q)=0.9856 KL(p||q)=0.1838
# H(p,q) - H(p) = 0.1838
# KL(q||p)=0.1920 (not equal: KL is asymmetric)
Forward versus reverse KL
Fit a simple q to a complex p. Which direction you minimize changes the answer a lot.
Forward KL(p‖q) = Ep[log p/q]. It explodes wherever p > 0 but q ≈ 0. So q must cover all of p. This is mass-covering or mean-seeking. MLE minimizes this direction. For a Gaussian q, the optimum matches the mean and covariance of p.
Reverse KL(q‖p) = Eq[log q/p]. It explodes wherever q > 0 but p ≈ 0. So q hides inside one region of high p. This is mode-seeking or zero-forcing. Variational inference minimizes this direction.
Forward KL spreads q across both modes and puts mass where p has almost none. Reverse KL locks onto one mode and ignores the other.
log p(x) does not depend on q. So raising the ELBO lowers the reverse KL to the true posterior. That is why mean-field VI tends to underestimate posterior variance.
Why it matters in practice
Every classifier and LLM trains on cross-entropy. Perplexity is exp of the per-token cross-entropy.
Distillation minimizes KL from a teacher’s soft labels to the student. Temperature T softens both, and the loss is scaled by T2.
Label smoothing is cross-entropy against a mix of the one-hot label and uniform. It caps over-confidence.
RLHF and PPO add a KL penalty to the reference model. It stops the policy drifting into reward hacks.
Drift monitoring compares feature histograms with KL or its symmetric cousin, PSI.
Interview check
Show that minimizing cross-entropy is MLE. Average log likelihood equals −H(p̂, qθ). H(p̂) is fixed, so this is also minimizing KL(p̂‖qθ).
Which KL does VI minimize, and what is the side effect? Reverse KL(q‖p). It is mode-seeking, so q underestimates spread and can miss modes.
Why is the softmax cross-entropy gradient so clean? Logsumexp is the log-partition of the categorical. Its gradient is softmax, so ∂L/∂z = s − y.
Mutual information
Mutual information. How much knowing one variable cuts your uncertainty about another. It is zero exactly when the two are independent.
Symmetric and non-negative. I(X;Y) = I(Y;X) ≥ 0, with zero only for independence.
Catches non-linear links. For Y = X2 with symmetric X, correlation is zero but MI is large.
Invariant to reparameterization. Any invertible map of X or Y leaves I unchanged.
Data processing inequality. If X → Y → Z is a Markov chain, then I(X;Z) ≤ I(X;Y). Processing cannot create information.
Gaussian case. For a jointly Gaussian pair with correlation ρ, I = −½ log(1 − ρ2).
Feature selection
Information gain in a decision tree is the mutual information between the split and the label. Filter methods rank features by I(Xj; Y). mRMR adds a penalty for redundancy. It picks features with high I(Xj; Y) and low average I(Xj; Xselected).
Estimating MI is the hard part. Binning works in one or two dimensions. The KSG estimator uses nearest neighbours for continuous data. In high dimensions, all estimators have high bias or variance.
InfoNCE and contrastive learning
Take a batch of K pairs (xi, yi) drawn from p(x, y). For each xi, the matching yi is the positive. The other K − 1 are negatives. Score pairs with a critic f, often a scaled cosine similarity.
The loss is cross-entropy for picking the positive out of K. A perfect critic drives the loss to zero. Then the bound reads I ≥ log K. So the bound can never exceed log K. That is one reason contrastive methods want huge batches or memory queues.
Why it matters in practice
CLIP, SimCLR and two-tower retrieval train with InfoNCE. In-batch negatives make it cheap.
Feature screening. MI flags non-linear signal that a correlation screen would drop.
Clustering evaluation. Normalized MI and adjusted MI compare a clustering to labels.
Privacy and fairness. MI between a representation and a sensitive attribute measures leakage.
Active learning. BALD picks points with the highest MI between the label and model weights.
Interview check
Why can InfoNCE not certify more than log K nats? The loss is at least zero, so the bound log K − L tops out at log K.
How does information gain in a tree relate to MI? It is exactly I(split; label) on the node’s data.
Correlation is 0.02 but MI is high. What does that tell you? The link is strong but non-linear or non-monotone. A tree or a spline will use it.
Jensen-Shannon and Wasserstein
Divergences between distributions. KL can be infinite and is lopsided. Jensen-Shannon fixes the symmetry. Wasserstein measures how far mass must move, so it stays useful when supports do not overlap.
Jensen-Shannon divergence
JS(p, q) = ½ KL(p ‖ m) + ½ KL(q ‖ m), m = ½(p + q)
0 ≤ JS ≤ log 2, symmetric, √JS is a true metric
Because m covers both p and q, JS is always finite. The original GAN links to it directly. With the optimal discriminator, the generator minimizes 2 · JS(pdata, pg) − log 4.
Here is the catch. If pdata and pg sit on disjoint low-dimensional manifolds, JS equals log 2 no matter how close they are. The gradient is zero. Early GAN training often failed this way.
Π(p, q) is the set of joint plans whose marginals are p and q. Think of p as piles of dirt and q as holes. W1 is the least total work to fill the holes.
One example shows the difference
Let p be a point mass at 0 and q a point mass at θ.
KL(p‖q) is ∞ for every θ ≠ 0.
JS(p, q) is log 2 for every θ ≠ 0. It is flat, so no gradient.
W1(p, q) is |θ|. It shrinks smoothly as q moves toward p.
WGAN uses the dual form. The critic f must be 1-Lipschitz. Weight clipping enforced this at first. A gradient penalty on ‖∇f‖ works better.
Why it matters in practice
GAN training. WGAN-GP gives a loss that tracks sample quality and rarely collapses.
FID. The standard image metric is the squared W2 between Gaussians fit to features. It equals ‖μ1 − μ2‖2 + Tr(Σ1 + Σ2 − 2(Σ1Σ2)1/2).
Drift detection. 1D Wasserstein on a feature or score is cheap and in the units of the feature. JS is bounded, which makes alert thresholds easy.
Domain adaptation. Optimal transport aligns source and target feature distributions.
Interview check
Why did WGAN help? JS is constant when supports do not overlap, so gradients vanish. W1 still shrinks as the generator moves closer.
What is the upper bound on JS? log 2, reached when p and q have disjoint supports.
How do you compute W1 between two 1D samples? Sort both. With equal sizes, average |x(i) − y(i)|. This is the area between the empirical CDFs.
Concentration and sketching
Sketch. A small summary of a huge stream that answers one kind of query with a provable error bound. Hashing makes the summary random. Concentration inequalities bound the error.
Count-min sketch: frequencies
Keep a d × w table of counters and d independent hash functions h1, …, hd. To add item x, increment cell (r, hr(x)) in every row r. To query x, take the minimum over the d cells.
w = ⌈ e / ε ⌉, d = ⌈ ln(1/δ) ⌉
true(x) ≤ estimate(x) ≤ true(x) + ε N with probability ≥ 1 − δ
The proof uses only Markov’s inequality. Fix a row. Other items that collide with x add an overcount Z ≥ 0. Each other item lands in x’s cell with probability 1/w. So E[Z] ≤ N/w = εN/e. Markov gives P(Z ≥ εN) ≤ 1/e. The minimum fails only if all d independent rows fail. That has probability at most e−d = δ.
Memory is O((1/ε) log(1/δ)), independent of the number of distinct items. The sketch never undercounts. It works best for heavy hitters, where εN is small next to the true count.
Hash each item to a uniform bit string. The chance that a hash starts with k zeros is 2−(k+1). Say the longest run of leading zeros so far is k. Then you have likely seen about 2k distinct items. Repeats do not matter, since the same item always gives the same hash.
One such estimate is very noisy. HyperLogLog splits the hash. The first b bits pick one of m = 2b registers. Each register keeps its own max leading-zero count. The estimate combines registers with a harmonic mean and a bias constant.
Error. The relative standard error is about 1.04 / √m. With m = 214 it is 0.81%.
Memory. Each register needs about 6 bits. So 214 registers take about 12 KB, for billions of distinct items.
Mergeable. The union of two sketches is the element-wise max of registers. So you can count distinct users across shards and days.
Two more you should know
Bloom filter. Set membership with no false negatives. With m bits, n items and k hashes, the false positive rate is about (1 − e−kn/m)k. The best k is (m/n) ln 2. About 9.6 bits per item give 1% false positives.
Johnson-Lindenstrauss. Take n points and project them at random to k = O(log n / ε2) dimensions. All pairwise distances stay within a factor 1 ± ε. The proof is a concentration bound on a sum of squared Gaussians.
The best q is proportional to |f(x)| p(x). So put samples where f p is large, not just where p is large. q needs heavier tails than p. Otherwise rare huge weights make the variance blow up. A tiny ESS is the warning sign.
import numpy as np
rng = np.random.default_rng(1)
n = 100_000
# Goal: P(X > 4) for X ~ N(0, 1). True value is about 3.167e-5.
hits = rng.standard_normal(n) > 4
naive = hits.mean()
# Importance sampling: draw from q = N(4, 1), reweight by w = p(x) / q(x).
x = rng.normal(4.0, 1.0, n)
w = np.exp(-4.0 * x + 8.0) # p(x) / q(x) for these two Gaussians
terms = (x > 4) * w
est = terms.mean()
print(f"naive : {naive:.3e} ({hits.sum()} hits)")
print(f"importance : {est:.3e} rel. std err {terms.std() / np.sqrt(n) / est:.2%}")
# naive : 2.000e-05 (2 hits)
# importance : 3.162e-05 rel. std err 0.67%
The naive estimate rests on two hits. Shifting the proposal into the tail gives a 0.67% relative error with the same budget.
The Gumbel-max trick
Draw Gk = −log(−log Uk) independently for each class. Then:
Why it works. The max of Gumbels shifted by log πk is again Gumbel. The chance that index k wins is exactly πk / ∑j πj. You never normalize, so it works straight from logits.
Gumbel-top-k. Take the top k of the perturbed logits. You get k draws without replacement.
Gumbel-softmax. Replace argmax with softmax((z + G)/τ). This gives a differentiable relaxed sample. As τ → 0 it approaches one-hot.
import numpy as np
rng = np.random.default_rng(0)
logits = np.array([2.0, 1.0, 0.0, -1.0])
probs = np.exp(logits - logits.max())
probs /= probs.sum()
# Gumbel-max: argmax(logits + Gumbel noise) is an exact draw from softmax(logits).
n = 200_000
g = -np.log(-np.log(rng.uniform(size=(n, logits.size))))
draws = np.argmax(logits + g, axis=1)
freq = np.bincount(draws, minlength=logits.size) / n
print("softmax :", np.round(probs, 4))
print("gumbel :", np.round(freq, 4))
# softmax : [0.6439 0.2369 0.0871 0.0321]
# gumbel : [0.6436 0.2374 0.0872 0.0318]
Why it matters in practice
Off-policy evaluation. Inverse propensity scoring is importance sampling. It estimates a new ranker’s value from logs: (1/n) ∑ ri π(ai|xi) / μ(ai|xi). Clip or self-normalize the weights to tame variance.
PPO. The ratio πθ/πold in its objective is an importance weight. Clipping keeps it near 1.
LLM decoding. Temperature and top-k sampling are Gumbel-max on scaled or truncated logits.
Negative sampling. Sampled softmax draws negatives from a proposal and corrects with importance weights, often log q.
Rare-event metrics. Fraud or outage rates need importance sampling or stratified oversampling to estimate well.
Interview check
Sample from an exponential with only a uniform RNG. X = −ln(1 − U)/λ, by the inverse CDF.
Your IPS estimate has huge variance. What do you do? Clip weights, self-normalize, or use a doubly robust estimator. Check the ESS first.
Why is Gumbel-max useful when softmax already exists? It samples straight from logits with no normalization. It extends to top-k without replacement and to a differentiable relaxation.
Recap
Independence and conditional independence are different. Conditioning on a common effect creates dependence.
Var(X) = E[Var(X|Y)] + Var(E[X|Y]). Averages of correlated parts keep a ρσ2 floor.
Markov and Chebyshev need little. Hoeffding gives exponential tails for bounded data. Jensen gives KL ≥ 0 and the ELBO.
The strong law needs a finite mean. The CLT needs finite variance. Berry-Esseen and the delta method make it usable.
Gaussian conditionals are linear regression with a Schur complement covariance. Precision zeros mean conditional independence.
Exponential families have convex log-partitions. Sigmoid and softmax are their mean maps. MLE is moment matching.
PageRank is the stationary distribution of a teleporting random walk. Power iteration converges like dt.
Cross-entropy loss is MLE is forward KL. Forward KL covers mass. Reverse KL seeks modes.
InfoNCE lower-bounds MI but caps at log K. W1 keeps gradients when JS goes flat.
Count-min, HyperLogLog and Bloom filters trade a bounded error for tiny memory. Importance sampling and Gumbel-max power OPE and decoding.