← Back to the Applied Science guide

Optimization

Every model you train is the answer to an optimization problem. Know the theory, and you can tell why a run diverged, stalled or overfit.

Interviewers rarely ask you to prove a rate. They ask why training blew up at step 200. They ask why the loss plateaued, or why a bigger batch hurt. The answers come from a small set of ideas. These are curvature, noise, conditioning and constraints. This page builds those ideas from the ground up. It goes deeper than the optimizer comparison in Deep Learning. The goal is to know where each rule comes from. Then you can reason about a case you have never seen.

Throughout, we minimize a loss f(θ) over parameters θ ∈ ℝ^d. The gradient is ∇f and the Hessian is ∇²f. The best value is f^* at a minimizer θ^*.

1. Convex sets and functions

The definitions

A set C is convex if for all x, y ∈ C and t ∈ [0, 1], the point tx + (1 − t)y is also in C. A function f is convex if its domain is convex and:

f(t x + (1 - t) y)  ≤  t f(x) + (1 - t) f(y)        for all t in [0, 1]

For smooth functions there are two easier tests. They say the same thing in different words.

• First-order test. f(y) ≥ f(x) + ∇f(x)^T(y − x). The tangent plane is a global lower bound. So ∇f(x) = 0 proves x is a global minimum.
• Second-order test. ∇²f(x) ≽ 0 everywhere. The Hessian is positive semidefinite. The surface curves up or stays flat in every direction.

Strong convexity and smoothness

Convexity alone says the bowl has no false bottoms. It says nothing about how steep or flat the bowl is. Two numbers pin that down.

• μ-strongly convex. f(y) ≥ f(x) + ∇f(x)^T(y − x) + (μ/2)‖y − x‖². Equivalently ∇²f ≽ μI. The bowl curves up at least as fast as a parabola. The minimizer is unique.
• L-smooth. ‖∇f(x) − ∇f(y)‖ ≤ L‖x − y‖. Equivalently ∇²f ≼ LI. The gradient cannot change too fast. This gives the quadratic upper bound f(y) ≤ f(x) + ∇f(x)^T(y − x) + (L/2)‖y − x‖².
• Condition number. κ = L/μ. It is the ratio of the steepest to the flattest curvature. For a quadratic it is λ_max/λ_min of the Hessian.

So f is trapped between two parabolas. One opens with width μ, and one opens with width L. Almost every rate on this page is a function of κ.

A useful relaxation: the PL condition

The Polyak–Łojasiewicz (PL) condition is weaker than strong convexity:

½ ‖∇f(θ)‖2  ≥  μ ( f(θ) - f* )

It says a small gradient means you are close to optimal in value. It does not need convexity. Gradient descent still converges linearly under PL. Wide, over-parameterized networks satisfy PL near their initialization. That is one theory for why they train so reliably.

Why it matters in practice

• Model choice. A convex model like logistic regression gives the same weights on every run. That makes it a stable baseline and easy to debug.
• Regularization fixes conditioning. Adding (λ/2)‖θ‖² makes the loss λ-strongly convex. It lowers κ and speeds up solvers.
• Feature scaling. Unscaled features stretch the Hessian. Standardizing them shrinks κ, often by orders of magnitude.

Interview check

• Is logistic regression strongly convex? Not in general. On separable data the loss keeps falling as ‖w‖ grows, so no minimizer exists. Add L2 and it becomes strongly convex.
• Is a neural net with one hidden layer convex in its weights? No. Swapping two hidden units gives the same function, so minima come in symmetric copies. The average of two copies is usually worse. That breaks convexity.
• What does κ = 1 mean? The bowl is perfectly round. Gradient descent with step 1/L reaches the minimum in one step.

2. Gradient descent and its rates

θk+1 = θk - η ∇f(θk)

The descent lemma

Plug the GD step into the L-smooth upper bound with y = θ_k+1:

f(θk+1)  ≤  f(θk) - η (1 - Lη/2) ‖∇f(θk)‖2

The loss must drop whenever 0 < η < 2/L. The guaranteed drop is largest at η = 1/L. There it is (1/2L)‖∇f‖². This one line is the root of the classic rule. Set the step to one over the largest curvature.

The three rates

• Smooth, non-convex. Sum the descent lemma over k steps. You get min_i≤k ‖∇f(θ_i)‖² ≤ 2L(f(θ₀) − f^*)/k. GD finds a near-stationary point at rate O(1/k). It may be a saddle.
• Smooth, convex. f(θ_k) − f^* ≤ L‖θ₀ − θ^*‖²/(2k). This is sublinear, O(1/k). To halve the error you double the steps.
• Smooth, μ-strongly convex. f(θ_k) − f^* ≤ (1 − 1/κ)^k(f(θ₀) − f^*). This is linear convergence. You need about κ log(1/ε) steps to reach error ε.

Why the condition number controls everything

Take a quadratic f(θ) = ½θ^TAθ. In the eigenbasis of A, each coordinate evolves on its own:

θi(k+1) = (1 - η λi) θi(k)

Each direction shrinks by the factor |1 − ηλ_i|. The steep direction needs η < 2/λ_max, or it blows up. But then the flat direction shrinks by only 1 − 2λ_min/λ_max per step. The best fixed step is η = 2/(L + μ). It gives the factor (κ − 1)/(κ + 1). With κ = 1000, that is 0.998 per step. You need thousands of steps to make progress along the flat valley.

Line search

You rarely know L. Backtracking line search finds a safe step on the fly. Start with a large η. Shrink it by a factor β until the Armijo condition holds:

f(θ - η∇f)  ≤  f(θ) - c η ‖∇f‖2        with c in (0, 1), often 1e-4

Full-batch solvers like L-BFGS use line search. Deep learning does not. Each extra loss evaluation costs a forward pass on noisy minibatches. So deep learning uses schedules instead.

Why it matters in practice

• Divergence diagnosis. If loss jumps to NaN in the first steps, η exceeded about 2/λ_max. Lower the rate or add warmup.
• Slow plateaus. A loss that falls fast then crawls often means high κ. Normalization layers and feature scaling attack κ directly.

Interview check

• Why is 2/L the stability limit? On the steepest eigen-direction, the update factor is 1 − ηL. Past η = 2/L its magnitude exceeds one, so errors grow each step.
• Does GD converge on a non-convex function? To a stationary point, yes, at rate O(1/k) in gradient norm. Nothing promises a global minimum.
• How does L2 regularization change the rate? It adds λ to every Hessian eigenvalue. So κ becomes (L + λ)/(μ + λ), which is smaller. Convergence gets faster.

3. Newton and quasi-Newton methods

The update

The second-order Taylor model is f(θ + p) ≈ f(θ) + ∇f^Tp + ½p^T∇²f p. Set its gradient in p to zero:

θk+1 = θk - [∇2f(θk)]-1 ∇f(θk)

• Quadratic convergence. Near a minimum with a well-behaved Hessian, ‖θ_k+1 − θ^*‖ ≤ C‖θ_k − θ^*‖². The number of correct digits doubles each step.
• Affine invariance. Rescale or rotate the parameters and Newton takes the same path. So κ does not matter to it. GD has no such property.
• Cost. Forming the Hessian is O(nd²). Solving the system is O(d³). Memory is O(d²). That is fine for d = 1,000 and hopeless for d = 10⁹.

Second-order intuition

In the Hessian eigenbasis, Newton divides each gradient component by its curvature λ_i. Steep directions get small steps. Flat directions get big steps. That is exactly the cure for the zigzag in the GD figure. Every adaptive optimizer is an attempt to get this effect cheaply.

Far from the optimum the quadratic model can be poor. Damped Newton takes θ − t[∇²f]^-1∇f with t from line search. For logistic regression, Newton is the classic IRLS algorithm (iteratively reweighted least squares).

Newton on regularized logistic regression

import numpy as np

rng = np.random.default_rng(0)
n, d = 500, 5
X = rng.normal(size=(n, d)); w_true = rng.normal(size=d)
y = (rng.random(n) < 1 / (1 + np.exp(-X @ w_true))).astype(float)
lam = 1e-2

def grad_hess(w):
    p = 1 / (1 + np.exp(-X @ w))
    g = X.T @ (p - y) / n + lam * w
    H = (X * (p * (1 - p))[:, None]).T @ X / n + lam * np.eye(d)
    return g, H

w = np.zeros(d)
for k in range(8):
    g, H = grad_hess(w)
    print(k, "%.2e" % np.linalg.norm(g))
    w = w - np.linalg.solve(H, g)       # Newton step
# gradient norm: 3.3e-01, 7.6e-02, 1.6e-02, 1.1e-03, 6.3e-06, 2.1e-10, 3.0e-17

Watch the exponent. It goes −3, −6, −10, −17. That is quadratic convergence. GD would need hundreds of steps for the same accuracy.

Quasi-Newton: BFGS

Quasi-Newton methods build a Hessian estimate from gradients alone. Define the step s_k = θ_k+1 − θ_k and the gradient change y_k = ∇f_k+1 − ∇f_k. A true Hessian roughly satisfies ∇²f s_k ≈ y_k. So we ask the new estimate B_k+1 to satisfy the secant condition B_k+1s_k = y_k. BFGS picks the closest such matrix to the old one. It updates the inverse H ≈ [∇²f]^-1 directly:

ρk = 1 / (ykT sk)
Hk+1 = (I - ρk sk ykT) Hk (I - ρk yk skT) + ρk sk skT

• The update is rank two and costs O(d²), not O(d³).
• H stays positive definite if y_k^Ts_k > 0. A Wolfe line search guarantees this curvature condition.
• Convergence is superlinear near the optimum. It is slower than Newton but far faster than GD.

L-BFGS

L-BFGS never stores H. It keeps the last m pairs (s, y), often m = 5 to 20. A "two-loop recursion" computes H∇f from those pairs in O(md) time and memory. That makes quasi-Newton work for millions of parameters.

Why it matters in practice

• Default for convex models. scikit-learn's LogisticRegression uses lbfgs by default. Full-batch L-BFGS often beats SGD on medium tabular data.
• Not for minibatches. Noisy gradients corrupt y_k, so the curvature pairs become garbage. That is why deep nets rarely use L-BFGS.
• Cheap curvature in deep nets. K-FAC and Shampoo approximate the Hessian or Fisher with Kronecker factors per layer. They aim for Newton-like steps at a manageable cost.

Interview check

• Why not use Newton for a 7B model? The Hessian has about 5 × 10¹⁹ entries. You cannot store it, let alone invert it.
• What is the secant condition? B_k+1s_k = y_k. The new curvature estimate must explain the gradient change seen on the last step.
• How do Gauss-Newton and natural gradient relate to Newton? Both replace the Hessian with a positive semidefinite matrix. Gauss-Newton uses J^TJ. Natural gradient uses the Fisher information. For log-loss models with canonical links, the two match.

4. Stochastic gradient descent theory

The setup

The loss is an average f(θ) = (1/n)Σ_i f_i(θ). Draw a minibatch S of size B. Use the minibatch gradient:

g = (1/B) Σi∈S ∇fi(θ)        E[g] = ∇f(θ)        Cov[g] ≈ Σ(θ) / B

Here Σ is the covariance of per-example gradients. The estimate is unbiased. Its variance falls as 1/B. Its standard error falls as 1/√B. So a 4× bigger batch only halves the noise.

What the noise does to convergence

Add noise to the descent lemma. With variance bound σ², one step gives:

E[f(θk+1)]  ≤  f(θk) - η(1 - Lη/2) ‖∇f‖2 + (Lη2/2) σ2/B

The last term never goes away with a fixed η. Progress stops once ‖∇f‖² is about Lησ²/B. So SGD with a constant step converges to a noise ball around the optimum. Its radius is proportional to ησ²/B.

• Robbins–Monro conditions. To converge exactly, the steps need Ση_t = ∞ and Ση_t² < ∞. The first lets you travel any distance. The second kills the noise.
• Convex rate. With η_t ∝ 1/√t, E[f(θ̄_T)] − f^* = O(1/√T). Here θ̄ is the averaged iterate.
• Strongly convex rate. With η_t = 1/(μt), the gap is O(1/(μT)). This is far slower than GD's linear rate. But each step costs 1/n as much.
• Polyak averaging. Averaging the iterates cancels noise. It reaches the optimal statistical rate with larger steps. EMA of weights is the deep learning version.

Noise scale and batch size

Model SGD as a stochastic differential equation. Then the noise "temperature" scales like η/B. Two runs with the same ratio of η to B behave alike. That gives the linear scaling rule. When you multiply the batch by k, multiply the learning rate by k.

• Goyal et al. trained ResNet-50 on ImageNet with batch 8,192 in one hour this way. They needed a 5-epoch warmup to survive the first steps.
• The rule breaks past the critical batch size. One estimate is B_noise = tr(Σ)/‖∇f‖². Below it, doubling B halves the steps needed. Above it, extra batch buys almost nothing.
• The critical batch grows as loss falls, because the signal shrinks faster than the noise. That is why LLM runs often ramp batch size up during training.
• Adam is closer to square-root scaling. Its step size is already normalized by gradient scale.

Why it matters in practice

• Scaling to more GPUs. Data parallelism raises B. Rescale η and add warmup, or the run diverges or underfits.
• Final LR decay. The last drop in loss at the end of a schedule is the noise ball shrinking. It is not new learning.
• Implicit regularization. SGD noise pushes toward wide, flat regions. Very large batches lose some of this, which can hurt test accuracy.

Interview check

• Why does SGD with a constant step not converge exactly? The gradient noise has non-zero variance at the optimum. Each step kicks the iterate by about ησ, so it jitters in a ball.
• You double the batch. What happens to the gradient variance and the step? Variance halves. Under the linear scaling rule you double η, keeping η/B fixed.
• When is the linear scaling rule wrong? Past the critical batch size, early in training without warmup, and for adaptive optimizers like Adam.

5. Momentum and Nesterov

Heavy ball (Polyak)

vk+1 = β vk - η ∇f(θk)
θk+1 = θk + vk+1

The physics view. Think of a ball with mass m and friction γ on the surface f. Its motion is mθ̈ + γθ̇ + ∇f(θ) = 0. Discretize it and you get heavy ball. Here β plays the role of 1 minus friction. Plain GD is the massless limit, where velocity equals force. A massless ball stops dead in a flat valley. A heavy ball keeps rolling. It also overshoots, which is the price of speed.

The quadratic analysis. On a quadratic with curvature between μ and L, choose:

η = 4 / (√L + √μ)2        β = ( (√κ - 1) / (√κ + 1) )2

The error then shrinks by (√κ − 1)/(√κ + 1) per step. That needs about √κ log(1/ε) steps instead of κ log(1/ε). With κ = 10,000 that is a 100× speedup.

Nesterov accelerated gradient

Nesterov looks ahead first. It evaluates the gradient where momentum is about to carry you:

yk     = θk + βk (θk - θk-1)
θk+1 = yk - (1/L) ∇f(yk)

• Convex, L-smooth. With β_k = (k − 1)/(k + 2), f(θ_k) − f^* = O(L‖θ₀ − θ^*‖²/k²). GD only gets O(1/k).
• Strongly convex. With β = (√κ − 1)/(√κ + 1), steps scale as √κ log(1/ε).
• Optimal. Nesterov proved that no first-order method can beat these rates on the worst case. Acceleration is as good as it gets with gradients alone.
• Why the look-ahead helps. If momentum is about to overshoot, the gradient at y already points back. Nesterov brakes one step earlier than heavy ball.

Seeing the √κ speedup

import numpy as np

def run(kappa, method, tol=1e-6, max_it=100000):
    # f(x) = 0.5 x^T A x with eigenvalues 1 and kappa
    A = np.diag([1.0, kappa]); L, mu = kappa, 1.0
    x = np.array([1.0, 1.0]); v = np.zeros(2); y = x.copy()
    for t in range(max_it):
        if np.linalg.norm(x) < tol:
            return t
        if method == "gd":
            x = x - (1 / L) * (A @ x)
        elif method == "heavy":
            a = 4 / (np.sqrt(L) + np.sqrt(mu)) ** 2
            b = ((np.sqrt(kappa) - 1) / (np.sqrt(kappa) + 1)) ** 2
            x_new = x - a * (A @ x) + b * v
            v = x_new - x; x = x_new
        elif method == "nesterov":
            b = (np.sqrt(kappa) - 1) / (np.sqrt(kappa) + 1)
            x_new = y - (1 / L) * (A @ y)
            y = x_new + b * (x_new - x); x = x_new
    return max_it

for k in [10, 100, 1000]:
    print(k, {m: run(k, m) for m in ["gd", "heavy", "nesterov"]})
# 10   {'gd': 132,   'heavy': 27,  'nesterov': 44}
# 100  {'gd': 1375,  'heavy': 95,  'nesterov': 158}
# 1000 {'gd': 13809, 'heavy': 321, 'nesterov': 519}

GD steps grow 10× each time κ grows 10×. Momentum steps grow about 3.2×, which is √10. On a pure quadratic, tuned heavy ball edges out Nesterov.

Momentum in deep learning

Deep learning uses a fixed β near 0.9, not the tuned values above. The steady-state step is η/(1 − β), or 10η at β = 0.9. So changing β changes the effective learning rate. With noisy gradients, momentum also acts as a low-pass filter. Its main gain is averaging out minibatch noise along consistent directions.

Why it matters in practice

• Tuning. If you raise β from 0.9 to 0.99, cut η by about 10×. Otherwise the effective step jumps.
• Ill-conditioned convex solves. Accelerated methods like FISTA solve large lasso and logistic problems in a fraction of the steps.

Interview check

• Why does momentum help on an ill-conditioned quadratic? It cuts the step count from order κ to order √κ. Along the flat axis velocity builds up. Across the steep axis it cancels.
• What is the difference between heavy ball and Nesterov? Nesterov takes the gradient at the look-ahead point. That gives a provable accelerated rate on all smooth convex functions. Heavy ball only has it on quadratics.
• What is the effective learning rate with momentum β? η/(1 − β) once velocity reaches steady state.

6. Adaptive methods and their caveats

AdaGrad

Gt = Gt-1 + gt ⊙ gt
θt+1 = θt - η gt / (√Gt + ε)

• AdaGrad comes from online convex optimization. Its regret bound is O(√T). It adapts to the geometry of the data with no tuning of per-feature rates.
• It shines with sparse features. A rare word or ad ID gets a big step on the few times it appears.
• Weakness. G only grows. So the step shrinks forever, roughly as 1/√t. In non-convex deep learning it stops learning too early.

RMSProp

RMSProp swaps the sum for an exponential moving average. It forgets old gradients:

vt = β2 vt-1 + (1 - β2) gt2
θt+1 = θt - η gt / (√vt + ε)

Adam

Adam is RMSProp plus momentum plus bias correction for the zero start:

mt = β1 mt-1 + (1 - β1) gt          m̂t = mt / (1 - β1t)
vt = β2 vt-1 + (1 - β2) gt2         v̂t = vt / (1 - β2t)
θt+1 = θt - η m̂t / (√v̂t + ε)

The update rules and AdamW are covered in Deep Learning, question 3. Here we focus on what the theory does and does not promise.

What Adam really does

• Scale invariance. Multiply every gradient by a constant c. Then m̂ and √v̂ both scale by c, so the step is unchanged. Adam is blind to loss scale. That is why η transfers well across models.
• Bounded steps. Each coordinate moves by at most about η per step when β₁² < β₂. So η is a trust region in parameter units.
• Sign descent in disguise. When the gradient signal dominates its noise, m̂/√v̂ ≈ sign(g). Adam then behaves like signSGD. When noise dominates, the ratio is small and steps shrink. This is a signal-to-noise filter.
• Diagonal preconditioner. √v̂ estimates the gradient scale per coordinate, not the curvature. It matches Newton only when the Hessian is diagonal and tied to gradient size. Neither holds in general.

The theory caveats

• Adam can diverge on convex problems. Reddi, Kale and Kumar built a simple online convex problem where Adam converges to the worst point. A rare large gradient gets washed out of v too fast. Their fix, AMSGrad, uses the running max of v̂. In practice AMSGrad rarely helps.
• Later results soften this. With β₂ close enough to 1 for the problem, Adam does converge to a stationary point. The failure needs β₂ too small. This is one reason LLM runs use β₂ = 0.95 to 0.999 with care.
• Generalization gap. Wilson et al. showed adaptive methods can find worse solutions than SGD on some tasks. On over-parameterized least squares, SGD finds the minimum-norm solution. Adam does not. For CNNs on vision, tuned SGD often wins by a little.
• ε is not cosmetic. When v is tiny, ε sets the step. Large ε turns Adam into momentum SGD. Too small can blow up low-precision runs.
• Memory. Two extra states per parameter. Adafactor factors v into row and column stats to save memory on large matrices.

Why it matters in practice

• Sparse recommender embeddings. Row-wise AdaGrad is still common for huge embedding tables. Rare IDs get large steps, and the state is one scalar per row.
• Transformers. AdamW wins clearly. Attention and embedding layers have very different gradient scales, and per-coordinate scaling handles that.
• Loss spikes in LLM training. When a rare large gradient hits after a quiet stretch, v is small, so the step is huge. Lower β₂, clip gradients, or skip the batch.

Interview check

• Why does AdaGrad stall in deep learning? Its denominator is a running sum, so every step is smaller than the last. Non-convex training needs to keep moving long after the early large gradients.
• Is Adam guaranteed to converge on convex problems? No. Reddi et al. gave a counterexample with fixed β₂. AMSGrad or a large enough β₂ restores convergence.
• Why is Adam's learning rate easier to set than SGD's? The step is scale invariant and roughly bounded by η per coordinate. So η means "how far each weight moves", independent of loss scale.

7. Constrained optimization and duality

The problem and the Lagrangian

minimize    f(x)
subject to  gi(x) ≤ 0,   i = 1..m
            hj(x) = 0,   j = 1..p

L(x, λ, ν) = f(x) + Σi λi gi(x) + Σj νj hj(x)        λi ≥ 0

Equality case first. At a constrained optimum, you cannot lower f by moving along the constraint surface. So ∇f must be normal to the surface. That gives ∇f = −ν∇h. The multiplier ν is the exchange rate between the two gradients.

KKT conditions

For a well-behaved problem, x^* is optimal only if there exist λ^*, ν^* with:

1. Stationarity. ∇f(x^*) + Σλ_i^*∇g_i(x^*) + Σν_j^*∇h_j(x^*) = 0.
2. Primal feasibility. g_i(x^*) ≤ 0 and h_j(x^*) = 0.
3. Dual feasibility. λ_i^* ≥ 0. An inequality can only push you back, never pull.
4. Complementary slackness. λ_i^*g_i(x^*) = 0. Either a constraint is tight, or its price is zero.

For convex problems that satisfy Slater's condition, KKT is also sufficient. Slater asks for one strictly feasible point.

Duality

The dual function is the best Lagrangian value for fixed prices:

q(λ, ν) = infx L(x, λ, ν)

• Always concave. It is a pointwise min of functions linear in (λ, ν). So the dual problem max q is convex, even when the primal is not.
• Weak duality. q(λ, ν) ≤ p^* for any λ ≥ 0. Every dual point gives a lower bound. The gap p^* − d^* is the duality gap.
• Strong duality. For convex problems with Slater, d^* = p^*. You can solve whichever side is easier.
• Shadow prices. If you relax g_i ≤ 0 to g_i ≤ u_i, then ∂p^*/∂u_i = −λ_i^*. The multiplier is the value of one more unit of slack.

Worked example: the SVM dual

The soft-margin SVM primal is:

minw,b,ξ  ½ ‖w‖2 + C Σi ξi
s.t.      yi(wTxi + b) ≥ 1 - ξi,    ξi ≥ 0

Give multiplier α_i to the margin rule and r_i to ξ_i ≥ 0. Set the Lagrangian's gradients to zero:

∂/∂w :  w = Σi αi yi xi
∂/∂b :  Σi αi yi = 0
∂/∂ξi:  C - αi - ri = 0   ⇒   0 ≤ αi ≤ C

Substitute back. The dual is a quadratic program in α:

maxα  Σi αi - ½ ΣiΣj αi αj yi yj xiTxj
s.t.   0 ≤ αi ≤ C,    Σi αi yi = 0

• Kernel trick. Data enters only through x_i^Tx_j. Replace it with k(x_i, x_j) and you get a nonlinear SVM for free.
• Support vectors from slackness. α_i = 0 for points safely outside the margin. 0 < α_i < C for points exactly on it. α_i = C for points inside the margin or misclassified.
• Sparsity. Only support vectors enter w. Prediction cost scales with their count, not with n.
• Solvers. SMO updates two α's at a time to keep Σαy = 0. It is coordinate ascent on the dual.

Why it matters in practice

• Budgeted allocation. Ads pacing and notification volume caps are constrained problems. The multiplier on a budget becomes a bid shading factor. You can update it online with dual ascent.
• Fairness and guardrail constraints. "Maximize engagement with no more than 1 percent drop in revenue" is a Lagrangian. Sweeping λ traces the trade-off curve.
• RLHF. Maximizing reward with a KL budget to the reference model gives the KL-penalized objective. The KL coefficient is a Lagrange multiplier.

Interview check

• What does complementary slackness say about an SVM? A point with α_i > 0 must sit on or inside the margin. Points outside the margin have zero weight and do not affect the model.
• Why is the dual always convex? q is an infimum of affine functions of the multipliers, so it is concave. Maximizing a concave function over a convex set is a convex problem.
• What does a large Lagrange multiplier on a budget mean? The budget is tight and costly. One more unit of budget would lower the loss by about λ.

8. Non-convex landscapes

Critical points and the Hessian

At a point with ∇f = 0, the Hessian eigenvalues tell you what you found:

• All positive. A local minimum.
• All negative. A local maximum. Rare and unstable.
• Mixed signs. A saddle. The index is the fraction of negative eigenvalues.
• Some zero. Degenerate. Flat directions are common in nets because of symmetries and over-parameterization.

Saddles dominate in high dimensions

Picture a random critical point in d dimensions. If each eigenvalue sign were a coin flip, the chance all d are positive is 2^−d. So most critical points are saddles. Dauphin et al. found that in real nets, the saddle index grows with loss. Critical points with high loss are almost all saddles. Local minima mostly sit near the global loss.

• Escaping is generic. Lee et al. showed GD from random init converges to a strict saddle with probability zero. A strict saddle has at least one clearly negative eigenvalue.
• But it can be slow. Plain GD may linger near a saddle for exponential time in bad cases. Perturbed GD adds small noise and escapes in polylog time. SGD noise does this for free.
• Plateaus look like convergence. A loss curve that flattens and then drops again was often near a saddle.

Flat vs sharp minima

Sharpness is often measured by λ_max(∇²f) or by the worst loss in a small ball. The intuition is in the figure. Test data shifts the loss surface a little. A flat minimum barely notices. A sharp one does. MDL and PAC-Bayes give the same story in theory terms. A flat minimum needs fewer bits to describe.

• Batch size link. Keskar et al. found large-batch training lands in sharper minima and generalizes worse. Less SGD noise means less push toward wide basins.
• The reparameterization catch. Dinh et al. showed you can rescale a ReLU net's layers to make any minimum look sharp. The function stays the same. So raw sharpness is not a complete story. Scale-aware measures fix part of this.
• SAM. Sharpness-aware minimization minimizes the worst loss in a radius ρ. It takes a gradient ascent step to θ + ρg/‖g‖, then descends using the gradient there. It costs two passes per step and often improves accuracy.
• SWA. Averaging weights along the end of a cyclic schedule lands in the center of a wide basin.

What the neural net loss surface looks like

• Symmetry. Permuting hidden units gives identical functions. So every minimum has a huge number of copies.
• Over-parameterization helps. Wide nets have many global minima that form connected sets. Gradient flow behaves close to a linear model near init (the neural tangent kernel regime).
• Mode connectivity. Two independently trained solutions are joined by a simple curved path of low loss. After aligning permutations, often a straight line works too.
• Hessian spectrum. A bulk of near-zero eigenvalues plus a few large outliers. The outliers number about the count of classes. Most directions are flat.
• Edge of stability. With full-batch GD, sharpness rises until λ_max ≈ 2/η. Then the loss keeps falling non-monotonically. Classic smooth theory does not predict this.

Why it matters in practice

• Restart vs wait. A long plateau is often a saddle, not a bad minimum. A higher LR or more noise usually moves you off it.
• Large-batch training. When you scale the batch, watch the generalization gap. LR warmup, longer schedules or SAM can recover it.
• Model soups and merging. Averaging fine-tuned checkpoints works because they share a basin. Merging models from different pre-trained inits usually fails.

Interview check

• Are local minima the main problem in deep nets? Mostly no. High-loss critical points are mostly saddles. Most local minima have loss close to the global one.
• Why might large batches generalize worse? Less gradient noise means less escape from sharp basins. Sharp minima are more sensitive to the train to test shift.
• Is sharpness a reliable generalization measure? Not in raw form. Rescaling layers changes it without changing the function. Normalized measures and controlled experiments are needed.

9. Learning rate schedules and warmup

Why decay at all

Recall the SGD noise ball. Its radius scales with η. Early on you want a large η to cover distance and explore. Late you want a small η to shrink the ball. A schedule trades between the two. It is the deep learning version of Robbins–Monro.

The common shapes

• Step decay. Multiply η by 0.1 at fixed epochs. Loss drops sharply right after each cut. This is the classic ResNet recipe.
• Cosine. η_t = η_min + ½(η_max − η_min)(1 + cos(πt/T)). Smooth, with one less knob than step decay. It spends a long time at high rates and decays gently at the end.
• Linear decay to zero. Simple and strong for fine-tuning and many LLM runs.
• Warmup-stable-decay (WSD). Warm up, hold a constant rate, then decay fast in the last 10 to 20 percent. You can branch a decay from any checkpoint. So one long run yields models for many budgets.
• One-cycle. Leslie Smith's recipe. Ramp η up to a high peak over about 30 to 45 percent of training. Then ramp it down to far below the start. Momentum moves the opposite way. It can train in far fewer epochs ("super-convergence").

Why warmup helps, especially for Adam

Warmup ramps η from near zero over the first few hundred to few thousand steps. Several reasons stack up:

• Noisy second moments. In the first steps, v̂ is an average of only a few squared gradients. Its variance is huge. For coordinates where it comes out small, the step is large. RAdam (Liu et al.) showed the variance of the adaptive rate is unbounded at the start. Warmup hides this until v̂ has enough samples.
• Bias correction makes step one big. At t = 1, m̂/√v̂ = g/|g| = ±1. Every weight moves by a full η on a single noisy gradient.
• High curvature at init. Early sharpness is often well above 2/η_max. Warmup lets training move to flatter regions before the rate gets large. Measured sharpness drops during warmup.
• Post-LN transformers. Xiong et al. showed gradients near the output layer are large at init with post-norm. Without warmup, they diverge. Pre-norm reduces but does not remove the need.
• Large batches. The linear scaling rule gives a large peak η. Early on, the loss is far from quadratic and the rule does not hold yet.

Why it matters in practice

• Finding the peak rate. Run an LR range test. Raise η exponentially over a few hundred steps and plot loss. Pick a value a bit below where loss starts to climb.
• Warmup length. One to five percent of steps is typical. A spike at the end of warmup means the peak is too high.
• Comparing runs. A cosine run's loss is only meaningful at its end. Comparing two cosine runs mid-way favors the shorter one. WSD makes intermediate checkpoints comparable.

Interview check

• Why does Adam need warmup when SGD often does not? Adam divides by √v̂, which is estimated from very few samples early. That estimate can be tiny, which makes steps huge. SGD has no such ratio.
• What does a cosine schedule give over a constant rate? The late decay shrinks the noise ball, so the final loss is lower. Most of the gain shows up in the last 10 to 20 percent.
• What is super-convergence? One-cycle with a very high peak trains some nets in a fraction of the usual epochs. The large rate acts as a regularizer.

10. Coordinate descent and proximal methods

Coordinate descent

Minimize over one coordinate at a time while holding the rest fixed. Cycle through coordinates.

• Each one-dimensional problem is often solvable in closed form. That makes each step very cheap.
• It converges for smooth convex f. It also converges for f(x) = g(x) + Σ_jh_j(x_j) with g smooth and h separable (Tseng). Lasso is exactly this form.
• It can stall if the non-smooth part couples coordinates. An example is the max of two linear functions.
• Its rate depends on per-coordinate smoothness L_j, not the global L. That can be much better when features have different scales.

Lasso and soft-thresholding

minw  (1/2n) ‖y - Xw‖2 + λ ‖w‖1

Fix all weights except w_j. Let r_−j be the residual without feature j. Define ρ_j = x_j^Tr_−j/n and z_j = ‖x_j‖²/n. The one-variable problem is a parabola plus λ|w_j|. Its subgradient condition gives:

wj = S(ρj, λ) / zj        S(z, t) = sign(z) · max(|z| - t, 0)

S is the soft-thresholding operator. If the correlation |ρ_j| is below λ, the weight is exactly zero. Otherwise it shrinks toward zero by λ. That is how lasso produces sparse models.

import numpy as np

def soft(z, t):
    return np.sign(z) * np.maximum(np.abs(z) - t, 0.0)

def lasso_cd(X, y, lam, n_iter=200):
    n, d = X.shape
    w = np.zeros(d)
    r = y - X @ w                     # residual
    col_sq = (X ** 2).sum(axis=0) / n
    for _ in range(n_iter):
        for j in range(d):
            r += X[:, j] * w[j]       # remove feature j
            rho = X[:, j] @ r / n
            w[j] = soft(rho, lam) / col_sq[j]
            r -= X[:, j] * w[j]       # add it back
    return w

rng = np.random.default_rng(0)
X = rng.normal(size=(200, 10))
w_true = np.array([3, -2, 0, 0, 1.5, 0, 0, 0, 0, 0.0])
y = X @ w_true + 0.5 * rng.normal(size=200)
print(np.round(lasso_cd(X, y, lam=0.1), 2) + 0.0)
# [ 2.93 -1.88  0.    0.    1.37 -0.01  0.    0.    0.    0.  ]

Lasso keeps the three true features and zeros almost all the rest. The kept weights shrink by about λ. That bias is why people often refit the kept features without the penalty.

The proximal operator

proxηh(v) = argminx  h(x) + (1/2η) ‖x - v‖2

It finds a point that keeps h small but stays close to v. Some well-known cases:

• h = λ‖x‖₁. The prox is soft-thresholding S(v, ηλ), per coordinate.
• h = indicator of a convex set C. The prox is projection onto C. So projected GD is a special case.
• h = (λ/2)‖x‖². The prox is v/(1 + ηλ). That is weight decay.
• h = nuclear norm. The prox soft-thresholds the singular values. It is used for low-rank matrix completion.

Proximal gradient (ISTA) and FISTA

For f = g + h with g L-smooth and h simple:

θk+1 = proxηh( θk - η ∇g(θk) )        η = 1/L

ISTA keeps the O(1/k) rate of GD, even though f is not smooth. Add Nesterov momentum and you get FISTA, with rate O(1/k²). The subgradient method, by contrast, only gets O(1/√k).

Why it matters in practice

• Feature selection at scale. glmnet uses coordinate descent with warm starts along a decreasing λ path. It fits a whole regularization path in about the time of one fit.
• Sparse CTR models. FTRL-Proximal adds L1 to online logistic regression. It keeps billion-feature ad models sparse enough to serve.
• Weight decay is a prox step. AdamW's decoupled decay is a proximal step for an L2 penalty, applied outside the adaptive scaling.

Interview check

• Why does L1 give exact zeros but L2 does not? The L1 subgradient at zero is the whole interval [−λ, λ]. Any small correlation fits in it, so zero is optimal. L2 has gradient zero at zero, so it only shrinks.
• What is the prox of the L1 norm? Soft-thresholding: sign(v)·max(|v| − ηλ, 0), per coordinate.
• Why not just run SGD on the lasso loss? Subgradient steps make weights hover near zero but rarely hit it exactly. You lose the sparsity. Prox steps set them to zero exactly.

11. EM as optimization

The lower bound

Observed data x, hidden z, parameters θ. For any distribution q(z):

log p(x; θ) = Eq[ log p(x, z; θ) ] - Eq[ log q(z) ]  +  KL( q(z) ‖ p(z | x; θ) )
              =            ELBO(q, θ)                      +  KL( q ‖ posterior )

KL is never negative. So the ELBO (evidence lower bound) is always below the log-likelihood. EM is coordinate ascent on the ELBO:

• E-step. Fix θ. Maximize over q. The best q is the posterior p(z | x; θ_old). Then KL = 0 and the bound touches the log-likelihood.
• M-step. Fix q. Maximize over θ. Only the first term depends on θ. So maximize Q(θ | θ_old) = E_q[log p(x, z; θ)], the expected complete-data log-likelihood.

Why the likelihood never drops. After the E-step, log p(x; θ_old) = ELBO(q, θ_old). The M-step gives ELBO(q, θ_new) ≥ ELBO(q, θ_old). The bound gives log p(x; θ_new) ≥ ELBO(q, θ_new). Chain the three and you are done.

Gaussian mixture example

Model p(x) = Σ_k π_k N(x; μ_k, σ_k²). The hidden z_i is the component for point i.

E-step:  rik = πk N(xi; μk, σk2) / Σj πj N(xi; μj, σj2)
M-step:  Nk = Σi rik,   πk = Nk/n
         μk = Σi rik xi / Nk,   σk2 = Σi rik (xi - μk)2 / Nk

The M-step is just weighted maximum likelihood. Each point counts toward component k with weight r_ik.

import numpy as np

def em_gmm_1d(x, k=2, n_iter=100, seed=0):
    rng = np.random.default_rng(seed)
    pi = np.full(k, 1 / k)
    mu = rng.choice(x, k, replace=False)
    var = np.full(k, x.var())
    prev = -np.inf
    for it in range(n_iter):
        # E-step: responsibilities r[i, j] = p(z = j | x_i), in log space
        logp = (-0.5 * np.log(2 * np.pi * var)
                - 0.5 * (x[:, None] - mu) ** 2 / var + np.log(pi))
        m = logp.max(axis=1, keepdims=True)
        log_lik = (m.ravel() + np.log(np.exp(logp - m).sum(axis=1))).sum()
        r = np.exp(logp - m); r /= r.sum(axis=1, keepdims=True)
        # M-step: weighted MLE
        nk = r.sum(axis=0)
        pi = nk / len(x)
        mu = (r * x[:, None]).sum(axis=0) / nk
        var = (r * (x[:, None] - mu) ** 2).sum(axis=0) / nk
        assert log_lik >= prev - 1e-9   # EM never lowers the likelihood
        if log_lik - prev < 1e-8:
            break
        prev = log_lik
    return pi, mu, var, it

rng = np.random.default_rng(1)
x = np.concatenate([rng.normal(-2, 0.5, 300), rng.normal(3, 1.0, 700)])
pi, mu, var, it = em_gmm_1d(x)
print(np.round(pi, 2), np.round(mu, 2), np.round(np.sqrt(var), 2), it)
# [0.3 0.7] [-2.05  2.96] [0.46 1.02] 21

EM recovers the true weights 0.3 and 0.7, means −2 and 3, and spreads 0.5 and 1.0. The assert checks the monotone guarantee on every step.

Properties and pitfalls

• Local optima. The likelihood is non-convex. Run several random starts, or start from k-means++.
• Degenerate solutions. If one component sits on a single point, σ² → 0 and the likelihood → ∞. Add a variance floor or a prior.
• Speed. Convergence is linear. The rate depends on the fraction of missing information. Heavily overlapping clusters converge slowly.
• k-means is hard EM. Shrink all variances to zero with equal weights. Responsibilities become 0 or 1, and you get Lloyd's algorithm.
• Generalized EM and variational EM. If the M-step is hard, any improving step keeps the guarantee. If the posterior is intractable, restrict q to a family. That is variational inference, and the VAE objective comes from the same ELBO.

Why it matters in practice

• Noisy labels. Treat the true label as hidden. EM estimates each annotator's confusion matrix (Dawid–Skene). It is a standard way to clean crowd labels.
• Position bias in ranking. Click models treat "examined" as hidden. EM separates position bias from true relevance in click logs.
• Missing data. EM gives maximum likelihood estimates with missing features, without ad hoc imputation.

Interview check

• Why does EM never decrease the likelihood? The E-step makes the ELBO tight at the current θ. The M-step raises the ELBO. The ELBO is a lower bound, so the likelihood rises at least as much.
• How is k-means related to EM? It is EM on a GMM with equal, shared, vanishing variances. Soft responsibilities become hard assignments.
• Can EM blow up on a GMM? Yes. A component can collapse onto one point, driving its variance to zero and the likelihood to infinity. A variance floor or prior prevents it.

12. Hyperparameter optimization

The problem

Minimize validation loss F(λ) over hyperparameters λ. F is a black box. It has no gradients, each call is costly, and it is noisy across seeds. Some dimensions are continuous, some are integers, and some are choices.

Grid search

Try every combination on a grid. With k values per dimension and d dimensions, that is k^d runs. It wastes budget badly. Suppose only one of three dimensions matters. A 5×5×5 grid runs 125 trials but tests only 5 distinct values of the one that matters.

Random search

Sample each hyperparameter from a distribution. Bergstra and Bengio showed it beats grid search for the same budget. Every trial tests a new value of every dimension. With 125 random trials, you test 125 values of the one that matters.

• Sample learning rate and weight decay log-uniformly. Their effect is multiplicative.
• Take 60 random trials. You then have a 95 percent chance to hit the top 5 percent of the space. That is because 1 − 0.95⁶⁰ ≈ 0.95.
• Random search parallelizes perfectly and is a strong baseline.

Bayesian optimization with a GP and expected improvement

BO fits a cheap surrogate model to the trials so far. It then picks the next trial by maximizing an acquisition function that balances exploration and exploitation.

The surrogate. A Gaussian process gives a posterior mean μ(λ) and standard deviation σ(λ) at every point. With kernel matrix K over observed points, k_* between a new point and observed points, and targets y:

μ(λ)  = k*T (K + σn2I)-1 y
σ2(λ) = k(λ, λ) - k*T (K + σn2I)-1 k*

A Matérn 5/2 kernel is the usual choice. It is smooth but not too smooth.

Expected improvement. Let f_best be the lowest loss so far. The improvement at λ is max(f_best − F(λ), 0). Under the GP it has a closed-form mean:

Z     = (fbest - μ(λ) - ξ) / σ(λ)
EI(λ) = (fbest - μ(λ) - ξ) Φ(Z) + σ(λ) φ(Z)

Here Φ and φ are the standard normal CDF and PDF. The first term rewards a low predicted mean, which is exploitation. The second rewards high uncertainty, which is exploration. The small ξ ≥ 0 tilts toward exploration.

import numpy as np
from math import erf, sqrt, pi

def rbf(a, b, ls=0.3):
    return np.exp(-0.5 * (a[:, None] - b[None, :]) ** 2 / ls ** 2)

def gp_posterior(X, y, Xs, noise=1e-6):
    K = rbf(X, X) + noise * np.eye(len(X))
    Ks = rbf(X, Xs)
    Lc = np.linalg.cholesky(K)
    alpha = np.linalg.solve(Lc.T, np.linalg.solve(Lc, y))
    mu = Ks.T @ alpha
    v = np.linalg.solve(Lc, Ks)
    sd = np.sqrt(np.maximum(1.0 - (v ** 2).sum(axis=0), 1e-12))
    return mu, sd

def expected_improvement(mu, sd, best, xi=0.01):
    z = (best - mu - xi) / sd                 # minimization
    cdf = 0.5 * (1 + np.vectorize(erf)(z / sqrt(2)))
    pdf = np.exp(-0.5 * z ** 2) / sqrt(2 * pi)
    return (best - mu - xi) * cdf + sd * pdf

f = lambda x: np.sin(3 * x) + 0.5 * x         # stand-in for validation loss
grid = np.linspace(0, 2, 401)
X = np.array([0.1, 1.0, 1.9]); y = f(X)
for step in range(6):
    mu, sd = gp_posterior(X, y - y.mean(), grid)
    ei = expected_improvement(mu + y.mean(), sd, y.min())
    x_next = grid[np.argmax(ei)]
    X, y = np.append(X, x_next), np.append(y, f(x_next))
print("best x = %.3f, best f = %.3f" % (X[y.argmin()], y.min()))
# best x = 1.525, best f = -0.228   (true minimum near x = 1.515)

Nine evaluations land within 0.01 of the true minimum. Random search would need many more.

• Limits. GP fitting costs O(n³) in the number of trials. It struggles past about 20 dimensions and with categorical choices.
• TPE. The tree-structured Parzen estimator models p(λ | good) and p(λ | bad) instead. It handles conditional and categorical spaces well. Optuna and Hyperopt use it.

Hyperband and successive halving

Most bad configs look bad early. So stop them early and give their budget to the promising ones.

Successive halving (SH). Start n configs with budget r each, like r epochs. Keep the top 1/η by validation loss. Multiply their budget by η. Repeat until one remains. With η = 3, 81 configs at 1 epoch become 27 at 3, 9 at 9, 3 at 27, and 1 at 81.

The catch. SH assumes early rankings predict final rankings. That fails for configs that start slow but finish strong, like a low learning rate. Choosing n versus r is a bet on how reliable early signals are.

Hyperband. It hedges that bet. It runs several SH brackets, from "many configs, tiny budget" to "few configs, full budget". The last bracket is plain random search. Hyperband is never much worse than the best bracket, and you do not have to guess.

• BOHB. Replaces random sampling in Hyperband with a TPE-style model. It gets early stopping plus learned proposals.
• ASHA. An asynchronous version of SH that scales to hundreds of workers with no idle waiting.
• Population-based training. Trains many models in parallel. It copies the weights of good ones into bad ones and perturbs their hyperparameters. This yields a schedule, not just a fixed value.

Why it matters in practice

• Tune what matters. Learning rate dominates for most nets. Tune it first, on a log scale. Then weight decay, warmup and batch size.
• Compute budgets. For a model that trains in hours, ASHA or Hyperband with 50 to 200 trials is a solid default. For one that trains in days, tune on a smaller proxy.
• Hyperparameter transfer. The maximal update parameterization (μP) makes the best learning rate stable across widths. You tune a small model and reuse the settings on the big one.
• Noise and overfitting. With many trials, the best validation score is biased upward. Confirm the winner on fresh seeds and a held-out test set.

Interview check

• Why does random search beat grid search? Usually only a few hyperparameters matter. Random search tests a new value of each one on every trial. Grid search repeats the same few values.
• What does expected improvement trade off? Low predicted loss against high uncertainty. The mean term exploits, and the σφ(Z) term explores.
• When does Hyperband fail? When early performance does not predict final performance. Low learning rates or long warmups look bad early but can win. Use a large minimum budget there.