Statistics & classic ML, for inference-performance work

Benchmarking statistics first, then the linear algebra, optimization, and classic-ML core that interviews and day-to-day perf work lean on. Every number on this page came from the code beside it.

1. Probability essentials

A random variable is a number produced by a process you cannot predict exactly; a distribution says how likely each value is. Everything in benchmarking is about estimating properties of a distribution (its mean, its 99th percentile) from a finite sample of it.

Expectation, variance, independence

• E[X] is the probability-weighted average. It is linear: E[aX + bY] = aE[X] + bE[Y], always, even for dependent variables.
• Var(X) = E[(X-E[X])²] = E[X²] - E[X]². Variance adds only for independent variables: Var(X+Y) = Var(X)+Var(Y) + 2Cov(X,Y). Correlated errors do not average out, which is why "run it 10 times on the same noisy host" helps less than you hope.
• Independent means knowing one tells you nothing about the other: P(A and B) = P(A)P(B). Back-to-back requests on one GPU are usually not independent (shared queue, thermal state, cache warmth).
• Conditional probability P(A|B) = P(A and B)/P(B). Bayes: P(A|B) = P(B|A)P(A)/P(B). Base rates dominate when the event is rare.

Demo: distributions and Bayes

Copyimport numpy as np
rng = np.random.default_rng(42)
# lognormal latency: median 40ms, sigma 0.6
x = rng.lognormal(mean=np.log(40), sigma=0.6, size=1_000_000)
print(f"mean={x.mean():.1f}  p50={np.percentile(x,50):.1f}  p90={np.percentile(x,90):.1f}  p99={np.percentile(x,99):.1f}  p99.9={np.percentile(x,99.9):.1f}")
print(f"theory mean = 40*exp(0.6^2/2) = {40*np.exp(0.18):.1f}")
mu, sd = x.mean(), x.std()
print(f"normal fit would say p99 = mean+2.326*sd = {mu+2.326*sd:.1f}; P(x<0) under that fit = {__import__('scipy.stats').stats.norm.cdf(0,mu,sd):.3f}")
print(f"skew={__import__('scipy.stats').stats.skew(x):.2f}")
# Poisson / exponential
lam=50
k = rng.poisson(lam, 100000)
print(f"Poisson(50): mean={k.mean():.2f} var={k.var():.2f}")
gaps = rng.exponential(1/lam, 100000)
print(f"Exp gaps: mean={gaps.mean()*1000:.2f} ms, P(gap>40ms)={np.mean(gaps>0.04):.4f} theory={np.exp(-lam*0.04):.4f}")
# memoryless: P(gap>60 | gap>40) == P(gap>20)
print(f"memoryless: P(gap>60|gap>40)={np.mean(gaps>0.06)/np.mean(gaps>0.04):.4f}  P(gap>20ms)={np.mean(gaps>0.02):.4f}")
# Bayes
p_bad, tpr, fpr = 0.01, 0.95, 0.05
post = tpr*p_bad/(tpr*p_bad+fpr*(1-p_bad))
print(f"Bayes: P(regression | alert) = {post:.4f}")
# Pareto heavy tail: sample mean unstable
pa = rng.pareto(1.1, 1_000_000)+1
for n in (1_000, 100_000, 1_000_000):
    print(f"Pareto(a=1.1) mean of first {n:>9,}: {pa[:n].mean():.2f}")
# variance of sum of independent vs perfectly correlated
a = rng.normal(0,1,100000); b = rng.normal(0,1,100000)
print(f"Var(a+b) indep={np.var(a+b):.2f}  Var(a+a)={np.var(a+a):.2f}")

mean=47.9  p50=40.0  p90=86.3  p99=161.7  p99.9=255.4
theory mean = 40*exp(0.6^2/2) = 47.9
normal fit would say p99 = mean+2.326*sd = 121.4; P(x<0) under that fit = 0.065
skew=2.30
Poisson(50): mean=50.03 var=49.78
Exp gaps: mean=19.98 ms, P(gap>40ms)=0.1344 theory=0.1353
memoryless: P(gap>60|gap>40)=0.3666  P(gap>20ms)=0.3675
Bayes: P(regression | alert) = 0.1610
Pareto(a=1.1) mean of first     1,000: 5.55
Pareto(a=1.1) mean of first   100,000: 7.37
Pareto(a=1.1) mean of first 1,000,000: 10.47
Var(a+b) indep=1.99  Var(a+a)=3.96

2. Latency statistics: tails, fan-out, queues

Why latency is not normal

Latency has a hard floor (the work must be done) and no ceiling (queueing, GC, a cold cache, a retry, a noisy neighbor each multiply the time). Multiplicative effects make a lognormal-ish, right-skewed shape; mixtures of fast and slow paths (cache hit vs miss, batch-full vs timeout) make it multimodal. Consequences: the mean is dragged up by the tail, the standard deviation is not a useful "±" band, and a normal-based "mean + 3σ" misstates the tail (section 1 output).

Percentiles vs the mean

pXX is the value below which XX% of requests fall. SLOs are written on percentiles because users experience individual requests, and the mean hides the bad ones. The mean is still useful for capacity math (throughput = concurrency / mean latency, Little's law below) because totals add and percentiles do not.

You cannot average percentiles

The p99 of a union is not a function of the p99s of the parts: it depends on how many requests each part had and the whole shape of each distribution. Averaging per-host or per-minute p99s is a very common dashboard bug. Correct options: keep the raw samples, or keep mergeable histograms (Prometheus histogram buckets, HdrHistogram, t-digest/DDSketch) and compute the percentile after merging. Even then, note that Prometheus histogram_quantile interpolates inside a bucket, so its accuracy is bounded by bucket width.

Tail amplification under fan-out

If one call exceeds its own p99 with probability p = 0.01, a request that waits on N independent calls is slowed by at least one of them with probability 1 - (1-p)^N. The request latency is the max of the N calls, so the request's median sits near the individual calls' high percentiles. This is exactly the situation when a ranking request with C candidates is split into ceil(C/100) Triton calls of up to 100 candidates each: at 1000 candidates the request makes 10 calls and 9.6% of requests see a "p99-slow" call.

Little's law and queueing

L = λW: the average number in the system equals the arrival rate times the average time each spends in it. It needs no distributional assumptions, only a stable system. Uses: 200 req/s at 45 ms means 9 requests in flight on average, so you need at least that much concurrency headroom; also how many in-flight requests a pod needs, how big a queue is, how big a connection pool must be.

M/M/1 (Poisson arrivals, exponential service, one server) has utilization ρ = λ/μ and mean time in system W = 1/(μ-λ) = (1/μ)/(1-ρ). The 1/(1-ρ) factor is the hockey stick: the closer you run to 100%, the faster latency explodes. In M/M/1 the time in system is itself exponential with rate (μ-λ), so p99 = ln(100)/(μ-λ) = 4.6 × the mean. Real servers (batching, multiple workers, non-exponential service) differ in constants but not in the shape.

Coordinated omission

A closed-loop load generator sends the next request only after the previous one returns. When the server stalls, the client stops sending, so it records one slow request instead of the hundreds that would have arrived during the stall in a real open-loop workload. The percentiles then look great. Fix: schedule requests on an intended timeline (open loop) and measure latency from the intended send time, as wrk2 and HdrHistogram's correction do.

Demo

Copyimport numpy as np
rng = np.random.default_rng(42)
# 1) averaging percentiles
fast = rng.lognormal(np.log(30), 0.3, 900_000)   # 900k reqs
slow = rng.lognormal(np.log(80), 0.8, 100_000)   # 100k reqs
allx = np.concatenate([fast, slow])
p99f, p99s = np.percentile(fast,99), np.percentile(slow,99)
print(f"host A p99={p99f:.1f} (900k reqs)  host B p99={p99s:.1f} (100k reqs)")
print(f"mean of p99s      = {(p99f+p99s)/2:.1f}")
print(f"traffic-weighted  = {0.9*p99f+0.1*p99s:.1f}")
print(f"true pooled p99   = {np.percentile(allx,99):.1f}")
# 2) fan-out
p=0.01; 
print("\nfan-out: P(at least one of N calls exceeds its own p99)")
for N in (1,5,10,30,100):
    sim = (rng.random((200_000,N))<p).any(axis=1).mean()
    print(f"N={N:>3}  formula={1-(1-p)**N:.4f}  sim={sim:.4f}")
import math
for c in (100,500,1000,3000):
    N = math.ceil(c/100)
    print(f"candidates={c:>5} -> N={N:>2} Triton calls, P(request hits a per-call-p99 straggler)={1-0.99**N:.3f}")
# request latency = max of N calls
N=30; calls = rng.lognormal(np.log(10),0.5,(200_000,N)); r = calls.max(axis=1)
print(f"\nN=30 calls, single-call p50={np.percentile(calls[:,0],50):.1f} p99={np.percentile(calls[:,0],99):.1f}; request(max) p50={np.percentile(r,50):.1f} p99={np.percentile(r,99):.1f}")
print(f"request p50 lands at single-call percentile {np.mean(calls[:,0]<np.percentile(r,50))*100:.1f}")
# 3) Little
lam, W = 200, 0.045
print(f"\nLittle: lambda=200 req/s, W=45ms -> L={lam*W:.1f} in flight")
# 4) M/M/1 sim via Lindley
def mm1(lam, mu, n=400_000):
    ia = rng.exponential(1/lam,n); sv = rng.exponential(1/mu,n)
    w = np.zeros(n)
    for i in range(1,n): w[i] = max(0, w[i-1]+sv[i-1]-ia[i])
    return (w+sv)[n//10:]
mu=100.0
print("\nM/M/1 mu=100/s: mean sojourn (ms) sim vs 1/(mu-lam), plus p99 sim vs ln(100)/(mu-lam)")
for rho in (0.5,0.7,0.8,0.9,0.95):
    lam=rho*mu; t=mm1(lam,mu)
    print(f"rho={rho:.2f}  sim={t.mean()*1000:7.1f}  formula={1000/(mu-lam):7.1f}   p99 sim={np.percentile(t,99)*1000:7.1f} formula={1000*np.log(100)/(mu-lam):7.1f}")
# 5) coordinated omission: intended rate 1 req/10ms for 60s, server takes 5ms, but stalls 1s once at t>=10s.
# Closed-loop client: one in flight, sends late if the previous response was late, and times from the actual send.
def run(correct):
    t=0.0; sched=0.0; lat=[]; stalled=False
    while sched<60.0:
        send=max(t,sched)
        svc=0.005
        if not stalled and send>=10.0: svc=1.0; stalled=True
        t=send+svc
        lat.append(t-(sched if correct else send))
        sched+=0.01
    return np.array(lat)*1000
naive=run(False); corr=run(True)
for nm,l in (("naive (time from send)",naive),("corrected (time from intended)",corr)):
    print(f"{nm:31s} n={len(l):>5} p50={np.percentile(l,50):6.1f} p99={np.percentile(l,99):7.1f} p99.9={np.percentile(l,99.9):7.1f} max={l.max():7.1f}  #>100ms={int((l>100).sum())}")

host A p99=60.3 (900k reqs)  host B p99=514.0 (100k reqs)
mean of p99s      = 287.2
traffic-weighted  = 105.7
true pooled p99   = 224.8

fan-out: P(at least one of N calls exceeds its own p99)
N=  1  formula=0.0100  sim=0.0097
N=  5  formula=0.0490  sim=0.0484
N= 10  formula=0.0956  sim=0.0968
N= 30  formula=0.2603  sim=0.2603
N=100  formula=0.6340  sim=0.6333
candidates=  100 -> N= 1 Triton calls, P(request hits a per-call-p99 straggler)=0.010
candidates=  500 -> N= 5 Triton calls, P(request hits a per-call-p99 straggler)=0.049
candidates= 1000 -> N=10 Triton calls, P(request hits a per-call-p99 straggler)=0.096
candidates= 3000 -> N=30 Triton calls, P(request hits a per-call-p99 straggler)=0.260

N=30 calls, single-call p50=10.0 p99=32.0; request(max) p50=27.2 p99=54.6
request p50 lands at single-call percentile 97.7

Little: lambda=200 req/s, W=45ms -> L=9.0 in flight

M/M/1 mu=100/s: mean sojourn (ms) sim vs 1/(mu-lam), plus p99 sim vs ln(100)/(mu-lam)
rho=0.50  sim=   20.0  formula=   20.0   p99 sim=   92.2 formula=   92.1
rho=0.70  sim=   33.9  formula=   33.3   p99 sim=  155.7 formula=  153.5
rho=0.80  sim=   49.7  formula=   50.0   p99 sim=  229.3 formula=  230.3
rho=0.90  sim=   99.1  formula=  100.0   p99 sim=  432.8 formula=  460.5
rho=0.95  sim=  195.1  formula=  200.0   p99 sim=  780.6 formula=  921.0
naive (time from send)          n= 6001 p50=   5.0 p99=    5.0 p99.9=    5.0 max= 1000.0  #>100ms=1
corrected (time from intended)  n= 6001 p50=   5.0 p99=  700.0 p99.9=  970.0 max= 1000.0  #>100ms=181

3. Estimation & uncertainty

A benchmark gives you a sample; you want a fact about the population (all the requests this config would ever serve). Any statistic you compute (mean, p99) is itself a random variable: rerun and you get a different number. The spread of that number across hypothetical reruns is its standard error.

• CLT: the mean of n independent draws is approximately normal with sd = σ/√n, almost regardless of the underlying shape, once n is moderate. SE = s/√n.
• 95% CI for a mean: x̄ ± 1.96 · SE. Meaning: if you repeated the experiment many times, ~95% of intervals built this way would contain the true mean. It is not "95% probability the truth is in this interval".
• Precision costs √n: halving the CI width needs 4x the samples. Required n for a CI half-width h: n = (1.96σ/h)².
• Percentiles have no neat formula, and the CLT does not apply to the mean of a heavy tail the same way. Use the bootstrap: resample your n observations with replacement B times, recompute p99 each time, take the 2.5th and 97.5th percentiles of those B values. It needs no distributional assumption.
• The p99 of n=100 samples is essentially the 1st or 2nd largest value. With ~1 observation in the tail, nothing can estimate it well. Rule of thumb: you want at least ~10 observations beyond the percentile you report: n≥1,000 for p99, n≥10,000 for p99.9. The bootstrap cannot manufacture information that is not in the sample; it fails for extremes like the max.

Copyimport numpy as np
rng = np.random.default_rng(42)
pop = rng.lognormal(np.log(40),0.6,2_000_000)
true99 = np.percentile(pop,99); print(f"true population p99 = {true99:.1f} ms (theory {40*np.exp(0.6*2.3263):.1f})")
def boot(x,B=2000):
    idx = rng.integers(0,len(x),(B,len(x)))
    return np.percentile(np.percentile(x[idx],99,axis=1),[2.5,97.5])
for n in (100,1000,10000):
    x = rng.choice(pop,n)
    lo,hi = boot(x, 2000 if n<10000 else 500)
    print(f"n={n:>6} sample p99={np.percentile(x,99):6.1f}  95% bootstrap CI=[{lo:6.1f}, {hi:6.1f}] width={hi-lo:6.1f}")
# coverage of the mean CI via CLT, and spread of p99 estimator
ests = [np.percentile(rng.choice(pop,100),99) for _ in range(2000)]
print(f"n=100: p99 estimate across 2000 repeats: mean={np.mean(ests):.1f} sd={np.std(ests):.1f} (biased low: p99 of 100 samples is ~ the 2nd largest)")
# mean SE
x=rng.choice(pop,1000); se=x.std(ddof=1)/np.sqrt(1000)
print(f"n=1000 mean={x.mean():.2f} SE={se:.2f} 95% CI=[{x.mean()-1.96*se:.2f},{x.mean()+1.96*se:.2f}]  true mean={pop.mean():.2f}")
# sample size for mean to +-1ms
sd=pop.std(); print(f"pop sd={sd:.1f}; n for +-1ms 95% on mean = (1.96*sd/1)^2 = {(1.96*sd)**2:,.0f}")
# how many samples needed so ~10 obs sit above p99
print("rule of thumb: need >= ~1000 samples for p99 to have 10 obs in tail; >= 10000 for p99.9")
# bootstrap failing on max
x=rng.choice(pop,100); m=x.max(); bm=x[rng.integers(0,100,(2000,100))].max(axis=1)
print(f"bootstrap of max: P(boot max == sample max)={np.mean(bm==m):.3f} (theory 1-(1-1/n)^n=0.634)")

true population p99 = 161.7 ms (theory 161.5)
n=   100 sample p99= 155.4  95% bootstrap CI=[ 121.6,  192.1] width=  70.5
n=  1000 sample p99= 145.6  95% bootstrap CI=[ 131.8,  160.4] width=  28.6
n= 10000 sample p99= 162.8  95% bootstrap CI=[ 155.4,  168.9] width=  13.5
n=100: p99 estimate across 2000 repeats: mean=149.4 sd=30.2 (biased low: p99 of 100 samples is ~ the 2nd largest)
n=1000 mean=47.80 SE=0.98 95% CI=[45.88,49.71]  true mean=47.90
pop sd=31.6; n for +-1ms 95% on mean = (1.96*sd/1)^2 = 3,826
rule of thumb: need >= ~1000 samples for p99 to have 10 obs in tail; >= 10000 for p99.9
bootstrap of max: P(boot max == sample max)=0.638 (theory 1-(1-1/n)^n=0.634)

4. Hypothesis testing & A/B for performance

A test asks: if there were no real difference, how surprising is what I saw? The p-value is P(data at least this extreme | no difference). It is not P(no difference | data), not the probability the result is a fluke, and says nothing about size. Alpha (usually 0.05) is the false-positive rate you accept. Power is P(detect a real effect of a given size) (usually target 80%). Effect size is how big the difference is relative to noise (Cohen's d = diff/sd; or simply "ms of p99").

Multiple comparisons: test 20 metrics at 0.05 and you expect one false positive; P(at least one) = 1-0.95²⁰ = 64%. Pre-register the primary metric, or correct (Bonferroni: divide alpha by the number of tests; Benjamini-Hochberg: control false discovery rate). Peeking (checking significance repeatedly and stopping when p<0.05) inflates false positives the same way.

Copyimport numpy as np
from scipy import stats
rng = np.random.default_rng(42)
# A: baseline, B: config where 90% of reqs the same but median shifts -3ms, plus rare tail changes
def lat(n, shift=1.0, tail=0.0):
    x = rng.lognormal(np.log(40*shift),0.6,n)
    k = rng.random(n)<tail; x[k]*=8
    return x
A=lat(500); B=lat(500,shift=0.93)
print(f"A mean={A.mean():.1f} med={np.median(A):.1f} | B mean={B.mean():.1f} med={np.median(B):.1f}")
print("Welch t p =", round(stats.ttest_ind(A,B,equal_var=False).pvalue,3), "| Mann-Whitney p =", round(stats.mannwhitneyu(A,B).pvalue,3))
# power sim: skewed data, true median shift 7%
def power(test,n,reps=2000):
    h=0
    for _ in range(reps):
        a=rng.lognormal(np.log(40),0.6,n); b=rng.lognormal(np.log(40*0.93),0.6,n)
        h+= test(a,b)<0.05
    return h/reps
t=lambda a,b: stats.ttest_ind(a,b,equal_var=False).pvalue
m=lambda a,b: stats.mannwhitneyu(a,b).pvalue
lt=lambda a,b: stats.ttest_ind(np.log(a),np.log(b)).pvalue
for n in (100,500,1000):
    print(f"n={n:>4}/arm  power: welch-t={power(t,n):.2f}  mann-whitney={power(m,n):.2f}  t-on-log={power(lt,n):.2f}")
# Type I: no difference
def alpha(test,n=200,reps=3000):
    return np.mean([test(rng.lognormal(0,1,n),rng.lognormal(0,1,n))<0.05 for _ in range(reps)])
print(f"false-positive rate with no real diff (n=200): welch-t={alpha(t):.3f}")
# Mann-Whitney is NOT a median test: same median, different tails
a=rng.lognormal(np.log(40),0.3,2000); b=rng.lognormal(np.log(40),0.9,2000)
print(f"same median ({np.median(a):.1f} vs {np.median(b):.1f}); MW p={stats.mannwhitneyu(a,b).pvalue:.3g}; P(A<B)={stats.mannwhitneyu(a,b).statistic/(2000*2000):.3f}")
# multiple comparisons
print("\nmultiple comparisons: 20 metrics, no true effect")
p=[stats.ttest_ind(rng.normal(0,1,100),rng.normal(0,1,100)).pvalue for _ in range(20)]
print(f"min p={min(p):.3f}; #p<0.05 = {sum(np.array(p)<0.05)}; P(at least one) analytic = {1-0.95**20:.3f}; Bonferroni threshold = {0.05/20:.4f}")
# power analysis: n per arm for d
from math import ceil
for d in (0.2,0.5,0.8):
    n=2*((1.96+0.8416)/d)**2; print(f"effect size d={d}: n/arm for 80% power ~ {ceil(n)}")
print(f"d for 3ms diff with sd=24ms: {3/24:.3f}")
# paired design: noisy machine drift
drift=rng.normal(0,10,30)   # per-run shared noise
A=40+drift+rng.normal(0,2,30); B=37+drift+rng.normal(0,2,30)
print(f"\npaired runs: unpaired p={stats.ttest_ind(A,B).pvalue:.3f}  paired p={stats.ttest_rel(A,B).pvalue:.4f}")
# bootstrap CI for diff in p99
x=rng.lognormal(np.log(40),0.6,5000); y=rng.lognormal(np.log(40),0.6,5000)*1.05
d=[np.percentile(rng.choice(y,5000),99)-np.percentile(rng.choice(x,5000),99) for _ in range(1000)]
print(f"p99 diff (B-A) 95% bootstrap CI: [{np.percentile(d,2.5):.1f},{np.percentile(d,97.5):.1f}] ms")

# tail-only regression: bulk identical, 2% of B requests are 8x slower
a=rng.lognormal(np.log(40),0.6,5000); b=rng.lognormal(np.log(40),0.6,5000); k=rng.random(5000)<0.02; b[k]*=8
print(f"\ntail-only regression: p50 {np.median(a):.1f} vs {np.median(b):.1f}; p99 {np.percentile(a,99):.1f} vs {np.percentile(b,99):.1f}; MW p={stats.mannwhitneyu(a,b).pvalue:.3f}  Welch p={stats.ttest_ind(a,b,equal_var=False).pvalue:.3f}")
# outliers: t loses power, MW does not
def pw2(test,n=100,reps=2000):
    h=0
    for _ in range(reps):
        a=rng.lognormal(np.log(40),0.4,n); b=rng.lognormal(np.log(40*0.9),0.4,n); o=rng.random(n)<0.02; b[o]*=60
        h+=test(a,b)<0.05
    return h/reps
print(f"10% median shift + 2% 60x outliers (n=100/arm): power welch-t={pw2(t):.2f} mann-whitney={pw2(m):.2f}")

A mean=47.0 med=40.1 | B mean=45.2 med=37.8
Welch t p = 0.352 | Mann-Whitney p = 0.254
n= 100/arm  power: welch-t=0.13  mann-whitney=0.13  t-on-log=0.14
n= 500/arm  power: welch-t=0.41  mann-whitney=0.45  t-on-log=0.46
n=1000/arm  power: welch-t=0.69  mann-whitney=0.75  t-on-log=0.79
false-positive rate with no real diff (n=200): welch-t=0.045
same median (40.1 vs 40.2); MW p=0.548; P(A<B)=0.505

multiple comparisons: 20 metrics, no true effect
min p=0.035; #p<0.05 = 2; P(at least one) analytic = 0.642; Bonferroni threshold = 0.0025
effect size d=0.2: n/arm for 80% power ~ 393
effect size d=0.5: n/arm for 80% power ~ 63
effect size d=0.8: n/arm for 80% power ~ 25
d for 3ms diff with sd=24ms: 0.125

paired runs: unpaired p=0.246  paired p=0.0000
p99 diff (B-A) 95% bootstrap CI: [-13.5,18.0] ms

tail-only regression: p50 40.0 vs 40.0; p99 162.3 vs 354.0; MW p=0.613  Welch p=0.000
10% median shift + 2% 60x outliers (n=100/arm): power welch-t=0.11 mann-whitney=0.34

Comparing two serving configs properly

1. Pick the one primary metric and the minimum effect that matters (e.g. p99 −5%) before running. Everything else is secondary/exploratory.
2. Use the same hardware, same model, same request mix (replay the same trace), fresh warm-up discarded, steady state reached.
3. Interleave A and B (ABABAB across the same time window) or run on identical hosts simultaneously; never "A in the morning, B in the afternoon".
4. Use an open-loop load at the production arrival rate, and also sweep load, since configs often cross over (B better at low load, worse near saturation).
5. Pre-compute n from the noise level and target effect; collect enough for the tail you report. Report the effect with a CI (ratio of medians, difference in p99), not a bare p-value.
6. Repeat the whole run several times; run-to-run variance is usually larger than within-run sampling error.

5. Linear algebra for ML

• Matmul shapes: (m×k)·(k×n) = (m×n); inner dimensions must match. FLOPs = 2mkn. A transformer linear layer is activations (tokens×d_in) · weights (d_in×d_out). Attention scores: (B,H,T,d)·(B,H,d,T) = (B,H,T,T), which is why attention memory grows with T².
• Rank = number of linearly independent rows (= columns) = dimension of the output space the matrix can reach. An outer product u vᵀ is rank 1.
• SVD: any matrix W = U Σ Vᵀ with orthonormal U, V and singular values s₁≥s₂≥…≥0. Keeping the top r terms gives the best rank-r approximation in Frobenius norm (Eckart-Young); the squared error equals the sum of the discarded s². Many trained weight updates are close to low rank, which is the bet behind LoRA: learn ΔW = B·A with B: d×r, A: r×d, r≪d, so you train 2dr parameters instead of d².
• Norms: L1 = sum|x| (promotes sparsity), L2 = √(sum x²) (smooth, weight decay), L∞ = max|x| (worst case; used in quantization range clipping).
• Eigen: A v = λ v: directions the matrix only stretches. For a symmetric matrix the eigenvectors are orthogonal; the covariance matrix's eigenvectors are PCA's directions. The condition number (largest/smallest singular value) says how badly a problem is stretched; large condition numbers make gradient descent slow and numerical solves unstable.

Copyimport numpy as np
rng = np.random.default_rng(42)
A=rng.normal(size=(4,3)); B=rng.normal(size=(3,5)); print("matmul shapes:", (A@B).shape)
# attention-style: (B,H,T,d) @ (B,H,d,T)
print("QK^T:", (np.zeros((2,8,128,64))@np.zeros((2,8,64,128))).shape)
# low-rank matrix + noise
m,n,r=512,512,16
W = rng.normal(size=(m,r))@rng.normal(size=(r,n)) + 0.05*rng.normal(size=(m,n))
U,S,Vt=np.linalg.svd(W,full_matrices=False)
print("top singular values:", np.round(S[:3],1), " ... s[15]=",round(S[15],1)," s[16]=",round(S[16],1)," s[100]=",round(S[100],1))
print(f"{'rank':>5} {'rel Frobenius err':>18} {'params':>9} {'vs full':>8}")
for k in (1,4,8,16,32,64):
    Wk=(U[:,:k]*S[:k])@Vt[:k]
    print(f"{k:>5} {np.linalg.norm(W-Wk)/np.linalg.norm(W):>18.4f} {k*(m+n):>9,} {k*(m+n)/(m*n):>8.1%}")
k=16; print(f"energy in top-{k}: {(S[:k]**2).sum()/(S**2).sum():.4f}; Eckart-Young check err^2 = sum of tail s^2: {np.linalg.norm(W-(U[:,:k]*S[:k])@Vt[:k])**2:.3f} vs {(S[k:]**2).sum():.3f}")
# LoRA param count for 4096x4096 r=16
d=4096;r=16; print(f"LoRA d=4096 r=16: full={d*d:,} lora={2*d*r:,} ({2*d*r/(d*d):.2%})")
# norms
v=np.array([3.,4.,12.]); print("L1,L2,Linf:", np.abs(v).sum(), np.linalg.norm(v), np.abs(v).max())
# eigen
C=np.array([[2,1],[1,2.]]); w,V=np.linalg.eigh(C); print("eig of [[2,1],[1,2]]:", w, " vectors:\n",np.round(V,3))
# power iteration
x=np.array([1.,0]);
for _ in range(20): x=C@x; x/=np.linalg.norm(x)
print("power iteration ->", np.round(x,3))
# rank
R=np.outer(rng.normal(size=5),rng.normal(size=5)); print("rank of outer product:", np.linalg.matrix_rank(R), " rank of random 5x5:", np.linalg.matrix_rank(rng.normal(size=(5,5))))
print("condition number of [[1,0],[0,1e-6]]:", np.linalg.cond(np.diag([1,1e-6])))

matmul shapes: (4, 5)
QK^T: (2, 8, 128, 128)
top singular values: [604.7 599.2 583.3]  ... s[15]= 406.0  s[16]= 2.2  s[100]= 1.6
 rank  rel Frobenius err    params  vs full
    1             0.9560     1,024     0.4%
    4             0.8197     4,096     1.6%
    8             0.6352     8,192     3.1%
   16             0.0120    16,384     6.2%
   32             0.0113    32,768    12.5%
   64             0.0101    65,536    25.0%
energy in top-16: 0.9999; Eckart-Young check err^2 = sum of tail s^2: 616.474 vs 616.474
LoRA d=4096 r=16: full=16,777,216 lora=131,072 (0.78%)
L1,L2,Linf: 19.0 13.0 12.0
eig of [[2,1],[1,2]]: [1. 3.]  vectors:
 [[-0.707  0.707]
 [ 0.707  0.707]]
power iteration -> [0.707 0.707]
rank of outer product: 1  rank of random 5x5: 5
condition number of [[1,0],[0,1e-6]]: 1000000.0

6. Calculus & optimization

• Gradient: vector of partial derivatives of the loss with respect to each parameter; it points uphill, so we step the other way: θ ← θ - lr · ∇L.
• Chain rule / backprop: if L = f(g(h(w))), dL/dw = f'(g)·g'(h)·h'(w). Backprop is this applied from the output backward, reusing intermediate values saved from the forward pass (that stored activation memory is what activation checkpointing trades for compute).
• SGD uses a noisy minibatch gradient. Momentum keeps an exponentially weighted velocity, which damps oscillation across steep directions and accelerates along shallow ones. Adam divides each parameter's step by a running RMS of its gradients (per-parameter step size) and bias-corrects the early averages; it is far less sensitive to gradient scale.
• Learning rate: on a quadratic with curvature λ_max, plain gradient descent converges only for lr < 2/λ_max. Too small: slow. Too big: diverges. The slow direction is set by the smallest curvature, hence ill-conditioning hurts.
• Convexity: a convex loss (linear/logistic regression) has one global minimum, so any local minimum is global. Neural nets are non-convex; the starting point and noise matter, and in practice many minima are about equally good.

Logistic gradient: derivation

p = σ(w·x), σ'(z) = σ(1-σ). Loss per example L = -[y ln p + (1-y) ln(1-p)]. dL/dp = -y/p + (1-y)/(1-p) = (p-y)/(p(1-p)). dp/dz = p(1-p). So dL/dz = p - y, and dL/dw = (p - y)·x. Over a batch: ∇L = Xᵀ(p - y)/n. The p(1-p) cancels, which is why cross-entropy with a sigmoid/softmax trains well where MSE would stall on saturated outputs.

Copyimport numpy as np
rng = np.random.default_rng(42)
sig=lambda z:1/(1+np.exp(-z))
X=rng.normal(size=(50,3)); y=(rng.random(50)<sig(X@[1.5,-2,0.5])).astype(float); w=rng.normal(size=3)*0.1
def loss(w): p=sig(X@w); return -np.mean(y*np.log(p)+(1-y)*np.log(1-p))
g=X.T@(sig(X@w)-y)/len(y)
num=np.array([(loss(w+e)-loss(w-e))/2e-5 for e in np.eye(3)*1e-5])
print("analytic grad:", np.round(g,6)); print("numeric  grad:", np.round(num,6)); print("max abs diff:", np.abs(g-num).max())
# tiny net backprop: y=w2*relu(w1*x), L=(yhat-t)^2/2
x,t,w1,w2=2.0,1.0,0.5,-0.3  # h = relu(1.0)=1.0 ; yhat=-0.3
h=max(0,w1*x); yh=w2*h; L=0.5*(yh-t)**2
dL=yh-t; gw2=dL*h; gw1=dL*w2*(1.0 if w1*x>0 else 0)*x
f=lambda a,b:0.5*(b*max(0,a*x)-t)**2
print(f"\ntiny net: L={L:.4f} dL/dw2={gw2:.4f} (num {(f(w1,w2+1e-6)-f(w1,w2-1e-6))/2e-6:.4f}) dL/dw1={gw1:.4f} (num {(f(w1+1e-6,w2)-f(w1-1e-6,w2))/2e-6:.4f})")
# optimizers on ill-conditioned quadratic f=0.5*(x^2+25 y^2)
H=np.array([1.,25.]); fq=lambda p:0.5*(H*p*p).sum()
def run(opt,lr,steps=100):
    p=np.array([5.,1.]); m=np.zeros(2); v=np.zeros(2); b1,b2=0.9,0.999
    for t in range(1,steps+1):
        g=H*p
        if opt=="sgd": p=p-lr*g
        elif opt=="mom": m=0.9*m+g; p=p-lr*m
        else:
            m=b1*m+(1-b1)*g; v=b2*v+(1-b2)*g*g
            p=p-lr*(m/(1-b1**t))/(np.sqrt(v/(1-b2**t))+1e-8)
        if not np.isfinite(fq(p)) or fq(p)>1e12: return float('inf')
    return fq(p)
print("\nf after 100 steps (start f=%.1f, Hessian eigs 1 and 25 -> SGD stable if lr<2/25=0.08):"%fq(np.array([5.,1.])))
for lr in (0.001,0.01,0.04,0.079,0.081,0.2):
    print(f"lr={lr:<6} sgd={run('sgd',lr):<12.3e} momentum(b=.9)={run('mom',lr):<12.3e}")
for lr in (0.01,0.1,0.5,1.0):
    print(f"adam lr={lr:<5} f={run('adam',lr):.3e}")
# convexity: non-convex multiple minima
f2=lambda x:x**4-3*x**2+x; g2=lambda x:4*x**3-6*x+1
for x0 in (-2,2):
    x=x0
    for _ in range(500): x-=0.01*g2(x)
    print(f"non-convex x^4-3x^2+x from {x0:+d} -> x={x:.3f}, f={f2(x):.3f}")

analytic grad: [-0.121981  0.278311  0.101886]
numeric  grad: [-0.121981  0.278311  0.101886]
max abs diff: 7.746497887595183e-12

tiny net: L=0.8450 dL/dw2=-1.3000 (num -1.3000) dL/dw1=0.7800 (num 0.7800)

f after 100 steps (start f=25.0, Hessian eigs 1 and 25 -> SGD stable if lr<2/25=0.08):
lr=0.001  sgd=1.031e+01    momentum(b=.9)=1.699e+00   
lr=0.01   sgd=1.675e+00    momentum(b=.9)=8.731e-05   
lr=0.04   sgd=3.558e-03    momentum(b=.9)=1.641e-04   
lr=0.079  sgd=7.904e-02    momentum(b=.9)=3.250e-04   
lr=0.081  sgd=1.745e+03    momentum(b=.9)=3.841e-04   
lr=0.2    sgd=inf          momentum(b=.9)=inf         
adam lr=0.01  f=8.774e+00
adam lr=0.1   f=8.685e-04
adam lr=0.5   f=3.439e-04
adam lr=1.0   f=2.802e-04
non-convex x^4-3x^2+x from -2 -> x=-1.301, f=-3.514
non-convex x^4-3x^2+x from +2 -> x=1.131, f=-1.070

7. Supervised learning basics

• Linear regression: ŷ = w·x + b; minimize mean squared error (MSE). Closed form w = (XᵀX)⁻¹Xᵀy. MSE is the MLE under Gaussian noise.
• Logistic regression: p = σ(w·x + b); minimize cross-entropy = average negative log-likelihood of the labels (so "cross-entropy loss" and "maximum likelihood" are the same thing for Bernoulli/categorical labels). Decision boundary is linear in x. Output is a probability only if the model is calibrated (section 9).
• Regularization adds a penalty to shrink weights: L2 / ridge (λ∑w²) shrinks all weights smoothly; L1 / lasso (λ∑|w|) drives many weights to exactly zero (feature selection). It is a Gaussian (L2) or Laplace (L1) prior on the weights in Bayesian terms.
• Bias-variance: expected test error = bias² (systematic error from a too-simple model) + variance (sensitivity to which training set you got) + irreducible noise. More capacity lowers bias and raises variance. Overfitting = low train error, high test error.
• Splits: train fits parameters; validation chooses hyperparameters/models; test is touched once at the end. Leakage = information from outside the training world sneaking in: preprocessing/feature selection/scaling fit on all data, duplicates across splits, features that encode the label, random splits of time-series (training on the future). Cross-validation rotates the validation fold for a more stable estimate; for time series use forward-chaining splits, and fit all preprocessing inside each fold.

Copyimport numpy as np
from sklearn.linear_model import Ridge, Lasso, LinearRegression, LogisticRegression
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import cross_val_score, KFold, train_test_split
rng = np.random.default_rng(42)
f=lambda x:np.sin(2*np.pi*x)
x=rng.random(30); y=f(x)+rng.normal(0,0.3,30); xt=rng.random(1000); yt=f(xt)+rng.normal(0,0.3,1000)
print("polynomial degree: train MSE vs test MSE (noise var = 0.09)")
for d in (1,3,5,9,15):
    m=make_pipeline(PolynomialFeatures(d),StandardScaler(),LinearRegression()).fit(x[:,None],y)
    print(f"deg={d:>2} train={np.mean((m.predict(x[:,None])-y)**2):.3f} test={np.mean((m.predict(xt[:,None])-yt)**2):.3f}")
# bias-variance decomposition at x0=0.25 by resampling
print("\nbias^2 / variance at x0=0.3 over 500 training sets (n=30):")
for d in (1,3,9):
    pr=[]
    for _ in range(500):
        xs=rng.random(30); ys=f(xs)+rng.normal(0,0.3,30)
        pr.append(make_pipeline(PolynomialFeatures(d),StandardScaler(),LinearRegression()).fit(xs[:,None],ys).predict([[0.3]])[0])
    pr=np.array(pr); print(f"deg={d} bias^2={(pr.mean()-f(0.3))**2:.4f} var={pr.var():.4f}")
# L1 vs L2
X=rng.normal(size=(100,20)); wt=np.zeros(20); wt[:3]=[3,-2,1.5]; yy=X@wt+rng.normal(0,1,100)
r=Ridge(alpha=10).fit(X,yy); l=Lasso(alpha=0.2).fit(X,yy)
print(f"\ntrue nonzero=3 | ridge nonzero coefs={np.sum(np.abs(r.coef_)>1e-6)}, lasso nonzero={np.sum(np.abs(l.coef_)>1e-6)}")
print("ridge coef[:5]", np.round(r.coef_[:5],2), " lasso coef[:5]", np.round(l.coef_[:5],2))
# CE = MLE
p=np.array([.9,.2,.7]); yv=np.array([1,0,1]); print(f"\nbinary CE = -mean log-lik = {-np.mean(yv*np.log(p)+(1-yv)*np.log(1-p)):.4f}")
# leakage: feature selection on all data with pure noise
Xn=rng.normal(size=(60,2000)); yn=rng.integers(0,2,60)
from sklearn.feature_selection import SelectKBest, f_classif
Xs=SelectKBest(f_classif,k=10).fit_transform(Xn,yn)
leak=cross_val_score(LogisticRegression(),Xs,yn,cv=5).mean()
from sklearn.pipeline import Pipeline
ok=cross_val_score(Pipeline([("s",SelectKBest(f_classif,k=10)),("m",LogisticRegression())]),Xn,yn,cv=5).mean()
print(f"leakage on pure-noise labels: select-then-CV acc={leak:.2f} (should be ~0.5); selection inside CV={ok:.2f}")
# time split leakage-ish: k-fold CV score
Xa=rng.normal(size=(200,5)); ya=Xa@[1,2,0,0,0]+rng.normal(size=200)
print("5-fold R2:", np.round(cross_val_score(LinearRegression(),Xa,ya,cv=KFold(5,shuffle=True,random_state=0)),3))

polynomial degree: train MSE vs test MSE (noise var = 0.09)
deg= 1 train=0.206 test=0.294
deg= 3 train=0.064 test=0.102
deg= 5 train=0.047 test=0.122
deg= 9 train=0.044 test=0.113
deg=15 train=0.031 test=1.885

bias^2 / variance at x0=0.3 over 500 training sets (n=30):
deg=1 bias^2=0.3179 var=0.0163
deg=3 bias^2=0.0036 var=0.0115
deg=9 bias^2=0.0000 var=0.0393

true nonzero=3 | ridge nonzero coefs=20, lasso nonzero=4
ridge coef[:5] [ 2.72 -1.66  1.23  0.12 -0.15]  lasso coef[:5] [ 2.88 -1.56  1.21  0.   -0.  ]

binary CE = -mean log-lik = 0.2284
leakage on pure-noise labels: select-then-CV acc=0.88 (should be ~0.5); selection inside CV=0.52
5-fold R2: [0.716 0.828 0.816 0.881 0.808]

8. Trees & ensembles

• Decision tree: repeatedly split on the feature/threshold that most reduces impurity. Gini = 1 - ∑p²; entropy = -∑p log₂p. Gain = parent impurity - weighted child impurity. An unrestricted tree memorizes the training set (high variance).
• Random forest: many deep trees, each on a bootstrap sample and a random subset of features per split, predictions averaged. Averaging decorrelated high-variance trees cuts variance. Free validation: out-of-bag (OOB) score.
• Gradient boosting: build shallow trees sequentially, each fit to the residuals (negative gradient of the loss) of the current ensemble, add it with a small learning rate. Boosting mostly reduces bias; too many trees overfit, so tune n_trees, learning rate, depth, and use early stopping on a validation set.
• XGBoost (and LightGBM) are optimized gradient boosting.

Inference angle

Tree models are cheap, branchy, CPU-friendly, and have no matmuls: latency scales with (number of trees × depth), so on the serving side the levers are tree count, depth, batching rows per call, and thread count, not GPU tensor cores. This is the opposite regime from LLM serving (compute-/memory-bandwidth-bound matmuls).

Copyimport numpy as np, time
from sklearn.tree import DecisionTreeClassifier, export_text
from sklearn.ensemble import RandomForestClassifier, GradientBoostingRegressor
from sklearn.datasets import make_classification, make_friedman1
from sklearn.model_selection import cross_val_score, train_test_split
rng=np.random.default_rng(42)
gini=lambda p:1-p**2-(1-p)**2
ent=lambda p:-(p*np.log2(p)+(1-p)*np.log2(1-p))
# split example: parent 10 pos/10 neg; left 8/2, right 2/8
print(f"parent gini={gini(.5):.3f}; children (8/2, 2/8) gini={gini(.8):.3f} each; gain={gini(.5)-gini(.8):.3f}")
print(f"entropy parent={ent(.5):.3f} child={ent(.8):.3f} info gain={ent(.5)-ent(.8):.3f}")
X,y=make_classification(2000,20,n_informative=5,flip_y=0.1,random_state=0)
Xtr,Xte,ytr,yte=train_test_split(X,y,random_state=0)
for d in (2,5,None):
    t=DecisionTreeClassifier(max_depth=d,random_state=0).fit(Xtr,ytr); print(f"tree depth={d}: train acc={t.score(Xtr,ytr):.3f} test acc={t.score(Xte,yte):.3f}")
rf=RandomForestClassifier(300,random_state=0,n_jobs=-1).fit(Xtr,ytr); print(f"random forest: train={rf.score(Xtr,ytr):.3f} test={rf.score(Xte,yte):.3f}")
print("OOB score:", round(RandomForestClassifier(300,oob_score=True,random_state=0,n_jobs=-1).fit(Xtr,ytr).oob_score_,3))
# boosting by hand: residual fitting with stumps on Friedman
from sklearn.tree import DecisionTreeRegressor
Xf,yf=make_friedman1(1000,random_state=0); Xa,Xb,ya,yb=Xf[:700],Xf[700:],yf[:700],yf[700:]
pred=np.full(len(ya),ya.mean()); predb=np.full(len(yb),ya.mean()); lr=0.1
print("\nhand-rolled gradient boosting (depth-3 trees on residuals, lr=0.1):")
for m in range(1,201):
    t=DecisionTreeRegressor(max_depth=3).fit(Xa,ya-pred); pred+=lr*t.predict(Xa); predb+=lr*t.predict(Xb)
    if m in (1,10,50,200): print(f"  trees={m:>3} train MSE={np.mean((ya-pred)**2):6.2f} test MSE={np.mean((yb-predb)**2):6.2f}")
g=GradientBoostingRegressor(n_estimators=200,max_depth=3,learning_rate=0.1,random_state=0).fit(Xa,ya)
print(f"sklearn GBR test MSE={np.mean((yb-g.predict(Xb))**2):.2f}")
try:
    import xgboost as xgb; print("xgboost", xgb.__version__)
except Exception as e: print("xgboost not installed here:", type(e).__name__)

parent gini=0.500; children (8/2, 2/8) gini=0.320 each; gain=0.180
entropy parent=1.000 child=0.722 info gain=0.278
tree depth=2: train acc=0.821 test acc=0.808
tree depth=5: train acc=0.892 test acc=0.856
tree depth=None: train acc=1.000 test acc=0.816
random forest: train=1.000 test=0.874
OOB score: 0.867

hand-rolled gradient boosting (depth-3 trees on residuals, lr=0.1):
  trees=  1 train MSE= 22.38 test MSE= 23.30
  trees= 10 train MSE=  8.81 test MSE= 10.07
  trees= 50 train MSE=  1.04 test MSE=  2.01
  trees=200 train MSE=  0.11 test MSE=  0.87
sklearn GBR test MSE=0.89
xgboost not installed here: ModuleNotFoundError

9. Evaluation metrics

• Precision = TP/(TP+FP): of what I flagged, how much is right. Recall = TP/(TP+FN): of what exists, how much I found. F1 = harmonic mean. Accuracy misleads under imbalance (always predicting "negative" is 99% accurate at 1% prevalence).
• ROC-AUC = P(random positive scores above random negative); plots TPR vs FPR. It is insensitive to prevalence, which makes it look good even when the positive class is rare and FP dominate. PR-AUC (average precision) tracks precision, so it is the better lens when positives are rare and false positives are costly. Baseline for PR-AUC is the prevalence, for ROC-AUC it is 0.5.
• Calibration: when the model says 0.7, is it right 70% of the time? Needed whenever the probability is used (bids, expected value, thresholds across segments). Measure with reliability curves and ECE; fix with Platt scaling or isotonic regression. Log-loss punishes confident mistakes and measures calibration plus ranking; AUC measures ranking only.
• Ranking (ads/search): Precision@k, MRR (mean of 1/rank of the first relevant item), NDCG@k = DCG/ideal DCG with DCG = ∑ gain_i/log₂(i+1), where gain = 2^rel-1 (common variant) or rel (scikit-learn's). Position-discounted: a hit at rank 1 matters more than at rank 10.

Copyimport numpy as np
from sklearn.metrics import roc_auc_score, average_precision_score, log_loss, precision_score, recall_score, f1_score
from sklearn.linear_model import LogisticRegression
from sklearn.calibration import calibration_curve
rng=np.random.default_rng(42)
# imbalance: 1% positives
n=100000; y=(rng.random(n)<0.01).astype(int)
s=np.where(y==1, rng.normal(1.5,1,n), rng.normal(0,1,n))
print(f"prevalence={y.mean():.4f} ROC-AUC={roc_auc_score(y,s):.3f} PR-AUC(AP)={average_precision_score(y,s):.3f} (random baseline AP~{y.mean():.3f})")
thr=1.5; pr=s>thr
print(f"thr={thr}: precision={precision_score(y,pr):.3f} recall={recall_score(y,pr):.3f} F1={f1_score(y,pr):.3f}; FPR={((pr==1)&(y==0)).sum()/(y==0).sum():.4f}")
# manual P/R
tp=((pr==1)&(y==1)).sum(); fp=((pr==1)&(y==0)).sum(); fn=((pr==0)&(y==1)).sum(); print(f"TP={tp} FP={fp} FN={fn}")
# accuracy trap
print(f"always-negative accuracy = {1-y.mean():.3f}")
# calibration: ranking identical, probabilities wrong
X=rng.normal(size=(20000,3)); yy=(rng.random(20000)<1/(1+np.exp(-(X@[1,1,1])))).astype(int)
m=LogisticRegression().fit(X,yy); p=m.predict_proba(X)[:,1]
pbad=1/(1+np.exp(-(3*np.log(p/(1-p)))))  # overconfident monotone transform
print(f"\ncalibrated: AUC={roc_auc_score(yy,p):.4f} logloss={log_loss(yy,p):.4f} | overconfident (same ranking): AUC={roc_auc_score(yy,pbad):.4f} logloss={log_loss(yy,pbad):.4f}")
fp_,mp=calibration_curve(yy,pbad,n_bins=5,strategy='quantile'); print("overconfident: mean predicted vs observed:", np.round(mp,2), np.round(fp_,2))
ece=lambda p,y,b=10:sum(abs(y[(p>=i/b)&(p<(i+1)/b)].mean()-p[(p>=i/b)&(p<(i+1)/b)].mean())*((p>=i/b)&(p<(i+1)/b)).mean() for i in range(b) if ((p>=i/b)&(p<(i+1)/b)).any())
print(f"ECE calibrated={ece(p,yy):.4f} overconfident={ece(pbad,yy):.4f}")
# ranking metrics
def dcg(r): r=np.asarray(r,float); return np.sum((2**r-1)/np.log2(np.arange(2,len(r)+2)))
rel=[3,2,0,1,0]; print(f"\nrels in ranked order {rel}: DCG={dcg(rel):.4f} IDCG={dcg(sorted(rel,reverse=True)):.4f} NDCG={dcg(rel)/dcg(sorted(rel,reverse=True)):.4f}")
rel2=[0,0,3,2,1]; print(f"rels {rel2}: NDCG={dcg(rel2)/dcg(sorted(rel2,reverse=True)):.4f}")
rr=[1/2,1,1/3]; print("MRR for first-relevant ranks 2,1,3 =", round(np.mean(rr),4))
from sklearn.metrics import ndcg_score
print("sklearn ndcg check:", round(ndcg_score([rel],[[5,4,3,2,1]]),4))
# Pr@k
lin=lambda r:np.sum(np.asarray(r,float)/np.log2(np.arange(2,len(r)+2)));print("linear-gain NDCG (what sklearn uses):",round(lin(rel)/lin(sorted(rel,reverse=True)),4))
print("P@3 for [1,0,1,0,0] =", round(2/3,4))
# AUC = P(score_pos > score_neg)
a=rng.normal(1,1,2000); b=rng.normal(0,1,2000); print("AUC pairwise:", round(np.mean(a[:,None]>b[None,:]),4), " sklearn:", round(roc_auc_score(np.r_[np.ones(2000),np.zeros(2000)],np.r_[a,b]),4))

prevalence=0.0105 ROC-AUC=0.858 PR-AUC(AP)=0.108 (random baseline AP~0.010)
thr=1.5: precision=0.071 recall=0.478 F1=0.123; FPR=0.0662
TP=500 FP=6551 FN=547
always-negative accuracy = 0.990

calibrated: AUC=0.8409 logloss=0.4920 | overconfident (same ranking): AUC=0.8409 logloss=0.7122
overconfident: mean predicted vs observed: [0.   0.07 0.48 0.92 1.  ] [0.09 0.28 0.49 0.72 0.89]
ECE calibrated=0.0077 overconfident=0.1409

rels in ranked order [3, 2, 0, 1, 0]: DCG=9.3235 IDCG=9.3928 NDCG=0.9926
rels [0, 0, 3, 2, 1]: NDCG=0.5514
MRR for first-relevant ranks 2,1,3 = 0.6111
sklearn ndcg check: 0.9854
linear-gain NDCG (what sklearn uses): 0.9854
P@3 for [1,0,1,0,0] = 0.6667
AUC pairwise: 0.7558  sklearn: 0.7558

10. Unsupervised learning & embeddings

• k-means (Lloyd's algorithm): alternately assign each point to the nearest centroid, then move each centroid to its cluster's mean; minimizes within-cluster squared distance (inertia). Needs k, is sensitive to initialization (use several), assumes roughly spherical clusters. Choose k with the elbow of inertia or silhouette.
• PCA: rotate to the orthogonal directions of maximum variance; = SVD of the mean-centered data. Explained-variance ratio = s²/∑s². Forgetting to center is a classic mistake.
• Cosine similarity = angle between vectors, ignoring length: cos = a·b/(|a||b|). For unit vectors, |a-b|² = 2 - 2cos, so cosine ranking = L2 ranking = dot-product ranking after normalization. If vectors are not normalized, dot product, cosine, and L2 can rank differently.
• Nearest-neighbor search: brute force is O(N·d) per query and memory-bandwidth-bound (you read the whole index). Approximate indexes (HNSW graph, IVF clusters, product quantization) trade a little recall for orders of magnitude less work. In high dimensions distances concentrate (everything is roughly equally far), which is why recall is measured, not assumed.

Copyimport numpy as np, time
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.datasets import make_blobs
rng=np.random.default_rng(42)
X,_=make_blobs(600,2,centers=3,cluster_std=1.0,random_state=1)
for k in (1,2,3,4,5,6):
    km=KMeans(k,n_init=10,random_state=0).fit(X); print(f"k={k} inertia={km.inertia_:8.1f}")
# lloyd by hand
c=X[rng.choice(len(X),3,replace=False)]
for it in range(10):
    a=np.argmin(((X[:,None]-c[None])**2).sum(-1),1); c=np.array([X[a==j].mean(0) for j in range(3)])
print("hand-rolled lloyd inertia (k=3):", round(sum(((X[a==j]-c[j])**2).sum() for j in range(3)),1))
# PCA via SVD
Z=rng.normal(size=(1000,3))@np.array([[3,0,0],[0,1,0],[0,0,0.2]])@np.linalg.qr(rng.normal(size=(3,3)))[0]
p=PCA().fit(Z); print("\nPCA explained var ratio:", np.round(p.explained_variance_ratio_,4))
Zc=Z-Z.mean(0); S=np.linalg.svd(Zc,compute_uv=False); print("via SVD  s^2/(n-1) / total:", np.round(S**2/(S**2).sum(),4))
# cosine vs L2
a=rng.normal(size=128); b=3*a; 
cos=lambda u,v:u@v/np.linalg.norm(u)/np.linalg.norm(v)
print(f"\ncos(a,3a)={cos(a,b):.3f} L2 dist={np.linalg.norm(a-b):.2f}")
u=rng.normal(size=(5,128)); u/=np.linalg.norm(u,axis=1,keepdims=True); q=u[0]
print("unit vectors: ||a-b||^2 = 2-2cos ->", round(np.linalg.norm(u[0]-u[1])**2,4), round(2-2*cos(u[0],u[1]),4))
# curse of dimensionality + brute-force NN cost
for d in (2,10,100,1000):
    P=rng.random((2000,d)); D=np.linalg.norm(P[1:]-P[0],axis=1); print(f"d={d:>4} (max-min)/min dist = {(D.max()-D.min())/D.min():.3f}")
N,d=200_000,128; DB=rng.normal(size=(N,d)).astype(np.float32); q=rng.normal(size=d).astype(np.float32)
t=time.perf_counter(); s=DB@q; top=np.argpartition(-s,10)[:10]; dt=time.perf_counter()-t
print(f"\nbrute force dot+top10 over {N:,}x{d} fp32: {dt*1000:.1f} ms; memory scanned = {DB.nbytes/1e6:.0f} MB -> {DB.nbytes/dt/1e9:.1f} GB/s effective (machine-dependent)")

k=1 inertia= 24777.4
k=2 inertia=  3790.3
k=3 inertia=  1161.4
k=4 inertia=  1005.5
k=5 inertia=   879.4
k=6 inertia=   750.7
hand-rolled lloyd inertia (k=3): 1161.4

PCA explained var ratio: [0.8993 0.0967 0.004 ]
via SVD  s^2/(n-1) / total: [0.8993 0.0967 0.004 ]

cos(a,3a)=1.000 L2 dist=20.15
unit vectors: ||a-b||^2 = 2-2cos -> 2.0982 2.0982
d=   2 (max-min)/min dist = 171.064
d=  10 (max-min)/min dist = 2.094
d= 100 (max-min)/min dist = 0.521
d=1000 (max-min)/min dist = 0.118

brute force dot+top10 over 200,000x128 fp32: 10.9 ms; memory scanned = 102 MB -> 9.4 GB/s effective (machine-dependent)

11. Information theory for LLMs

• Entropy H(p) = -∑p log p: the average surprise / minimum average bits to encode samples. Max for uniform (log₂8 = 3 bits for 8 outcomes), 0 when deterministic.
• Cross-entropy H(p,q) = -∑p log q: the average cost of coding data from p using a code built for q. The LLM training loss is the cross-entropy between the data distribution (one-hot next token) and the model: just the negative log-probability of the true next token.
• KL divergence KL(p||q) = H(p,q) - H(p) ≥ 0: the extra cost of using q. Not symmetric. Used in distillation, RLHF/PPO penalties, speculative decoding acceptance analysis.
• Perplexity = exp(mean NLL in nats) = 2^(bits per token): the effective number of equally likely choices per token. A uniform model over V tokens has perplexity V. Compare perplexities only with the same tokenizer and data.
• Softmax with temperature: p_i = exp(z_i/T)/∑exp(z_j/T). T→0: argmax (greedy); T=1: model distribution; T↑: flatter, higher entropy. Compute it as exp(z - max z) for numerical stability. Gradient of softmax + cross-entropy w.r.t. logits is simply p - onehot(y).

Copyimport numpy as np
H=lambda p:-np.sum(p[p>0]*np.log2(p[p>0]))
print("entropy bits: fair coin",H(np.array([.5,.5])),"| biased .9/.1",round(H(np.array([.9,.1])),4),"| uniform over 8",H(np.ones(8)/8),"| deterministic",H(np.array([1.,0])))
p=np.array([.7,.2,.1]); q=np.array([.5,.3,.2])
ce=-np.sum(p*np.log2(q)); kl=np.sum(p*np.log2(p/q)); klr=np.sum(q*np.log2(q/p))
print(f"H(p)={H(p):.4f} CE(p,q)={ce:.4f} KL(p||q)={kl:.4f} -> H+KL={H(p)+kl:.4f}; KL(q||p)={klr:.4f} (asymmetric)")
# perplexity
nll=np.array([2.1,0.3,1.5,0.9,3.0]) # nats per token
print(f"mean NLL={nll.mean():.3f} nats -> perplexity={np.exp(nll.mean()):.3f}; uniform over V=50000 has ppl 50000, ln V={np.log(50000):.3f} nats")
print("bits per token:", round(nll.mean()/np.log(2),3))
sm=lambda z,T=1:np.exp(z/T-np.max(z/T))/np.sum(np.exp(z/T-np.max(z/T)))
z=np.array([2.0,1.0,0.1])
for T in (0.1,0.5,1,2,5): 
    s=sm(z,T); print(f"T={T:<4} probs={np.round(s,4)} entropy={H(s):.3f} bits")
print("overflow-safe softmax of [1000,1001]:", sm(np.array([1000.,1001.])))
with np.errstate(all='ignore'): print("naive softmax of [1000,1001]:", np.exp([1000.,1001.])/np.exp([1000.,1001.]).sum())
# softmax + CE gradient = p - onehot
z=np.array([1.,2.,0.5]); y=np.array([0,1,0]); p=sm(z); g=p-y
L=lambda z:-np.log(sm(z)[1]); num=[(L(z+e)-L(z-e))/2e-6 for e in np.eye(3)*1e-6]
print("softmax-CE grad analytic",np.round(g,5),"numeric",np.round(num,5))
# logsumexp
print("logsumexp stable:", np.max(z)+np.log(np.sum(np.exp(z-np.max(z)))))

entropy bits: fair coin 1.0 | biased .9/.1 0.469 | uniform over 8 3.0 | deterministic -0.0
H(p)=1.1568 CE(p,q)=1.2796 KL(p||q)=0.1228 -> H+KL=1.2796; KL(q||p)=0.1328 (asymmetric)
mean NLL=1.560 nats -> perplexity=4.759; uniform over V=50000 has ppl 50000, ln V=10.820 nats
bits per token: 2.251
T=0.1  probs=[1. 0. 0.] entropy=0.001 bits
T=0.5  probs=[0.8638 0.1169 0.0193] entropy=0.655 bits
T=1    probs=[0.659  0.2424 0.0986] entropy=1.222 bits
T=2    probs=[0.5017 0.3043 0.194 ] entropy=1.481 bits
T=5    probs=[0.3996 0.3272 0.2733] entropy=1.568 bits
overflow-safe softmax of [1000,1001]: [0.26894142 0.73105858]
naive softmax of [1000,1001]: [nan nan]
softmax-CE grad analytic [ 0.23122 -0.37147  0.14024] numeric [ 0.23122 -0.37147  0.14024]
logsumexp stable: 2.464368784107945

12. Interview prep: questions with model answers

Try answering aloud before opening. Two extra demos back the first and third questions.

Copyimport numpy as np
rng = np.random.default_rng(42)
# Q: "p99 went 80 -> 95 ms: real?"  Same true distribution (p99=80), 1000 requests per run, 2000 pairs of runs.
sigma=0.5; med=80/np.exp(sigma*2.3263)
def p99(n): return np.percentile(rng.lognormal(np.log(med),sigma,n),99)
for n in (500,1000,5000):
    a=np.array([p99(n) for _ in range(2000)]); b=np.array([p99(n) for _ in range(2000)])
    print(f"n={n:>5}/run, NO real change: sd of p99 across runs={a.std():5.1f}ms; P(run2 - run1 >= 15ms)={np.mean(b-a>=15):.3f}; P(|diff|>=15)={np.mean(np.abs(b-a)>=15):.3f}")

n=  500/run, NO real change: sd of p99 across runs=  6.5ms; P(run2 - run1 >= 15ms)=0.044; P(|diff|>=15)=0.092
n= 1000/run, NO real change: sd of p99 across runs=  4.5ms; P(run2 - run1 >= 15ms)=0.014; P(|diff|>=15)=0.021
n= 5000/run, NO real change: sd of p99 across runs=  2.1ms; P(run2 - run1 >= 15ms)=0.000; P(|diff|>=15)=0.000

Copyimport numpy as np
rng = np.random.default_rng(42)
# dynamic batching: Poisson arrivals, close batch after 10ms window or at max 8
lam, window, maxb = 300.0, 0.010, 8
t = np.cumsum(rng.exponential(1/lam, 200_000)); i=0; sizes=[]
while i < len(t):
    j=i
    while j+1<len(t) and j+1-i<maxb and t[j+1]-t[i] < window: j+=1
    sizes.append(j-i+1); i=j+1
s=np.array(sizes)
print(f"batches={len(s):,}  mean batch size (per batch) = {s.mean():.2f}")
print(f"mean batch size seen by a REQUEST (size-weighted) = {np.sum(s*s)/np.sum(s):.2f}")
print("batch size histogram:", {int(k):int(v) for k,v in zip(*np.unique(s,return_counts=True))})
# toy cost model: step time = 6ms fixed + 1.5ms per item (weights read once; compute scales)
for b in (1,2,4,8,16,32):
    st=6+1.5*b; print(f"batch={b:>2} step={st:5.1f}ms  throughput={b/st*1000:6.0f} items/s  per-item cost={st/b:5.2f}ms")

batches=50,224  mean batch size (per batch) = 3.98
mean batch size seen by a REQUEST (size-weighted) = 4.70
batch size histogram: {1: 2465, 2: 7638, 3: 11235, 4: 11140, 5: 8448, 6: 5079, 7: 2472, 8: 1747}
batch= 1 step=  7.5ms  throughput=   133 items/s  per-item cost= 7.50ms
batch= 2 step=  9.0ms  throughput=   222 items/s  per-item cost= 4.50ms
batch= 4 step= 12.0ms  throughput=   333 items/s  per-item cost= 3.00ms
batch= 8 step= 18.0ms  throughput=   444 items/s  per-item cost= 2.25ms
batch=16 step= 30.0ms  throughput=   533 items/s  per-item cost= 1.88ms
batch=32 step= 54.0ms  throughput=   593 items/s  per-item cost= 1.69ms