Optimization Theory

A function f is convex if for all x, y in its domain and t in [0,1]: f(tx + (1-t)y) ≤ tf(x) + (1-t)f(y). Geometrically, the chord between any two points lies above (or on) the function. Strict convexity holds with strict inequality for x ≠ y. For convex functions, any local minimum is also a global minimum: if x* is a local minimum and there existed a better point y, the convexity inequality along the segment from x* to y would contradict x* being a local minimum. Strictly convex functions have a unique global minimum. This guarantee makes optimization tractable—gradient descent converges to the global optimum for convex loss functions.

First-order necessary condition (FONC): at a local minimum x*, the gradient ∇f(x*) = 0 (stationary point). Second-order necessary condition (SONC): the Hessian H(x*) is positive semi-definite (all eigenvalues ≥ 0). Second-order sufficient condition (SOSC): if ∇f(x*) = 0 and H(x*) is positive definite (all eigenvalues > 0), then x* is a strict local minimum. A point satisfying FONC but with indefinite Hessian is a saddle point. For convex functions, FONC is both necessary and sufficient for a global minimum. In deep learning, most critical points in practice are saddle points rather than local minima, especially in high-dimensional parameter spaces.

Gradient descent minimizes f(x) by iterating x_{t+1} = x_t - η ∇f(x_t), where η > 0 is the learning rate (step size). Derivation: a first-order Taylor expansion gives f(x-ε∇f) ≈ f(x) - ε||∇f||² < f(x) for small ε, so moving against the gradient always decreases f locally. The learning rate controls step size: too large causes divergence (can overshoot the minimum); too small causes slow convergence. For L-smooth functions (||∇f(x)-∇f(y)||≤L||x-y||), the optimal constant learning rate is 1/L, giving convergence rate O(1/T) for convex f and O(1/T²) with Nesterov acceleration.

import numpy as np
def gradient_descent(f, grad_f, x0, lr=0.01, iters=1000):
    x = x0.copy()
    for _ in range(iters):
        x -= lr * grad_f(x)
    return x

SGD replaces the full gradient (computed over all n samples) with a noisy estimate from a mini-batch of size B: ∇f_B(x) = (1/B)Σ∇f_i(x). Full-batch (B=n): exact gradient, O(n) per step, slow convergence rate. SGD (B=1): cheapest step, highest noise, does not converge to exact minimum (oscillates). Mini-batch (1<B<n) balances variance and compute. Large B reduces noise but loses data augmentation's regularization benefit and requires larger learning rate. The linear scaling rule: scale learning rate by B/B_ref when increasing batch size. Gradient variance: O(1/B). Hardware efficiency peaks at B using all GPU memory without thrashing.

import numpy as np
def sgd_step(params, grads, lr=1e-3):
    for p, g in zip(params, grads):
        p -= lr * g
    return params

Momentum augments gradient descent with a velocity term: v_{t+1} = βv_t + ∇f(x_t); x_{t+1} = x_t - ηv_{t+1}. The momentum coefficient β (typically 0.9) accumulates past gradients, dampening oscillations in high-curvature directions and accelerating movement along consistent gradient directions. For a quadratic with condition number κ = L/μ, gradient descent converges in O(κ) steps; momentum reduces this to O(√κ). Physically, the iterate gains 'inertia'—noise in stochastic settings is averaged out by the accumulated velocity, and the effective step in consistent directions is ≈η/(1-β) times larger. Heavy ball (Polyak momentum) and Nesterov momentum are the two classical variants.

def sgd_momentum(params, grads, velocity, lr=1e-3, beta=0.9):
    for i in range(len(params)):
        velocity[i] = beta*velocity[i] + grads[i]
        params[i] -= lr * velocity[i]
    return params, velocity

Standard (Polyak) momentum computes the gradient at the current position x_t, then updates: v_{t+1} = βv_t + η∇f(x_t); x_{t+1} = x_t - v_{t+1}. NAG computes the gradient at the 'look-ahead' position x_t - βv_t (where momentum would take the iterate): v_{t+1} = βv_t + η∇f(x_t - βv_t); x_{t+1} = x_t - v_{t+1}. This look-ahead gradient is more informative: it corrects the trajectory before overshooting. For convex functions, NAG achieves the optimal O(1/T²) convergence rate matching the lower bound for first-order methods, compared to O(1/T) for gradient descent without momentum. PyTorch implements NAG via nesterov=True in SGD.

AdaGrad adapts the learning rate for each parameter individually by dividing by the accumulated sum of squared past gradients: G_t = G_{t-1} + g_t²; θ_{t+1} = θ_t - η·g_t/√(G_t + ε). Parameters with large historical gradients get smaller effective learning rates; rarely-updated parameters (sparse features) get larger rates. This is beneficial for sparse data (NLP, recommendation). The limitation: G_t only grows monotonically, causing the effective learning rate to shrink to zero, stopping learning in long runs. This motivated RMSprop and Adam, which use exponential moving averages of squared gradients instead of the cumulative sum, preventing the learning rate from decaying to zero.

import numpy as np
def adagrad_step(param, grad, G, lr=0.01, eps=1e-8):
    G += grad**2
    param -= lr * grad / (np.sqrt(G) + eps)
    return param, G

Adam (Adaptive Moment Estimation) maintains exponential moving averages of gradients (first moment m_t) and squared gradients (second moment v_t). Update: m_t = β₁m_{t-1} + (1-β₁)g_t; v_t = β₂v_{t-1} + (1-β₂)g_t². Bias correction: m̂_t = m_t/(1-β₁^t); v̂_t = v_t/(1-β₂^t). Parameter update: θ_{t+1} = θ_t - η·m̂_t/(√v̂_t + ε). Bias correction compensates for initialization at zero: without it, early estimates of m and v are biased toward zero, causing undersized steps. Defaults: β₁=0.9, β₂=0.999, ε=1e-8. Adam is the de-facto optimizer for deep learning due to robustness to hyperparameter choice, but may generalize slightly worse than tuned SGD with momentum.

import numpy as np
def adam_step(p, g, m, v, t, lr=1e-3, b1=0.9, b2=0.999, eps=1e-8):
    m = b1*m + (1-b1)*g
    v = b2*v + (1-b2)*g**2
    mh = m/(1-b1**t)
    vh = v/(1-b2**t)
    p -= lr * mh/(np.sqrt(vh)+eps)
    return p, m, v

L2 regularization in Adam adds λθ to the gradient, but Adam's per-parameter adaptive scaling divides this regularization term by √v̂+ε, making the effective weight decay smaller for parameters with large gradient variance. This decouples weight decay from the gradient magnitude adaptation, breaking the intended regularization. AdamW (Adam with decoupled Weight decay) applies weight decay directly to the parameter update: θ_{t+1} = (1-ηλ)θ_t - η·m̂_t/(√v̂_t+ε). This restores the clean regularization effect. AdamW is the standard optimizer in modern transformers (GPT, BERT, ViT) and has been shown to generalize better than L2-regularized Adam across a wide range of architectures.

Learning rate warmup linearly increases the learning rate from near zero to the target value over the first W steps before applying the main schedule. At the start of training, parameters are randomly initialized and gradients are high-variance. A large initial learning rate causes instability and divergence, especially for large models and batches where the second moment estimates in Adam are poorly initialized (denominator is small, causing large steps). Warmup stabilizes early training. Common schedule: linear warmup for W steps, then cosine decay or polynomial decay. Transformers typically use W=4000 steps (original Attention is All You Need) or W proportional to 1% of total steps.

def warmup_cosine_schedule(step, warmup_steps, total_steps, lr_max):
    if step < warmup_steps:
        return lr_max * step / warmup_steps
    progress = (step - warmup_steps) / (total_steps - warmup_steps)
    import math
    return lr_max * 0.5 * (1 + math.cos(math.pi * progress))

KKT conditions are necessary (and sufficient under constraint qualifications) for a point x* to be a local minimum of min f(x) subject to g_i(x) ≤ 0 and h_j(x) = 0. The conditions are: (1) Stationarity: ∇f(x*) + Σλ_i∇g_i(x*) + Σμ_j∇h_j(x*) = 0. (2) Primal feasibility: g_i(x*) ≤ 0, h_j(x*) = 0. (3) Dual feasibility: λ_i ≥ 0. (4) Complementary slackness: λ_i·g_i(x*) = 0 for all i (either the constraint is active or its multiplier is zero). For convex optimization problems, KKT conditions are necessary and sufficient. Used in SVM derivation (support vectors satisfy complementary slackness).

For minimize f(x) subject to h(x) = 0, the Lagrangian is L(x, λ) = f(x) + λh(x). At a constrained minimum x*, the gradient of the Lagrangian must be zero: ∇_x f(x*) + λ∇h(x*) = 0, meaning the gradient of f must be parallel to the gradient of h. The scalar λ (Lagrange multiplier) adjusts the constraint gradient to cancel f's gradient. For inequality constraints, the Lagrangian becomes L(x, λ) = f(x) + Σλ_i g_i(x) with λ_i ≥ 0. This transforms a constrained problem into finding saddle points of L—the basis of dual methods. Lagrange multipliers are the shadow prices of constraints in economics.

The Lagrangian dual problem maximizes the dual function g(λ,μ) = min_x L(x,λ,μ) over dual variables. Weak duality: the dual objective g(λ,μ) ≤ f(x) for any primal feasible x and dual feasible (λ,μ)—the dual always provides a lower bound on the primal. Strong duality: the optimal dual value equals the optimal primal value (duality gap is zero). Strong duality holds under Slater's condition for convex problems: if there exists a strictly feasible point (satisfying all inequality constraints strictly), strong duality holds. Strong duality justifies solving the (often easier) dual problem instead. SVMs exploit this: the dual has only support vectors as variables.

Newton's method uses the quadratic approximation of f at x_t: minimize the second-order Taylor expansion. The update is x_{t+1} = x_t - H_t^{-1} ∇f(x_t), where H_t is the Hessian at x_t. Convergence is quadratic near the optimum: ||x_{t+1} - x*|| ≤ C||x_t - x*||². Starting close enough to x*, the number of correct digits roughly doubles per iteration. The cost per step is O(n²) to compute the Hessian and O(n³) to solve the linear system, making it impractical for large n (deep learning has millions of parameters). Modified Newton methods add damping or use line search for global convergence. Newton's method is scale-invariant—unaffected by ill-conditioning.

import numpy as np
def newton_step(x, grad_f, hessian_f):
    g = grad_f(x)
    H = hessian_f(x)
    return x - np.linalg.solve(H, g)  # O(n^3)

BFGS (Broyden-Fletcher-Goldfarb-Shanno) is a quasi-Newton method that approximates the inverse Hessian H_t^{-1} using only gradient evaluations, avoiding the O(n³) cost of computing and inverting the exact Hessian. It maintains a positive definite approximation B_t that is updated at each step using the secant condition: B_{t+1}(s_t) = y_t, where s_t = x_{t+1}-x_t and y_t = ∇f(x_{t+1})-∇f(x_t). The BFGS rank-2 update preserves positive definiteness. Convergence is superlinear. Memory cost is O(n²) for the dense approximation. L-BFGS uses the last m gradient difference pairs (m≈20) for O(mn) memory—the standard for large-scale ML problems like logistic regression and shallow neural nets.

A line search finds an acceptable step size η along a descent direction p_t: x_{t+1} = x_t + ηp_t. Exact line search minimizes f(x_t + ηp_t) exactly (expensive). Inexact line search uses conditions to accept sufficient decrease. Wolfe conditions: (1) Sufficient decrease (Armijo): f(x_t+ηp_t) ≤ f(x_t) + c₁η∇f_t·p_t, c₁=10⁻⁴. (2) Curvature condition: ∇f(x_t+ηp_t)·p_t ≥ c₂∇f_t·p_t, c₂=0.9. The curvature condition ensures the accepted step is not too small (that the slope hasn't reduced enough). Strong Wolfe conditions additionally bound the slope from above. Wolfe conditions are required for BFGS to maintain positive definiteness of the Hessian approximation.

Conjugate Gradient (CG) solves symmetric positive definite linear systems Ax=b in at most n steps (exact arithmetic) by generating search directions that are A-conjugate (p_i^T A p_j = 0 for i≠j). Applied to minimizing a quadratic f(x) = ½x^TAx - b^Tx, CG achieves the exact minimum in n steps. For nonlinear optimization, nonlinear CG (Fletcher-Reeves, Polak-Ribière) extends the conjugacy idea with periodic restarts. CG is optimal for quadratic problems among all first-order methods: it cannot be improved without using curvature information. Convergence rate depends on condition number κ: O(√κ) iterations to reduce error by a constant factor. Used for large-scale linear systems in ML (Gaussian processes, kernel methods).

Proximal gradient applies to problems f(x) + g(x) where f is smooth (has Lipschitz gradient) and g is convex but possibly non-smooth (e.g., L1 norm). Gradient descent would require f+g to be smooth. Proximal gradient alternates: gradient step on f: y = x - η∇f(x); proximal step on g: x' = prox_{ηg}(y) = argmin_z (g(z) + ||z-y||²/(2η)). For g = λ||·||₁, the proximal operator is soft-thresholding: prox(y)_i = sign(y_i)·max(|y_i|-ηλ, 0). ISTA is the basic version; FISTA adds Nesterov momentum achieving O(1/T²). Used for LASSO, elastic net, and sparse deep learning applications.

import numpy as np
def soft_threshold(y, lam):
    return np.sign(y) * np.maximum(np.abs(y) - lam, 0)

def ista_step(x, grad_f_x, lr, lam):
    y = x - lr * grad_f_x
    return soft_threshold(y, lr * lam)

Frank-Wolfe minimizes f(x) over a convex set C by: at each step, finding s_t = argmin_{s∈C} ∇f(x_t)·s (a linear minimization oracle, often easier than projection), then taking x_{t+1} = (1-γ_t)x_t + γ_t·s_t with γ_t = 2/(t+2). The iterate stays in C by convex combination. Convergence is O(1/T) for convex f, same as projected gradient. Advantage: avoids projections, which for complex convex sets (positive semidefinite matrices, nuclear norm balls, polytopes) can be O(n³) while the linear oracle is O(n²) or O(n log n). Additionally, Frank-Wolfe produces sparse solutions naturally when C is a polytope (vertices are sparse). Used in matrix completion and sparse recovery.

RMSprop (Root Mean Square Propagation) replaces AdaGrad's cumulative sum of squared gradients with an exponential moving average: E[g²]_t = ρ·E[g²]_{t-1} + (1-ρ)g_t². Parameter update: θ_{t+1} = θ_t - η·g_t/√(E[g²]_t + ε). The exponential moving average (typically ρ=0.9) discounts older gradients, keeping the effective learning rate from decaying to zero in long training runs. This is suitable for non-stationary problems and RNNs where gradient magnitudes change over time. RMSprop was proposed by Hinton in 2012 (unpublished Coursera lecture); Adam can be viewed as RMSprop with bias-corrected first-moment estimation (momentum).

import numpy as np
def rmsprop_step(p, g, cache, lr=1e-3, rho=0.9, eps=1e-8):
    cache = rho*cache + (1-rho)*g**2
    p -= lr * g / (np.sqrt(cache) + eps)
    return p, cache

Gradient clipping prevents exploding gradients by rescaling the gradient vector when its norm exceeds a threshold τ: if ||g|| > τ, set g = τ·g/||g||. This is norm clipping (clip the entire gradient vector uniformly). Alternatively, value clipping clips each component independently. RNNs suffer from exploding gradients because BPTT multiplies gradients through T time steps, and if the recurrent weight matrix W has spectral radius > 1, the gradient grows exponentially. Without clipping, parameter updates become huge, causing NaN losses. The threshold τ is typically 1.0 or 5.0. LSTM and GRU architectures mitigate vanishing gradients but not exploding; clipping remains necessary for long sequences.

import numpy as np
def clip_grad_norm(grads, max_norm=1.0):
    total_norm = np.sqrt(sum((g**2).sum() for g in grads))
    clip_coef = max_norm / (total_norm + 1e-6)
    if clip_coef < 1.0:
        grads = [g * clip_coef for g in grads]
    return grads

Cosine annealing decays the learning rate following a cosine curve from η_max to η_min: η_t = η_min + ½(η_max - η_min)(1 + cos(πt/T)). At t=0, η=η_max; at t=T, η=η_min. The smooth, non-monotonic-in-gradient nature avoids abrupt learning rate changes. SGDR (with restarts) resets t to 0 periodically, allowing the optimizer to escape local minima and explore new regions. Cosine annealing tends to find flatter minima that generalize better than step decay (sharp drops), because it keeps the learning rate high for longer before a smooth decay. Used universally in vision and NLP training, often combined with warmup.

import math
def cosine_anneal(t, T, lr_max, lr_min=0.):
    return lr_min + 0.5*(lr_max-lr_min)*(1+math.cos(math.pi*t/T))

The condition number κ(A) = σ_max/σ_min (ratio of largest to smallest singular value, or eigenvalue for symmetric positive definite A) measures sensitivity of the linear system Ax=b to perturbations. For gradient descent on a quadratic f(x) = ½x^TAx - b^Tx, convergence rate is ((κ-1)/(κ+1))^t. High κ means gradient descent converges slowly because the loss landscape is a elongated ellipse—the gradient direction is nearly perpendicular to the optimal direction. Preconditioning (transforming to a better-conditioned system) mitigates this. In neural networks, batch normalization and weight normalization reduce effective condition numbers, explaining their acceleration of training.

import numpy as np
A = np.random.randn(10,10)
A = A.T @ A  # symmetric positive definite
kappa = np.linalg.cond(A)
print(f'Condition number: {kappa:.1f}')
# convergence rate per step approx:
rho = (kappa-1)/(kappa+1)
print(f'Steps for 1e-6 accuracy: {int(6/(-np.log10(rho)))+1}')

Convex optimization: all local minima are global; gradient methods converge to the global optimum. Non-convex optimization (typical in deep learning): multiple local minima, saddle points, and flat regions; gradient methods converge to a stationary point not guaranteed to be globally optimal. Empirically, deep learning loss landscapes have many saddle points but few poor local minima—most local minima are similarly good (Choromanska et al.). Quasiconvex: every sublevel set {x: f(x) ≤ c} is convex; local minima are global but the gradient condition ∇f=0 is not sufficient for minimality. Quasiconvex functions include log-sum-exp, cumulative distribution functions, and ratio of convex functions.

An LP minimizes c^Tx subject to Ax ≤ b, x ≥ 0. The feasible set is a convex polytope; the optimal solution always lies at a vertex (corner). The simplex method moves from vertex to vertex along edges, improving the objective. It terminates at the optimal vertex because the polytope has finitely many vertices. The simplex method is not polynomial in the worst case (exponential adversarial examples exist) but is extremely fast in practice (polynomial average case). Interior-point methods solve LPs in polynomial time O(n^3.5 L) where L is input bit length. LP arises in ML for SVM training, basis pursuit (L1 minimization), and resource allocation problems.

Coordinate descent minimizes f by cycling through coordinates (or features), minimizing exactly along each coordinate while fixing all others. Each coordinate update is often closed-form and cheap. Convergence per epoch is O(n·d) same as gradient descent, but inner loops are much cheaper when each coordinate minimization has a closed form. Preferred when: the full gradient is expensive but per-coordinate updates are cheap (LASSO with cyclic coordinate descent), features are very high-dimensional and sparse, or the problem decouples into independent subproblems per coordinate. Stochastic coordinate descent (randomly sample coordinate) works well for sparse problems. Not ideal for strongly coupled variables.

import numpy as np
def lasso_coord_descent(X, y, lam, max_iter=1000):
    n,d = X.shape
    w = np.zeros(d)
    for _ in range(max_iter):
        for j in range(d):
            r = y - X @ w + X[:,j]*w[j]
            rho = X[:,j] @ r / n
            w[j] = np.sign(rho)*max(abs(rho)-lam,0) / ((X[:,j]**2).mean())
    return w

The subgradient method extends gradient descent to non-smooth convex functions. A subgradient g at x satisfies f(y) ≥ f(x) + g·(y-x) for all y—a generalized gradient. For smooth points, the subgradient equals the gradient; at non-smooth points (e.g., the kink in |x| at x=0), multiple subgradients exist. Update: x_{t+1} = x_t - η_t·g_t, where g_t is any subgradient. Unlike gradient descent, subgradient iterates do not necessarily decrease f at each step; the method tracks the best iterate seen. Convergence for convex f with step η_t = c/√t: O(1/√T). Slower than gradient descent on smooth functions but applicable to L1 regularization, hinge loss, and max functions.

Trust region methods build a local quadratic model m_t(p) = f(x_t) + ∇f_t·p + ½p^TH_tp of f and minimize it subject to ||p|| ≤ Δ (the trust region radius). If the ratio ρ = (f(x_t)-f(x_t+p))/(m_t(0)-m_t(p)) is high (model was accurate), accept the step and expand Δ. If low, reject and shrink Δ. Unlike line search (which fixes direction and searches along it), trust region jointly determines both direction and step size. Trust region methods are more robust to negative curvature and can escape saddle points more reliably. The Cauchy point gives a gradient descent step constrained to the trust region; the Dogleg method refines this.

Mirror descent generalizes gradient descent by replacing the Euclidean projection onto the feasible set with a Bregman projection: x_{t+1} = argmin_x {η∇f(x_t)·(x-x_t) + D_φ(x, x_t)}, where D_φ(x,y) = φ(x) - φ(y) - ∇φ(y)·(x-y) is the Bregman divergence for a strictly convex mirror function φ. With φ(x) = ½||x||², D_φ is squared Euclidean distance, recovering gradient descent. With φ(x) = Σx_i log x_i (negative entropy) on the simplex, the update is multiplicative weights (exponentiated gradient): x_{t+1} ∝ x_t · exp(-η∇f_t). Mirror descent exploits the geometry of the constraint set, achieving tighter bounds for problems on the simplex or nuclear norm ball.

Online learning processes examples sequentially: at each round t, predict x_t, observe loss f_t(x_t), update. Regret after T rounds measures performance vs. the best fixed decision in hindsight: R_T = Σf_t(x_t) - min_x Σf_t(x). Online gradient descent achieves R_T = O(√T) for convex losses. Algorithms with sublinear regret (R_T/T → 0) are competitive with the best fixed strategy. Follow The Regularized Leader (FTRL) and AdaGrad are online algorithms with strong regret bounds. The online-to-batch conversion shows that the average of online iterates converges to the offline ERM minimizer at rate O(1/√T). Online learning models adversarial settings where data distributions can shift.

The Polyak step size η_t = (f(x_t) - f*)/(||∇f(x_t)||² + ε) uses knowledge of the optimal value f* to set the step size adaptively. It is theoretically optimal for subgradient methods on convex functions when f* is known, achieving O(1/T) convergence without needing to tune learning rates. The step size is large when far from the optimum (large f-f*) and small near convergence (small numerator). In practice, f* is unknown; heuristics use f_best (best value seen so far) or f* ≈ 0 for unnormalized losses. Variants appear in SVM training (SMO) and Polyak averaging. Rarely used in deep learning but important in theoretical analysis.

SVRG reduces gradient variance in stochastic optimization. Every m steps, compute a full gradient ∇f(x̃) at a snapshot x̃. At each stochastic step, use the variance-reduced estimator: g_t = ∇f_i(x_t) - ∇f_i(x̃) + ∇f(x̃). The correction term reduces variance because the stochastic gradient difference ∇f_i(x_t) - ∇f_i(x̃) vanishes when x_t approaches x̃. SVRG achieves linear convergence O(log(1/ε)) for strongly convex functions—exponentially faster than SGD's O(1/ε). Full gradient evaluation every m=2n steps adds 2n/m=2 full gradient equivalents per stochastic step, giving a constant-factor overhead vs. SGD but dramatically faster convergence to high accuracy.

Natural gradient (Amari) uses the Fisher Information Matrix F (expected outer product of score function gradients) as a Riemannian metric on the space of probability distributions. The natural gradient step θ_{t+1} = θ_t - η·F^{-1}∇L(θ) is invariant to reparameterization—unlike gradient descent, it does not depend on how parameters are encoded. F captures the geometry of the model distribution, making the update correspond to the direction of steepest descent on the KL divergence manifold. Computing F^{-1} is O(n²) for storage and O(n³) for inversion. KFAC approximates F as a Kronecker product of smaller matrices for practical use in deep networks.

Batch gradient descent converges at O(1/T) for convex functions (linear convergence for strongly convex). SGD converges at O(1/√T) for convex and O(1/T) for strongly convex (with noise floor, no true convergence). Variance-reduced methods (SVRG, SARAH, SPIDER) achieve linear convergence O(log(1/ε)) for strongly convex—same as batch GD but with O(n + κ) cost per epoch vs O(nκ) for batch. For non-convex: batch GD finds ε-stationary point in O(1/ε) gradient evaluations; SGD in O(1/ε²). The total oracle complexity (gradient evaluations to reach ε accuracy) is the key metric: SVRG requires O((n+κ)log(1/ε)) vs SGD O(1/(με)) and batch O(n/ε).

Polyak-Ruppert averaging maintains the running average of SGD iterates: x̄_T = (1/T)Σx_t. While the last iterate of SGD with decaying step size oscillates and may have high variance, the average converges to the optimal solution at the minimax optimal rate O(1/(μT)) for strongly convex objectives. The averaged iterate has lower variance than the final iterate without slowing convergence—it smooths out the noise in the trajectory. Stochastic weight averaging (SWA) is a modern variant: average model weights from the last portion of training to find flatter, better-generalizing minima. SWA is essentially Polyak averaging applied to neural network training.

import numpy as np
def sgd_with_averaging(grad_f, x0, n_iter, lr):
    x = x0.copy(); x_avg = x0.copy()
    for t in range(1, n_iter+1):
        x = x - lr / np.sqrt(t) * grad_f(x)
        x_avg = x_avg + (x - x_avg) / t  # running average
    return x_avg

The proximal point algorithm iterates x_{t+1} = prox_{ηf}(x_t) = argmin_y {f(y) + ||y-x_t||²/(2η)}, computing an implicit gradient step. Unlike explicit gradient descent (x_{t+1} = x_t - η∇f(x_t) using the gradient at x_t), the proximal step uses the gradient at the new point x_{t+1}. This implicit method is unconditionally convergent for any η > 0, unlike explicit gradient descent which requires η ≤ 1/L. The proximal point algorithm is theoretically optimal and motivates many practical algorithms: ADMM is a proximal method for separable problems; prox operators for common regularizers (L1, L2, nuclear norm) have closed forms, enabling efficient implementation.

Lookahead (Zhang et al. 2019) wraps any base optimizer (Adam, SGD). It maintains two sets of weights: 'fast weights' updated by the base optimizer for k steps, and 'slow weights' updated by interpolating toward the fast weights: θ_slow = θ_slow + α(θ_fast - θ_slow), then resetting fast weights to slow weights. This outer loop (k=5, α=0.5 typically) reduces variance: the fast weights explore locally while slow weights track the average direction. Lookahead with SGD matches or exceeds Adam on many tasks with less hyperparameter sensitivity. The slow weight trajectory avoids oscillation in poorly-conditioned landscapes and tends to find flatter minima that generalize better.

class Lookahead:
    def __init__(self, base_opt, k=5, alpha=0.5):
        self.base=base_opt; self.k=k; self.alpha=alpha
        self.slow={n:p.clone() for n,p in base_opt.param_groups[0]['params']}
        self.step_count=0
    def step(self):
        self.base.step(); self.step_count+=1
        if self.step_count%self.k==0:
            for p,sp in zip(self.base.param_groups[0]['params'],self.slow.values()):
                sp.data.add_(self.alpha*(p.data-sp.data))
                p.data.copy_(sp.data)

Online gradient descent (OGD) on a sequence of convex loss functions f_1,...,f_T achieves regret bound R_T ≤ D²/(2η) + η·Σ||g_t||²/2, where D is the diameter of the feasible set and g_t = ∇f_t(x_t). With constant step η = D/(G√T) (where G bounds gradient norms), regret is O(DG√T). This is optimal for general convex losses—no online algorithm achieves o(√T) regret against adversarial loss sequences. For strongly convex losses, OGD with η_t = 1/(μt) achieves O(log T) regret. These bounds are achieved without any distributional assumptions on the loss sequence, making online gradient descent robust to distribution shift and adversarial perturbations.

L-smoothness means ||∇f(x)-∇f(y)|| ≤ L||x-y|| (bounded gradient change), implying f(y) ≤ f(x)+∇f(x)·(y-x)+L||y-x||²/2. μ-strong convexity means f(y) ≥ f(x)+∇f(x)·(y-x)+μ||y-x||²/2. For gradient descent on a μ-strongly convex, L-smooth function, convergence is linear: ||x_t - x*||² ≤ (1-μ/L)^t ||x_0-x*||². The condition number κ=L/μ determines the rate—high κ (ill-conditioned) slows convergence. Convex but not strongly convex (μ=0): O(1/T) convergence. Strongly convex: O(exp(-t/κ)). These rates justify preconditioning (reducing κ) and Nesterov acceleration (achieving O(exp(-t/√κ))).

import numpy as np
# Check L-smoothness empirically via gradient Lipschitz constant
def estimate_L(grad_f, x, n_samples=1000, eps=1e-4):
    L_est = 0
    for _ in range(n_samples):
        d = np.random.randn(*x.shape)
        y = x + eps*d/np.linalg.norm(d)
        L_est = max(L_est, np.linalg.norm(grad_f(y)-grad_f(x))/eps)
    return L_est

Label smoothing replaces one-hot targets y_i ∈ {0,1} with soft targets ỹ_i = y_i(1-ε) + ε/K for K classes. Instead of assigning full probability 1.0 to the correct class, the model is trained toward ε/K for incorrect classes. This prevents overconfidence: without smoothing, cross-entropy loss pushes logits toward ±∞ to achieve probability 1, causing gradient saturation and poor calibration. Label smoothing reduces this over-optimization, acts as a regularizer by reducing the effective margin between logit differences, and improves calibration. With ε=0.1 (Szegedy et al.), it improved ImageNet top-1 accuracy on Inception models. Interestingly, label smoothing hurts knowledge distillation.

import numpy as np
def label_smooth_loss(logits, targets, num_classes, eps=0.1):
    n = logits.shape[0]
    log_probs = logits - np.log(np.sum(np.exp(logits - logits.max(1,keepdims=True)), axis=1, keepdims=True)) - logits.max(1,keepdims=True)
    soft_targets = np.full_like(log_probs, eps/num_classes)
    soft_targets[np.arange(n), targets] += 1.0 - eps
    return -np.sum(soft_targets * log_probs) / n

In high-dimensional non-convex optimization (deep networks), saddle points (∇f=0, indefinite Hessian) are far more common than local minima. Standard gradient descent slows dramatically near saddle points because the gradient is near zero in all directions. The saddle-free Newton method modifies Newton's step by using |H| (absolute values of eigenvalues) instead of H: x_{t+1} = x_t - |H|^{-1}∇f. This escapes saddle points along negative-curvature directions. In practice, stochastic noise in SGD provides natural escape from saddle points—gradient estimates have noise that perturbs the trajectory off the saddle manifold. Adaptive methods (Adam) also help escape saddle points faster than vanilla SGD.

Armijo backtracking starts with a large step size η=1, then halves it until the sufficient decrease condition is met: f(x_t - η∇f(x_t)) ≤ f(x_t) - c₁·η·||∇f(x_t)||², where c₁=10⁻⁴. This ensures the function value decreases by at least a fraction c₁ of the predicted decrease from the gradient. Backtracking terminates in O(log(1/η_min)) iterations. Unlike Wolfe conditions, Armijo does not require the curvature condition, making it simpler but potentially choosing steps that are too small in nonlinear regions. Armijo backtracking guarantees global convergence for gradient descent on smooth objectives without requiring knowledge of the Lipschitz constant L.

def armijo_backtrack(f, grad_f, x, p, c1=1e-4, rho=0.5):
    eta = 1.0
    f0 = f(x)
    g0 = grad_f(x)
    while f(x + eta*p) > f0 + c1*eta*(g0 @ p):
        eta *= rho
    return eta

Mixed precision training (Micikevicius et al. 2018) stores weights and activations in FP16 (16-bit float), halving memory and doubling bandwidth. Forward and backward passes use FP16 arithmetic; NVIDIA Tensor Cores accelerate FP16 matrix multiplications by 8× over FP32. Three techniques maintain accuracy: (1) FP32 master weights updated by the optimizer, copied to FP16 for each forward pass (avoids precision loss in small weight updates). (2) Loss scaling: multiply loss by a large factor (e.g., 2¹³) before backward pass to prevent FP16 underflow in gradients; unscale before optimizer step. (3) Keep batch normalization and loss computation in FP32. Enables training larger models or larger batches within the same GPU memory.

ReduceLROnPlateau monitors a validation metric (e.g., validation loss) and reduces the learning rate by a factor (e.g., 0.1) when no improvement is seen for 'patience' epochs (e.g., 10). It is adaptive: the schedule reacts to training progress rather than a fixed time. Use it when: training duration is not predetermined, convergence is unpredictable, or the dataset/task is new. Prefer cosine annealing when: training duration is fixed, you want warm restarts (SGDR), or reproducibility is critical. ReduceLROnPlateau is problematic for large batch training (validation metrics are noisier) and when the plateau is temporary (premature reduction). It is the most common schedule for fine-tuning pre-trained models.

Lagrangian relaxation dualizes difficult constraints by incorporating them into the objective with multipliers λ: L(x,λ) = f(x) + λ·g(x). The relaxed problem (minimize L over x for fixed λ) is often easier to solve than the original. The Lagrangian dual g(λ) = min_x L(x,λ) provides a lower bound for any λ ≥ 0. The dual problem max_{λ≥0} g(λ) finds the tightest lower bound. For integer programs, the LP relaxation solved at each dual iteration gives tractable subproblems. In ML: constrained fair learning (fairness constraints relaxed with multipliers), domain adaptation (distribution matching via moment constraints), and rate-distortion optimization in compression. Subgradient method updates λ using the constraint violation.

In non-convex optimization (e.g., deep learning), gradient methods converge to stationary points, not necessarily the global optimum. However, the neural loss landscape empirically has saddle points more than poor local minima, and wide local minima often generalize better than sharp ones. Heuristics to avoid poor local optima: (1) Random restarts from multiple initializations. (2) SGD noise provides implicit annealing, escaping shallow basins. (3) Warm restarts (SGDR) cyclically raise and lower the learning rate. (4) Stochastic weight averaging (SWA) finds wider minima. (5) Sharpness-Aware Minimization (SAM) explicitly seeks flat minima by perturbing weights to find the worst-case loss in a neighborhood.

First-order methods use only gradient information: gradient descent, SGD, Adam, RMSprop. Cost per step: O(d) for computing the gradient (with autograd). They require O(1) per step beyond gradient computation. Second-order methods use curvature (Hessian) information: Newton's method, BFGS, L-BFGS, natural gradient. They converge faster (quadratic for Newton's) but each step costs O(d²) to O(d³). In deep learning (d ~ 10⁸ parameters), second-order methods are impractical directly. Approximate second-order methods (K-FAC, Shampoo, Sophia) reduce cost while retaining curvature benefits. The fundamental tradeoff: first-order = cheap steps, more steps; second-order = expensive steps, fewer steps.

SWA (Izmailov et al. 2018) averages model weights collected at regular intervals during SGD training (typically the last portion with a cyclically varying or constant learning rate). The averaged weights lie in the middle of a flat region (basin) that multiple individually-good solutions all share. Flat minima generalize better because small perturbations to the weights cause small changes in loss; they correspond to a larger effective posterior in Bayesian terms. SWA with a cyclic learning rate collects 'diverse' weight snapshots near the basin, and their average sits at the wide center. SWA improves generalization across vision, NLP, and reinforcement learning with negligible computational overhead over standard SGD.

SAM (Foret et al. 2021) minimizes the maximum loss within a neighborhood of the current parameters: min_w max_{||ε||≤ρ} L(w+ε). This seeks flat minima by penalizing sharp regions where the worst-case perturbation still has high loss. Implementation: (1) Compute the ascent direction ε̂ = ρ·∇L(w)/||∇L(w)||. (2) Evaluate gradient at w+ε̂. (3) Update w using this gradient. Requires two forward-backward passes per step, doubling cost. SAM significantly reduces the sharpness metric and improves test accuracy, especially on larger models and datasets. ASAM (adaptive SAM) normalizes the perturbation per-parameter, accounting for parameter scale differences.

import numpy as np
def sam_step(w, loss_fn, grad_fn, lr, rho=0.05):
    g = grad_fn(w)  # gradient at w
    eps = rho * g / (np.linalg.norm(g) + 1e-12)  # ascent step
    g_sam = grad_fn(w + eps)  # gradient at perturbed point
    w -= lr * g_sam
    return w

Adafactor (Shazeer & Stern 2018) reduces Adam's memory from O(2d) (storing m and v) to O(d) by factorizing the second moment matrix. For a 2D parameter matrix W (m×n), instead of storing a full m×n second moment matrix, Adafactor stores only row factor r (m,) and column factor c (n,), approximating v ≈ (r·c^T)/sum(c). This reduces second-moment memory by a factor of min(m,n). Additionally, Adafactor can eliminate the first moment (no m), saving another O(d). For a model with d=10¹² parameters, this saves terabytes. Adafactor also uses relative learning rate scaling (no fixed lr needed in default config), reducing hyperparameter tuning. Used in T5, PaLM training.