Optimization: how models learn

Training a model is an optimisation problem: pick the parameters that make the loss smallest, and everything else you touch — learning rate, batch size, $\lambda$, warmup, architecture — is a knob on that one problem. This module hands you the map of the machine: what each knob changes mathematically, which knobs are safe to turn, and which ones people turn by ritual without knowing what they do.

Cold questions — answer from memory, in writing

Model answers

Only open these after your own answers are written down.

1 · The problem statement

Strip a training script to its mathematics and one line remains. The parameters you search for are $\theta$; the loss of a single example is $L$; the model is $f(x;\theta)$; the penalty $\Omega$ discourages silly solutions.

Read it as two factions. The sum is the empirical risk: the average error on the $n$ samples you actually have. It pulls $\theta$ toward whatever reproduces the training data, including its noise. The second term is the regulariser: a penalty whose only job is to keep the solution tame — small weights, sparse weights, bounded complexity. $\lambda$ is the dial between the two factions, and it is set by you, not by the gradient.

That distinction is the whole parameter/hyperparameter split. Parameters $\theta$ are learned by the optimiser — they move every step. Hyperparameters are outside the gradient: $\lambda$, the learning rate $\eta$, the batch size, the number of layers, the number of epochs, $\beta$. Nothing in $J$ can tell you the right learning rate, because the learning rate is not a variable of $J$. Choosing them is a second, outer optimisation — usually by cross-validation or a schedule you write by hand (Module 05 covers the honest version of that loop).

2 · Convexity: the property that makes everything easy

A set is convex if the straight segment between any two of its points stays inside it. A function is convex if the same is true one level up: the segment joining two points on the graph stays above the graph. In words: the bowl holds water. Formally, for every $\alpha$ between 0 and 1:

The picture to sketch: draw a bowl, mark two points on its rim, draw the straight chord between them, and note the chord is always the pessimistic one — it over-estimates the function in between. The equivalent "epigraph" reading: if you shade everything above the graph, that shaded region is a convex set. Convexity of the function is convexity of a set, just lifted one dimension.

Why convex ⇒ local minima are global. Suppose $\theta^\star$ is a local minimum but not global, so some point $\theta'$ has $J(\theta') \lt J(\theta^\star)$. Take the segment between them: convexity says every point on it satisfies $J(\alpha\theta^\star + (1-\alpha)\theta') \le \alpha J(\theta^\star) + (1-\alpha)J(\theta')$, and the right side is strictly less than $J(\theta^\star)$ for any $\alpha \in (0,1)$. But those points are inside the neighbourhood of $\theta^\star$ — contradicting it being a local minimum. No dip, no hidden basin, no local trap.

Examples that matter in ML. Linear regression is convex: $J(w)=\frac{1}{n}\|Xw-y\|^2$ has Hessian $\frac{2}{n}X^\top X$, which is positive semidefinite because for any vector $v$, $v^\top X^\top X v = \|Xv\|^2 \ge 0$ — a Module 02 fact doing a Module 07 job. Logistic regression is convex too (its Hessian is $\frac{1}{n}X^\top W X$ with $W$ a diagonal matrix of positive weights — positive semidefinite for the same reason). A neural network is not convex: composing non-linearities with the parameters, not just with the input, destroys the property. Neither is matrix factorisation, low-rank approximation, or any problem where the unknown sits inside a product of matrices.

Convex combinations show up again immediately. A convex combination is a weighted average of points with non-negative weights summing to 1 — the $\alpha x + (1-\alpha)y$ of the definition. Momentum's velocity is a convex combination of past gradients; Adam's first moment $m_t$ is exactly a convex combination of all past gradients; the mixture interpretation of a regulariser interpolates between the data's answer and zero. When you see $\sum_k c_k v_k$ with $c_k \ge 0$ and $\sum c_k = 1$, read it as "a weighted blend of the $v_k$", and the result stays inside the hull of them.

3 · Gradient descent, derived and understood

Module 03 gave you the update rule; here it is derived, so it stops being a memorised line. Start from the first-order Taylor expansion — the local linear model of $J$ near $\theta$:

The second term is $\|\nabla J\|\|\Delta\|\cos\phi$, maximally negative when the step points opposite the gradient. So the best linear decrease is $\Delta = -\eta\nabla J$ with $\eta$ the step length — and substituting back gives the update rule by construction, not by fiat:

The learning-rate regimes, all four of which you will see in a career. Too small: the loss falls, but the linear model is re-evaluated far slower than necessary — you crawl, and a budget of steps gets you nowhere. Just right: each step buys a healthy fraction of the remaining error, and the loss falls geometrically. Too large but still stable: the iterates oscillate — you overshoot the valley floor, land on the other wall, get flung back. Past the stability limit: divergence — the loss climbs, doubles, hits overflow, and prints nan.

The complete, checkable result. Take the simplest possible convex function, $J(\theta)=\frac{1}{2}\theta^2$, whose gradient is $\theta$. The update is a multiplication by a constant:

There is nothing approximate about that statement — it is an equality, and you can check it by hand. $\eta = 0.9$ shrinks the error by $0.1$ per step; $\eta = 1.5$ flips the sign and shrinks by $0.5$ (zig-zagging but converging); $\eta = 2$ flips the sign with no shrinkage — a perfect oscillation forever; $\eta = 2.1$ multiplies the error by $1.1$ each step and runs away. Pass 3 sweeps $\eta$ and draws exactly this curve. The general statement replaces $1$ with the curvature: on a quadratic with Hessian eigenvalues in $[\lambda_{\min},\lambda_{\max}]$, the same one-number argument gives the stability threshold $\eta \lt 2/\lambda_{\max}$ — and Section 9's $\kappa=\lambda_{\max}/\lambda_{\min}$ is the ratio between the step you can afford and the step the flat directions want.

Convergence rates, stated not proved. On a strongly convex, smooth loss (Hessian eigenvalues in $[\mu, L]$ with $0 \lt \mu \le L$) and $\eta \le 1/L$, the error shrinks by a constant factor each step: the distance to the minimum is multiplied by at most $\left(1-\eta\mu\right)$ — linear (geometric) convergence, which means "halve the error every $c$ steps" forever. With convexity but no strong convexity you only get sublinear convergence: the objective gap falls like $O(1/k)$, so ten times more steps buys one extra digit. Non-convex losses get neither statement — only "the gradient norm gets small". And $\kappa = L/\mu$ is baked into the constant: a worse-conditioned problem converges more slowly no matter how carefully you pick $\eta$.

4 · Batch, mini-batch, stochastic

The gradient of the average is the average of the gradients. If $J(\theta)=\frac{1}{n}\sum_i L_i(\theta)$, then $\nabla J = \frac{1}{n}\sum_i \nabla L_i$ — which means you may trade exactness for cost:

Batch uses all $n$ samples: an exact gradient, zero noise, and a cost proportional to $n$ per step. Stochastic uses $B=1$: a single sample's gradient — biased in the large, but an unbiased estimator at every point, and $n$ times cheaper per step. Mini-batch is everything real training does: $B$ between about 32 and 512, mostly to make the estimate less noisy and to use the GPU, which is a matrix-multiply engine that idles below a few thousand elements.

The unit of progress is the epoch — one full pass over the $n$ samples. Compare methods by epochs as well as by steps, or you will rank a method that saw 40 times less data as faster. Shuffling at every epoch is not decoration: if batches are not i.i.d. samples, consecutive gradients are correlated, the effective batch size collapses, and the optimizer can chase a sequence of single-class or single-source batches. This bites hardest exactly where data arrives in order — time series, per-user records, files sorted by directory.

Practical sizes. Start at 32–128 for anything you iterate on quickly; go to 256–512 when the GPU is idle (utilisation below ~80%) and memory allows; go larger only when you have a schedule that scales with it, because the generalisation behaviour changes with $B$. Doubling the batch size means you can roughly double the learning rate and halve the epochs — the interesting part is the "roughly", which is empirical and problem-dependent, not a law.

5 · Momentum, Nesterov, RMSProp, Adam

Plain gradient descent reacts to each gradient as if it were the whole truth. All four optimisers below are answers to the same complaint: that estimate is noisy, badly scaled, and cheap to improve.

Momentum. Keep a velocity — an exponential moving average of past gradients — and step along it instead of the raw gradient:

The physical analogy is literal: $v$ is a velocity and $\beta$ is friction. On a narrow curved valley, plain GD zig-zags because the across-the-valley component flips sign every step while the along-the-valley component is tiny. Momentum averages away the flipping component (they cancel in the running average) and accumulates the consistent one (they add). The "look-back window" follows from the weights: $v_t = -\eta\sum_{k\ge0}\beta^k g_{t-k}$, and $\sum_{k\ge0}(1-\beta)\beta^k = 1$ normalises those weights — so the average lives over about $1/(1-\beta)$ steps. $\beta=0.9$ ≈ 10 steps, $\beta=0.95$ ≈ 20, $\beta=0.99$ ≈ 100. Two consequences worth knowing: a constant gradient produces a steady-state velocity $\eta/(1-\beta)\cdot\nabla J$ — so momentum does scale your effective step by $1/(1-\beta)$; and a changing gradient takes $1/(1-\beta)$ steps to be fully forgotten, which is why a large $\beta$ makes an optimizer sluggish to correct course (and why warmup exists, below).

Nesterov is momentum with the gradient evaluated where you are about to land: compute $g$ at $\theta+v$ instead of $\theta$. It is the same averaging with better lookahead, and it typically converges a few steps sooner. One sentence is enough — you will meet it as nesterov=True.

AdaGrad and RMSProp: per-parameter step sizes. The second complaint is scaling: parameters with rare, large gradients (the vocabulary of a sparse feature, say) need small steps, while frequent small-gradient parameters need large ones. AdaGrad divides each coordinate's step by the running root of squared gradients. RMSProp fixes AdaGrad's fatal flaw — its accumulator only ever grows, so the effective learning rate decays to zero and learning stops — by making the accumulation exponential rather than cumulative:

Adam = momentum (first moment $m$) + RMSProp (second moment $v$) + one correction. Both running averages start at zero, and with decay $\beta_1 = 0.9$ the estimate is heavily biased toward zero in early steps — the update would be far too small exactly when you know least. Dividing by the unbiased normaliser $1-\beta^t$ removes exactly that bias, because $\mathbb{E}[\beta^t m_0] = \beta^t$:

Defaults are $\beta_1=0.9$, $\beta_2=0.999$, $\varepsilon=10^{-8}$, $\eta=10^{-3}$ — and they are good defaults, which is most of why Adam won.

• AdamW — decoupled weight decay: apply $w \leftarrow (1-\eta\lambda)w$ to the weights directly, not by folding $\lambda$ into the gradient. With Adam, folded-in L2 gets divided by $\sqrt{s_t}$, so "the same $\lambda$" shrinks different parameters by different amounts. Decoupling makes weight decay mean what it says. This is the correct default for networks.
• Lion (2023) — update by the sign of the momentum, nothing else: one hyperparameter-free rule that often matches Adam at lower memory cost. Sign updates are a natural fit with quantisation.
• Sophia (2023) — a diagonal second-order (Hessian-based) preconditioner with clipping, which reacts to curvature rather than to gradient magnitude, and reports better sample efficiency on language models.
• 8-bit optimisers — store the state in 8-bit integers with block-wise quantisation, cutting optimizer memory by ~4× at almost no quality cost: the state is three vectors the size of the model, which matters when the model is large.

Schedules. The learning rate is a function of time, not a constant. Step decay: drop $\eta$ by 10× every fixed number of epochs — crude and effective. Cosine: decay $\eta$ from a peak to near zero along half a cosine wave — the modern default for fine-tuning, because it spends most of the run at a small step where fine detail gets settled. Warmup: ramp $\eta$ linearly from ~0 to the peak over the first 500–2000 steps. Warmup matters because of Section 9: early in training the gradients are uninformative and, crucially, the second moments $v_t$ have not been estimated yet — a few large early gradients can poison $v_t$ for the rest of the run (the running average with $\beta_2=0.999$ forgets on a 1000-step timescale), so the safest way to survive estimation error is to be small while it is worst.

6 · Second-order methods, honestly

Newton's method uses all three terms of the second-order Taylor expansion — value, slope, and curvature — by solving the local quadratic model exactly:

Why it is beautiful, in two sentences. First, scale invariance: $H^{-1}$ rescales each direction by its own curvature, so a valley that is $10^6$ times steeper in one direction than another gets equalised instead of demanding a compromise step size — it is a preconditioner computed from the true local geometry. Second, quadratic convergence near a minimum: once you are close enough that the quadratic model is accurate, one step lands on its minimum, which is the minimum; the error squares every iteration, so three or four steps buy full double precision.

Why it is not the default at scale, in two numbers. $H$ has $d^2$ entries — for $d = 10^9$ parameters that is $10^{18}$ numbers, far beyond any memory on Earth — and inverting or solving it costs $O(d^3)$, worse still. Even forming a matrix-vector product with $H$ requires second-order derivatives you would rather not store. So Newton's method lives where $d$ is small and the objective is smooth: the inner loops of classical scientific computing, small convex models, and — the part you will actually meet — as the idea behind every method below.

• Quasi-Newton (L-BFGS). Never form $H$. Instead keep the last $m$ (say 10) gradient and position changes and build a low-rank approximate inverse Hessian from them, implicitly, in $O(md)$ memory. L-BFGS is what scikit-learn reaches for on small convex problems, it needs no learning-rate tuning, and Pass 3 Block 4 shows it hitting the closed-form ridge solution to $10^{-9}$ in 24 iterations.
• Natural gradient. The clever reframe: rather than the geometry of the parameter space, use the geometry of the model distribution. Two parameter settings that barely change the predicted probabilities are close, however far apart the raw parameters are. The metric that says so is the Fisher information matrix $F = \mathbb{E}[\nabla\log p\,\nabla\log p^\top]$, and the natural gradient step is $F^{-1}\nabla J$ instead of $H^{-1}\nabla J$. The idea — precondition by a statistic of the model, not by the raw curvature — is the direct ancestor of modern variance-reduced and diagonal preconditioning methods (K-FAC, Shampoo, SOAP), and it explains why Adam's $1/\sqrt{v}$ feels like it works: it is a crude diagonal stand-in for exactly this.

7 · Constrained optimisation and Lagrange multipliers

So far the constraint was invisible: minimise $J(\theta)$ over all of $\mathbb{R}^d$. Real problems add requirements — weights bounded, budget spent, outputs summing to one, margin at least 1. The idea that makes constraints tractable is geometric, and it is the picture to keep:

Turn the picture into algebra with the Lagrangian: the objective with the constraint folded in as a term whose strength is a new unknown, the multiplier $\mu$. For an inequality constraint $g(\theta) \le 0$:

The KKT conditions are the full statement of optimality for a constrained problem. Stated plainly, at the solution $\theta^\star$:

Complementary slackness is the one people forget and the one that does the work: $\mu g = 0$ with both factors non-negative forces one of them to vanish. Either the constraint is active — you are sitting exactly on the boundary — and then $\mu$ can be positive, meaning the constraint is genuinely costing you; or it is inactive — the unconstrained optimum was already feasible — and then $\mu = 0$, and the constraint may be deleted without changing the answer.

The primal–dual story, one paragraph. For every constrained ("primal") problem there is a dual problem, obtained by asking: what is the best lower bound on $J(\theta^\star)$ that I can certify by choosing a $\mu$? Minimising over $\theta$ with $\mu$ fixed gives a bound; the dual problem maximises that bound over $\mu$. Weak duality always holds — the dual value never exceeds the primal optimum, so any dual-feasible $\mu$ certifies that you cannot do better than its bound. Strong duality — no gap at all, dual optimum = primal optimum — holds for convex problems under mild conditions (Slater's: a strictly feasible point exists). That is why convexity buys more than convergence: it guarantees a certificate. And it is what makes the kernel trick possible, next.

Payoff (a): ridge regression, with the algebra. Constrained least squares with a budget on the weight norm, and its unconstrained twin:

Write the Lagrangian, differentiate, and solve. Four steps, all mechanical — but they are the correspondence between "$\lambda$ in scikit-learn" and "a budget you could actually state to a stakeholder":

Compare with the ridge objective $\frac{1}{n}\|Xw-y\|^2 + \lambda\|w\|^2$, whose stationarity condition is $(X^\top X + n\lambda I)w = X^\top y$ — identical, with $\lambda = \mu$. So the constrained and penalised forms trace out the same curve of solutions; $\lambda$ and the radius $t$ are two dials on one control. Read off the correspondence: if $\|w_{\text{OLS}}\|^2 \le t$, the constraint is slack, $\mu = 0$ and $\lambda = 0$ — no shrinkage at all; as $t$ shrinks, $\mu$ (and hence $\lambda$) must grow to push the solution in; at $t \to 0$ you shrink toward zero. Complementary slackness $\mu(\|w\|^2 - t) = 0$ says exactly this: for every $t$ below the unconstrained norm, the budget binds and the solution sits on the sphere of radius $\sqrt{t}$.

Payoff (b): the SVM and why its dual buys kernels. The maximum-margin classifier is a constrained quadratic program: maximise the margin $\|w\|^2$ subject to $y_i(w^\top x_i + b) \ge 1$ for every sample. Its dual has one multiplier $\mu_i$ per training sample, and two features of that dual explain everything about SVMs. First: samples with $\mu_i = 0$ contribute nothing to the solution — only those with $\mu_i \gt 0$ matter, and they are exactly the points on or inside the margin: the support vectors. The model is defined by a handful of points, not by the dataset. Second: the dual objective and the solution depend on the data only through inner products $x_i^\top x_j$ — so you may replace them by a kernel $k(x_i,x_j) = \phi(x_i)^\top\phi(x_j)$ and implicitly optimise in a feature space you never wrote down. That is the kernel trick, and it is available because the problem was convex and its dual was tractable.

Payoff (c): projection and resource allocation. When the feasible set is simple, the cheapest solver is projected gradient descent: take the step, then project back onto the set — clip a parameter to $[-c,c]$, renormalise a weight vector to unit norm, mask a budget. When the set is not simple, the multiplier $\mu$ becomes a price: allocate resource where $\nabla J$ is largest per unit of constraint, which is how capacity, rate limits and budgets get optimised in systems too. The vocabulary transfers directly — "dual", "multiplier", "price", "shadow cost" — from optimisation into the systems work you already know.

8 · Regularisation as an optimisation choice

Section 1 introduced $\lambda\Omega(\theta)$ as a term. Here are the three versions you will actually use, and each one is a different answer to the same question: which solution among the many that fit the data equally well do you prefer?

L2 pulls every coefficient toward zero in proportion to its size, so big coefficients shrink hardest. Because the gradient $2\lambda w$ vanishes as $w \to 0$, a coefficient gets close to zero but never exactly there — the update is $w \leftarrow (1-2\eta\lambda)w$ plus the data term, an exponential decay that only reaches zero in the limit. (Pass 3 Block 3 prints the smallest surviving L2 coefficient: $9 \times 10^{-4}$ — small, and not zero.)

L1 produces sparse solutions: exact zeros, and therefore feature selection. Two arguments, both worth owning. The subgradient argument: $|w|$ is not differentiable at 0 — its subdifferential there is the whole interval $[-1,1]$ — so the data gradient must exceed $\lambda$ in magnitude before the coefficient can leave zero. Weak signals never clear that hurdle and are set exactly to zero. The geometric argument: draw the contours of the data-fitting term and the constraint set $\|w\|_1 \le t$, which is a diamond with its corners on the axes. When a tilted ellipse first touches the diamond, it very often touches at a corner — and a corner is a point where all coordinates but one are zero. The $\ell_2$ ball is a circle with no corners, so contact is generically at a non-axis point and no coefficient is exactly zero.

Getting the gradient wrong here is the standard bug: $\operatorname{sign}(w)$ is a subgradient, and plain subgradient steps on L1 oscillate around zero rather than settling on it. The correct simple algorithm is proximal gradient (ISTA): take the gradient step for the smooth part, then apply the proximal operator of the $\ell_1$ term, which is a soft threshold:

• Elastic net — both penalties: $\lambda_1\|w\|_1 + \lambda_2\|w\|^2$. Sparsity with stable selection among correlated features, where pure L1 picks one of a correlated group at random and L2 keeps all of them.
• Early stopping — regularise by not finishing. The trajectory itself limits effective capacity: after $k$ steps an iterative solver has mostly fitted the directions with large eigenvalues (the ones it can reach fastest), and stopping before convergence acts like an $\ell_2$ penalty with $\lambda$ decreasing in $k$. It is implicit regularisation — the optimiser is doing the constraint without writing it down.
• Dropout — regularise by randomly deleting units each step: the surviving sub-network must fit without any single unit, which pushes against brittle co-adaptation. Data augmentation — regularise by making the training distribution fatter: every rotated or cropped image is a new constraint on the solution. Same idea, different faces: shrink the hypothesis space, or fill the data space.

9 · Numerical reality: conditioning, and why your model might be untrainable

Sometimes the model is fine, the loss is right, the code is correct — and training still crawls. The reason is usually arithmetic, not ideas.

$\kappa = 1$ means every direction is equally curved — a circular bowl, one learning rate fits all. Large $\kappa$ means the loss is a long narrow valley: the gradient in the steep direction is $\kappa$ times the one in the flat direction, so any $\eta$ that survives the steep wall crawls along the valley floor. The number to remember is the square law: $\kappa(X^\top X) = \kappa(X)^2$, because the eigenvalues of $X^\top X$ are the squares of the singular values of $X$. Two features on scales $1$ and $100$ give $\kappa(X) \approx 100$ and $\kappa(X^\top X) \approx 10^4$ — your optimisation problem is ten thousand times worse-scaled than your data looked. The second symptom is numerical: solving $X^\top X w = X^\top y$ loses roughly $\log_{10}\kappa$ decimal digits to rounding, which is why the normal equations are avoided in serious linear algebra code (Module 02).

The cheapest fix in all of ML: scale the features. Standardising each column to zero mean and unit variance makes the singular values cluster near 1, $\kappa \approx 1$, and one learning rate serves every coordinate. Pass 3's exercise measures a 5000× speed-up from scaling alone — no algorithm change, no tuning, one line of preprocessing.

Deep networks add a second failure of the same family: gradients along the chain. Backpropagation multiplies one factor per layer (Module 03, Section 7), and the product of many numbers that are each a bit below 1 is a number that is essentially 0 — vanishing gradients, where early layers never learn. The product of many a bit above 1 is overflow — exploding gradients, where one step destroys the weights and prints nan. Same mechanism, opposite sign, and the mitigations are all attempts to keep those per-layer factors near 1:

• Initialisation. Weights must be drawn with variance matched to the layer's width, so signals neither blow up nor die in the forward pass and the gradients behave in the backward one. Xavier/Glorot: variance $\propto 1/\text{fan-in}$, derived for linear and sigmoid/tanh layers. He: variance $\propto 2/\text{fan-in}$ — twice as large, to compensate for the ReLU units that it zeroes out.
• Normalisation layers (BatchNorm, LayerNorm): re-centre and re-scale activations to a controlled distribution each step, so the next layer always sees inputs of sane magnitude — which keeps gradients in the non-saturating range.
• Residual connections: $y = F(x) + x$ gives the gradient a clean additive path — $\partial L/\partial x = \partial L/\partial y + \ldots$ — so the signal can flow backwards through hundreds of layers without being multiplied by anything.
• Gradient clipping: rescale the whole gradient to norm at most $c$ when it spikes. It treats the symptom, deliberately, because sometimes the symptom is all you need.

10 · The map

Every knob you turn in training, translated into what it changes mathematically. Keep this table: it is the answer to "what will halving this number actually do?"

Block 1 · GD vs SGD vs momentum on one anisotropic quadratic

Same loss, three optimisers. The quadratic is rotated so the valley is not axis-aligned — that is what makes plain GD zig-zag, and what makes the difference between the methods visible.

import numpy as np
import matplotlib.pyplot as plt

# anisotropic quadratic, rotated so plain GD visibly zig-zags
th = np.deg2rad(30)
R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
A = R @ np.diag([1.0, 50.0]) @ R.T          # eigenvalues 1 and 50 -> kappa = 50
J  = lambda p: 0.5 * p @ A @ p
gJ = lambda p: A @ p
print("eigs:", np.round(np.linalg.eigvalsh(A), 6), " kappa:", np.linalg.cond(A))

# SGD data: n per-sample quadratics whose MEAN is exactly J
rng = np.random.default_rng(0)
n = 40
a = rng.uniform(0.6, 1.4, n); a *= 1.0 / a.mean()   # rescale so the means hit exactly
b = rng.uniform(40.0, 60.0, n); b *= 50.0 / b.mean()
Ai = np.array([R @ np.diag([a[i], b[i]]) @ R.T for i in range(n)])
print("mean per-sample Hessian == A:", np.allclose(Ai.mean(0), A))

def run(kind, eta, steps=20000, tol=1e-6, seed=0):
    p, v = np.array([2.0, -1.5]), np.zeros(2)
    r = np.random.default_rng(seed)
    curve, hit = [J(p)], None
    for k in range(1, steps + 1):
        if kind == "momentum":
            v = 0.9 * v - eta * gJ(p); p = p + v
        elif kind == "gd":
            p = p - eta * gJ(p)
        else:                                   # sgd: one sample's gradient
            p = p - eta * (Ai[r.integers(n)] @ p)
        curve.append(J(p))
        if hit is None and J(p) < tol: hit = k
    return np.array(curve), hit

for kind, eta in [("gd", 0.02), ("sgd", 0.005), ("momentum", 0.02)]:
    curve, hit = run(kind, eta)
    print(f"{kind:9s} eta={eta:<6} iterations to J<1e-6: {hit}")
    # gd 324 | sgd 1293 (~32 epochs of 40 samples) | momentum 141

Three honest readings of that line. Momentum needs 141 steps where plain GD needs 324 — the same $\eta$, roughly 2.3× fewer iterations, because the consistent valley-floor direction accumulates while the across-valley zig-zag cancels. SGD needs 1293 steps — but each of its steps used one sample where the others used 40, so in data terms SGD used 1293 samples to GD's 12960: the trade of Section 4, quantified. And note what the tolerance means: SGD's noise keeps it from settling far below $10^{-6}$, which is the noise floor — the same feature Section 4 calls implicit regularisation, seen as a floor on the achievable training loss.

Block 2 · A learning-rate sweep, against the exact threshold

import numpy as np
import matplotlib.pyplot as plt

# 1-D quadratic: J(theta) = 0.5 theta^2  ->  theta_{k+1} = (1-eta) theta_k
etas = np.logspace(-3, np.log10(4.0), 400)
final1 = np.abs(0.5 * (1 - etas) ** 100)          # closed form after 100 steps
thr1 = etas[final1 < 0.5].max()
print("1-D : largest grid eta still converging:", round(thr1, 4), "(theory: 2.0)")

# same experiment on a 2-D anisotropic quadratic: threshold set by the curvature
th = np.deg2rad(30); R = np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
A = R @ np.diag([1.0, 50.0]) @ R.T
lam_max = np.linalg.eigvalsh(A).max()
p0 = np.array([2.0, -1.5])

def gd2d(eta, steps=100):
    p = p0.copy()
    with np.errstate(over='ignore', invalid='ignore'):   # diverging eta overflows on purpose
        for _ in range(steps): p = p - eta * (A @ p)
        return 0.5 * p @ p

final2 = np.array([gd2d(e) for e in etas])
thr2 = etas[final2 < 0.5 * p0 @ A @ p0].max()
print("2-D : largest grid eta still converging:", round(thr2, 4),
      "(theory 2/lambda_max =", round(2 / lam_max, 4), ")")

plt.semilogx(etas, np.maximum(final1, 1e-20), label="1-D: J = 0.5 theta^2")
plt.semilogx(etas, np.maximum(final2, 1e-20), label="2-D: lambda_max = 50")
plt.axvline(2.0, ls='--', c='k'); plt.axvline(2 / lam_max, ls=':', c='r')
plt.xlabel("learning rate eta"); plt.ylabel("J after 100 steps")
plt.legend(); plt.show()

for e in [0.1, 1.0, 1.5, 1.9, 2.1]:                # the exact factor, five times
    print(f"eta={e:<4} |1-eta| = {abs(1-e):.2f}  -> {'converges' if abs(1-e) < 1 else 'DIVERGES'}")

The measured thresholds land where the theory puts them: $1.973$ against a theoretical $2.0$ (the grid's largest point below 2), and $0.0396$ against $2/\lambda_{\max} = 0.04$. The second number is the important one — the same method on the same kind of problem has a different threshold once the curvature changes. That is why "the learning rate" is not a property of an algorithm: it is a property of the algorithm and the loss, joined by $\lambda_{\max}$. It is also why Adam's per-parameter division by $\sqrt{v_t}$ helps so much: it removes the dependence on scale, direction by direction.

Block 3 · L1 vs L2: what "sparsity" actually means

A synthetic regression with a genuinely sparse truth: 5 of 60 coefficients are non-zero. Fit it twice, once with each penalty, and count exact zeros.

import numpy as np
from sklearn.linear_model import Lasso, Ridge

rng = np.random.default_rng(1)
n, d = 200, 60
X = rng.normal(size=(n, d))
X = (X - X.mean(0)) / X.std(0)                       # standardise: one step size for all
w_true = np.zeros(d)
support = rng.choice(d, 5, replace=False)
w_true[support] = [2.0, -3.0, 1.5, -1.0, 2.5]
y = X @ w_true + 0.3 * rng.normal(size=n)

L = np.linalg.eigvalsh(2 / n * X.T @ X).max()        # Lipschitz constant of the gradient
step = 1.0 / L
grad = lambda w: (2 / n) * X.T @ (X @ w - y)
soft = lambda z, t: np.sign(z) * np.maximum(np.abs(z) - t, 0.0)
lam = 0.15

w2 = np.zeros(d)                                     # L2: plain gradient steps
for _ in range(4000):
    w2 -= step * (grad(w2) + 2 * lam * w2)

w1 = np.zeros(d)                                     # L1: proximal gradient (ISTA)
for _ in range(4000):
    w1 = soft(w1 - step * grad(w1), step * lam)

print("L1 zeros:", d - np.count_nonzero(w1), " support found:", np.sort(np.flatnonzero(w1)))
print("L2 zeros:", d - np.count_nonzero(w2), " smallest |w|:", f"{np.abs(w2).min():.4f}")
lasso = Lasso(alpha=lam / 2, fit_intercept=False, max_iter=50000, tol=1e-10).fit(X, y)
ridge = Ridge(alpha=n * lam, fit_intercept=False).fit(X, y)
print("max |w_L1 - Lasso|:", f"{np.abs(w1 - lasso.coef_).max():.2e}",
      " max |w_L2 - Ridge|:", f"{np.abs(w2 - ridge.coef_).max():.2e}")

Printed: 55 exact zeros from L1, and the five survivors are exactly the true support $\{4, 32, 37, 40, 56\}$ — a clean recovery. Zero exact zeros from L2: the smallest coefficient is $9\times10^{-4}$, small and not null, exactly as Section 8's argument predicts ($2\lambda w$ vanishes as $w$ does). Both match scikit-learn to $5\times10^{-12}$ or better, which also verifies the scaling correspondence between the two objectives: Lasso minimises $\frac{1}{2n}\|Xw-y\|^2 + \alpha\|w\|_1$, so our $\lambda$ is its $2\alpha$; Ridge minimises $\|Xw-y\|^2 + \alpha\|w\|^2$, so its $\alpha$ is our $n\lambda$. Get those factors wrong and your "same" model silently isn't — a real and common bug.

Block 4 · When a quasi-Newton method is the right tool

import numpy as np
from scipy.optimize import minimize

# the ridge problem of Block 3: min_w (1/n)||Xw-y||^2 + lam||w||^2
rng = np.random.default_rng(1)
n, d = 200, 60
X = rng.normal(size=(n, d)); X = (X - X.mean(0)) / X.std(0)
w_true = np.zeros(d); support = rng.choice(d, 5, replace=False)
w_true[support] = [2.0, -3.0, 1.5, -1.0, 2.5]
y = X @ w_true + 0.3 * rng.normal(size=n)
lam = 0.15

f  = lambda w: np.mean((X @ w - y) ** 2) + lam * w @ w
gf = lambda w: (2 / n) * X.T @ (X @ w - y) + 2 * lam * w

w_cf = np.linalg.solve(X.T @ X + n * lam * np.eye(d), X.T @ y)   # closed form

res = minimize(f, np.zeros(d), jac=gf, method='L-BFGS-B',
               options=dict(maxiter=2000, ftol=1e-16, gtol=1e-12))
print(f"L-BFGS: {res.nit} iterations, {res.nfev} gradients,"
      f"  max |w - closed form| = {np.abs(res.x - w_cf).max():.2e}")

L = np.linalg.eigvalsh(2 / n * X.T @ X + 2 * lam * np.eye(d)).max()
w = np.zeros(d); k = 0
while np.abs(w - w_cf).max() > 1e-8 and k < 100000:   # plain GD to the same accuracy
    w -= (1 / L) * gf(w); k += 1
print(f"plain GD: {k} iterations,  max |w - closed form| = {np.abs(w - w_cf).max():.2e}")

L-BFGS reaches the closed-form answer to $2.7\times10^{-9}$ in 24 iterations and 35 gradient evaluations, where tuned plain gradient descent needs 83 — about 3.5× more, on a problem whose condition number is only 5.74. Scale that gap up to a badly conditioned problem and the ratio explodes. The lesson is not "always use L-BFGS"; it is the ordering: on a small convex problem with a closed form, the closed form wins outright (one solve, no iterations); when there is no closed form and the problem is small enough to hold its curvature implicitly, quasi-Newton wins; at billion-parameter scale neither is even representable, which is where Sections 5's Adam and Section 6's diagonal preconditioners come back in.

Multiple choice — instant feedback

Open questions — answer aloud, without notes

Re-answer these in 1 / 3 / 7 days — cold, no notes. If you cannot re-derive the ridge correspondence or the $0 \lt \eta \lt 2$ condition on the spot, they are not yours yet.