Limited-Memory L-BFGS and Inexact (Truncated) Newton-CG

Intuition Beginner

Full BFGS keeps a complete curvature table — a square array the size of the problem squared. For a problem in a million variables that array has a trillion entries, far too many to store. Limited-memory BFGS notices that the table was never used on its own: it only ever multiplied the gradient to produce a search direction. So instead of keeping the table, L-BFGS keeps just the raw measurements that built it — the last handful of step and slope-change pairs — and rebuilds the one multiplication it needs on the fly. Keeping the last five or ten pairs is usually enough, so the storage drops from a trillion numbers to a few million.

There is a second large-scale idea in this unit. Newton's method aims each step by solving a system of equations involving the curvature. Solving that system exactly is expensive, and early on it is wasted effort — the curvature you are solving against is only a local snapshot, so an exact answer to an approximate question buys you little. A truncated Newton method solves the system only roughly, using a cheap inner loop that you stop early, and tightens the accuracy as you close in on the answer. Far away you take quick sloppy steps; near the solution you solve carefully and recover Newton's full speed.

Picture packing for a trip with one small bag (L-BFGS): you cannot carry your whole closet, so you bring the few recent outfits that actually tell you what the weather has been like, and reconstruct a sensible choice from those. The truncated idea is like sketching a route: a rough sketch is fine while you are far away, and you fill in the fine detail only as you near your destination.

Visual Beginner

Picture the curvature information as a stack of note cards. Full BFGS would keep one enormous card holding the entire table. L-BFGS instead keeps a short stack of the most recent cards, each recording one step and the slope change across it; when an older card is pushed off the bottom, its information is simply forgotten. To get a search direction, the method runs down the stack and back up it once — a quick two-pass sweep — which is exactly the two-loop recursion.

method	what it stores	cost per step	best when
full BFGS	whole curvature table	medium	a few hundred variables
L-BFGS	last $m$ measurement pairs	cheap, fixed	millions of variables
truncated Newton	nothing (rebuilds curvature action)	varies (inner loop)	curvature action is cheap and accurate

The takeaway: both methods refuse to store a full curvature table. L-BFGS keeps a short stack of recent measurements and sweeps it twice; truncated Newton keeps no table at all and instead solves the Newton equations approximately with an early-stopped inner loop.

Worked example Beginner

Let us run the L-BFGS direction sweep with a memory of just one pair, on a two-variable problem, so every number is concrete. Suppose the most recent step was $s=(1,0)$ and the slope changed by $y=(2,0)$ across it. The current gradient is $g=(3,1)$, and we want the search direction.

First compute the helper number $\rho =1/(y\cdot s)$. Here $y\cdot s=2\cdot 1+0\cdot 0=2$, so $\rho =1/2=0.5$.

Backward pass. Set $q=g=(3,1)$. Compute $a=\rho (s\cdot q)=0.5\cdot (1\cdot 3+0\cdot 1)=0.5\cdot 3=1.5$. Update $q\leftarrow q-ay=(3,1)-1.5(2,0)=(3-3,1)=(0,1)$.

Scaling. The initial guess scales the identity by $\gamma =(s\cdot y)/(y\cdot y)=2/(2^2+0^2)=2/4=0.5$. Apply it: $r=\gamma q=0.5(0,1)=(0,0.5)$.

Forward pass. Compute $b=\rho (y\cdot r)=0.5(2\cdot 0+0\cdot 0.5)=0$. Update $r\leftarrow r+(a-b)s=(0,0.5)+(1.5-0)(1,0)=(1.5,0.5)$.

The direction is $-r=(-1.5,-0.5)$.

What this tells us: from one stored pair and the current gradient, two short passes produce a search direction that already bends the plain downhill direction $-(3,1)$ toward the measured curvature — it stretched the first coordinate (where curvature was large) less than the second. We never built or stored a curvature matrix; we only swept the stored measurements.

Check your understanding Beginner

Formal definition Intermediate+

Let $f:{\mathbb{R}}^n\to \mathbb{R}$ be twice continuously differentiable with gradient $g(x)=\nabla f(x)$ and Hessian ${\nabla }^{2}f(x)$. This unit treats two large-scale strategies that both refuse to store an $n\times n$ matrix.

Limited-memory BFGS. Recall from 44.03.04 the BFGS inverse update $H_{k+1}=V_k^{\top }{H}_k{V}_k+{\rho }_k{s}_k{s}_k^{\top }$ with $s_k=x_{k+1}-x_k$, $y_k=g_{k+1}-g_k$, ${\rho }_k=1/(y_k^{\top }{s}_k)$, and $V_k=I-{\rho }_k{y}_k{s}_k^{\top }$. L-BFGS never forms $H_k$. It stores only the most recent $m$ correction pairs $\{(s_i,y_i){\}}_{i=k-m}^{k-1}$ and defines the implicit inverse-Hessian approximation by applying $m$ BFGS updates to a scaled initial matrix $H_k^0={\gamma }_{k}I$:

$$H_k=(V_{k-1}^{\top }\cdots V_{k-m}^{\top })H_k^0(V_{k-m}\cdots V_{k-1})+\sum _{i=k-m}^{k-1}(V_{k-1}^{\top }\cdots V_{i+1}^{\top }){\rho }_i{s}_i{s}_i^{\top }(V_{i+1}\cdots V_{k-1}),$$

with the standard scaling ${\gamma }_k=(s_{k-1}^{\top }{y}_{k-1})/(y_{k-1}^{\top }{y}_{k-1})$, which sizes the unit initial matrix to the most recent curvature. The product $H_k{g}_k$ is computed by the two-loop recursion: set $q\leftarrow g_k$; for $i=k-1,\dots ,k-m$ set ${\alpha }_i\leftarrow {\rho }_i{s}_i^{\top }q$ and $q\leftarrow q-{\alpha }_i{y}_i$; set $r\leftarrow {\gamma }_{k}q$; for $i=k-m,\dots ,k-1$ set $\beta \leftarrow {\rho }_i{y}_i^{\top }r$ and $r\leftarrow r+({\alpha }_i-\beta )s_i$; return $H_k{g}_k=r$. The cost is $4mn$ multiply-adds and $O(mn)$ storage, with no matrix formed and no linear solve. The search direction is $p_k=-r$, and a line search satisfying the Wolfe conditions 44.03.01 supplies the step, guaranteeing $s_k^{\top }{y}_k>0$ so each stored pair is curvature-admissible.

Inexact (truncated) Newton. An inexact Newton method computes a step $p_k$ that solves the Newton system only approximately: with residual $r_k={\nabla }^{2}f(x_k)p_k+g_k$, the step is accepted when

$$\Vert r_k\Vert \le {\eta }_k\Vert g_k\Vert ,$$

where $\{{\eta }_k\}\subset [0,1)$ is the forcing sequence. The truncated Newton-CG method realises this by applying the linear conjugate gradient method 43.07.04 to ${\nabla }^{2}f(x_k)p=-g_k$, terminating CG as soon as the residual test holds. Because CG accesses its matrix only through matrix-vector products, the Hessian is never assembled: each inner iteration needs one Hessian-vector product ${\nabla }^{2}f(x_k)v$, supplied either by automatic differentiation or by the finite-difference approximation ${\nabla }^{2}f(x_k)v\approx (g(x_k+hv)-g(x_k))/h$ of 44.03.08, making the method Hessian-free. When CG generates a direction $d_j$ with nonpositive curvature, $d_j^{\top }{\nabla }^{2}f(x_k)d_j\le 0$, the iteration is truncated: on the first inner iteration it returns $-g_k$ (steepest descent), and thereafter it returns the current inner iterate, which is a descent direction. The new symbols $H_k^0$, ${\gamma }_k$, the two-loop scalars ${\alpha }_i,\beta$, the forcing sequence ${\eta }_k$, the inner residual $r_k$, and the Hessian-vector product are recorded in _meta/NOTATION.md.

Counterexamples to common slips Intermediate+

• "L-BFGS is BFGS with a smaller matrix." L-BFGS stores no matrix at all. It stores $2m$ vectors and reconstructs only the action $H_k{g}_k$. The implicit $H_k$ it represents changes its initial matrix ${\gamma }_{k}I$ every iteration, so it is not a truncation of a single evolving matrix but a fresh $m$-update product each step.

• "A more accurate inner solve is always better." Solving the Newton system to high accuracy far from the solution wastes inner iterations against a Hessian that is only a local model; the forcing sequence deliberately keeps ${\eta }_k$ large early. Over-solving raises cost with no gain in the outer rate, which is capped by how good the local quadratic model is.

• "CG can be run to convergence on any Newton system." Away from a minimiser the Hessian can be indefinite, and CG breaks down on indefinite systems. Truncated Newton must detect nonpositive curvature and stop, returning a guaranteed descent direction; ignoring this produces an ascent step or a CG division by a nonpositive curvature.

• "Hessian-free means derivative-free." The Hessian-vector product still needs first derivatives: automatic differentiation differentiates the gradient, and the finite-difference form evaluates the gradient at a perturbed point. The saving is the $n^2$ storage and assembly of the Hessian, not the use of derivatives.

Key theorem with proof Intermediate+

The signature result is that the two-loop recursion is not an approximation of the BFGS direction but computes it exactly from the stored pairs: the $O(mn)$ sweep returns precisely $H_k{g}_k$ for the limited-memory $H_k$ defined by $m$ BFGS updates applied to ${\gamma }_{k}I$. This is what makes L-BFGS a genuine quasi-Newton method 44.03.04 rather than a heuristic, and it is the reason the line-search descent theory carries over unchanged.

Theorem (two-loop recursion computes the L-BFGS direction). Let $H_k^0={\gamma }_{k}I\succ 0$ and let $\{(s_i,y_i){\}}_{i=k-m}^{k-1}$ satisfy $s_i^{\top }{y}_i>0$. Define $H_k$ by applying the BFGS inverse update successively to $H_k^0$ with the pairs $(s_{k-m},y_{k-m}),\dots ,(s_{k-1},y_{k-1})$. Then the two-loop recursion returns $r=H_k{g}_k$, and $H_k\succ 0$, so $p_k=-H_k{g}_k$ is a descent direction.

Proof. Write $V_i=I-{\rho }_i{y}_i{s}_i^{\top }$, ${\rho }_i=1/(y_i^{\top }{s}_i)$. One BFGS inverse update is $H^{+}=V_i^{\top }H{V}_i+{\rho }_i{s}_i{s}_i^{\top }$. Applying $H^{+}$ to a vector $u$ gives

$$H^{+}u=V_i^{\top }H(V_{i}u)+{\rho }_i{s}_i(s_i^{\top }u)=V_i^{\top }H(u-{\rho }_i{y}_i(s_i^{\top }u))+{\rho }_i(s_i^{\top }u)s_i.$$

Set ${\alpha }_i={\rho }_i{s}_i^{\top }u$. The vector $u-{\alpha }_i{y}_i$ is precisely the backward-loop update $q\leftarrow q-{\alpha }_i{y}_i$, and $V_i^{\top }w=w-{\rho }_i{s}_i(y_i^{\top }w)$ acting on $w=H(u-{\alpha }_i{y}_i)$ produces $H(u-{\alpha }_i{y}_i)-{\rho }_i{s}_i(y_i^{\top }H(u-{\alpha }_i{y}_i))$. Adding the term ${\alpha }_i{s}_i$ gives

$$H^{+}u=H(u-{\alpha }_i{y}_i)+({\alpha }_i-{\rho }_i{y}_i^{\top }H(u-{\alpha }_i{y}_i))s_i.$$

Reading this recursively: the innermost operation is $H_k^0$ applied to the fully-reduced $q$, and each layer outward restores one $s_i$ with coefficient ${\alpha }_i-{\beta }_i$ where ${\beta }_i={\rho }_i{y}_i^{\top }r$ and $r$ is the partial result. The backward loop computes all ${\alpha }_i$ and reduces $g_k$ to $q=(V_{k-m}\cdots V_{k-1})g_k$; the scaling computes $r=H_k^{0}q$; the forward loop applies the $V_i^{\top }$ factors and adds the ${\rho }_i{s}_i{s}_i^{\top }$ terms in the correct order, yielding $r=H_k{g}_k$. Each ${\alpha }_i$ stored in the backward loop is exactly the ${\alpha }_i$ reused in the forward loop because $s_i^{\top }q$ at the backward visit equals $s_i^{\top }u$ at that layer.

For positive-definiteness, hereditary positive-definiteness of the BFGS inverse update 44.03.04 applies to each of the $m$ updates: $H_k^0={\gamma }_{k}I\succ 0$ since ${\gamma }_k=(s_{k-1}^{\top }{y}_{k-1})/(y_{k-1}^{\top }{y}_{k-1})>0$ by the curvature condition, and each update preserves positive-definiteness given $s_i^{\top }{y}_i>0$. Hence $H_k\succ 0$, so $g_k^{\top }{p}_k=-g_k^{\top }{H}_k{g}_k<0$ whenever $g_k\ne 0$. $\square$

Bridge. This identity is the foundational reason L-BFGS inherits the quasi-Newton theory of 44.03.04 wholesale: because the two-loop output is the exact BFGS direction for the implicit $H_k$, hereditary positive-definiteness makes $p_k$ a descent direction, and the Wolfe line search of 44.03.01 then supplies the next admissible pair $s_k^{\top }{y}_k>0$, so the invariant propagates and the globalisation theory needs no change — this is exactly the same self-perpetuating curvature condition that drove full BFGS. The argument builds toward the compact-representation theorem of the Advanced results, where the same $m$-update product is rewritten as a single low-rank correction ${\gamma }_{k}I+UM{U}^{\top }$, and it appears again in the truncated-Newton analysis: there, positive curvature of the Hessian along the CG directions plays the role that $s_i^{\top }{y}_i>0$ plays here, both being the condition under which the inner model stays a genuine bowl. The central insight is that limited memory changes only the representation of the quasi-Newton operator, not its mathematics; putting these together, L-BFGS generalises BFGS to large $n$ by storing the data that built the matrix instead of the matrix, and the bridge is the recursion identity that makes the two provably the same direction.

Exercises Intermediate+

Advanced results Master

Large-scale unconstrained optimization is organised here around two refusals to store an $n\times n$ matrix — the limited-memory secant representation and the matrix-free inexact Newton solve — and the theory below settles the compact form of the L-BFGS operator, the exact convergence dichotomy of inexact Newton, the negative-curvature handling that makes Newton-CG robust, and the practical comparison that governs which method a large problem wants.

Theorem 1 (compact representation of L-BFGS). With $S_k=[s_{k-m}\cdots s_{k-1}]$, $Y_k=[y_{k-m}\cdots y_{k-1}]\in {\mathbb{R}}^{n\times m}$, let $L_k$ be the strictly-lower-triangular part of $S_k^{\top }{Y}_k$ and $D_k=\mathrm{diag}(s_i^{\top }{y}_i)$. The limited-memory inverse-Hessian admits the compact form

$$H_k={\gamma }_{k}I+\left[\begin{array}{cc}{S}_k& {\gamma }_k{Y}_k\end{array}\right]\left[\begin{array}{cc}{R}_k^{-\top }(D_k+{\gamma }_k{Y}_k^{\top }{Y}_k)R_k^{-1}& -R_k^{-\top }\\ -R_k^{-1}& 0\end{array}\right]\left[\begin{array}{c}{S}_k^{\top }\\ {\gamma }_k{Y}_k^{\top }\end{array}\right],$$

where $R_k$ is the upper-triangular part of $S_k^{\top }{Y}_k$. This expresses $H_k$ as a scaled identity plus a rank-$2m$ correction, so $H_{k}v$ costs $O(mn)$ via two tall-skinny matrix products and a small $2m\times 2m$ solve, the matrix-form equivalent of the two-loop recursion. The direct-Hessian approximation $B_k=H_k^{-1}$ has the dual compact form $B_k=\frac{1}{{\gamma }_k}I-[\frac{1}{{\gamma }_k}{S}_k{Y}_k]N_k^{-1}[\cdots {]}^{\top }$, which is what trust-region L-BFGS uses to apply $B_k$ and solve the subproblem ^{[Nocedal, J. & Wright, S. J. — Numerical Optimization (2nd ed.)]}.

Theorem 2 (Dembo-Eisenstat-Steihaug local convergence of inexact Newton). Let $F=\nabla f$ be continuously differentiable with $F^{\prime }(x^{\ast })={\nabla }^{2}f(x^{\ast })$ nonsingular at a root $x^{\ast }$, and run the inexact Newton iteration accepting $p_k$ whenever $\Vert F(x_k)+F^{\prime }(x_k)p_k\Vert \le {\eta }_k\Vert F(x_k)\Vert$. If ${\eta }_k\le {\eta }_{\max }<1$ then for $x_0$ near $x^{\ast }$ the iterates converge to $x^{\ast }$, q-linearly in the norm $\Vert v{\Vert }_{\ast }=\Vert F^{\prime }(x^{\ast })v\Vert$ with asymptotic factor ${\eta }_{\max }$. Moreover ${\eta }_k\to 0$ forces q-superlinear convergence, and ${\eta }_k=O(\Vert F(x_k){\Vert }^p)$ for $p\in (0,1]$ forces q-order $1+p$; in particular ${\eta }_k=O(\Vert F(x_k)\Vert )$ recovers Newton's quadratic rate. The forcing sequence is thus a precise dial trading inner-solve accuracy against outer convergence order, and the practical choice ${\eta }_k=\min \{\frac{1}{2},\sqrt{\Vert F(x_k)\Vert }\}$ (Eisenstat-Walker) buys superlinear convergence at a controlled inner cost ^{[Dembo, R. S., Eisenstat, S. C. & Steihaug, T. — Inexact Newton methods]}.

Theorem 3 (Newton-CG truncation and Steihaug-Toint trust-region CG). Applying CG 43.07.04 to ${\nabla }^{2}f(x_k)p=-g_k$ from $p^{(0)}=0$: while the curvature $d_j^{\top }{\nabla }^{2}f(x_k)d_j>0$ along every generated direction, the iterates $p^{(j)}$ monotonically decrease the quadratic model and increase in norm, $\Vert p^{(j)}\Vert$ strictly increasing. If a nonpositive-curvature direction appears, the line-search Newton-CG returns $p^{(j)}$ (or $-g_k$ if $j=0$), a descent direction; the trust-region variant (Steihaug-Toint) instead follows the CG path to the point where it crosses the trust-region boundary, exploiting the negative-curvature direction to reach the boundary. The monotone-norm-increase property is what makes the trust-region truncation well defined: the CG path leaves the trust region exactly once. This furnishes a globally convergent Hessian-free method whose inner cost adapts to the local difficulty of the Newton system ^{[Nocedal, J. & Wright, S. J. — Numerical Optimization (2nd ed.)]}.

Theorem 4 (global convergence of L-BFGS on convex problems). On a twice-differentiable $f$ with ${\nabla }^{2}f$ bounded, $mI⪯{\nabla }^{2}f⪯MI$ on the level set $\{f\le f(x_0)\}$ (strong convexity), L-BFGS with a Wolfe line search 44.03.01 and any fixed memory $m\ge 1$ converges to the unique minimiser at an R-linear rate, $f(x_k)-f^{\ast }\le c{\rho }^k$ for some $\rho \in (0,1)$. The proof bounds the eigenvalues of the implicit $H_k$ above and below uniformly in $k$ (the trace and determinant of $B_k=H_k^{-1}$ stay controlled because every stored pair satisfies $m\le y_i^{\top }{y}_i/(s_i^{\top }{y}_i)$ and $s_i^{\top }{y}_i/(s_i^{\top }{s}_i)\le M$), which keeps the search direction angle-bounded and feeds Zoutendijk's condition 44.03.01. Unlike full BFGS the rate is only R-linear, not superlinear, because the finite memory cannot accumulate the full curvature; superlinearity is recovered only as $m\to \infty$ ^{[Liu, D. C. & Nocedal, J. — On the limited memory BFGS method for large scale optimization]}.

Theorem 5 (L-BFGS versus truncated Newton: the cost dichotomy). For a problem with $n$ variables, L-BFGS costs a fixed $O(mn)$ per outer iteration with $O(mn)$ storage and a single gradient evaluation, and needs no Hessian access; truncated Newton-CG costs $O(n)$ storage but a variable number $j_k$ of Hessian-vector products per outer iteration, total inner cost $O(j_k\cdot \mathrm{cost}({\nabla }^{2}fv))$, where $j_k$ grows as the forcing sequence tightens. The asymptotic comparison: truncated Newton attains superlinear or quadratic outer convergence (Theorem 2) and wins when Hessian-vector products are cheap and accurate (e.g. supplied by automatic differentiation 44.03.08) and the Hessian is well-conditioned, so CG converges in few inner steps; L-BFGS attains only R-linear convergence (Theorem 4) but with a fixed, predictable, low per-iteration cost and no Hessian access, so it wins when Hessian-vector products are unavailable or expensive, when the problem is noisy, or when robustness and simplicity outweigh asymptotic rate. This is the standard division in large-scale smooth optimization and machine learning ^{[Nocedal, J. & Wright, S. J. — Numerical Optimization (2nd ed.)]}.

Synthesis. Large-scale unconstrained optimization is two answers to one constraint — do not store an $n\times n$ matrix — and the foundational reason both work is that a quasi-Newton or Newton step needs only the action of a curvature operator on a vector, never the operator itself. L-BFGS realises this by storing the last $m$ secant pairs and replaying them through the two-loop recursion, which computes the exact BFGS direction of 44.03.04 for an implicit $H_k={\gamma }_{k}I+UM{U}^{\top }$; the central insight is that limited memory changes the representation, not the mathematics, so hereditary positive-definiteness and the Wolfe-line-search globalisation of 44.03.01 carry over unchanged, and this is exactly the same self-perpetuating curvature condition $s_i^{\top }{y}_i>0$ that drove full BFGS. Truncated Newton realises the same matrix-free principle differently: it solves the Newton system of 44.03.03 approximately by CG 43.07.04 through Hessian-vector products of 44.03.08, and the Dembo-Eisenstat-Steihaug theorem reveals that the forcing sequence ${\eta }_k$ is dual to the convergence order, so inner accuracy and outer rate are two faces of one dial. Putting these together, L-BFGS is BFGS with a memory budget and truncated Newton is Newton with an accuracy budget, and the bridge between them is the conjugate-gradient method 43.07.04, at once the inner solver of truncated Newton and the limiting behaviour of BFGS, so the two families meet exactly where the Broyden-class conjugate-direction property of 44.03.04 said they would.

The bridge forward is explicit. The negative-curvature truncation that makes Newton-CG robust generalises to the trust-region setting 44.03.02, where the Steihaug-Toint inner CG follows its path to the trust-region boundary and the indefinite Hessian is handled by the subproblem rather than forbidden; the compact representation of L-BFGS is dual to the compact form of the direct approximation $B_k$ used inside limited-memory trust-region methods; and the cost dichotomy between fixed-budget L-BFGS and accuracy-budget truncated Newton is the same trade between memory and asymptotic rate that separates nonlinear conjugate gradient 44.03.06 from L-BFGS one rung lower. What full BFGS 44.03.04 does for medium problems, these two methods do for the largest — one by forgetting all but the recent curvature, the other by never forming the curvature at all.

Full proof set Master

Proposition 1 (two-loop recursion computes $H_k{g}_k$). Let $H_k$ be defined by $m$ BFGS inverse updates with pairs $(s_i,y_i)$, $s_i^{\top }{y}_i>0$, applied to $H_k^0={\gamma }_{k}I$. The two-loop recursion returns $r=H_k{g}_k$.

Proof. As in the Key theorem. The single-update action $H^{+}u=H(u-{\alpha }_i{y}_i)+({\alpha }_i-{\rho }_i{y}_i^{\top }H(u-{\alpha }_i{y}_i))s_i$ with ${\alpha }_i={\rho }_i{s}_i^{\top }u$ is precisely one nesting of the backward reduction $q\leftarrow q-{\alpha }_i{y}_i$ and the forward restoration $r\leftarrow r+({\alpha }_i-{\beta }_i)s_i$ with ${\beta }_i={\rho }_i{y}_i^{\top }r$. The backward loop visits $i=k-1,\dots ,k-m$ peeling off the $V_i$ factors from $g_k$; the scaling applies $H_k^0$; the forward loop visits $i=k-m,\dots ,k-1$ restoring the $V_i^{\top }$ factors and the rank-one terms in the reverse order, so the composed operation is exactly $H_k{g}_k$. The stored ${\alpha }_i$ are reusable because $s_i^{\top }q$ at the backward visit to index $i$ equals $s_i^{\top }{u}_i$ where $u_i$ is the argument seen by the $i$-th nested update. $\square$

Proposition 2 (positive-definiteness of the implicit $H_k$). If ${\gamma }_k>0$ and $s_i^{\top }{y}_i>0$ for all stored $i$, then $H_k\succ 0$ and $p_k=-H_k{g}_k$ is a descent direction.

Proof. $H_k^0={\gamma }_{k}I\succ 0$. Each BFGS inverse update preserves positive-definiteness under $s_i^{\top }{y}_i>0$ by hereditary positive-definiteness 44.03.04: for the update $H^{+}=V_i^{\top }H{V}_i+{\rho }_i{s}_i{s}_i^{\top }$ and any $z\ne 0$, $z^{\top }{H}^{+}z=(V_{i}z{)}^{\top }H(V_{i}z)+{\rho }_i(s_i^{\top }z{)}^2\ge 0$, vanishing only if $V_{i}z=0$ and $s_i^{\top }z=0$; but $V_{i}z=0$ means $z={\rho }_i{y}_i(s_i^{\top }z)=0$, a contradiction. Inducting over the $m$ updates gives $H_k\succ 0$, so $g_k^{\top }{p}_k=-g_k^{\top }{H}_k{g}_k<0$ for $g_k\ne 0$. $\square$

Proposition 3 (first-order accuracy of the finite-difference Hessian-vector product). For $f\in C^3$ near $x$, $\Vert (g(x+hv)-g(x))/h-{\nabla }^{2}f(x)v\Vert =O(h)$.

Proof. Taylor's theorem 02.05.05 applied to $g=\nabla f$ gives $g(x+hv)=g(x)+h{\nabla }^{2}f(x)v+\frac{h^2}{2}{D}^{3}f(\xi )[v,v]$ for some $\xi$ on the segment. Subtract $g(x)$, divide by $h$: the difference quotient equals ${\nabla }^{2}f(x)v+\frac{h}{2}{D}^{3}f(\xi )[v,v]$, and $\Vert D^{3}f\Vert$ bounded near $x$ makes the remainder $O(h)$. $\square$

Proposition 4 (inexact Newton local convergence). Under $x_k\to x^{\ast }$, $H^{\ast }={\nabla }^{2}f(x^{\ast })\succ 0$, ${\nabla }^{2}f$ Lipschitz near $x^{\ast }$, unit steps with $\Vert H_k{p}_k+g_k\Vert \le {\eta }_k\Vert g_k\Vert$: ${\eta }_k\le {\eta }_{\max }<1$ gives q-linear convergence in $\Vert \cdot {\Vert }_{\ast }=\Vert H^{\ast }\cdot \Vert$ with factor ${\eta }_{\max }$; ${\eta }_k\to 0$ gives q-superlinear; ${\eta }_k=O(\Vert g_k\Vert )$ gives q-quadratic.

Proof. Write $e_k=x_k-x^{\ast }$, $r_k=H_k{p}_k+g_k$. The exact Newton estimate under the Lipschitz Hessian gives $\Vert e_k-H_k^{-1}{g}_k{\Vert }_{\ast }\le C\Vert e_k{\Vert }_{\ast }^2$. The inexact step is $p_k=-H_k^{-1}(g_k-r_k)$, so $e_{k+1}=(e_k-H_k^{-1}{g}_k)+H_k^{-1}{r}_k$, whence

$$\Vert e_{k+1}{\Vert }_{\ast }\le C\Vert e_k{\Vert }_{\ast }^2+\Vert H^{\ast }{H}_k^{-1}{r}_k\Vert \le C\Vert e_k{\Vert }_{\ast }^2+(1+o(1))\Vert r_k\Vert .$$

With $\Vert r_k\Vert \le {\eta }_k\Vert g_k\Vert$ and $\Vert g_k\Vert =\Vert H_k^{\text{aste}}\Vert (1+o(1))=\Vert e_k{\Vert }_{\ast }(1+o(1))$, the linear term is ${\eta }_k\Vert e_k{\Vert }_{\ast }(1+o(1))$. Thus $\Vert e_{k+1}{\Vert }_{\ast }\le ({\eta }_k+o(1))\Vert e_k{\Vert }_{\ast }+C\Vert e_k{\Vert }_{\ast }^2$. For $x_0$ close enough the quadratic term is dominated and $\Vert e_{k+1}{\Vert }_{\ast }/\Vert e_k{\Vert }_{\ast }\le {\eta }_{\max }+o(1)<1$: q-linear with factor ${\eta }_{\max }$. If ${\eta }_k\to 0$ the ratio tends to $0$: q-superlinear. If ${\eta }_k\le K\Vert g_k\Vert =O(\Vert e_k{\Vert }_{\ast })$, the linear term becomes $O(\Vert e_k{\Vert }_{\ast }^2)$ and $\Vert e_{k+1}{\Vert }_{\ast }=O(\Vert e_k{\Vert }_{\ast }^2)$: q-quadratic. $\square$

Proposition 5 (Newton-CG returns a descent direction). CG on ${\nabla }^{2}f(x_k)p=-g_k$ from $p^{(0)}=0$, truncated at the first nonpositive-curvature direction (returning $-g_k$ if this is the first inner step), returns $p_k$ with $g_k^{\top }{p}_k<0$.

Proof. CG minimises $\varphi (p)=\frac{1}{2}{p}^{\top }{\nabla }^{2}f(x_k)p+g_k^{\top }p$ over the Krylov subspace generated, taking only positive-curvature steps before truncation. Since $\varphi (0)=0$ and each accepted step strictly decreases $\varphi$, the returned $p^{(j)}$ has $\varphi (p^{(j)})<0$, i.e. $\frac{1}{2}(p^{(j)}{)}^{\top }{\nabla }^{2}f(x_k)p^{(j)}+g_k^{\top }{p}^{(j)}<0$. The quadratic term is $\ge 0$ on the explored positive-curvature subspace, so $g_k^{\top }{p}^{(j)}<0$. If truncation occurs at $j=0$ the fallback $p_k=-g_k$ gives $g_k^{\top }{p}_k=-\Vert g_k{\Vert }^2<0$. $\square$

Proposition 6 (monotone norm increase of the CG path). While curvature stays positive, the CG iterates from $p^{(0)}=0$ satisfy $\Vert p^{(j+1)}\Vert >\Vert p^{(j)}\Vert$.

Proof. CG iterates expand as $p^{(j)}={\sum }_{i=0}^{j-1}{\xi }_i{d}_i$ over ${\nabla }^{2}f(x_k)$-conjugate directions $d_i$ with step coefficients ${\xi }_i=(r^{(i)\top }{r}^{(i)})/(d_i^{\top }{\nabla }^{2}f(x_k)d_i)>0$ (positive curvature). Conjugacy gives $d_i^{\top }{d}_l$ not necessarily zero, but the standard CG identity $\Vert p^{(j+1)}{\Vert }^2=\Vert p^{(j)}{\Vert }^2+2{\xi }_j{p}^{(j)\top }{d}_j+{\xi }_j^2\Vert d_j{\Vert }^2$ with $p^{(j)\top }{d}_j>0$ (each $p^{(j)}$ has nonnegative components along the running directions, provable by induction since ${\xi }_i>0$ and $d_j=r^{(j)}+{\tau }_j{d}_{j-1}$ with ${\tau }_j\ge 0$) makes the increment strictly positive. Hence $\Vert p^{(j)}\Vert$ strictly increases, so the path crosses any trust-region boundary exactly once, justifying the Steihaug-Toint truncation. $\square$

Connections Master

• Quasi-Newton methods (BFGS/DFP/SR1) 44.03.04 are the parent construction: L-BFGS is BFGS with the inverse-Hessian approximation stored as its last $m$ generating pairs rather than as a dense matrix, and the two-loop recursion computes the exact BFGS direction for the implicit $H_k={\gamma }_{k}I+UM{U}^{\top }$. Hereditary positive-definiteness, the curvature condition $s_i^{\top }{y}_i>0$, and the secant interpretation all transfer verbatim; what changes is only that finite memory degrades the asymptotic rate from superlinear to R-linear, recovered as $m\to \infty$. The Broyden-class conjugate-direction property of that unit is the structural reason L-BFGS behaves like a preconditioned nonlinear conjugate gradient on large problems.

• The linear conjugate gradient method 43.07.04 is the inner engine of truncated Newton: CG solves the symmetric positive-definite Newton system ${\nabla }^{2}f(x_k)p=-g_k$ using only matrix-vector products, terminated early by the forcing-sequence residual test. The same CG, run on a quadratic with exact line search, coincides with full BFGS 44.03.04, so conjugate gradient is the common ancestor of both large-scale families — the inner solver of inexact Newton and the limiting behaviour of the limited-memory secant method. The negative-curvature breakdown of CG on indefinite systems is exactly what forces the Newton-CG truncation.

• Finite-difference and automatic differentiation 44.03.08 supply the Hessian-vector products ${\nabla }^{2}f(x)v$ that make truncated Newton Hessian-free: automatic differentiation differentiates the gradient to give the product to working precision, while the finite difference $(g(x+hv)-g(x))/h$ gives it to $O(h)$ accuracy with a cancellation trade-off optimised at $h\approx \sqrt{{\epsilon }_{\mathrm{mach}}}$. Without cheap Hessian-vector products truncated Newton loses its advantage, which is precisely the regime where L-BFGS — needing no Hessian access — is preferred.

• Newton's method 44.03.03, trust-region methods 44.03.02, and nonlinear conjugate gradient 44.03.06 complete the network: truncated Newton is the inexact, matrix-free realisation of Newton's method 44.03.03 with the forcing sequence dialing the convergence order; the negative-curvature truncation generalises to the Steihaug-Toint trust-region CG of 44.03.02, where the indefinite Hessian is absorbed by the subproblem; and L-BFGS sits one rung above nonlinear conjugate gradient 44.03.06 in the memory hierarchy, carrying $m$ directions of curvature memory against CG's single direction, the same memory-versus-rate trade made one level down.

Historical & philosophical context Master

The limited-memory idea is due to Jorge Nocedal, whose 1980 paper showed that storing only the last few BFGS correction pairs and replaying them through a recursion preserves most of the method's speed at a fraction of the storage ^{[Nocedal 1980]}. The method was placed on a firm convergence footing by Dong C. Liu and Jorge Nocedal in 1989, who proved R-linear convergence on strongly convex problems and demonstrated that small memory ($m$ between 3 and 20) is competitive with full BFGS on large problems while scaling to dimensions where a dense matrix is impossible ^{[Liu-Nocedal 1989]}. The two-loop recursion they popularised has become the standard direction computation in large-scale optimization libraries.

The inexact-Newton framework was formalised by Ron Dembo, Stanley Eisenstat, and Trond Steihaug in 1982, who introduced the forcing sequence ${\eta }_k$ and proved the local-convergence dichotomy relating it to linear, superlinear, and quadratic rates ^{[Dembo-Eisenstat-Steihaug 1982]}. Trond Steihaug in 1983 and, independently, Philippe Toint combined conjugate gradient with trust regions, giving the negative-curvature-tolerant inner solver now called Steihaug-Toint CG ^{[Steihaug 1983]}. The practical forcing-sequence choice ${\eta }_k=\min \{\frac{1}{2},\sqrt{\Vert \nabla f_k\Vert }\}$ that secures superlinear convergence at controlled cost is due to Stanley Eisenstat and Homer Walker in 1996 ^{[Eisenstat-Walker 1996]}. The convergence theory of inexact Newton descends from the classical analysis of Newton's method for operator equations by L. V. Kantorovich, whose 1948 fixed-point estimates underlie the quadratic-rate baseline these methods relax.

Bibliography Master

@article{nocedal1980lbfgs,
  author  = {Nocedal, Jorge},
  title   = {Updating quasi-Newton matrices with limited storage},
  journal = {Mathematics of Computation},
  volume  = {35},
  number  = {151},
  year    = {1980},
  pages   = {773--782}
}

@article{liunocedal1989,
  author  = {Liu, Dong C. and Nocedal, Jorge},
  title   = {On the limited memory BFGS method for large scale optimization},
  journal = {Mathematical Programming},
  volume  = {45},
  number  = {1--3},
  year    = {1989},
  pages   = {503--528}
}

@article{dembo1982inexact,
  author  = {Dembo, Ron S. and Eisenstat, Stanley C. and Steihaug, Trond},
  title   = {Inexact {N}ewton methods},
  journal = {SIAM Journal on Numerical Analysis},
  volume  = {19},
  number  = {2},
  year    = {1982},
  pages   = {400--408}
}

@article{steihaug1983,
  author  = {Steihaug, Trond},
  title   = {The conjugate gradient method and trust regions in large scale optimization},
  journal = {SIAM Journal on Numerical Analysis},
  volume  = {20},
  number  = {3},
  year    = {1983},
  pages   = {626--637}
}

@article{eisenstatwalker1996,
  author  = {Eisenstat, Stanley C. and Walker, Homer F.},
  title   = {Choosing the forcing terms in an inexact {N}ewton method},
  journal = {SIAM Journal on Scientific Computing},
  volume  = {17},
  number  = {1},
  year    = {1996},
  pages   = {16--32}
}

@article{byrd1994compact,
  author  = {Byrd, Richard H. and Nocedal, Jorge and Schnabel, Robert B.},
  title   = {Representations of quasi-Newton matrices and their use in limited memory methods},
  journal = {Mathematical Programming},
  volume  = {63},
  number  = {1--3},
  year    = {1994},
  pages   = {129--156}
}

@book{nocedalwright2006ln,
  author    = {Nocedal, Jorge and Wright, Stephen J.},
  title     = {Numerical Optimization},
  edition   = {2},
  series    = {Springer Series in Operations Research and Financial Engineering},
  publisher = {Springer},
  year      = {2006}
}