# theory-demos Five small, self-contained numerical demos of ML theory. Each one is a single script that prints its key numbers and saves a PNG to `out/`. They run on a laptop CPU with numpy, scipy and matplotlib, and the neural networks are written in plain numpy. | script | what it shows | |---|---| | `double_descent.py` | test error vs model size peaks at the interpolation threshold p = n, then falls again | | `high_dim.py` | random vectors are nearly orthogonal; ball volume sits near the surface; exponentially many near-orthogonal directions | | `diffusion_memorization.py` | a diffusion model with the *exact* score regurgitates its training set; imperfect scores generalise | | `randomized_response.py` | Warner's randomized response: deniable answers with unbiased population estimates, ε-DP, plus a poll-analysis CLI | | `certified_robustness.py` | Lipschitz-margin certificates and randomized smoothing, checked against a PGD attack | | `run_all.py` | runs everything (quick settings by default) | ## Install and run ```bash pip install -r requirements.txt python run_all.py # quick settings: ~70 s wall-clock python run_all.py --full # full settings: ~3 min (171 s measured) python double_descent.py # any single demo; most accept --quick python -m pytest -q # 18 tests, ~40 s ``` The timings are wall-clock on a 2-CPU machine that was shared with other jobs, using one process. Full-mode per-demo times were: double descent ≈ 48 s, high-dim ≈ 32 s, diffusion ≈ 37 s, randomized response < 1 s, certified robustness ≈ 50–80 s. --- ## 1. Double descent (`double_descent.py`) **Idea.** 12 noisy samples of sin(2πx) (inputs on a jittered grid, so they are random but never clumped), fit with p random Fourier features `cos(w·x + b)`. For p < n we use least squares. For p ≥ n we use the **minimum-norm interpolant** (the pseudo-inverse solution, which is also where gradient descent from zero converges). * At **p = n** there is exactly one interpolant. The feature matrix is square and ill-conditioned, so the fit has to wiggle wildly to pass through the noise, and the test error spikes. * For **p ≫ n** there are many interpolants, and the minimum-norm one is smooth, so the error falls again (Belkin, Hsu, Ma & Mandal, 2019). * A tiny ridge penalty (λ = 1e-3) removes the peak. **Output (full run; median over 100 feature draws):** ``` best under-parameterised: p=8 test MSE 0.0138 interpolation peak: p=12 test MSE 0.3165 (23x worse) far over-parameterised: p=2000 test MSE 0.0108 ridge 1e-3 near p = n: max test MSE 0.0165 (no peak) ``` Caveat: whether the second descent ends *below* the first minimum depends on the setup (noise level, feature bandwidth). With this configuration it does (0.0108 vs 0.0138), but this is not a universal law. Fixed parameters are `SCALE = 15`, noise sd 0.1, data seed 0. ## 2. High-dimensional geometry (`high_dim.py`) * For random unit vectors in R^d, `cos(u, v)` has mean 0 and **standard deviation exactly 1/√d**: | d | 10 | 100 | 1000 | 10000 | |---|---|---|---|---| | measured std | 0.3175 | 0.1010 | 0.0319 | 0.0101 | | 1/√d | 0.3162 | 0.1000 | 0.0316 | 0.0100 | * The fraction of a d-ball within ε of its surface is `1 − (1 − ε)^d`. This was checked by Monte Carlo for d ≤ 8; for example, d = 3, ε = 0.1 gives 0.2710 from the formula and 0.2716 from simulation. At d = 100, 99.997 % of the volume lies in the outer 10 % of the radius. * **Near-orthogonal packing** (all pairwise |cos| < 0.1). Only d vectors can be exactly orthogonal. For two random vectors, `cos² ~ Beta(½, (d−1)/2)`, so the chance that a pair violates the threshold, p_d, decays exponentially in d. By the birthday argument, a purely random set survives up to about `√(2 ln 2 / p_d)` vectors. That is about 30 at d = 1000 (simulated: 25), but about 7.7·10⁴ at d = 4000 and 10⁶ at d = 5000. A greedy search with 20 000 random tries keeps 972 vectors at d = 800, which is already more than d. The growth is exponential, but it starts late. ## 3. Diffusion memorisation (`diffusion_memorization.py`) **Idea.** For a training set {x_i}, the noised data distribution is the Gaussian mixture `p_σ = (1/N) Σ N(x_i, σ²I)`, and its score has a closed form: a softmax-weighted pull toward the training points. A variance-exploding reverse-diffusion sampler with this *exact* score samples from `p_{σ_min}`, which means it reproduces the training points. This is the "perfect model memorises" phenomenon behind training-data extraction from real diffusion models (Carlini et al., 2023). Two imperfect scores generalise instead: * a **blurred** score (evaluated at noise level `√(σ² + τ²)`), which is a kernel density estimate; * a **small MLP** (2×64 SiLU) trained by denoising score matching in numpy. Its smoothness bias spreads samples along the circle. **Output (20 training points on the unit circle, 2000 samples each; "copy" = within 0.03 of a training point):** ``` source median NN dist % copies median |r-1| fresh samples from true circle 0.1235 11.3% 0.0136 exact empirical score 0.0115 99.0% 0.0121 blurred score (tau=0.15) 0.1482 3.1% 0.1035 MLP score (DSM, numpy) 0.0845 17.5% 0.0154 ``` The MLP's samples stay on the circle (|r − 1| is about the same as for real data) and are mostly *not* copies. The blur spreads samples off the circle (|r − 1| ≈ 0.10): generalisation by smoothing is not the same thing as learning the manifold. Caveat: the MLP result depends on its capacity and training length. A bigger net trained much longer would move toward the exact (memorising) score. ## 4. Randomized response (`randomized_response.py`) **Idea (Warner, 1965).** Each respondent privately answers either "Am I in group A?" (with probability p) or the complement question. Then `P(yes) = λ = pπ + (1−p)(1−π)`, and the estimator `π̂ = (λ̂ − (1−p))/(2p−1)` is unbiased with variance `π(1−π)/n + p(1−p)/(n(2p−1)²)`. Every single answer is ε-locally differentially private with `ε = ln(p/(1−p))`. **Output (π = 0.3, 5000 simulated surveys per n):** ``` design eps n mean est sd (sim) sd (theory) n-inflation Warner p=0.6 0.41 1000 0.2989 0.0807 0.0788 29.6x Warner p=0.75 1.10 1000 0.2995 0.0310 0.0310 4.6x Warner p=0.9 2.20 1000 0.2999 0.0187 0.0187 1.7x honest answers sd = sqrt(0.21/1000) = 0.0145 ``` "n-inflation" is how many more respondents you need to match an honest survey's precision. That is the price of deniability. **Poll CLI:** ```bash python randomized_response.py make-example out/example_poll.csv --n 800 --pi 0.3 --p 0.75 python randomized_response.py poll out/example_poll.csv --p 0.75 --column answer ``` ``` answers used: 800 (yes 322, no 478; skipped 0 unreadable rows) design: truthful-question probability p = 0.75 -> local DP epsilon = 1.099 raw 'yes' rate: 0.4025 (NOT the answer: it is mixed with the complement question) debiased estimate pi_hat: 0.3050 (unbiased; may fall outside [0,1] by chance) 95% CI, exact (Clopper-Pearson mapped): [0.2366, 0.3748] 95% CI, normal approx: [0.2370, 0.3730] privacy cost: this survey is about as precise as ~200 direct answers (at pi ~ 1/2) ``` (The true sample fraction in that example file is 0.2825.) The exact interval maps a Clopper–Pearson interval for λ through the linear Warner transform, so it is conservative and valid. The CSV can be headerless or have a header (use `--column`), and yes/no, 1/0 and true/false are all accepted. ## 5. Certified robustness (`certified_robustness.py`) **Idea.** For a classifier `sign f(x)` where f is L-Lipschitz in ℓ2, no perturbation smaller than `|f(x)|/L` can flip the decision. That is a *certificate*, not an empirical claim. We train two numpy MLPs on two moons: a **wiggly** unconstrained one, and a **spectrally normalised** one (every weight matrix projected to spectral norm ≤ 1 after each Adam step, so that with ReLU, `L ≤ Π‖W_i‖₂ ≤ 1`). Each is then attacked with ℓ2 PGD. We also certify the wiggly net with **randomized smoothing** (Cohen, Rosenfeld & Kolter, 2019): `g(x) = majority vote of sign f(x + σZ)`, with certified radius `σ Φ⁻¹(p_A,lower)`. Here `p_A,lower` is a one-sided Clopper–Pearson bound from n = 10 000 Monte Carlo samples (α = 0.001). **Output (200 test points):** ``` model train acc test acc Lipschitz bound wiggly 1.000 0.960 2.01e+03 Lipschitz 0.955 0.985 1 PGD inside certified radii of the Lipschitz net: 0 violations out of 197 certified points (must be 0) smoothed wiggly net (sigma=0.25, n=10000, alpha=0.001): accuracy 0.980, abstain 0.000 radius Lip cert Lip PGD wig cert wig PGD smooth cert 0.00 0.985 0.985 0.960 0.960 0.980 0.10 0.830 0.945 0.000 0.890 0.925 0.20 0.235 0.900 0.000 0.735 0.805 0.30 0.000 0.815 0.000 0.575 0.620 0.50 0.000 0.500 0.000 0.190 0.170 ``` "cert" is a guaranteed lower bound on robust accuracy. "PGD" is an attack-based upper bound, and the truth lies in between. The wiggly net is empirically fairly robust, but its Lipschitz bound (~2000) certifies nothing. Smoothing turns the same network into a certifiable classifier. Note that it certifies the *smoothed* classifier g, not f. The product-of-spectral-norms bound is loose, which is why "Lip cert" falls well below "Lip PGD". ## Ideas for turning this into a tool * **Privacy-preserving survey kit**: a web form that runs the randomization in the respondent's browser (only the randomized bit is ever sent), plus the `poll` analysis as the backend. Useful for HR or climate surveys and for sensitive-behaviour questions. * **Memorisation auditor**: generalise the nearest-training-point metric into a script that takes generated samples and a training set (as embeddings) and flags near-copies, with the exact-score sampler as a calibrated "fully memorising" reference. * **Robustness certificate service**: given a small model's weights, compute Lipschitz upper bounds and randomized-smoothing certificates for incoming inputs, and return "certified radius" alongside each prediction. * **Model-size planner**: turn the double-descent script into a quick check that shows how close a regression or classification setup is to the interpolation threshold, and what ridge strength removes the peak. * **Embedding-capacity calculator**: the packing numbers in `high_dim.py` answer "how many near-orthogonal concepts fit in d dimensions at tolerance τ?" ## File layout ``` theory-demos/ ├── README.md ├── requirements.txt ├── pytest.ini ├── run_all.py ├── double_descent.py ├── high_dim.py ├── diffusion_memorization.py ├── randomized_response.py ├── certified_robustness.py ├── tests/ │ └── test_demos.py # cos std ≈ 1/√d, shell formula, DD peak/ridge, RR unbiasedness/CI coverage, │ # exact score vs finite differences, memorisation, Lipschitz bound validity, │ # PGD never succeeds inside certified radii, smoothing certificate check └── out/ # generated PNGs (+ example_poll.csv) ```