---
title: Second-Order & Approximate Methods
module: Optimization
moduleNumber: 3
lessonNumber: 5
order: 305
summary: >
  Newton's method reads the curvature of the loss off its Hessian and jumps to the
  minimum of the local quadratic in a single step, rescaling away the
  ill-conditioning that slows first-order descent. We derive it, then explain the
  three obstacles that keep it out of deep learning: a $d \times d$ Hessian for $d$ in
  the billions, an attraction to saddle points, and minibatch noise. The alternative is
  approximation (conjugate gradients, BFGS and L-BFGS, the natural gradient and
  Hessian-free methods), each buying some of Newton's curvature information without
  ever forming or inverting $H$.
topics: [Optimization]
sources:
  - book: Goodfellow
    ref: "§8.6 — Approximate Second-Order Methods"
  - book: Goodfellow
    ref: "§4.3.1 — Beyond the Gradient: Jacobian and Hessian Matrices; §8.2 Challenges in Optimization"
  - book: Chollet
    ref: "§2.4 — Gradient-Based Optimization"
---

Every optimizer so far ([SGD](/deep-learning/optimization/gradient-descent-and-sgd),
[momentum, and the adaptive methods](/deep-learning/optimization/momentum-and-adaptive-methods))
is **first-order**: it sees only the gradient $g = \nabla\mathcal{L}(\theta)$ and steps
$\theta \gets \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](/deep-learning/optimization/the-optimization-landscape) 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.[^gf-secondorder]

## Newton's method

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

$$
\mathcal{L}(\theta + \Delta) \;\approx\; \mathcal{L}(\theta) \;+\; g^\top \Delta \;+\;
\tfrac12\, \Delta^\top H\, \Delta,
\qquad g = \nabla\mathcal{L}(\theta),\quad H = \nabla^2\mathcal{L}(\theta).
$$

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:

$$
\nabla_\Delta\!\brackets{ g^\top\Delta + \tfrac12 \Delta^\top H \Delta }
= g + H\Delta = 0
\;\Longrightarrow\;
\Delta^\star = -\,H^{-1} g.
$$

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 $\tfrac12\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.

$$
% caption: First order fits a tangent plane (no bottom, needs a learning rate); second
% order fits a bowl $\mathcal{L}+g^\top\Delta+\tfrac12\Delta^\top H\Delta$ whose minimum
% $\Delta^\star=-H^{-1}g$ sets both the direction and the length of the step.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  \definecolor{red}{HTML}{C0392B}
  % the true loss curve (a bowl)
  \draw[acc, very thick] plot[domain=-2.6:3.4, samples=120]
    (\x, {0.35*(\x-0.4)*(\x-0.4) + 0.25});
  \node[acc, anchor=west, font=\footnotesize] at (2.7,3.15) {\texttt{true loss}};
  % current iterate on the left flank
  \coordinate (c) at (-2.0,2.06);
  \fill[black] (c) circle (2.4pt);
  \node[anchor=east, font=\footnotesize] at (-2.2,2.06) {\texttt{here}};
  % first-order tangent plane (a line): slope of the bowl at c is 0.7*(-2.4) = -1.68
  \draw[red, thick] (-3.0,3.74) -- (0.55,-2.15);
  \node[red, anchor=east, font=\footnotesize] at (-2.75,3.55) {\texttt{tangent plane}};
  \node[red, anchor=west, font=\footnotesize] at (0.65,-2.1) {\texttt{no bottom}};
  % second-order model: a wider quadratic tangent at c, with its own minimum
  \draw[green, thick] plot[domain=-3.2:1.2, samples=80]
    (\x, {0.55*(\x+0.85)*(\x+0.85) + 0.42});
  \node[green, anchor=north east, font=\footnotesize] at (-2.7,2.55) {\texttt{quadratic model}};
  % minimum of the second-order model
  \fill[green] (-0.85,0.42) circle (2.2pt);
  \draw[green, densely dashed] (-0.85,0.42) -- (-1.9,-0.55);
  \node[green, anchor=east, font=\footnotesize] at (-1.9,-0.55) {\texttt{model min}};
  % the Newton step from here to model min
  \draw[black, ->, thick] (-2.0,1.55) -- (-0.9,1.55);
  \node[anchor=south, font=\footnotesize] at (-1.45,1.55) {\texttt{step}};
\end{tikzpicture}
$$

> **Definition (Newton's method).** The second-order update
> $$\theta \;\gets\; \theta \;-\; H^{-1} g,$$
> which moves directly to the minimum of the second-order Taylor model of $\mathcal{L}$
> at $\theta$. When $\mathcal{L}$ is _exactly_ quadratic and $H \succ 0$, it reaches the
> true minimum in a single step from any starting point.

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.

$$
% caption: On an ill-conditioned bowl, gradient descent (red) zig-zags across the
% contours while the Newton step (green) rescales by curvature and lands in one move.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  \definecolor{red}{HTML}{C0392B}
  % elongated elliptical contours (ill-conditioned bowl)
  \foreach \r in {0.5,1.0,1.5,2.0,2.5}
    \draw[black] (0,0) ellipse ({\r*2.0} and {\r*0.72});
  % minimum
  \fill[green] (0,0) circle (2.6pt);
  \node[green, anchor=north] at (1.5,-1.9) {\texttt{minimum}};
  % start point on the long axis, high up
  \coordinate (s) at (-3.7,1.35);
  \fill[black] (s) circle (2.4pt);
  \node[anchor=east, font=\footnotesize] at (-3.85,1.35) {\texttt{start}};
  % gradient descent: zig-zag across the narrow valley
  \draw[red, very thick, ->] (s) -- (-2.45,-0.55);
  \draw[red, very thick, ->] (-2.45,-0.55) -- (-1.55,0.78);
  \draw[red, very thick, ->] (-1.55,0.78) -- (-0.95,-0.42);
  \draw[red, very thick, ->] (-0.95,-0.42) -- (-0.5,0.32);
  \node[red, anchor=north, font=\footnotesize] at (-1.5,-2.0) {\texttt{gradient descent zig-zags}};
  % newton: straight to the minimum
  \draw[green, very thick, ->] (s) -- (0,0);
  \node[green, anchor=south, font=\footnotesize] at (-1.5,2.0) {\texttt{Newton: one step}};
\end{tikzpicture}
$$

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.[^gf-hessian] To see the one-step claim algebraically, take an exact quadratic
$\mathcal{L}(\theta) = \tfrac12(\theta - \theta^\star)^\top H (\theta - \theta^\star)$,
for which $g = H(\theta - \theta^\star)$:

$$
\theta - H^{-1} g
= \theta - H^{-1} H (\theta - \theta^\star)
= \theta - (\theta - \theta^\star)
= \theta^\star.
$$

> **Theorem (Local quadratic convergence).** If $\mathcal{L}$ is twice continuously
> differentiable, $H(\theta^\star) \succ 0$ at a minimum $\theta^\star$, and the iterate
> starts close enough to $\theta^\star$, Newton's method converges _quadratically_:
> $\norm{\theta_{k+1} - \theta^\star} \le C\,\norm{\theta_k - \theta^\star}^2$.
> The number of correct digits roughly doubles each step.

> **Proof sketch.** The Newton update is $\theta_{k+1} = \theta_k - H(\theta_k)^{-1} g(\theta_k)$.
> Taylor-expand $g$ about $\theta^\star$: $g(\theta_k) = H(\theta^\star)(\theta_k - \theta^\star) + O(\norm{\theta_k - \theta^\star}^2)$,
> and $H(\theta_k) = H(\theta^\star) + O(\norm{\theta_k - \theta^\star})$. Substituting,
> the linear error term cancels exactly — that is what $H^{-1}$ buys — leaving only the
> quadratic remainder, so $\theta_{k+1} - \theta^\star = O(\norm{\theta_k - \theta^\star}^2)$. $\qed$

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

$$
-H_\phi^{-1} g_\phi
= -\parens{A^{-\top} H A^{-1}}^{-1} A^{-\top} g
= -A\,H^{-1} A^\top A^{-\top} g
= A\parens{-H^{-1} g},
$$

the original step carried through $A$. Gradient descent has no such property — the
direction $-g_\phi = -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 — $\binom{d+1}{2} \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.

$$
% caption: The Hessian grows quadratically in $d$ while the gradient grows only
% linearly, so at $d=10^9$ it already holds $10^{18}$ entries.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{acclo}{HTML}{A7B5FB}
  % small gradient column (linear in d)
  \node[font=\footnotesize, anchor=south] at (0.35,2.5) {\texttt{gradient} $g$};
  \draw[acc, thick, fill=acclo!40] (0.05,0) rectangle (0.65,2.3);
  \node[font=\scriptsize, anchor=north] at (0.35,-0.1) {$d$};
  % three growing Hessian squares (quadratic in d)
  \node[font=\footnotesize, anchor=south] at (2.5,2.5) {\texttt{small} $d$};
  \draw[acc, thick, fill=acclo!30] (2.0,1.3) rectangle (3.0,2.3);
  \node[font=\footnotesize, anchor=south] at (5.2,2.5) {\texttt{medium }$d$};
  \draw[acc, thick, fill=acclo!30] (4.2,0.5) rectangle (6.2,2.5);
  \node[font=\footnotesize, anchor=south] at (9.1,2.85) {\texttt{large} $d$};
  \draw[acc, thick, fill=acclo!30] (7.4,-0.6) rectangle (10.8,2.8);
  \node[font=\footnotesize, anchor=center, text=acc] at (9.1,1.1) {$d$ \texttt{by} $d$};
  % the explosion annotation
  \node[font=\footnotesize, text=black, anchor=north] at (5.4,-0.95) {\texttt{storage grows as} $d^2$, \texttt{inversion as} $d^3$};
\end{tikzpicture}
$$

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.

| $d$ (parameters) | gradient entries | Hessian entries | inversion cost $\sim d^3$ |
| --- | --- | --- | --- |
| $10^3$ | $10^3$ | $10^6$ | $10^9$ |
| $10^6$ | $10^6$ | $10^{12}$ | $10^{18}$ |
| $10^9$ | $10^9$ | $10^{18}$ | $10^{27}$ |

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](/deep-learning/optimization/the-optimization-landscape).
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

$$
\Delta_i = -\frac{g_i}{\lambda_i}.
$$

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}$.[^gf-saddlefree]

| Eigendirection | $\lambda_i$ | first-order step $-g_i$ | Newton step $-g_i/\lambda_i$ | effect |
| --- | --- | --- | --- | --- |
| up-curving | $\lambda_i > 0$ | downhill | downhill, rescaled | toward minimum (correct) |
| down-curving | $\lambda_i < 0$ | downhill | _uphill_ (sign flips) | toward saddle (wrong) |
| flat | $\lambda_i \approx 0$ | small | blows up ($\div \approx 0$) | catastrophic step |

$$
% caption: Near a saddle, the negative-curvature axis flips the sign of its Newton
% component, so the Newton step (red) climbs back toward the saddle.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  \definecolor{red}{HTML}{C0392B}
  % hyperbolic contours of a saddle u1^2 - u2^2
  \foreach \c in {0.4,1.1,2.0} {
    \draw[black] plot[domain=-2.0:2.0, samples=50] ({sqrt(\c + \x*\x)}, {\x});
    \draw[black] plot[domain=-2.0:2.0, samples=50] ({-sqrt(\c + \x*\x)}, {\x});
    \draw[black] plot[domain=-2.0:2.0, samples=50] ({\x}, {sqrt(\c + \x*\x)});
    \draw[black] plot[domain=-2.0:2.0, samples=50] ({\x}, {-sqrt(\c + \x*\x)});
  }
  % the saddle
  \fill[black] (0,0) circle (2.6pt);
  \node[anchor=south west, font=\footnotesize] at (0.12,0.08) {\texttt{saddle}};
  % current iterate, below the saddle on the escape (down-curving) axis
  \coordinate (p) at (0,-2.2);
  \fill[acc] (p) circle (2.4pt);
  \node[acc, anchor=west, font=\footnotesize] at (0.15,-2.2) {\texttt{here}};
  % gradient step: continues away from the saddle (downhill, correct)
  \draw[green, very thick, ->] (p) -- (0,-3.3);
  \node[green, anchor=west, font=\footnotesize] at (0.15,-3.1) {\texttt{gradient leaves}};
  % newton step: sign-flipped, climbs back toward the saddle
  \draw[red, very thick, ->] (p) -- (0,-0.55);
  \node[red, anchor=center, font=\footnotesize] at (0,2.45) {\texttt{Newton climbs in}};
\end{tikzpicture}
$$

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.

| Assumption Newton needs | Reality in deep learning | Consequence |
| --- | --- | --- |
| $H$ small and storable | $d$ up to $10^9$, $H$ is $d \times d$ | cannot form or invert $H$ |
| $H \succ 0$ (a bowl) | mostly saddles, mixed-sign $H$ | attracted to saddles |
| $g, H$ exact | minibatch estimates, noisy | division amplifies the noise |
| problem deterministic | stochastic objective | step direction unstable |

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.

> **Definition (Conjugate directions).** Two directions $d_i, d_j$ are **conjugate**
> with respect to $H$ if $d_i^\top H d_j = 0$ for $i \ne j$. Minimizing a quadratic
> exactly along $d_j$ leaves its minimum along every earlier conjugate $d_i$ undisturbed,
> so $n$ conjugate line-searches reach the exact minimum of an $n$-dimensional quadratic
> — with no zig-zag and no Hessian inverse.

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.

$$
% caption: Steepest descent (red) takes orthogonal steps that zig-zag across the
% stretched valley; conjugate gradients (green) bends each direction to be
% $H$-conjugate to the last and reaches the minimum without retracing.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  \definecolor{red}{HTML}{C0392B}
  % elongated contours
  \foreach \r in {0.5,1.0,1.5,2.0,2.6}
    \draw[black] (0,0) ellipse ({\r*1.95} and {\r*0.7});
  \fill[acc] (0,0) circle (2.4pt);
  \node[acc, anchor=north] at (1.6,-1.9) {\texttt{minimum}};
  % start
  \coordinate (s) at (-3.6,1.45);
  \fill[black] (s) circle (2.2pt);
  \node[anchor=south, font=\footnotesize] at (-3.6,1.6) {\texttt{start}};
  % steepest descent: many orthogonal zig-zag steps
  \draw[red, very thick, ->] (s) -- (-2.3,-0.55);
  \draw[red, very thick, ->] (-2.3,-0.55) -- (-1.45,0.72);
  \draw[red, very thick, ->] (-1.45,0.72) -- (-0.85,-0.4);
  \draw[red, very thick, ->] (-0.85,-0.4) -- (-0.42,0.3);
  \node[red, anchor=south, font=\footnotesize] at (-1.6,2.05) {\texttt{steepest descent}};
  % conjugate gradients: two steps to the minimum
  \draw[green, very thick, ->] (s) -- (-1.6,-1.35);
  \draw[green, very thick, ->] (-1.6,-1.35) -- (0,0);
  \node[green, anchor=north, font=\footnotesize] at (-1.6,-2.05) {\texttt{conjugate gradients (2 steps)}};
\end{tikzpicture}
$$

```algorithm
caption: $\textsc{ConjugateGradients}(\mathcal{L}, \theta_0)$ — minimize without forming $H$
$g_0 \gets \nabla\mathcal{L}(\theta_0)$
$d_0 \gets -g_0$ // first direction is steepest descent
for $t \gets 0, 1, 2, \dots$ do
  $a_t \gets \arg\min_a \mathcal{L}(\theta_t + a\,d_t)$ // line search along $d_t$
  $\theta_{t+1} \gets \theta_t + a_t\,d_t$
  $g_{t+1} \gets \nabla\mathcal{L}(\theta_{t+1})$
  if $\norm{g_{t+1}} \le \epsilon$ then // converged
    return $\theta_{t+1}$
  $b_t \gets (g_{t+1}^\top g_{t+1}) / (g_t^\top g_t)$ // Fletcher-Reeves coefficient
  $d_{t+1} \gets -g_{t+1} + b_t\,d_t$ // bend to be conjugate to $d_t$
```

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.[^gf-cg]

## 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

$$
s_t = \theta_{t+1} - \theta_t,
\qquad
y_t = g_{t+1} - g_t,
$$

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):

$$
B_{t+1}
= \parens{I - \frac{s_t y_t^\top}{y_t^\top s_t}} B_t
  \parens{I - \frac{y_t s_t^\top}{y_t^\top s_t}}
  + \frac{s_t s_t^\top}{y_t^\top s_t}.
$$

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.

$$
% caption: BFGS reads curvature off successive gradients: the secant relation
% $y_t \approx H s_t$ feeds a rank-two update of the inverse-Hessian estimate.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  \definecolor{red}{HTML}{C0392B}
  % a 1D loss curve (a bowl) with two iterates
  \draw[->, thick] (-0.3,0) -- (8.4,0) node[right, font=\footnotesize] {\texttt{parameter}};
  \draw[->, thick] (0,-0.3) -- (0,3.6) node[above, font=\footnotesize] {\texttt{loss}};
  \draw[acc, very thick] plot[domain=0.5:7.8, samples=120]
    (\x, {0.45 + 0.28*(\x-5.3)*(\x-5.3)});
  % iterate t at left, steep tangent
  \fill[black] (1.7,3.06) circle (2.4pt);
  \node[anchor=north west, font=\footnotesize] at (1.85,2.95) {\texttt{here at} $t$};
  \draw[red, thick] (0.7,3.62) -- (2.7,1.5);
  \node[red, anchor=south, font=\footnotesize] at (1.7,3.7) {\texttt{steep gradient}};
  % iterate t+1 at right of center, gentle tangent
  \fill[black] (4.0,0.94) circle (2.4pt);
  \node[anchor=south, font=\footnotesize] at (3.5,1.85) {\texttt{here at} $t{+}1$};
  \draw[green, thick] (3.0,1.66) -- (5.0,0.5);
  \node[green, anchor=center, font=\footnotesize] at (5.5,1.35) {\texttt{gentler gradient}};
  % the step s_t between them
  \draw[black, <->, thick] (1.7,-0.35) -- (4.0,-0.35);
  \node[black, anchor=north, font=\footnotesize] at (2.85,-0.4) {\texttt{step} $s$ \texttt{(change in weight)}};
  % annotation: change in slope = curvature
  \node[anchor=west, font=\footnotesize, text=black] at (5.2,2.7) {\texttt{slope change gives curvature}};
\end{tikzpicture}
$$

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).[^gf-bfgs]

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.

| Method | Hessian use | Memory | Per-step cost | Reaches quadratic min in |
| --- | --- | --- | --- | --- |
| Newton | exact $H$, inverts it | $O(d^2)$ | $O(d^3)$ | $1$ step |
| Conjugate gradients | implicit (via $Hd$ products) | $O(d)$ | $O(d)$ + line search | $\le d$ steps |
| BFGS | approximates $H^{-1}$ from $(s, y)$ | $O(d^2)$ | $O(d^2)$ | $\le d$ steps |
| L-BFGS | last $m$ curvature pairs | $O(md)$ | $O(md)$ | superlinear, $\le d$ steps |

## 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} = \tfrac12\sum_i r_i(\theta)^2$,
write the Jacobian of the residual vector as $J = \partial r/\partial\theta$. The exact
Hessian is

$$
H = J^\top J \;+\; \sum_i r_i\, \nabla^2 r_i.
$$

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 \succeq 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,

$$
F = \mathbb{E}_{x}\,\mathbb{E}_{y\sim p_\theta}\!\brackets{ \nabla_\theta \log p_\theta(y\mid x)\,
\nabla_\theta \log p_\theta(y\mid x)^\top }.
$$

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:

$$
\theta \;\gets\; \theta \;-\; \eta\, F^{-1} g.
$$

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 \tfrac12\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.

| Method | Curvature matrix | Trick that makes it affordable |
| --- | --- | --- |
| Natural gradient | Fisher information $F$ | curvature in output-distribution space, always $\succeq 0$ |
| K-FAC | block-diagonal $F$ | factor each layer's block as a Kronecker product $A \otimes G$ |
| Hessian-free | Hessian $H$ (implicit) | inner CG using $Hv$ products via [autodiff](/deep-learning/neural-networks/backpropagation) |

> **Definition (Hessian-free optimization).** A truncated-Newton method that solves
> $H\Delta = -g$ approximately by an inner [conjugate-gradient](#conjugate-gradients)
> loop, supplying each Hessian–vector product $Hv$ through a single extra
> backpropagation pass — never forming $H$. Curvature is used; the $d \times d$ matrix
> is never materialized.

> **Remark (K-FAC).** Kronecker-Factored Approximate Curvature approximates each layer's
> Fisher block as a Kronecker product $A \otimes G$ of two small factors — one from the
> layer's inputs, one from its output gradients. Inverting a Kronecker product is just
> inverting the two small factors, $(A \otimes G)^{-1} = A^{-1} \otimes G^{-1}$,
> collapsing a per-layer $O(d^3)$ inverse to two cheap ones.

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$,

$$
Hv = \nabla_\theta\parens{ g(\theta)^\top v } = \nabla_\theta\parens{ (\nabla_\theta\mathcal{L})^\top 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.

$$
% caption: Pearlmutter's trick computes $Hv$ as the gradient of the scalar $g^\top v$:
% forward pass, backward pass for $g$, contract with $v$, backward again. The matrix
% $H$ is never formed; storage stays $O(d)$.
\begin{tikzpicture}[>=stealth, font=\footnotesize, node distance=6mm]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  \tikzset{
    bx/.style={draw=acc, thick, fill=acc!8, minimum width=20mm, minimum height=8mm, align=center},
    op/.style={draw=black, thick, fill=black!4, minimum width=17mm, minimum height=8mm, align=center},
  }
  \node[bx] (fwd) {\texttt{forward}\\\texttt{pass}};
  \node[bx, right=13mm of fwd] (bwd) {\texttt{backward}\\\texttt{pass}};
  \node[op, minimum width=22mm] (dot) [right=13mm of bwd] {\texttt{dot} $g,v$};
  \node[bx, right=13mm of dot] (bwd2) {\texttt{backward}\\\texttt{again}};
  \node[op, right=13mm of bwd2] (out) {$Hv$};
  \draw[->, thick, acc] (fwd) -- (bwd);
  \draw[->, thick, acc] (bwd) -- node[above, font=\scriptsize, black]{$g$} (dot);
  \draw[->, thick, acc] (dot) -- (bwd2);
  \draw[->, thick, green] (bwd2) -- (out);
  % v injected into the dot product
  \node[op, above=9mm of dot] (v) {\texttt{vector} $v$};
  \draw[->, thick, black] (v) -- (dot);
  % annotation
  \node[font=\scriptsize, black, anchor=north, text width=76mm, align=center]
    at (dot|-fwd.south) {\texttt{one forward pass, two backward passes; H never stored}};
\end{tikzpicture}
$$

Curvature is _accessible_ — just never as a stored matrix.[^gf-natural]

## 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.

$$
% caption: The cost/accuracy ladder. Climbing captures more of the true curvature $H$
% (diagonal $\to$ block $\to$ low-rank $\to$ full) at rising per-step cost
% ($O(d)\to O(d)\to O(md)\to O(d^3)$). Deep learning lives on the lowest rungs.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  % axes
  \draw[->, thick, black] (0,0) -- (0,6.2) node[above, font=\footnotesize, black, anchor=south west]{\texttt{curvature captured}};
  \draw[->, thick, black] (0,0) -- (10.6,0) node[below, font=\footnotesize, black]{\texttt{cost per step}};
  \tikzset{rung/.style={draw=acc, thick, fill=acc!8, minimum width=40mm, minimum height=8.5mm, align=center, font=\scriptsize}}
  \node[rung, anchor=west] (a) at (0.4,0.85) {\texttt{gradient descent}};
  \node[rung, anchor=west] (b) at (1.9,2.05) {\texttt{Adam (diagonal)}};
  \node[rung, anchor=west] (c) at (3.4,3.25) {\texttt{K-FAC (block)}};
  \node[rung, anchor=west] (d) at (4.9,4.45) {\texttt{L-BFGS (low-rank)}};
  \node[rung, anchor=west, fill=green!10, draw=green] (e) at (6.4,5.65) {\texttt{full Newton}};
  % connecting climb in the horizontal overlap band, pointing up into the next rung
  \draw[->, green, thick] (3.15,1.28) -- (3.15,1.62);
  \draw[->, green, thick] (4.65,2.48) -- (4.65,2.82);
  \draw[->, green, thick] (6.15,3.68) -- (6.15,4.02);
  \draw[->, green, thick] (7.65,4.88) -- (7.65,5.22);
  % cost labels to the right of each rung, clear of boxes and arrows
  \node[font=\scriptsize, black, anchor=west] at (a.east) {\;$O(d)$};
  \node[font=\scriptsize, black, anchor=west] at (b.east) {\;$O(d)$};
  \node[font=\scriptsize, black, anchor=west] at (c.east) {\;$O(d)$};
  \node[font=\scriptsize, black, anchor=west] at (d.east) {\;$O(md)$};
  \node[font=\scriptsize, black, anchor=west] at (e.east) {\;$O(d^3)$};
\end{tikzpicture}
$$

## 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.

| | First-order (SGD, Adam) | Second-order (Newton, quasi-Newton) |
| --- | --- | --- |
| Information used | gradient $g$ | gradient $g$ + curvature $H$ |
| Per-step cost | $O(d)$ | $O(d^2)$–$O(d^3)$ (exact), $O(md)$ (L-BFGS) |
| Steps to converge | many | few (quadratic / superlinear) |
| Robust to conditioning | no (needs momentum / adaptive rates) | yes (rescales by curvature) |
| Robust to minibatch noise | yes (noise averages out) | no (division amplifies it) |
| Saddle behavior | escapes (with noise / momentum) | _attracted_ unless modified |
| Deep-learning suitability | the default | rare; full-batch or approximate only |

> **Remark (Why Adam, not Newton).** [Adam and the adaptive methods](/deep-learning/optimization/momentum-and-adaptive-methods)
> can be read as the cheapest possible second-order approximation: they precondition by
> a _diagonal_ estimate of curvature (one scale per parameter, $O(d)$ memory) instead
> of the full $d \times d$ Hessian. They keep the conditioning benefit of second-order
> methods at first-order cost, and that trade is why they, not
> Newton, dominate deep learning.

## 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.[^kfac]
- **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.[^shampoo]
- **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.[^sophia]

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 \gets \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.

[^gf-secondorder]: **Goodfellow**, _Deep Learning_, §8.6 — Approximate Second-Order Methods: why curvature converges in fewer steps and the price that keeps exact second-order methods out of deep learning.
[^gf-hessian]: **Goodfellow**, _Deep Learning_, §4.3.1 — Beyond the Gradient: Jacobian and Hessian Matrices: the second-order Taylor model, the Newton step $-H^{-1}g$, and condition-number invariance.
[^gf-saddlefree]: **Goodfellow**, _Deep Learning_, §8.2.3, §8.6 — Newton's attraction to saddle points (dividing by a negative eigenvalue) and the saddle-free / damped fixes that force positive curvature.
[^gf-cg]: **Goodfellow**, _Deep Learning_, §8.6 — Conjugate Gradients: $H$-conjugate search directions that eliminate the zig-zag of steepest descent without forming the Hessian.
[^gf-bfgs]: **Goodfellow**, _Deep Learning_, §8.6 — BFGS and L-BFGS: learning $B \approx H^{-1}$ from the secant condition, and the limited-memory recursion that drops storage to $O(md)$.
[^gf-natural]: **Goodfellow**, _Deep Learning_, §8.6, §12.1 — natural gradient via the Fisher information matrix, K-FAC's Kronecker factorization, and Hessian-free optimization through $Hv$ products.
[^kfac]: **Martens & Grosse**, _Optimizing Neural Networks with Kronecker-Factored Approximate Curvature_, ICML 2015 — the K-FAC approximation of the Fisher block as a Kronecker product of two small factors.
[^shampoo]: **Gupta, Koren & Singer**, _Shampoo: Preconditioned Stochastic Tensor Optimization_, ICML 2018 — per-dimension Kronecker-factored preconditioning, later scaled to production language models.
[^sophia]: **Liu et al.**, _Sophia: A Scalable Stochastic Second-Order Optimizer for Language Model Pre-training_, ICLR 2024 — a clipped diagonal-Hessian preconditioner that roughly halves pre-training steps.