Optimization: convexity, duality, and the algorithms behind every trainer

Why this subject matters now

A decade ago a practitioner could treat optimization as a library call. Today the field's largest artifacts are produced by an optimizer running for months on tens of thousands of accelerators, and the difference between a well-tuned and a badly-tuned update rule is measured in millions of dollars of compute. The update rule inside every large training run, AdamW under a warmup-plus-decay schedule, is a stack of results from this subject. It combines per-coordinate adaptive scaling from the AdaGrad line of work (Duchi, Hazan, and Singer, 2011), exponential moving averages with bias correction (Kingma and Ba, 2015), decoupled weight decay (Loshchilov and Hutter, 2019), and beneath all of it the stochastic approximation framework that goes back to Robbins and Monro in 1951. A practitioner who cannot derive these updates cannot reason about why a run diverges at a given learning rate, why the loss plateaus at a level set by gradient noise, or what the newest optimizer papers (Shampoo, Lion, Muon, schedule-free methods) are actually claiming.

The continuous story is half of it. The other half is the constrained and discrete side that runs the world outside the trainer. Linear and integer programming schedule airline fleets, clear electricity markets, allocate ads, and route packages. Submodular maximization places sensors and selects data. Interior-point methods solve the convex subproblems inside model-predictive controllers at kilohertz rates. The two halves share one spine, duality. The same Lagrangian that produces the SVM dual produces the LP dual that certifies a branch-and-bound node can be pruned, and the same KKT conditions that give water-filling its closed form define the central path that interior-point methods follow. This page treats the spine as one object. Complexity-theoretic aspects of the discrete side, max-flow as an LP, NP-completeness, and the vertex-cover rounding argument, are developed on the algorithms page. The practical craft of tuning a neural-network run, schedules, warmup, and diagnosis, lives on the deep learning engineering page. Here the goal is the theory both of those stand on, derived rather than asserted.

Core theory

Convex sets and convex functions

A set \( C \subseteq \R^n \) is convex if it contains the segment between any two of its points, meaning \( x, y \in C \) and \( \theta \in [0,1] \) imply \( \theta x + (1-\theta) y \in C \). A function \( f : \R^n \to \R \) is convex if its domain is convex and for all \( x, y \) and \( \theta \in [0,1] \),

the chord lies above the graph. It is strictly convex if the inequality is strict for \( x \ne y \) and \( \theta \in (0,1) \), and \( \mu \)-strongly convex if \( f(x) - \tfrac{\mu}{2}\|x\|^2 \) is convex, which quantifies curvature from below. The single most important consequence is global. For a convex \( f \), every local minimum is a global minimum. If \( x \) were locally minimal but \( f(y) < f(x) \) somewhere, then points \( \theta y + (1-\theta) x \) arbitrarily close to \( x \) (small \( \theta \)) would satisfy \( f(\theta y + (1-\theta)x) \le \theta f(y) + (1-\theta) f(x) < f(x) \), contradicting local minimality. This is the license that lets every algorithm in the first half of this page stop when the gradient is small.

First-order characterization. Suppose \( f \) is differentiable. Then \( f \) is convex if and only if for all \( x, y \),

The tangent plane at any point is a global underestimator. Both directions deserve proof. (\( \Rightarrow \)) Convexity gives \( f(x + \theta(y-x)) \le (1-\theta) f(x) + \theta f(y) \), so

where the limit is the directional derivative. (\( \Leftarrow \)) Let \( z = \theta x + (1-\theta) y \) and apply the inequality twice, at \( z \) toward \( x \) and toward \( y \).

Multiply the first by \( \theta \), the second by \( 1-\theta \), and add. The gradient terms combine to \( \nabla f(z)\T (\theta x + (1-\theta) y - z) = 0 \), leaving \( \theta f(x) + (1-\theta) f(y) \ge f(z) \), which is the definition. As a corollary, \( \nabla f(x^\star) = 0 \) at a point of a differentiable convex function makes the underestimator constant, so \( f(y) \ge f(x^\star) \) for all \( y \). Stationarity is sufficient for global optimality, not merely necessary.

Second-order characterization. If \( f \) is twice differentiable, \( f \) is convex if and only if \( \nabla^2 f(x) \succeq 0 \) everywhere. The proof reduces to one dimension. Fix \( x \) and a direction \( d \), and let \( g(t) = f(x + t d) \), so \( g''(t) = d\T \nabla^2 f(x + t d)\, d \). If the Hessian is positive semidefinite everywhere then \( g'' \ge 0 \), and Taylor's theorem with the integral form of the remainder gives

which is the first-order condition, hence convexity. Conversely, if \( d\T \nabla^2 f(x)\, d = -2c < 0 \) for some \( x, d \), then \( f(x + t d) = f(x) + t\, \nabla f(x)\T d - c t^2 + o(t^2) \) drops below the tangent plane for small \( t \), violating the first-order condition. The same reduction shows \( \mu \)-strong convexity is equivalent to \( \nabla^2 f \succeq \mu I \), and \( L \)-smoothness (defined in the gradient-descent section) to \( \nabla^2 f \preceq L I \) for convex twice-differentiable \( f \). For a quadratic \( f(x) = \tfrac12 x\T Q x \) with symmetric \( Q \succ 0 \), the constants are exactly the extreme eigenvalues, \( \mu = \lambda_{\min}(Q) \) and \( L = \lambda_{\max}(Q) \), and their ratio \( \kappa = L/\mu \) is the condition number that will control every convergence rate below.

Operations that preserve convexity. Convexity survives exactly the operations that preserve the chord inequality, and recognizing them is how one verifies convexity of anything complicated without touching a Hessian. Nonnegative weighted sums survive. If \( f_1, f_2 \) are convex and \( \alpha, \beta \ge 0 \), then \( \alpha f_1 + \beta f_2 \) is convex, term by term. Composition with an affine map survives, since \( g(x) = f(Ax + b) \) is convex when \( f \) is, because affine maps carry segments to segments. The pointwise maximum survives too. If each \( f_i \) is convex then \( g(x) = \max_i f_i(x) \) is convex, since

and the same argument works for a supremum over an arbitrary index set, which is why the dual function derived next is always concave. Scalar composition \( h(g(x)) \) is convex when \( h \) is convex and nondecreasing with \( g \) convex, or \( h \) convex and nonincreasing with \( g \) concave. Partial minimization \( g(x) = \inf_y f(x, y) \) of a jointly convex \( f \) over a convex set in \( y \) also preserves convexity. Two consequences are used constantly in machine learning. Log-sum-exp, \( f(x) = \log \sum_i e^{x_i} \), is convex because its Hessian is \( \diag(p) - p p\T \) with \( p = \softmax(x) \), a covariance matrix and hence positive semidefinite, so softmax cross-entropy is convex in the logits. The hinge loss \( \max(0,\, 1 - y\, w\T x) \) is a pointwise max of two affine functions of \( w \), hence convex, so the SVM objective below is convex. What is not preserved is general composition. The function \( e^{-x^2} \) composes two convex functions and is not convex.

Lagrangian duality, derived from the primal

Take the standard-form problem over \( x \in \R^n \),

Nothing yet is assumed convex. The Lagrangian attaches a price to each constraint.

The construction is not arbitrary. For any fixed \( x \),

because a violated inequality \( f_i(x) > 0 \) lets \( \lambda_i \to \infty \) blow the supremum up, a violated equality does the same through \( \nu_j \) with the appropriate sign, and at a feasible point the supremum is attained by zeroing every product \( \lambda_i f_i(x) \). So the primal is exactly \( p^\star = \min_x \sup_{\lambda \ge 0, \nu} L \). The dual swaps the order. Define the dual function

a pointwise infimum of functions affine in \( (\lambda, \nu) \), hence concave regardless of any convexity in the primal, and the dual problem \( d^\star = \max_{\lambda \ge 0, \nu} g(\lambda, \nu) \), which is therefore always a convex problem even when the primal is not.

Weak duality, \( d^\star \le p^\star \), is two lines. Let \( \tilde x \) be primal feasible and \( \lambda \ge 0 \). Then

the last step because \( \lambda_i \ge 0 \) and \( f_i(\tilde x) \le 0 \) make each inequality term nonpositive while every equality term vanishes. Taking the supremum on the left and the infimum on the right gives \( d^\star \le p^\star \). Every dual feasible point is a certificate, a lower bound on the primal optimum that can be exhibited and checked without trusting the solver. Branch and bound, at the end of this page, is this observation industrialized.

Strong duality and Slater's condition. The gap \( p^\star - d^\star \ge 0 \) can be strictly positive for nonconvex problems, and even for convex ones in degenerate cases. Slater's condition rules the degeneracy out. If the primal is convex (each \( f_i \) convex, each \( h_j \) affine) and some strictly feasible \( \bar x \) exists with \( f_i(\bar x) < 0 \) for every non-affine \( f_i \), then \( d^\star = p^\star \) and the dual optimum is attained. The proof runs through a separating-hyperplane argument on the epigraph-like set \( \mathcal{A} = \{ (u, t) : \exists x \text{ with } f_i(x) \le u_i \text{ for all } i, f_0(x) \le t \} \), which is convex when the problem is. The point \( (0, p^\star) \) lies on its boundary. A supporting hyperplane there has normal \( (\lambda, \lambda_0) \) with \( \lambda \ge 0 \), and strict feasibility forces \( \lambda_0 > 0 \) (the hyperplane cannot be vertical, because a vertical hyperplane through the origin would be violated by the strictly feasible point pushing \( u \) into the negative orthant), so the normal can be rescaled into a dual solution achieving \( p^\star \). The full argument is section 5.3.2 of Boyd and Vandenberghe. What this page needs is the statement and the intuition that strict feasibility keeps the geometry from collapsing. For linear programs, strong duality holds whenever the primal is feasible and bounded, with no Slater point required.

The KKT conditions, derived

Assume strong duality holds with both optima attained, at \( x^\star \) and \( (\lambda^\star, \nu^\star) \), and all functions differentiable. Chain the definitions.

The two ends are equal, so every inequality in the chain is an equality, and each forced equality is one KKT condition. The first says \( x^\star \) attains \( \inf_x L(x, \lambda^\star, \nu^\star) \). Since the Lagrangian is differentiable, its gradient in \( x \) vanishes there, which is stationarity.

The second forced equality says \( \sum_i \lambda_i^\star f_i(x^\star) = 0 \). Each term is nonpositive, so each term is individually zero, which is complementary slackness, \( \lambda_i^\star f_i(x^\star) = 0 \) for every \( i \). A constraint is either active, \( f_i(x^\star) = 0 \), or its price is zero. Together with primal feasibility and dual feasibility \( \lambda^\star \ge 0 \), these four families are the Karush-Kuhn-Tucker conditions (Karush 1939, Kuhn and Tucker 1951). For a convex problem satisfying Slater's condition they are necessary and sufficient. Any pair satisfying them makes the chain above hold with equality, exhibiting a matched primal-dual pair. Sufficiency is what turns KKT from a certificate checker into a solution method. The two problems below are solved by guessing the active set and letting the conditions do the rest.

Gradient descent and the descent lemma

A differentiable \( f \) is \( L \)-smooth if its gradient is \( L \)-Lipschitz, \( \|\nabla f(x) - \nabla f(y)\| \le L \|x - y\| \). Smoothness buys the one inequality that every first-order convergence proof starts from, the descent lemma. For all \( x, y \),

Derivation. Write the difference along the segment \( x + t(y-x) \) with the fundamental theorem of calculus, then add and subtract \( \nabla f(x) \).

Bound the integrand with Cauchy-Schwarz and then Lipschitzness, \( \big[\nabla f(x + t(y-x)) - \nabla f(x)\big]\T (y-x) \le L t \|y-x\|^2 \), and \( \int_0^1 L t \|y-x\|^2 dt = \tfrac{L}{2}\|y-x\|^2 \). Convexity was never used. The descent lemma is why gradient descent makes progress even on nonconvex losses. Now take a gradient step \( y = x - \eta \nabla f(x) \).

Any \( \eta < 2/L \) decreases \( f \) wherever the gradient is nonzero. The right-hand side is minimized at \( \eta = 1/L \), giving guaranteed progress \( f(x_{k+1}) \le f(x_k) - \tfrac{1}{2L}\|\nabla f(x_k)\|^2 \). This already proves something useful without convexity. Summing over \( k \) iterations, \( \min_{i \le k} \|\nabla f(x_i)\|^2 \le \tfrac{2L (f(x_0) - f^\star)}{k} \), so gradient descent finds an \( \varepsilon \)-stationary point of a smooth nonconvex function in \( O(1/\varepsilon^2) \) steps. That is the honest guarantee under which all of deep learning operates.

The \( O(1/k) \) rate for convex functions. Add convexity and run \( x_{k+1} = x_k - \eta \nabla f(x_k) \) with \( \eta = 1/L \). Track the distance to a minimizer \( x^\star \).

Convexity gives \( \nabla f(x_k)\T (x_k - x^\star) \ge f(x_k) - f^\star \), and the guaranteed progress above gives \( \|\nabla f(x_k)\|^2 \le 2L\,(f(x_k) - f(x_{k+1})) \). Substituting both with \( \eta = 1/L \) gives

Rearrange and sum from \( 0 \) to \( k-1 \). The distances telescope, leaving

Since the iterates only ever decrease \( f \), the last one is at most the average, and

This is sublinear. Reaching accuracy \( \varepsilon \) needs \( O(L R^2 / \varepsilon) \) iterations, where \( R = \|x_0 - x^\star\| \). Each extra digit of accuracy costs ten times more work.

The linear rate under strong convexity, with the condition number explicit. Strong convexity converts the sublinear rate into a geometric one. First derive the Polyak-Łojasiewicz inequality from \( \mu \)-strong convexity. The strong convexity lower bound \( f(y) \ge f(x) + \nabla f(x)\T (y-x) + \tfrac{\mu}{2}\|y-x\|^2 \) holds for all \( y \), so minimize both sides over \( y \). The left side becomes \( f^\star \), and the right side is a quadratic in \( y \) minimized at \( y = x - \nabla f(x)/\mu \) with value \( f(x) - \tfrac{1}{2\mu}\|\nabla f(x)\|^2 \). Therefore

Combine with the guaranteed progress of a \( \eta = 1/L \) step.

with \( \kappa = L/\mu \) the condition number. After \( k \) steps the gap is at most \( (1 - 1/\kappa)^k (f(x_0) - f^\star) \le e^{-k/\kappa}(f(x_0) - f^\star) \), so accuracy \( \varepsilon \) costs \( k \approx \kappa \ln\big((f(x_0)-f^\star)/\varepsilon\big) \) iterations. The iteration count scales linearly in \( \kappa \). On iterate distances one can do slightly better with the step \( \eta = 2/(\mu + L) \), which balances the contraction of the smallest and largest curvature directions and yields \( \|x_{k+1} - x^\star\| \le \frac{\kappa - 1}{\kappa + 1} \|x_k - x^\star\| \) (Nesterov, Lectures, Theorem 2.1.15). Note that the PL argument never used convexity beyond the inequality itself, so any function satisfying PL, including some nonconvex ones, inherits the linear rate. This observation is the seed of a large modern literature on why overparameterized networks train faster than the worst case suggests.

Subgradients: first-order methods without differentiability

Convex functions need not be differentiable (the hinge loss, the \( \ell_1 \) norm, any pointwise max), but they always support linear underestimators. A vector \( g \) is a subgradient of \( f \) at \( x \) if

and the set of them, \( \partial f(x) \), is the subdifferential, a nonempty, closed, convex set wherever a convex \( f \) is finite in the interior of its domain. At a point of differentiability \( \partial f(x) = \{\nabla f(x)\} \), while at a kink it fattens. Take \( f(x) = |x| \) on \( \R \). For \( x > 0 \) the subdifferential is \( \{1\} \), for \( x < 0 \) it is \( \{-1\} \), and at \( 0 \) it is the whole interval \( [-1, 1] \), since \( |y| \ge g\, y \) holds precisely when \( |g| \le 1 \). Optimality generalizes cleanly, in that \( x^\star \) minimizes \( f \) iff \( 0 \in \partial f(x^\star) \). The subgradient method iterates \( x_{k+1} = x_k - \eta_k g_k \) with any \( g_k \in \partial f(x_k) \). It is not a descent method (a subgradient step can increase \( f \)), and the descent lemma is unavailable, so the analysis tracks distances only. With \( \|g_k\| \le G \) and \( R = \|x_0 - x^\star\| \), the same expansion as before gives

using the subgradient inequality in place of convexity of a gradient. Summing and rearranging,

which with the optimized constant step \( \eta = R/(G\sqrt{k}) \) yields \( \min_i f(x_i) - f^\star \le R G / \sqrt{k} \). The \( O(1/\sqrt{k}) \) rate is tight for this method, an order slower than smooth gradient descent, which is the theoretical reason smoothing and proximal operators (below) matter. They recover the smooth rate while keeping the nonsmooth term exact.

Momentum and Nesterov acceleration

Gradient descent's \( O(1/k) \) and \( O(\kappa \log(1/\varepsilon)) \) are not the best a first-order method can do. Polyak's heavy-ball method (1964) adds momentum,

and on quadratics, with \( \eta = 4/(\sqrt{L} + \sqrt{\mu})^2 \) and \( \beta = \big((\sqrt\kappa - 1)/(\sqrt\kappa + 1)\big)^2 \), contracts at rate \( (\sqrt\kappa - 1)/(\sqrt\kappa + 1) \) per step. The dependence improves from \( \kappa \) to \( \sqrt\kappa \), which for \( \kappa = 10^4 \) is the difference between ten thousand iterations per digit and one hundred. The analysis is eigenvalue-by-eigenvalue. Momentum turns each mode into a second-order recurrence whose characteristic roots are complex with modulus \( \sqrt\beta \) across the whole spectrum, so all modes decay uniformly. Heavy ball is not globally convergent on general smooth strongly convex functions. Lessard, Recht, and Packard (2016) exhibit a three-piece-quadratic counterexample where it cycles. Nesterov's accelerated gradient (1983) fixes this with a subtle change, evaluating the gradient at an extrapolated point.

With \( \beta_k = (t_k - 1)/t_{k+1} \) where \( t_1 = 1 \) and \( t_{k+1} = \tfrac12\big(1 + \sqrt{1 + 4 t_k^2}\,\big) \), the guarantee for convex \( L \)-smooth \( f \) is

The proof device is Nesterov's estimate sequence, and the idea is worth stating honestly even though the full bookkeeping is long. One maintains a sequence of simple quadratic functions \( \phi_k \) and weights \( \lambda_k \downarrow 0 \) such that two invariants hold at every step, (i) \( \phi_k(x) \le (1 - \lambda_k) f(x) + \lambda_k \phi_0(x) \), so \( \phi_k \) is squeezed toward \( f \) from above as \( \lambda_k \to 0 \), and (ii) \( f(x_k) \le \min_x \phi_k(x) \), so the algorithm's actual value stays below the moving quadratic's minimum. Together these give \( f(x_k) - f^\star \le \lambda_k (\phi_0(x^\star) - f^\star) \), and the whole game is choosing the extrapolation \( \beta_k \) so that invariant (ii) survives the update of \( \phi_k \) by mixing in the tangent lower bound at \( y_k \). The recurrence for \( t_k \) is exactly what makes \( \lambda_k = O(1/k^2) \) achievable. The complete proof is Nesterov's Lectures on Convex Optimization, section 2.2. Beck and Teboulle (2009) give a shorter induction for the composite version reproduced in the proximal section below. The momentum coefficient \( \beta_k = (t_k - 1)/t_{k+1} \to 1 \), which matches the practical observation that momentum near 0.9 to 0.99 helps. Acceleration is momentum, scheduled correctly.

The lower bound, and why acceleration is optimal. Nemirovski and Yudin (1983) initiated oracle complexity, asking how many gradient evaluations any algorithm needs, in the worst case, if all it can see are function values and gradients. For the class of algorithms whose iterates stay in \( x_0 + \mathrm{span}\{\nabla f(x_0), \dots, \nabla f(x_{k-1})\} \), Nesterov's worst-case function (Lectures, Theorem 2.1.7), a tridiagonal quadratic resembling a discretized string, satisfies

The mechanism is information propagation. The function couples adjacent coordinates in a chain, the starting gradient touches only the first coordinate, and each oracle call can light up at most one more coordinate, so after \( k \) calls the iterate lives in a \( k \)-dimensional slice whose best value is still \( \Omega(1/k^2) \) away. Matching upper and lower bounds means \( O(1/k^2) \) is the exact complexity of smooth convex first-order optimization. Nesterov's method cannot be beaten by a better idea, only by a stronger oracle (second-order information) or a smaller function class (strong convexity, where the accelerated rate becomes \( (1 - 1/\sqrt\kappa)^k \) and is again tight).

Proximal methods, ISTA, and FISTA

Composite objectives \( F(x) = f(x) + h(x) \), with \( f \) smooth and \( h \) convex but possibly nonsmooth (for the lasso, \( f = \tfrac12\|Ax - b\|^2 \) and \( h = \lambda \|x\|_1 \)), sit between the smooth and nonsmooth worlds. The subgradient method would drop to \( O(1/\sqrt k) \). The proximal view recovers \( O(1/k) \) and then accelerates. Define the proximal operator

the unique minimizer since the objective is strongly convex. The proximal gradient step replaces \( f \) by its descent-lemma upper bound around \( x_k \) and minimizes the bound plus \( h \) exactly.

where the second equality is completing the square. For \( h = 0 \) this is gradient descent. For \( h \) the indicator of a convex set it is projected gradient descent. For \( h = \lambda\|\cdot\|_1 \) it is ISTA, iterative soft thresholding, and the convergence guarantee \( F(x_k) - F^\star \le L\|x_0 - x^\star\|^2 / (2k) \) carries over with essentially the same telescoping proof as gradient descent because the prox step satisfies the same key inequality. FISTA (Beck and Teboulle, 2009) applies Nesterov's extrapolation to the prox step, with exactly the \( t_k \) recurrence above, and achieves

the accelerated rate with a nonsmooth term handled exactly. On a small lasso instance (\( A \in \R^{20 \times 5} \), \( \lambda = 1 \), \( L = \lambda_{\max}(A\T A) = 38.188 \)) the gap is visible early. At iteration 10, ISTA's objective is \( 3.69401 \) while FISTA's is \( 3.43920 \), already within \( 7 \times 10^{-3} \) of the optimum \( 3.4323500 \) that ISTA needs roughly 50 iterations to reach. Both recover exactly the two truly nonzero coordinates, with every other coordinate identically zero, the exactness being the point of prox over subgradients.

Stochastic gradient descent: noise, rates, and batching

Machine learning objectives are finite sums or expectations, \( f(x) = \tfrac1n \sum_{i=1}^n f_i(x) = \E_i[f_i(x)] \), and computing \( \nabla f \) exactly costs a full pass over the data. Robbins and Monro (1951) observed that an unbiased estimate suffices. Sample \( i_k \) uniformly and step

The price of cheapness is variance. Decompose the second moment around the mean,

and note the crucial asymmetry, that near a minimizer the signal \( \|\nabla f\|^2 \) vanishes but the noise \( \sigma^2 \) generally does not, because individual \( \nabla f_i(x^\star) \ne 0 \) even though their average is zero. A constant step size therefore converges only to a noise ball whose radius scales like \( \eta \sigma^2 \). Decaying steps are forced by the noise, not by the analysis. The convex rate comes from the same distance telescoping as the subgradient method, now in expectation. Assume \( f \) convex, \( \E\|g_k\|^2 \le G^2 \), \( R = \|x_0 - x^\star\| \), and a constant step \( \eta \). Then

where the middle term used the tower property, \( \E[g_k\T (x_k - x^\star)] = \E[\nabla f(x_k)\T (x_k - x^\star)] \ge \E[f(x_k) - f^\star] \). Unbiasedness is exactly what lets the conditional expectation slide through the inner product. Sum over \( k \) steps, telescope, divide, and apply Jensen to the average iterate \( \bar x_k = \tfrac1k \sum_{i<k} x_i \).

The two terms are the bias-variance tradeoff of the step size in one line. Big steps forget the initialization fast but rattle in the noise, small steps do the reverse, and the balance point gives the \( O(1/\sqrt k) \) rate, which Nemirovski-Yudin lower bounds show is optimal for stochastic convex optimization. No acceleration is possible against noise of constant variance. Under strong convexity the rate improves to \( O(1/(\mu k)) \) with \( \eta_k \propto 1/(\mu k) \), still sublinear. Noise costs an exponential-to-polynomial downgrade compared with the deterministic linear rate.

Mini-batching averages \( B \) independent samples per step, cutting the variance term to \( \sigma^2 / B \) while multiplying the per-step cost by \( B \). Total sample complexity is unchanged in the noise-dominated regime, and the win is parallelism, not statistics. The practical consequences, the linear-scaling rule, critical batch size, and gradient-noise-scale measurement, are treated on the deep learning engineering page. The theory here explains why they exist at all. Variance reduction is the observation that for finite sums the noise floor can be removed entirely. SVRG (Johnson and Zhang, 2013) keeps a snapshot \( \tilde x \) with its full gradient \( \nabla f(\tilde x) \), computed once per epoch, and uses the control-variate estimator

which is unbiased (the last two terms have equal expectation) and whose variance is bounded by a constant times \( f(x_k) - f^\star + f(\tilde x) - f^\star \). As both the iterate and the snapshot approach the optimum, the noise vanishes with the signal instead of persisting. The result is a linear rate for strongly convex finite sums at cost \( O\big((n + \kappa)\log(1/\varepsilon)\big) \) component gradients, versus \( O(n \kappa \log(1/\varepsilon)) \) for full-batch gradient descent. SVRG and its relatives (SAG, SAGA) dominate on convex problems like logistic regression. They have largely failed to help deep learning, where the snapshot goes stale too fast relative to the epoch cost, a negative result documented carefully in the Bottou, Curtis, and Nocedal (2018) survey and subsequent empirical studies.

Adaptive methods: AdaGrad to Adam to AdamW

SGD scales every coordinate by the same \( \eta \). AdaGrad (Duchi, Hazan, and Singer, 2011) instead divides each coordinate by the history of its own gradient magnitudes.

The motivation is a regret bound in online convex optimization. Against any comparator \( x^\star \), projected online gradient descent with a per-coordinate step \( \eta_j \) accumulates regret roughly \( \sum_j \big( R_j^2/\eta_j + \eta_j \sum_t g_{t,j}^2 \big) \). Optimizing \( \eta_j \) in hindsight gives \( \sum_j R_j \sqrt{\sum_t g_{t,j}^2} \), and AdaGrad's running normalization achieves this up to a factor of 2 without hindsight,

which is much smaller than the dimension-free \( G R \sqrt{T} \) of plain online gradient descent when gradients are sparse. Rare features keep large effective step sizes while frequent ones are damped. The weakness for deep learning is that \( v \) only grows, so the effective step decays monotonically and the method starves. RMSProp replaces the sum with an exponential moving average, and Adam (Kingma and Ba, 2015) applies EMAs to both first and second moments.

The hats are the bias correction, and it is a short derivation rather than a heuristic. Unrolling the EMA from \( m_0 = 0 \) gives \( m_k = (1-\beta_1) \sum_{t=1}^{k} \beta_1^{k-t} g_t \). If the gradient distribution were stationary with mean \( \bar g \), then

by the geometric series, and identically \( \E[v_k] = (1-\beta_2^k)\, \E[g^2] \). Early in training the EMAs are shrunk toward their zero initialization by the factor \( 1 - \beta^k \). Dividing by it, \( \hat m_k = m_k/(1-\beta_1^k) \) and \( \hat v_k = v_k/(1-\beta_2^k) \), removes the bias exactly under stationarity. The correction matters most for \( v \). With \( \beta_2 = 0.999 \), at step 10 the uncorrected \( v \) is about \( 1\% \) of its stationary size, so \( \sqrt{v} \) underestimates by \( 10\times \) and the uncorrected update would be \( 10\times \) too large, exactly when the landscape is least trusted. The practical companion to this fact, why Adam still wants learning-rate warmup, is discussed on the deep learning engineering page.

What is actually known about Adam's convergence. The original paper's regret proof contained an error, and the gap is not a technicality. Reddi, Kale, and Kumar (2018) exhibit a one-dimensional online convex problem on \( [-1, 1] \), with losses cycling as \( f_t(x) = C x \) every third step and \( f_t(x) = -x \) otherwise for a constant \( C > 2 \), on which Adam with any fixed \( \beta_1, \beta_2 \) suffers linear regret. The rare large gradient \( C \) is squashed by the EMA in \( v \) faster than it can pull the iterate toward the optimum at \( x = -1 \), and the frequent small \( -1 \) gradients win, driving \( x \) to \( +1 \). The fix in that paper, AMSGrad, enforces a nondecreasing \( \hat v \) (take the running max), restoring the regret guarantee at some empirical cost. Later analyses established convergence of Adam-family methods to stationary points in the smooth nonconvex setting under parameter conditions (for instance \( \beta_1 \) small relative to \( \beta_2 \), or increasing \( \beta_2 \)), so the honest summary is that Adam as commonly configured has no online-convex regret guarantee, does have nonconvex stationarity guarantees under conditions practitioners do not check, and is used because it works, a state of affairs the theory community regards as unfinished business.

AdamW and the decoupling argument. Adding L2 regularization \( \tfrac\lambda2\|x\|^2 \) to the loss feeds \( \lambda x \) into the gradient, where Adam divides it by \( \sqrt{\hat v} + \epsilon \) like everything else. Parameters with large recent gradients get less regularization, coupling the decay strength to gradient noise in a way L2's derivation never intended. Loshchilov and Hutter (2019) decouple it, applying decay directly to the weights outside the adaptive machinery.

For SGD the two formulations are identical up to a reparameterization of \( \lambda \). For Adam they are not, and decoupling both restores the intended geometry (every weight decays at the same relative rate) and, empirically, makes the best \( \lambda \) independent of the learning-rate schedule, which is why AdamW rather than Adam-with-L2 is the default in essentially every modern pretraining codebase. The implementation section verifies a from-scratch AdamW against torch.optim.AdamW to the last float32 digit.

Second-order methods: Newton, BFGS, L-BFGS

Newton's method minimizes the local quadratic model exactly. Around \( x_k \),

whose minimizer over \( p \) (when the Hessian is positive definite) solves \( \nabla^2 f(x_k)\, p = -\nabla f(x_k) \), giving \( x_{k+1} = x_k + p_k \). Near a minimizer with Lipschitz Hessian, convergence is quadratic, \( \|x_{k+1} - x^\star\| \le C \|x_k - x^\star\|^2 \). The number of correct digits doubles per iteration, and crucially the rate has no condition-number dependence, because the method rescales the space by the true curvature at every step. The barrier-method table in the interior-point section shows this concretely. The gradient magnitude falls from \( 2 \) to \( 0.667 \) to \( 0.075 \) to \( 0.0011 \) to \( 2.3 \times 10^{-7} \) in four steps, each exponent roughly doubling. The costs are a Hessian (\( O(n^2) \) memory) and a linear solve (\( O(n^3) \) time), which is why pure Newton is reserved for small-to-medium smooth problems and for the inner loop of interior-point methods where the Hessian has exploitable structure.

BFGS from the secant condition. Quasi-Newton methods maintain an approximation \( B_k \approx \nabla^2 f(x_k) \) built only from gradient differences. After stepping \( s_k = x_{k+1} - x_k \) and observing \( y_k = \nabla f(x_{k+1}) - \nabla f(x_k) \), the fundamental theorem of calculus gives \( y_k = \big(\int_0^1 \nabla^2 f(x_k + t s_k)\, dt\big)\, s_k \), so the averaged Hessian maps \( s_k \) to \( y_k \). Demand the same of the model, the secant condition.

One equation cannot determine a matrix, so make the smallest symmetric correction with rank two, using the natural basis for the discrepancy, \( y_k \) and \( B_k s_k \).

Imposing the secant condition gives \( B_k s_k + a\, y_k (y_k\T s_k) + b\, B_k s_k\, (s_k\T B_k s_k) = y_k \). Matching the \( y_k \) and \( B_k s_k \) components separately forces \( a = 1/(y_k\T s_k) \) and \( b = -1/(s_k\T B_k s_k) \), which is the BFGS update (Broyden, Fletcher, Goldfarb, Shanno, 1970, four independent derivations in the same year).

As a check, \( B_{k+1} s_k = B_k s_k - B_k s_k + y_k = y_k \). The update preserves positive definiteness exactly when the curvature condition \( y_k\T s_k > 0 \) holds, which a Wolfe-condition line search guarantees. This is why BFGS is always stated with Wolfe line search. Applying the Sherman-Morrison-Woodbury identity twice converts it into a direct update of the inverse \( H_k = B_k^{-1} \),

so the search direction \( -H_k \nabla f_k \) costs a matrix-vector product, no solve. BFGS converges superlinearly under standard conditions (Nocedal and Wright, chapter 6). It can also be derived, more elegantly but less transparently, as the unique minimizer of a weighted Frobenius distance to \( H_k \) subject to symmetry and the secant equation, which is the treatment in Nocedal and Wright section 6.1.

L-BFGS and the memory arithmetic. Dense \( H_k \) needs \( n^2 \) floats, which for \( n = 10^9 \) parameters is \( 4 \times 10^{18} \) bytes in fp32, four exabytes, out of the question. L-BFGS (Liu and Nocedal, 1989) never forms \( H_k \). It stores the last \( m \) pairs \( (s_i, y_i) \) and reconstructs \( H_k \nabla f \) with the two-loop recursion in \( O(mn) \) time, implicitly applying the update formula \( m \) times to a diagonal seed. The arithmetic that decides whether it fits starts from \( m \) pairs of two \( n \)-vectors each, \( 2mn \) floats. For a 7-billion parameter model with the customary \( m = 20 \), that is \( 2 \times 20 \times 7 \times 10^9 \times 4 \) bytes \( = 1.12 \) TB, versus Adam's two EMA states at \( 2 \times 7 \times 10^9 \times 4 = 56 \) GB. L-BFGS costs \( 20\times \) Adam's optimizer memory, needs coherent full-batch (or very large batch) gradients for its curvature pairs to mean anything, and its line search multiplies the number of forward-backward passes. This arithmetic, not any subtlety of theory, is why deep learning abandoned it after early experiments and why it remains the default for medium-scale smooth problems in scipy.optimize.minimize(method="L-BFGS-B"), classical MAP estimation, and full-batch scientific fitting, where it is hard to beat.

Interior-point methods: the log barrier and the central path

Constrained problems \( \min f_0(x) \) s.t. \( f_i(x) \le 0 \) resist Newton's method directly because the feasible region has corners. The log barrier replaces the hard constraints with a smooth repulsion that blows up at the boundary.

for a parameter \( t > 0 \) that sets the price of objective relative to boundary repulsion. The curve \( \{x^\star(t) : t > 0\} \) is the central path, and its defining property drops out of stationarity, \( t\, \nabla f_0(x^\star(t)) + \sum_i \frac{1}{-f_i(x^\star(t))} \nabla f_i(x^\star(t)) = 0 \). Setting \( \lambda_i(t) = \frac{1}{-t f_i(x^\star(t))} > 0 \), this is exactly KKT stationarity for the original problem with multipliers \( \lambda(t) \), and the deviation from complementary slackness is precisely \( \lambda_i(t) \cdot (-f_i(x^\star(t))) = 1/t \) per constraint. Evaluating the dual function at \( \lambda(t) \) therefore gives

so by weak duality \( f_0(x^\star(t)) - p^\star \le m/t \). Every point on the central path comes with a computable optimality certificate, and driving \( t \to \infty \) drives the gap to zero. The barrier method alternates two moves, recentering with a few Newton steps at the current \( t \), then multiplying \( t \) by a factor \( \mu \) (typically 10 to 50), starting the next recentering warm. Nesterov and Nemirovski's self-concordance theory (1994) explains why this is fast. The log barrier's third derivative is controlled by its second in the exact way that makes Newton's quadratic convergence region large enough that \( O(\sqrt m) \) recenterings per constant-factor gap reduction suffice, giving \( O(\sqrt m \log(m/(t_0 \varepsilon))) \) Newton steps overall, in practice a few dozen regardless of problem size. This line completed the arc begun by Karmarkar (1984), whose projective method was the first practical polynomial-time LP algorithm, following the Grötschel-Lovász-Schrijver development of the ellipsoid method (1981) which had settled polynomiality in principle but not in practice.

A barrier recentering is worth doing once by hand. Take the toy LP \( \min x \) s.t. \( 0 \le x \le 1 \), whose barrier objective at parameter \( t \) is \( \psi_t(x) = t x - \log x - \log(1 - x) \), with \( \psi_t'(x) = t - \tfrac1x + \tfrac{1}{1-x} \) and \( \psi_t''(x) = \tfrac{1}{x^2} + \tfrac{1}{(1-x)^2} \). At \( t = 2 \) the central point solves \( 2 - \tfrac1x + \tfrac{1}{1-x} = 0 \), a quadratic. Multiplying by \( x(1-x) \) gives \( 2x^2 - 4x + 1 = 0 \), so \( x^\star(2) = 1 - \tfrac{\sqrt2}{2} = 0.2928932 \), with certified gap \( m/t = 2/2 = 1 \). Newton from the analytic center \( x = 0.5 \) runs as follows.

The first row is pure hand arithmetic, \( \psi'(0.5) = 2 - 2 + 2 = 2 \), \( \psi''(0.5) = 4 + 4 = 8 \), step \( -2/8 = -0.25 \). From step 1 the gradient magnitudes go \( 0.667 \to 0.0748 \to 0.00109 \to 2.3 \times 10^{-7} \). The exponent doubles each iteration, quadratic convergence made visible. Five steps deliver \( x^\star(2) \) to ten digits. A production barrier method would now set \( t \leftarrow 20 \), warm-start from \( 0.2928932 \), and repeat, walking the central path toward \( x = 0 \) with the duality gap shrinking by \( 10\times \) per outer iteration.

Linear programming duality, on real numbers

Linear programming is where duality stops being abstract. Specialize the Lagrangian machinery to \( \max\, c\T x \) s.t. \( Ax \le b \), \( x \ge 0 \). Dualizing the row constraints with prices \( y \ge 0 \) produces \( \min\, b\T y \) s.t. \( A\T y \ge c \), \( y \ge 0 \), and weak duality \( c\T x \le y\T A x \le b\T y \) is the same two-line sandwich as before. Strong duality holds whenever the primal is feasible and bounded (the general derivation, and its equivalence to max-flow min-cut and Farkas' lemma, is on the algorithms page). The concrete instance used throughout the rest of this page is

a two-product planning problem, with profits 5 and 4 per unit, one machine with 24 hours consuming 6 and 4 hours per unit, and one with 6 hours consuming 1 and 2. Solving both problems with scipy.optimize.linprog (the HiGHS simplex backend) returns primal optimum \( x^\star = (3.0,\, 1.5) \) with value \( 21.0 \) and dual optimum \( y^\star = (0.75,\, 0.5) \) with value \( 21.0 \), and HiGHS reports the primal constraint marginals as \( (-0.75, -0.5) \), the dual solution with the sign convention flipped, no gap to fourteen digits. Problem 5 below reproduces both by hand and verifies complementary slackness. The dual prices are management information. An extra hour on machine 1 is worth \( 0.75 \), on machine 2 \( 0.5 \), and the simplex method's stopping certificate is exactly a feasible price vector proving no production plan can beat the current one. Dantzig's simplex method (1947, with the standard account in his 1963 book) walks vertices of the feasible polytope and remains, seventy-five years later, the workhorse inside every branch-and-bound solver because it re-solves a perturbed LP from a warm basis in a handful of pivots, something interior-point methods cannot match even though their worst-case complexity is polynomial and simplex's is exponential (Klee-Minty cubes).

Submodularity and the greedy 1 − 1/e guarantee

Some discrete problems have enough structure that a trivial algorithm carries a nontrivial guarantee. A set function \( f : 2^V \to \R \) is submodular if it has diminishing returns. For all \( S \subseteq T \) and \( e \notin T \),

and monotone if \( f(S) \le f(T) \) whenever \( S \subseteq T \). Coverage functions (how many elements do these sets cover), influence spread, entropy of a subset of random variables, and log-determinants of submatrices are all submodular. Maximizing a monotone submodular function under a cardinality constraint \( |S| \le k \) is NP-hard (it contains max coverage), and the greedy algorithm, add the element with the largest marginal gain \( k \) times, is the classical answer. Nemhauser, Wolsey, and Fisher (1978) proved it is a \( (1 - 1/e) \)-approximation, and the proof is short enough to give in full.

First comes a telescoping consequence of submodularity. For any sets \( S \) and \( T = \{t_1, \dots, t_r\} \),

each term bounded by diminishing returns, since \( S \subseteq S \cup \{t_1,\dots,t_{j-1}\} \). Now let \( S_i \) be the greedy set after \( i \) picks, \( O \) the optimal \( k \)-set, and use monotonicity, then the telescoped bound with \( S = S_i \), \( T = O \).

the last step because greedy's next gain is the largest of the at most \( k \) candidate gains in the sum. Rearrange into a contraction of the deficit \( \Delta_i = f(O) - f(S_i) \).

using \( 1 - u \le e^{-u} \) and \( f(S_0) = f(\varnothing) \ge 0 \). Therefore \( f(S_k) \ge (1 - 1/e)\, f(O) \approx 0.632\, f(O) \). The constant is tight. Nemhauser and Wolsey showed no polynomial number of value queries beats it, and Feige (1998) showed the same threshold is NP-hard to beat even for explicit max coverage.

A worked instance, verified by brute force. Universe \( \{1,\dots,6\} \), sets \( A = \{1,2,3,4\} \), \( B = \{1,2,5\} \), \( C = \{3,4,6\} \), budget \( k = 2 \). Greedy first compares gains \( 4, 3, 3 \) and picks \( A \). The remaining marginal gains are \( B: |\{5\}| = 1 \) and \( C: |\{6\}| = 1 \), so the second pick covers one more element, total \( f = 5 \). The optimum is \( \{B, C\} \), which covers all \( 6 \) elements. Greedy's first pick, locally best, forecloses the complementary pair. The exhaustive check over all \( \binom{3}{2} = 3 \) pairs confirms \( f(\{A,B\}) = 5 \), \( f(\{A,C\}) = 5 \), \( f(\{B,C\}) = 6 \). The ratio \( 5/6 = 0.833 \) respects the guarantee \( 1 - (1 - \tfrac12)^2 = 0.75 \) for \( k = 2 \) (and a fortiori the limiting \( 1 - 1/e = 0.632 \)). Constructions pushing the ratio to exactly \( 1 - 1/e \) require the universe to grow with \( k \). This pattern, greedy plus a certificate that near-optimality is the best polynomial-time outcome, is the discrete analogue of a dual bound, and it is why greedy data selection, sensor placement, and facility location ship in production without anyone solving the NP-hard problem exactly.

Integer programming: relaxation, branch and bound, cutting planes

Add integrality to the LP and the geometry breaks. The problem \( \max c\T x \) s.t. \( Ax \le b \), \( x \in \mathbb{Z}^n_{\ge 0} \) is NP-hard in general. Dropping integrality gives the LP relaxation, whose optimum bounds the integer optimum from above (for maximization) because the feasible set only grew, and the entire practice of exact integer programming is the disciplined use of that bound. Branch and bound maintains a tree of LPs. At a node, solve the relaxation. If it is infeasible, prune. If its value is at most the best integer solution found so far (the incumbent), prune. No descendant can beat what is already in hand, because descendants only add constraints. If its solution is integral, update the incumbent. Otherwise pick a fractional variable \( x_j = v \) and branch on \( x_j \le \lfloor v \rfloor \) versus \( x_j \ge \lceil v \rceil \), two child LPs that exclude the fractional point but no integer point. When the objective coefficients are integers there is a free strengthening. A node with LP value \( 20.667 \) cannot lead to an integer solution better than \( \lfloor 20.667 \rfloor = 20 \). Problem 6 traces the full tree for the running LP with integrality added. It closes after five node solves with optimum \( (4, 0) \) and value \( 20 \), and the trace shows both prune types firing.

Cutting planes attack from the other side. Instead of splitting, add a valid inequality, one satisfied by every integer point but violated by the current fractional LP optimum, and re-solve. Gomory cuts read such inequalities directly out of the simplex tableau. If a basic row says \( x_B + \sum_j \bar a_j x_j = \bar b \) with \( \bar b \) fractional, then because every feasible integer solution makes the left side integral, the fractional parts must satisfy \( \sum_j \mathrm{frac}(\bar a_j)\, x_j \ge \mathrm{frac}(\bar b) \), a new constraint the current vertex violates. In the running example the node solution \( (10/3, 1) \) is cut off by exactly such an inequality without any branching. Pure cutting-plane methods converge finitely in theory and stall in practice. Modern solvers (HiGHS, CBC, SCIP, Gurobi, CPLEX) run branch and cut, generating cuts at nodes to tighten bounds and branching when cuts stop paying, plus presolve, heuristics for incumbents, and conflict analysis. The measure of all of it is the integrality gap between the relaxation and the integer optimum, here \( 21 / 20 = 1.05 \), a 5% gap closed by one level of branching. When a relaxation's gap is provably bounded, rounding its solution yields an approximation algorithm. The canonical example, rounding the vertex-cover LP at threshold \( \tfrac12 \) for a factor-2 guarantee, is derived in full on the algorithms page, and the same recipe generalizes. In a covering IP where every constraint row \( \sum_{j \in R_i} x_j \ge 1 \) touches at most \( \varphi \) variables, at least one variable in each row has LP value \( \ge 1/\varphi \), so rounding up every \( x_j \ge 1/\varphi \) is feasible and costs at most \( \varphi \) times the LP optimum, hence at most \( \varphi \cdot \mathrm{OPT} \). That is a \( \varphi \)-approximation with the LP optimum as the certificate, vertex cover being the case \( \varphi = 2 \).

Worked problems

N0: root LP            x = (3, 1.5)      z = 21      fractional x2 -> branch
 |
 +-- N1: x2 <= 1       x = (10/3, 1)     z = 20.667  bound floor(20.667) = 20; fractional x1 -> branch
 |    |
 |    +-- N3: x1 <= 3  x = (3, 1)        z = 19      integral -> incumbent 19
 |    +-- N4: x1 >= 4  x = (4, 0)        z = 20      integral -> incumbent 20
 |
 +-- N2: x2 >= 2       x = (2, 2)        z = 18      18 < incumbent 20 -> pruned by bound

Implementation

Three programs, all of which were run to produce the numbers quoted in this page. The first solves the constrained problems, the LP pair, water-filling, and the SVM dual, with scipy, so every closed-form answer derived above has a machine check next to it.

import numpy as np
from scipy.optimize import linprog, minimize

# ---- LP primal and dual (HiGHS) -------------------------------------
# (P) max 5x1 + 4x2   s.t. 6x1 + 4x2 <= 24,  x1 + 2x2 <= 6,  x >= 0
# linprog minimizes, so negate c.
p = linprog(c=[-5, -4], A_ub=[[6, 4], [1, 2]], b_ub=[24, 6],
            bounds=[(0, None)] * 2, method="highs")
print(p.x, -p.fun)                 # [3.  1.5] 21.0
print(p.ineqlin.marginals)         # [-0.75 -0.5 ]  (dual prices, sign-flipped)

# (D) min 24y1 + 6y2  s.t. 6y1 + y2 >= 5,  4y1 + 2y2 >= 4,  y >= 0
d = linprog(c=[24, 6], A_ub=[[-6, -1], [-4, -2]], b_ub=[-5, -4],
            bounds=[(0, None)] * 2, method="highs")
print(d.x, d.fun)                  # [0.75 0.5 ] 21.0   -- zero duality gap

# ---- water-filling: max sum log(alpha_i + x_i), sum x = 1, x >= 0 ----
alpha = np.array([0.5, 1.0, 2.0])                       # noise floors, (3,)
res = minimize(lambda x: -np.sum(np.log(alpha + x)),
               x0=np.full(3, 1 / 3),
               constraints=[{"type": "eq", "fun": lambda x: x.sum() - 1}],
               bounds=[(0, None)] * 3, method="SLSQP", tol=1e-12)
print(res.x)                       # [0.75 0.25 0.  ]  water level w = 1.25
print(-res.fun)                    # 1.1394342832 = 2 log 1.25 + log 2

# ---- hard-margin SVM dual as a QP -----------------------------------
X = np.array([[1., 1.], [-1., -1.], [2., 2.]])          # (3, 2)
y = np.array([1., -1., 1.])                             # (3,)
K = (y[:, None] * y[None, :]) * (X @ X.T)               # (3, 3) y_i y_j x_i.x_j
res = minimize(lambda a: -(a.sum() - 0.5 * a @ K @ a),
               x0=np.full(3, 0.1),
               constraints=[{"type": "eq", "fun": lambda a: a @ y}],
               bounds=[(0, None)] * 3, method="SLSQP", tol=1e-14)
alpha_sv = res.x                   # [0.25 0.25 0.  ]  x3 is not a support vector
w = (alpha_sv * y) @ X             # [0.5 0.5],  b = 1 - w @ X[0] = 0
print(alpha_sv, w, -res.fun)       # dual value 0.25 = primal 0.5 ||w||^2

The second implements SGD with momentum and AdamW from scratch, in PyTorch and JAX side by side. The PyTorch tab runs its hand-rolled AdamW against torch.optim.AdamW on a least-squares problem for 200 steps and compares every iterate. The maximum absolute parameter difference across the whole trajectory came out to \( 1.49 \times 10^{-7} \), one float32 ulp at this scale, meaning the two are the same algorithm to machine precision, including the order of operations (decay is applied to the weights before the adaptive step, matching the PyTorch source). The JAX tab is the identical update algebra under lax.scan.

import torch

torch.manual_seed(0)
A = torch.randn(32, 8)                       # (n, d) design
b = torch.randn(32)                          # (n,) targets
lr, b1, b2, eps, wd = 1e-2, 0.9, 0.999, 1e-8, 0.01

def grad(w):                                 # d/dw mean((Aw - b)^2), (d,)
    return 2.0 / A.shape[0] * A.T @ (A @ w - b)

# --- SGD with (heavy-ball) momentum, for reference -------------------
w, buf = torch.zeros(8), torch.zeros(8)
for t in range(200):
    g = grad(w)
    buf = 0.9 * buf + g                      # PyTorch convention: v = mu*v + g
    w = w - lr * buf

# --- AdamW from scratch ----------------------------------------------
w = torch.zeros(8)
m, v = torch.zeros(8), torch.zeros(8)
for t in range(1, 201):
    g = grad(w)
    w = w - lr * wd * w                      # decoupled decay, BEFORE the step
    m = b1 * m + (1 - b1) * g                # EMA of gradient
    v = b2 * v + (1 - b2) * g * g            # EMA of squared gradient
    mhat = m / (1 - b1 ** t)                 # bias corrections
    vhat = v / (1 - b2 ** t)
    w = w - lr * mhat / (vhat.sqrt() + eps)

# --- verify against torch.optim.AdamW step by step -------------------
w_ref = torch.zeros(8, requires_grad=True)
opt = torch.optim.AdamW([w_ref], lr=lr, betas=(b1, b2), eps=eps,
                        weight_decay=wd)
w2, m2, v2, max_diff = torch.zeros(8), torch.zeros(8), torch.zeros(8), 0.0
for t in range(1, 201):
    opt.zero_grad()
    ((A @ w_ref - b) ** 2).mean().backward()
    opt.step()
    g = grad(w2)
    w2 = w2 - lr * wd * w2
    m2 = b1 * m2 + (1 - b1) * g
    v2 = b2 * v2 + (1 - b2) * g * g
    w2 = w2 - lr * (m2 / (1 - b1 ** t)) / ((v2 / (1 - b2 ** t)).sqrt() + eps)
    max_diff = max(max_diff, (w2 - w_ref.detach()).abs().max().item())
print(max_diff)   # 1.4901161e-07: trajectories match to a float32 ulp

import jax
import jax.numpy as jnp

kA, kb = jax.random.split(jax.random.PRNGKey(0))
A = jax.random.normal(kA, (32, 8))           # (n, d)
b = jax.random.normal(kb, (32,))             # (n,)
lr, b1, b2, eps, wd = 1e-2, 0.9, 0.999, 1e-8, 0.01

def loss(w):                                 # scalar
    r = A @ w - b
    return jnp.mean(r * r)

grad = jax.jit(jax.grad(loss))

# --- SGD with momentum as a scanned pure function --------------------
def sgd_step(carry, _):
    w, buf = carry
    buf = 0.9 * buf + grad(w)
    w = w - lr * buf
    return (w, buf), loss(w)

(w, _), losses = jax.lax.scan(sgd_step, (jnp.zeros(8), jnp.zeros(8)),
                              None, length=200)

# --- AdamW: same algebra as the PyTorch tab, same operation order ----
def adamw_step(carry, t):
    w, m, v = carry
    g = grad(w)
    w = w - lr * wd * w                      # decoupled decay first
    m = b1 * m + (1 - b1) * g
    v = b2 * v + (1 - b2) * g * g
    mhat = m / (1 - b1 ** t)                 # t is 1-indexed float
    vhat = v / (1 - b2 ** t)
    w = w - lr * mhat / (jnp.sqrt(vhat) + eps)
    return (w, m, v), loss(w)

init = (jnp.zeros(8), jnp.zeros(8), jnp.zeros(8))
(w, m, v), losses = jax.lax.scan(adamw_step, init,
                                 jnp.arange(1, 201, dtype=jnp.float32))
print(losses[-1])   # 0.9346 on this random instance; monotone after warm start
# optax.adamw implements the same update; writing it out once is the point.

The third generates the gradient-descent table of Problem 3 and the ISTA/FISTA comparison quoted in the proximal section, with the \( t_k \) recurrence exactly as derived.

import numpy as np

# ---- GD on f = 1/2 x'Qx, Q = diag(2, 20): the Problem 3 table --------
Q = np.diag([2.0, 20.0])
mu, L = 2.0, 20.0
x = np.array([1.0, 1.0])
eta = 2 / (mu + L)                            # = 1/11, optimal for distances
for k in range(11):
    print(k, x, 0.5 * x @ Q @ x)              # contracts by 9/11 per step
    x = x - eta * (Q @ x)

# ---- ISTA and FISTA on a small lasso ---------------------------------
rng = np.random.default_rng(0)
A = rng.standard_normal((20, 5))              # (n, d)
x_true = np.array([1.5, 0.0, -2.0, 0.0, 0.0])
b = A @ x_true + 0.01 * rng.standard_normal(20)
lam = 1.0
L = np.linalg.eigvalsh(A.T @ A).max()         # 38.188 = Lipschitz const of grad f

def soft(v, t):                               # prox of t*||.||_1, elementwise
    return np.sign(v) * np.maximum(np.abs(v) - t, 0.0)

def F(x):                                     # full composite objective
    return 0.5 * np.sum((A @ x - b) ** 2) + lam * np.abs(x).sum()

x = np.zeros(5)                               # ISTA
for k in range(200):
    x = soft(x - (1 / L) * A.T @ (A @ x - b), lam / L)
print(F(x), x)      # 3.4323500, support {0, 2} exact zeros elsewhere

x = np.zeros(5); y = x.copy(); t = 1.0        # FISTA
for k in range(200):
    x_new = soft(y - (1 / L) * A.T @ (A @ y - b), lam / L)
    t_new = (1 + np.sqrt(1 + 4 * t * t)) / 2  # Nesterov t_k recurrence
    y = x_new + ((t - 1) / t_new) * (x_new - x)
    x, t = x_new, t_new
print(F(x))         # 3.4323500; hits 3.43920 by k=10 vs ISTA's 3.69401

How it is done in practice

Convex solvers: modeling layers and the engines underneath

Almost nobody hand-derives KKT systems in production. The modeling layer CVXPY applies the convexity calculus from the first section mechanically, disciplined convex programming, to verify a problem is convex by construction and compile it to a cone program handed to a backend, Clarabel or SCS for general cones, OSQP for QPs, HiGHS for LPs and MIPs. The compilation is the same reduction this page performed on the SVM by hand. Introduce slack variables, express everything as linear maps plus cone membership, dualize. Below the modeling layer, the engineering picture has been stable for a decade. Simplex and interior-point implementations both survive because they fail differently. Interior-point wins large sparse LPs from cold starts (predictable \( O(\sqrt m) \) iterations, heavy per-iteration linear algebra amenable to sparse Cholesky). Simplex wins the re-solve workload inside branch and bound, where each node's LP differs from its parent by one bound and a dual-feasible warm basis finishes in a few pivots. Every serious MIP solver therefore runs interior-point (or nowadays a first-order method) once at the root and dual simplex everywhere else in the tree. For the largest LPs of the last few years, Google's PDLP, a restarted primal-dual hybrid gradient method that needs only matrix-vector products, solved instances with over \( 10^{10} \) nonzeros that do not fit any factorization in memory, trading iteration count for bandwidth, an explicit return of first-order methods to a domain interior-point had owned since Karmarkar.

Integer programming at industrial scale

Production MIP is branch and cut plus twenty years of engineering, presolve that removes a large fraction of variables before the root LP, cut families (Gomory mixed-integer, knapsack covers, flow covers) generated and aged out adaptively, primal heuristics (feasibility pump, RINS) that supply incumbents so pruning can fire early, and portfolio branching rules learned from the tree so far. Airlines re-solve crew schedules with millions of binaries overnight. Kidney-exchange programs clear national pools with IP solvers under legal audit requirements, the dual bound serving as the proof of optimality that policy demands. The open-source tier, HiGHS and OR-Tools CP-SAT, now wins independent benchmarks regularly and is what scipy.optimize.milp ships. The craft knowledge that matters day to day is to tighten formulations before buying hardware, because the integrality gap, not the node throughput, decides tractability. Big-M constraints with loose M are the classic self-inflicted wound, exactly because they make the LP relaxation nearly meaningless and every prune late.

Training-scale optimization: what survives contact with a trillion tokens

Inside large training runs the update rule is AdamW essentially everywhere, with the engineering questions shifted one level up, to state memory (two fp32 EMAs double the parameter memory, motivating 8-bit optimizer states and sharding of \( m \) and \( v \) across data-parallel ranks), numerics (bf16 gradients accumulated in fp32, and \( \epsilon \) placement inside versus outside the square root differing across codebases and mattering at small \( v \)), and the schedule and warmup questions covered on the deep learning engineering page. The theory in this page maps onto that practice with one honest caveat. The loss is nonconvex, so what the proofs actually license is stationarity at rate \( O(1/\varepsilon^2) \) for SGD under smoothness, not the convex rates. The convex analysis earns its keep as the design language, every piece of AdamW traces to a theorem about a convex surrogate, and as the exact theory of the last layer, of linear probes, and of the deterministic full-batch limits used to debug. When a training run misbehaves, the differential diagnosis is precisely the terms of the SGD bound. Is it the \( R^2/(2\eta k) \) term (learning rate too small, warmup too long), the \( \eta G^2/2 \) term (noise ball, batch too small or learning rate too high), or a violated smoothness assumption (loss spikes, gradient clipping engaged)?

The current research frontier

Optimizers beyond AdamW

The empirical frontier is a contest over whether structure-aware preconditioning beats coordinate-wise scaling at scale. Shampoo (Gupta, Koren, and Singer, 2018) preconditions each weight matrix by Kronecker-factored second-moment roots, a practical descendant of full-matrix AdaGrad, and a distributed implementation won the external tuning track of the AlgoPerf benchmark run by MLCommons. Muon (Jordan et al., 2024) orthogonalizes each layer's momentum matrix with a few Newton-Schulz iterations, spiritually a steepest-descent step under the spectral norm rather than the Euclidean one. It set nanoGPT speedrun records, and Moonshot's Kimi team reported it scaling to multi-billion-parameter pretraining with about half the training FLOPs of a tuned AdamW baseline (Liu et al., 2025). Lion (Chen et al., 2023), found by symbolic program search at Google, keeps only a momentum sign, halving optimizer memory. Sophia (Liu et al., 2023) uses a cheap diagonal Hessian estimate with clipping. Meta's schedule-free methods (Defazio et al., 2024) remove the decay schedule by averaging, won AlgoPerf's self-tuning track, and connect back to Polyak-Ruppert averaging from the classical stochastic approximation literature. The sober reading across independent evaluations is that coordinate-wise methods are hard to beat by more than tens of percent, the gains are real but contested, and every few years the community re-learns that fair baselines (tuned AdamW) erase many claimed victories.

Theory catching up to practice

On the analysis side, the assumptions are being renegotiated to match observed training. The \( (L_0, L_1) \)-smoothness model of Zhang, He, Sra, and Jadbabaie (2020), where the local smoothness constant grows with the gradient norm, was introduced precisely because clipped SGD provably beats plain SGD under it, matching what every LSTM and transformer practitioner already did. Work spanning several groups has pinned down conditions for Adam's nonconvex convergence and the role of \( \beta_2 \), the edge-of-stability phenomenon (Cohen et al., 2021) documented full-batch training living where classical descent analysis forbids, with the sharpness \( \lambda_{\max}(\nabla^2 f) \) hovering at \( 2/\eta \), and a parallel line studies implicit bias, meaning which minimizer, among the many, the optimizer selects. On the discrete side, machine learning is entering the solvers themselves. Learned branching policies (Gasse et al., 2019, imitation of strong branching with GNNs), learned cut selection, and neural primal heuristics now appear inside commercial solver releases, while first-order LP (PDLP and successors, including GPU implementations) is redrawing the boundary of what counts as a solvable linear program. Even interpretability research has looped back to this material. The sparse autoencoders used by Anthropic and Google DeepMind to decompose model activations are descendants of the sparse-coding objective, an \( \ell_1 \)-penalized reconstruction of exactly the form the soft-thresholding derivation above solves.

Open source to read

• cvxpy/cvxpy implements disciplined convex programming. Read cvxpy/atoms/ to see the convexity calculus of this page encoded as a ruleset, each atom declaring curvature and monotonicity.
• scipy/scipy offers scipy/optimize/_linprog_highs.py for the LP interface used here, and _optimize.py for the BFGS and L-BFGS-B drivers with Wolfe line search.
• ERGO-Code/HiGHS is the strongest open LP/MIP engine. Start with src/mip/HighsMipSolver.cpp to watch branch and cut orchestrated for real.
• google/or-tools contains CP-SAT, a SAT-based integer solver that routinely beats MIP on scheduling. The examples/python/ directory is the fastest education in modeling.
• In pytorch/pytorch, torch/optim/adamw.py is the reference the implementation above matches digit for digit. Note the single-tensor versus fused code paths and the decay-before-step order.
• google-deepmind/optax treats optimizers as composable gradient transformations. optax/_src/transform.py shows AdamW factored into scale_by_adam plus add_decayed_weights, the decoupling argument as an API.
• google/jaxopt provides differentiable solvers, implicit differentiation through KKT systems, the modern use of the KKT conditions as a layer.
• casadi/casadi is the standard tool for optimal control and MPC. It generates the sparse KKT systems interior-point methods consume, wired to Ipopt.
• coin-or/Ipopt is the canonical open interior-point NLP solver. src/Algorithm/IpIpoptAlg.cpp is the barrier loop of this page with a filter line search around it.
• osqp/osqp is an ADMM-based QP solver small enough to read in one sitting and deployed on embedded controllers.

Common misconceptions

• "Convex means easy, nonconvex means hopeless." Both halves are wrong. Convex problems can be practically hard (semidefinite programs with \( 10^6 \) constraints are convex and routinely unsolvable), and structured nonconvex problems are solved to global optimality daily by branch and bound, or provably locally at scale under PL-type conditions. Convexity buys certificates and composable theory, not a small running time by itself.

• "The dual is just a trick for deriving algorithms." The dual is a certificate machine. Any feasible dual point proves a bound on the primal that a third party can check without re-solving. Branch-and-bound prunes, auction and market prices, and the SVM's support-vector sparsity are all the same theorem wearing different clothes.

• "Gradient descent with a smaller learning rate always converges more reliably." Below \( 2/L \) smaller steps only slow the same convergence, and in the stochastic setting the step size also sets the noise-ball radius. Too small means the \( R^2/(2\eta k) \) term dominates and progress stalls for reasons that look like, but are not, a bad model. The rate-optimal choices balance the two terms and are not minimal.

• "Momentum works because it averages out mini-batch noise." Momentum's provable role is geometric, curing curvature anisotropy. It improves the deterministic \( \kappa \) dependence to \( \sqrt\kappa \), and Nesterov's variant is optimal for exactly that oracle. Against pure noise with constant variance, no first-order method beats \( O(1/\sqrt k) \), momentum included. Its stochastic benefits are real but second-order (implicit larger effective step, variance shaping), not a variance-reduction guarantee.

• "Adam converges, and the proof is in the paper." The original regret proof was wrong, and Reddi, Kale, and Kumar exhibited a convex counterexample on which Adam provably fails. What exists are corrected variants (AMSGrad) and nonconvex stationarity results under side conditions. Adam's dominance is an empirical fact with partial theory, and repeating the original claim signals not having read the 2018 paper.

• "Newton-type methods are obsolete in the deep learning era." The opposite is true in two directions. Interior-point methods, which are Newton's method on the central path, solve the convex programs running power grids and MPC controllers, and the most credible AdamW challengers (Shampoo, Muon, Sophia) are all curvature-aware methods rationing second-order information under a memory budget. What died was dense BFGS at \( n^2 \) memory, and the L-BFGS arithmetic above shows exactly where.

• "The LP relaxation is a heuristic without guarantees." The relaxation is a rigorous bound in one direction, and that single direction is what makes exact MIP possible. Every prune in Problem 6 is a theorem. When the gap is bounded, rounding converts the relaxation into an approximation algorithm with a proof. What is heuristic is naive rounding without the gap analysis, which Problem 6 shows losing 5% on a two-variable problem.

• "Greedy is a heuristic, and you can always do better with more compute." For monotone submodular maximization, greedy's \( 1 - 1/e \) is optimal for any polynomial algorithm unless P = NP (Feige, 1998). More compute inside the query model provably does not help. Doing better requires changing the problem, extra structure like curvature bounds, or exponential time. Some "heuristics" are optimal algorithms wearing modest clothing.

Self-check

References

1. Boyd, Vandenberghe. Convex Optimization. Cambridge University Press, 2004. web.stanford.edu/~boyd/cvxbook
2. Nocedal, Wright. Numerical Optimization, 2nd ed. Springer, 2006.
3. Nesterov. Lectures on Convex Optimization, 2nd ed. Springer, 2018.
4. Nemirovski, Yudin. Problem Complexity and Method Efficiency in Optimization. Wiley, 1983.
5. Dantzig. Linear Programming and Extensions. Princeton University Press, 1963.
6. Nesterov, Nemirovskii. Interior-Point Polynomial Algorithms in Convex Programming. SIAM, 1994.
7. Kuhn, Tucker. Nonlinear programming. Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, 1951. (Karush's 1939 thesis contains the conditions independently.)
8. Robbins, Monro. A stochastic approximation method. Annals of Mathematical Statistics 22(3), 1951. doi:10.1214/aoms/1177729586
9. Polyak. Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics 4(5), 1964.
10. Broyden, Fletcher, Goldfarb, and Shanno. The BFGS update, published independently in four 1970 papers (Journal of the Institute of Mathematics and Its Applications 6, The Computer Journal 13, and Mathematics of Computation 24, twice).
11. Nemhauser, Wolsey, Fisher. An analysis of approximations for maximizing submodular set functions I. Mathematical Programming 14, 1978. doi:10.1007/BF01588971
12. Grötschel, Lovász, Schrijver. The ellipsoid method and its consequences in combinatorial optimization. Combinatorica 1, 1981.
13. Nesterov. A method of solving a convex programming problem with convergence rate O(1/k²). Soviet Mathematics Doklady 27, 1983.
14. Karmarkar. A new polynomial-time algorithm for linear programming. Combinatorica 4, 1984. doi:10.1007/BF02579150
15. Liu, Nocedal. On the limited memory BFGS method for large scale optimization. Mathematical Programming 45, 1989.
16. Beck, Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences 2(1), 2009. doi:10.1137/080716542
17. Duchi, Hazan, Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of Machine Learning Research 12, 2011. jmlr.org/papers/v12/duchi11a
18. Johnson, Zhang. Accelerating stochastic gradient descent using predictive variance reduction. NeurIPS 2013.
19. Kingma, Ba. Adam: a method for stochastic optimization. ICLR 2015. arXiv:1412.6980
20. Bottou, Curtis, Nocedal. Optimization methods for large-scale machine learning. SIAM Review 60(2), 2018. arXiv:1606.04838
21. Reddi, Kale, Kumar. On the convergence of Adam and beyond. ICLR 2018. arXiv:1904.09237
22. Gupta, Koren, Singer. Shampoo: preconditioned stochastic tensor optimization. ICML 2018. arXiv:1802.09568
23. Loshchilov, Hutter. Decoupled weight decay regularization. ICLR 2019. arXiv:1711.05101
24. Gasse, Chételat, Ferroni, Charlin, Lodi. Exact combinatorial optimization with graph convolutional neural networks. NeurIPS 2019. arXiv:1906.01629
25. Liu, Su, Yao, et al. Muon is scalable for LLM training. 2025. arXiv:2502.16982

Key takeaway