Second-Order & Approximate Methods

Every optimizer so far (SGD, momentum, and the adaptive methods) is first-order: it sees only the gradient $g=\nabla \mathcal{L}(\theta )$ and steps $\theta \leftarrow \theta -\eta g$. The gradient is a flat plane tangent to the loss; it gives the downhill direction but nothing about how the surface curves. A second-order method also reads the curvature (the Hessian $H={\nabla }^2\mathcal{L}$) and uses it to choose both direction and step length at once. On the ill-conditioned ravines that make first-order descent crawl, that extra information is decisive. The problem is its cost, and most of this lesson is about approximations that capture a fraction of it cheaply.¹

Newton's method

Expand the loss to second order around the current iterate $\theta$. Writing the step as $\Delta ={\theta }^{\prime }-\theta$, Taylor's theorem gives the local quadratic model

This is a paraboloid in $\Delta$. If $H\succ 0$ it has a unique minimum, found by setting the gradient of the right-hand side to zero:

The minimizer of the quadratic is one step away, and that step is the Newton step. Taking it as the update gives Newton's method.

The model being minimized is a bowl, not a plane. First-order descent fits the loss with its tangent plane $\mathcal{L}(\theta )+g^{\top }\Delta$ and can only report a direction of decrease; the plane has no bottom, so the step length has to come from a separately tuned learning rate. The second-order model adds the quadratic term $\frac{1}{2}{\Delta }^{\top }H\Delta$, which closes the surface into a paraboloid with a definite minimum. Minimizing that bowl fixes direction and length together, and $-H^{-1}g$ is the vector from the current point to its floor.

The geometry explains the advantage. Gradient descent steps along $-g$, perpendicular to the contours of $\mathcal{L}$ as drawn in raw coordinates; on a stretched bowl that direction does not point at the minimum. Newton multiplies by $H^{-1}$, which rescales each direction by its curvature (stretching the gently curved axes and shrinking the steep ones), so that in the rescaled coordinates the bowl is round and $-H^{-1}g$ aims straight at the bottom.

Because the rescaling cancels the eigenvalue spread of $H$, Newton's method is invariant to the condition number $\kappa (H)={\lambda }_{\max }/{\lambda }_{\min }$: a bowl stretched a thousand-to-one is solved as fast as a round one. This removes exactly the problem that momentum and adaptive rates only mitigate.² To see the one-step claim algebraically, take an exact quadratic $\mathcal{L}(\theta )=\frac{1}{2}(\theta -{\theta }^{\star }{)}^{\top }H(\theta -{\theta }^{\star })$, for which $g=H(\theta -{\theta }^{\star })$:

The invariance runs deeper than the condition number. Newton's method is affine invariant: reparametrize by any invertible linear map $\varphi =A\theta$, and the iterates it produces are the exact images under $A$ of the iterates on the original problem. A quick check: under $\varphi =A\theta$ the gradient transforms as $g_{\varphi }=A^{-\top }g$ and the Hessian as $H_{\varphi }=A^{-\top }H{A}^{-1}$, so the Newton direction is

the original step carried through $A$. Gradient descent has no such property — the direction $-g_{\varphi }=-A^{-\top }g$ does not equal $A(-g)$, so a bad choice of units alone can make first-order descent crawl. Newton reads geometry, not coordinates, which is why the same step handles a bowl stretched a thousand-to-one and a round one alike.

The cost is high. Forming $H$ needs the second derivative with respect to every pair of parameters — $\left(\genfrac{}{}{0ex}{}{d+1}{2}\right)\approx d^2/2$ distinct entries — and solving the linear system $H\Delta =-g$ (whether by inverting $H$ or by Cholesky factorization) costs $O(d^3)$ arithmetic and $O(d^2)$ storage. For a $d=10^7$ network, $H$ alone is $10^{14}$ numbers; the factorization is $10^{21}$ operations, per step. The entire remainder of the lesson is a catalog of ways to keep the curvature information while retreating from that $O(d^2)$/$O(d^3)$ cost.

Why second-order methods are rare in deep learning

Newton's method works well on small, smooth, deterministic problems. In deep learning all three properties fail at once, and each failure alone rules the method out.

The Hessian is $d\times d$

The first wall is sheer size. The Hessian of a network with $d$ parameters is a $d\times d$ symmetric matrix (one entry per pair of parameters), so its storage grows as $d^2$ and forming its inverse costs $O(d^3)$ by Gaussian elimination.

The numbers are far from feasible. The table quantifies the gap: a gradient that fits in memory pairs with a Hessian that cannot be stored on any hardware.

Even storing $H$ once for a modest $10^6$-parameter network needs a terabyte at single precision; inverting it is hopeless. And the optimizer must do this every step.

Newton is attracted to saddle points

The second wall is qualitative, and worse: where it does run, Newton's method walks toward the very saddle points that dominate high-dimensional loss surfaces. Diagonalize $H=Q\Lambda Q^{\top }$ and resolve the Newton step into the eigenbasis. Along eigendirection $i$, with eigenvalue ${\lambda }_i$ and gradient component $g_i$, the step is

Here is the defect. The first-order step $-g_i$ always moves downhill along axis $i$. The Newton step divides by ${\lambda }_i$, and dividing by a negative eigenvalue flips the sign — turning a downhill move into an uphill one, straight toward the critical point. Newton does not distinguish whether a critical point is a minimum, a maximum, or a saddle; it converges to all of them, because it solves $g=0$ rather than minimizing $\mathcal{L}$.³

Since high-loss critical points are overwhelmingly saddles, an unmodified Newton step is not merely useless but actively harmful. The fix is to force the curvature positive: saddle-free Newton replaces $H^{-1}$ with $|H{|}^{-1}$ (take the absolute values of the eigenvalues), and classic damping (Levenberg–Marquardt) solves $(H+\mu I{)}^{-1}g$, adding $\mu I$ until the matrix is positive definite. Both restore descent, but both still need $H$.

Minibatch noise corrupts the estimate

The third wall is statistical. Deep nets are trained on minibatches, so both $g$ and $H$ are noisy estimates from a sample. First-order steps tolerate this (noise averages out over many steps), but the Newton step divides by the Hessian estimate, and a ratio amplifies error in the denominator. A small eigenvalue ${\lambda }_i$ estimated with even modest relative error produces a huge swing in $-g_i/{\lambda }_i$. The very curvature information Newton depends on is the part the noise corrupts most.

The rest of the lesson is the response: methods that capture some curvature without ever forming, storing, or inverting $H$.

Conjugate gradients

Steepest descent on a quadratic wastes effort because successive steps are orthogonal, and on a stretched bowl orthogonal steps zig-zag, each step partly undoing the progress of the last along the previous direction. Conjugate gradients (CG) fixes this by choosing search directions that are conjugate with respect to $H$, so that minimizing along a new direction never disturbs the minimization already achieved along the old ones.

The new direction reuses the previous one, bent just enough to be conjugate to it. The combining coefficient ${\beta }_t$ (the Fletcher–Reeves form ${\beta }_t=g_t^{\top }{g}_t/g_{t-1}^{\top }{g}_{t-1}$) is computed from gradients alone, so CG never touches $H$ explicitly. Each step needs one gradient and one line search.

On a strictly quadratic objective CG terminates in at most $d$ steps; on a general nonlinear loss it is run with periodic restarts and remains a strong batch optimizer, though minibatch noise in the line search limits its use in pure deep learning.⁴

Quasi-Newton: BFGS and L-BFGS

Newton needs $H^{-1}$ but the cost is prohibitive. Quasi-Newton methods learn an approximation $B\approx H^{-1}$ as they go, reading curvature off the change in the gradient. The secant condition anchors the construction: over a step, the gradient's change encodes the curvature along that step. Writing

a first-order expansion of $g$ gives $y_t\approx H{s}_t$, so $H^{-1}$ must satisfy $H^{-1}{y}_t\approx s_t$. BFGS updates $B_t$ at each step so it reproduces this relation on the latest $(s_t,y_t)$ pair while changing $B$ as little as possible (a rank-two correction):

The step is then ${\theta }_{t+1}={\theta }_t-{\alpha }_t{B}_t{g}_t$, a Newton-like move that never forms or inverts $H$, only matrix–vector products and a line search.

The rank-two form is engineered to satisfy three demands at once: it reproduces the secant relation $B_{t+1}{y}_t=s_t$, it keeps $B_{t+1}$ symmetric and positive definite whenever $y_t^{\top }{s}_t>0$ (the curvature condition, guaranteed by a proper line search), and among all matrices meeting the secant relation it is the closest to $B_t$ in a weighted Frobenius norm. Positive definiteness matters: it forces every step to be a descent direction, so unlike raw Newton, BFGS cannot sign-flip toward a saddle. The curvature is learned, one $(s_t,y_t)$ pair at a time, and it accumulates the effect of many past pairs into a single matrix.

The remaining problem is memory: $B$ is still $d\times d$. L-BFGS (limited-memory BFGS) never stores $B$ at all. It keeps only the last $m$ pairs $(s_t,y_t)$ (typically $m=5$ to $20$) and reconstructs the action $B_t{g}_t$ on the fly by a two-loop recursion over those vectors: a backward sweep peels off the recent curvature corrections into scalars ${\alpha }_i$, a scaling by $\gamma =(s_{t-1}^{\top }{y}_{t-1})/(y_{t-1}^{\top }{y}_{t-1})$ approximates the bulk of $H^{-1}$, and a forward sweep folds the corrections back in. Each sweep is $2m$ inner products and vector updates, so one matrix–vector action $B_t{g}_t$ costs $O(md)$ rather than $O(d^2)$. Memory drops from $O(d^2)$ to $O(md)$, the same order as the gradient. That is what makes a quasi-Newton method usable at deep-learning scale (in the full-batch regime).⁵

The reason L-BFGS is still uncommon on deep nets is the same noise wall Newton hit, now in a subtler place. The secant pair $(s_t,y_t)$ is built from the difference $y_t=g_{t+1}-g_t$ of two gradients. Under minibatching those gradients are computed on different samples, so $y_t$ mixes the true curvature signal $H{s}_t$ with the difference of two independent noise terms. That noise does not average away inside the BFGS update the way it does across SGD steps — it corrupts the very curvature estimate the method accumulates, and a stale or noisy $B$ can point the step in a wrong direction. L-BFGS works best with large or full batches, where the gradients are clean enough; on noisy minibatch gradients the first-order methods are preferable.

Gauss-Newton and the Fisher matrix

Before the natural gradient there is a simpler positive-semidefinite stand-in for $H$. For a loss that is a sum of squared residuals $\mathcal{L}=\frac{1}{2}{\sum }_i{r}_i(\theta {)}^2$, write the Jacobian of the residual vector as $J=\partial r/\partial \theta$. The exact Hessian is

The Gauss-Newton approximation drops the second term and keeps only $G=J^{\top }J$. That term is the part carrying the residual curvature ${\nabla }^2{r}_i$; near a good fit the residuals $r_i$ are small, so dropping it costs little, and what remains is an outer product $J^{\top }J⪰0$ that can never be indefinite. Gauss-Newton gets the same immunity to saddles that damping enforces by hand, and it needs only first derivatives of the residuals.

The Fisher information matrix generalizes this to probabilistic models. When the loss is the negative log-likelihood of the model's output distribution $p_{\theta }(y\mid x)$, the Fisher is the expected outer product of the score,

Two facts make $F$ attractive. It equals the expected Hessian of the negative log-likelihood (so it is curvature, not a crude substitute), and being an expectation of outer products it is positive semidefinite by construction — no sign-flips toward saddles. For the exponential-family output layers used in practice, $F$ coincides with the Gauss-Newton matrix $G$, tying the two views together.

Natural gradient, K-FAC, and Hessian-free

The natural gradient preconditions by $F$ instead of $H$, measuring distance in the space of the model's output distributions rather than its raw parameters:

The motivation is invariance. Ordinary gradient descent moves fastest in the direction of steepest descent measured by Euclidean distance in parameter space, but a parameter is an arbitrary label — reparametrize the network and the same distance means something different. The natural gradient instead takes the steepest-descent direction under the Kullback–Leibler divergence between output distributions, a quantity that does not depend on how the parameters are named. To second order that divergence is $\mathrm{KL}(p_{\theta }\Vert p_{\theta +\Delta })\approx \frac{1}{2}{\Delta }^{\top }F\Delta$, so $F$ plays exactly the role $H$ played in the Taylor bowl, and $F^{-1}g$ is the step that descends the loss per unit of distributional change. Because $F$ is positive semidefinite, that step never climbs toward a saddle the way $-H^{-1}g$ can.

$F$ is still $d\times d$, so like $H$ it is used only through approximations.

For a layer with $n_{\text{in}}$ inputs and $n_{\text{out}}$ outputs, its weight has $d_{\ell }=n_{\text{in}}{n}_{\text{out}}$ parameters and its Fisher block is $d_{\ell }\times d_{\ell }$. K-FAC models that block as $A\otimes G$, where $A$ is the $n_{\text{in}}\times n_{\text{in}}$ covariance of the layer's inputs and $G$ is the $n_{\text{out}}\times n_{\text{out}}$ covariance of its output gradients. Storing and inverting two factors of size $n_{\text{in}}$ and $n_{\text{out}}$ costs $O(n_{\text{in}}^3+n_{\text{out}}^3)$ instead of the $O(d_{\ell }^3)=O(n_{\text{in}}^3{n}_{\text{out}}^3)$ of the full block — a difference of many orders of magnitude for a wide layer. The Kronecker structure is an assumption (it holds exactly only if inputs and output-gradients are statistically independent), but it captures between-parameter curvature that a diagonal preconditioner throws away.

Hessian-vector products without forming $H$

The trick under Hessian-free optimization is that the whole $d\times d$ matrix is never needed — only its action $Hv$ on chosen vectors, which inner CG consumes one at a time. And $Hv$ is itself a gradient. For any fixed $v$,

because differentiating $g^{\top }v$ (a scalar, with $v$ held constant) with respect to $\theta$ contracts the second derivative ${\nabla }_{\theta }g=H$ against $v$. Pearlmutter's trick computes this at the cost of one extra backward pass: run the ordinary backward pass to get $g$, take the scalar $g^{\top }v$, and backpropagate that. No entry of $H$ is ever materialized, storage stays $O(d)$, and each CG iteration inside the inner solve costs one forward–backward pair.

Curvature is accessible — just never as a stored matrix.⁶

A ladder of approximations

Every method here occupies one rung of a single ladder that trades curvature accuracy against cost per step. At the bottom is plain gradient descent, which uses no curvature and costs $O(d)$. Each rung up costs more and captures more of $H$: a diagonal estimate (one number per parameter, which is what Adam does), then a per-layer Kronecker block (K-FAC), then a low-rank running estimate from curvature pairs (L-BFGS), and at the top the full Hessian solve of exact Newton at $O(d^3)$. Deep learning almost always sits on the lowest two rungs, because the higher ones cannot tolerate minibatch noise and their per-step cost is not worthwhile when the gradient itself is only an estimate.

First-order versus second-order

The whole chapter resolves into one comparison. First-order methods are cheap, noise- tolerant, and scale to billions of parameters but suffer under ill-conditioning; second-order methods fix conditioning and converge in far fewer steps but cost too much per step and break under saddles and noise. Deep learning sits in the top row, borrowing pieces of the bottom row (momentum, adaptive rates, occasionally L-BFGS) when it can.

Modern second-order optimizers

Goodfellow's §8.6 stops at the classical toolkit; the past decade pushed approximate second-order optimization back toward large-scale training.

• K-FAC in full. Martens and Grosse's paper is the source of the Kronecker-factored Fisher approximation above, and the first to make natural gradient practical on real networks.⁷
• Shampoo. Gupta, Koren, and Singer keep a preconditioner per tensor dimension (Kronecker-factored full-matrix AdaGrad); a distributed version has trained production language models faster than Adam.⁸
• Sophia. Liu et al. use a clipped diagonal Hessian estimate to stay saddle-safe and roughly halve language-model pre-training steps — the ladder's diagonal rung, but reading true curvature rather than Adam's gradient second moment.⁹

The trend runs up the ladder: as models grow, a little real curvature increasingly becomes worth its cost, and the saddle and noise walls are handled by the same positive-definite surrogates and clipping this lesson derived.

Takeaways

• Newton's method $\theta \leftarrow \theta -H^{-1}g$ jumps to the minimum of the local quadratic in one step and is invariant to the condition number — it rescales each direction by its curvature, curing the ill-conditioning that makes first-order descent zig-zag.
• Three walls keep it out of deep learning: the Hessian is $d\times d$ (storage $\propto d^2$, inversion $\propto d^3$, infeasible for $d$ up to billions); Newton is attracted to saddles because dividing by a negative eigenvalue flips a downhill step uphill; and minibatch noise in the denominator is amplified by the division.
• Conjugate gradients steps in $H$-conjugate directions to eliminate the zig-zag of steepest descent and reach a quadratic's minimum in $\le d$ steps using gradients only — no Hessian inverse.
• Quasi-Newton methods build $B\approx H^{-1}$ from gradient differences via the secant condition $y_t\approx H{s}_t$; BFGS stores the full $O(d^2)$ matrix, L-BFGS keeps only the last $m$ curvature pairs at $O(md)$ memory.
• The natural gradient preconditions by the always-positive Fisher matrix $F$; K-FAC factors $F$ per layer as a Kronecker product; Hessian-free optimization solves $H\Delta =-g$ by inner CG using $Hv$ products from one extra backward pass — curvature without ever storing $H$.
• Deep learning stays first-order: Adam is the diagonal, $O(d)$ approximation to second-order preconditioning, keeping the conditioning benefit at first-order cost.