Optimization Theory

We jumped from x=5 to x=1.8 — a big move because the gradient was steep (32). The true minimum is at x*=−1/3 ≈ −0.333, so we're getting closer but still far.

After 5 steps we're at x ≈ −0.279, close to the true minimum x* = −0.333. The linear recurrence x_k+1 = 0.4x_k − 0.2 converges because |0.4| < 1 (the contraction factor is 1 − αL = 1 − 0.6 = 0.4).

A condition number of 50 means the function is 50 times steeper in one direction than another. GD zigzags along the narrow valley, making ~59 steps for each order of magnitude of improvement. This is why preconditioning (rescaling coordinates) and Adam (per-parameter learning rates) help — they effectively reduce κ.

The contraction factor is |1 − αL| = |1 − 0.5 × 6| = |1 − 3| = 2 > 1. Since the contraction factor exceeds 1, each step amplifies the error. The iterates 1, −3, 5, −11, ... grow exponentially. The maximum stable learning rate is α < 2/L = 1/3 ≈ 0.333. This is why α = 0.5 diverges.

javascript
function gradDescent(gradf, x0, alpha, steps) {
  let x = x0;
  for (let i = 0; i < steps; i++) {
    x = x - alpha * gradf(x);
  }
  return x;
}

With the optimal rate, GD converges in exactly one step for a pure quadratic. This is because the Newton step for a quadratic is exact. In practice, loss functions aren't pure quadratics, so the optimal rate changes at every point — motivating adaptive methods like Adam.

The first step with momentum is identical to vanilla GD because v₀ = 0. The momentum effect kicks in starting from step 2, where v carries history from step 1.

Notice momentum overshoots — x went past the minimum at −0.333 to −2.36, then even further to −4.888. The accumulated velocity is too high. This is typical early behavior: momentum takes a few steps to stabilize. With vanilla GD, x₃ would be 0.008 — closer but slower to eventually converge. Momentum trades early overshoot for faster asymptotic convergence.

Standard momentum computes the gradient at where you are. Nesterov computes the gradient at where you're about to be. If momentum is carrying you past the minimum, the lookahead gradient points backward — so Nesterov applies the brakes before overshooting, not after. This gives O(1/k²) convergence on convex problems vs O(1/k) for standard momentum. That's provably optimal for first-order methods (Nesterov, 1983).

With β = 0.99, the velocity remembers roughly 100 past gradients. This is great for smoothing out minibatch noise (SGD on large datasets) but terrible for rapidly changing landscapes. Most practitioners use β = 0.9 (window ~10) as a balance. Adam uses β₁ = 0.9 for momentum and β₂ = 0.999 (window ~1000) for the second moment — the second moment changes more slowly than the mean.

Correct order: Evaluate f → Compute gradient → Update velocity → Update parameters. You must evaluate the function before computing its gradient (backprop requires the forward pass). Then velocity accumulates the gradient with history, and finally the parameter update uses the velocity.

Line 6 is the bug. It uses x + alpha * v instead of x - alpha * v. Gradient descent moves opposite to the gradient (downhill). Adding the gradient moves uphill — maximizing the function instead of minimizing it. The sign error is the most common optimization bug in practice.

The closest point on the line x+y=1 to the origin is (1/2, 1/2). By symmetry, x=y was expected. The Lagrange multiplier λ = −1 means: if we changed the constraint to x+y=1.01, the optimal value would decrease by about 0.01 (from 0.5 to 0.4999...).

From the derivation above: λ* = −1. This tells us that relaxing the constraint from x+y=1 to x+y=1+ε would decrease the objective by approximately ε. Tightening it to x+y=1−ε would increase the objective by ε.

The maximum entropy distribution is the uniform distribution — maximum uncertainty. This result generalizes: for n outcomes, max entropy is ln(n) achieved by the uniform distribution. This is the principle behind cross-entropy loss — it's the gap between your model's distribution and the data distribution.

javascript
function solveConstrained(a, b, c) {
  // -lambda/(2a) - lambda/(2b) = c
  // -lambda * (1/(2a) + 1/(2b)) = c
  // -lambda * (a+b)/(2ab) = c
  const lam = -c * 2 * a * b / (a + b);
  const x = -lam / (2 * a);
  const y = -lam / (2 * b);
  return [x, y];
}

With constraint x+y=1.1, the solution is x*=y*=0.55 and f*=2(0.55)² = 0.605. The linear approximation using λ gives f* ≈ 0.5 + |λ|×0.1 = 0.6, which is close. The extra 0.005 is the second-order correction. The sensitivity interpretation: λ gives the first-order rate of change of f* with respect to the constraint's right-hand side.

The unconstrained minimizer of x² is x=0, but 0 < 2 violates g(x) ≤ 0. So the constraint must be active.

All KKT conditions satisfied. The constraint pushes x away from the unconstrained minimum (x=0) to x=2. The multiplier μ=4 measures this push. In ML: this is like weight clipping pushing your weights to the constraint boundary.

μ* = 4 > 0, confirming the constraint is active. If we changed the constraint to x ≥ 1.9, the optimal value would improve (decrease) by approximately μ × 0.1 = 0.4. Actual change: 2² − 1.9² = 4 − 3.61 = 0.39 ≈ 0.4.

The unconstrained minimum x=5 satisfies x ≥ 2 (since 5 > 2), so the constraint is inactive (slack). By complementary slackness, μ=0. The constraint doesn't affect the solution at all — removing it wouldn't change x*. This is the "you don't feel the wall" case.

For convex problems (convex objective, convex feasible set), KKT conditions are both necessary AND sufficient. This is because convex problems have no spurious local minima — any point satisfying the first-order conditions must be globally optimal. This is the foundation of convex optimization: if you can prove your problem is convex and find a KKT point, you're done. SVMs, linear regression with L2 regularization, and logistic regression are all convex.

Both constraints are inactive. The unconstrained minimum lies safely inside the feasible region [−3, 1]. Neither constraint affects the solution. This is common in practice — many constraints in an optimization problem may be inactive at the solution.

Strong duality holds — the gap is zero. This always happens for LPs (linear programs). The dual variable μ* = 1.5 is the shadow price: increasing the right-hand side of 2x ≥ 4 to 2x ≥ 5 would increase the optimal cost by μ × 1 = 1.5 (from 6 to 7.5).

Duality gap = p* − d* = 5 − 3 = 2 > 0. The dual provides a certified lower bound on the primal optimum. We know the true optimal value is between 3 and 5. For non-convex problems, the gap can be large. For convex problems satisfying Slater's condition (a strictly feasible point exists), the gap is zero. This is why convexity is so valuable — it means the dual bound is tight.

Option 1: x²+y² is convex, x+y ≤ 5 is a linear (hence convex) constraint, and Slater's condition holds (e.g., x=0, y=0 is strictly feasible since 0 < 5). Strong duality guaranteed. Options 2 and 3 have non-convex objectives (sin is non-convex, x³ is non-convex over all R though it IS convex on [1,∞)) — no guarantee of strong duality.

Strong duality holds because the problem is convex. The dual computation independently confirms p* = 4 and μ* = 4, matching our KKT solution from Chapter 4.

The value 3 gets shrunk by λ=1 toward zero, becoming 2. It doesn't become zero because |3| > λ. Only values with |x| ≤ λ get zeroed out.

Since |0.5| < λ = 1, the value gets thresholded to exactly zero. This is the sparsity-inducing property of L1. Any weight smaller than the threshold gets killed. This is exactly what happens in LASSO regression.

The sign is preserved, and the magnitude is shrunk by λ = 0.5. The soft thresholding operator is symmetric: it shrinks positive and negative values equally toward zero.

javascript
function softThreshold(x, lambda) {
  return x.map(xi => {
    const s = Math.sign(xi);
    const m = Math.max(Math.abs(xi) - lambda, 0);
    return s * m;
  });
}

prox(6) = 6/(1 + 1×2) = 6/3 = 2. L2 regularization shrinks values toward zero by a multiplicative factor but never reaches zero. L1 regularization thresholds values to zero if they're small enough. This is the fundamental difference: L1 produces sparsity, L2 produces small weights. In neural networks, weight decay (L2) keeps all weights small; L1 pruning kills unimportant weights entirely.

Line 5 is the bug. It computes |x[i]| − lam without the max(..., 0) clamp. When |x[i]| < lam, this produces a negative value, and multiplied by sign gives a flipped-sign result instead of zero. For example, x=0.3, lam=1: should give 0, but gives sign(0.3)×(0.3−1) = −0.7. The fix: const m = Math.max(Math.abs(x[i]) - lam, 0);

One million iterations — painfully slow. This is the O(1/k) sublinear rate. Doubling accuracy requires doubling iterations. This motivates using Nesterov acceleration or exploiting strong convexity for faster convergence.

1,000 vs 1,000,000 — that's a 1000x speedup just from the Nesterov trick. The O(1/k²) rate is provably optimal for first-order methods on convex functions (Nesterov, 1983). You cannot do better with only gradient information — this is the information-theoretic lower bound.

Strong convexity gives linear convergence — each step removes a constant fraction (1/κ) of the remaining error. This is exponentially faster than the O(1/k) rate. The cost is proportional to κ: ill-conditioned problems (large κ) converge slowly. Most deep learning loss surfaces are extremely ill-conditioned (κ > 10⁶), which is why raw GD is impractical.

GD takes ~1382 steps, Nesterov takes ~138 steps (10x fewer), Newton takes ~5 steps (276x fewer). But Newton's method requires computing and inverting the Hessian at each step — O(n³) work per step where n is the problem dimension. For a neural network with n = 10⁸ parameters, one Newton step would require 10²⁴ operations. That's why deep learning uses first-order methods despite their slower convergence rate — the per-step cost dominates.

GD wins! 100 vs 500 seconds. Newton takes fewer steps but each step is enormously expensive. At 10M parameters, Newton is already 5x slower despite needing 200x fewer iterations. At 1B parameters, Newton would take 5 × 10⁶ seconds per step — completely infeasible. This is the fundamental reason deep learning uses first-order methods.

One step to the exact minimum. Since f is quadratic, Newton solves it perfectly. Compare to GD with optimal α = 1/L = 1/20: after one step, x₁ = 10 − 20/20 = 9 (y would reach 0, but x barely moved). GD would need many iterations because of the 10:1 condition number. Newton doesn't care about conditioning — it normalizes each direction by its curvature.

After one GD step: (9, 0). The y-direction converges fast (curvature 20, step size matches) but x-direction crawls (curvature 2, but step size is tuned for curvature 20). This is the zigzag problem caused by κ = 20/2 = 10. Newton would need 1 step; GD needs ~10 × ln(10) ≈ 23 steps for the same accuracy.

200 exabytes — that's ~200 million terabytes. The world's total data storage is ~120 exabytes. You literally cannot store the Hessian of a 7B model. This is why exact Newton's method is impossible for deep learning, and why approximations (Adam's per-parameter second moments, K-FAC's block-diagonal approximation, L-BFGS's low-rank updates) exist.

All options are real concerns, but the computational cost is the dealbreaker. Non-convexity and saddle points can be handled with modifications (regularized Newton, trust regions). Stochastic Hessians exist. But there's no getting around O(n²) storage and O(n³) per-step cost. Even approximate methods like L-BFGS (which stores a low-rank Hessian approximation) struggle beyond ~1M parameters because the inner products still scale with n. Adam is the practical compromise: it maintains per-parameter second moment estimates (diagonal Hessian approximation) at O(n) cost.

javascript
function newtonStep(df, ddf, x) {
  const h = ddf(x);
  if (Math.abs(h) < 1e-12) return x; // avoid div by zero
  return x - df(x) / h;
}

Correct order: Compute gradient → Compute Hessian → Solve linear system → Update x → Check convergence. The key insight: you solve HΔx = −∇f instead of computing H⁻¹ explicitly. Solving a linear system is O(n³) but inverting a matrix is also O(n³) with a worse constant. In practice, Cholesky decomposition (for PD Hessians) or conjugate gradient (for approximate solves) is used.

The eigenvalue −2 means the surface curves downward in one direction — you can decrease f by moving in that direction. This makes it a saddle point, not a minimum. In high-dimensional spaces, saddle points vastly outnumber local minima. For a random critical point in n dimensions, the probability that ALL eigenvalues are positive (true minimum) is roughly (1/2)ⁿ. For n = 10⁶ parameters, saddle points are overwhelmingly more common than minima.

Minimum B is 200x flatter and can tolerate 200x larger learning rates. Flat minima generalize better because small perturbations to the weights (which happen between train and test distributions) cause small changes in loss. In sharp minima, tiny weight perturbations can dramatically increase the loss — the model has memorized the training data at a fragile parameter configuration.

At random initialization, the gradient magnitudes are large and unpredictable. The effective Lipschitz constant L is very high early in training. A large learning rate violates α < 2/L and causes the loss to spike or diverge. Warmup keeps α small until the gradients stabilize, then ramps up to the target rate. Empirically, transformers are especially sensitive to this — without warmup, training often diverges in the first few hundred steps.

Correct order: Random init → LR warmup → Steady at peak LR → Cosine decay → Final low LR. Initialize randomly, warm up to avoid early divergence, train at peak rate for most of the budget (this is where most learning happens), then decay the LR to settle into a flat minimum. The final LR is typically 10% of peak. This is the standard "warmup + cosine decay" schedule used by GPT-3, LLaMA, and most modern LLMs.

AdamW is the standard choice for transformer training. Newton is infeasible (O(n²) memory). L-BFGS requires full-batch gradients and scales poorly beyond ~1M parameters. SGD works but requires extensive learning rate tuning and converges slower on transformers. AdamW provides per-parameter learning rate adaptation (diagonal Hessian approximation) at O(n) memory cost, with decoupled weight decay for regularization. The cost: 3x the memory of the model parameters (weights + m + v buffers).

84 GB — more than one A100 (80 GB). The optimizer states alone (56 GB) are 4x the model weights (14 GB). This is why model parallelism and ZeRO-style optimizer state sharding are essential for training large models. The optimizer is the memory bottleneck, not the model.

Line 8 has the sign wrong for weight decay. In AdamW (decoupled weight decay), the weight decay term should be added to the update, not subtracted: params[i] -= lr * (mHat / (Math.sqrt(vHat) + 1e-8)) + lr * wd * params[i]. The current code subtracts wd*params[i], which means it's growing the weights instead of shrinking them. Equivalently, write it as: params[i] = params[i] * (1 - lr * wd) - lr * mHat / (Math.sqrt(vHat) + 1e-8);