BFGS & Quasi-Newton Methods

You want to find the minimum of a function. Gradient descent is the obvious first tool: measure the slope, step downhill, repeat. But the slope alone is a blunt instrument. If the landscape curves steeply in one direction and gently in another, a step that works in one direction overshoots in the other. You waste iterations zigzagging when you could be gliding.

The fix is to model the curvature — the matrix of second derivatives known as the Hessian. If you know the Hessian, you can correct each step so that it lands precisely at the bottom of the local quadratic bowl. This is Newton's method, and it converges in very few iterations.

The catch is the Hessian itself. For a function of $n$ variables it is an $n \times n$ matrix, which costs $O(n^{2})$ to store and $O(n^{3})$ to invert on every step — completely impractical for modern machine learning where $n$ might be in the millions.

Quasi-Newton methods split the difference. Instead of computing the exact Hessian, they estimate it from the gradient differences that accumulate as the algorithm moves. They ask: given that the gradient changed by $\mathbf{y}_{k}$ when the position changed by $\mathbf{s}_{k}$, what does that tell us about curvature? The answer shapes a rank-2 update to the current Hessian approximation.

The most successful version of this idea is BFGS, named after Broyden, Fletcher, Goldfarb and Shanno, who independently derived it in 1970. It is the default optimizer in scipy, the backbone of many statistical fitting routines, and the intellectual parent of L-BFGS — the optimizer behind much of modern deep learning.

Below, both gradient descent (blue) and BFGS (orange) start at the same point on a curved 2-D surface and race to the minimum. The surface is deliberately ill-conditioned — stretched much further in one direction than the other — so gradient descent zigzags while BFGS adapts.

<!-- {{c_html_intro}} -->
<div class="controls">
  <button id="btn-step" type="button">{{btn_step}}</button>
  <button id="btn-run"  type="button">{{btn_run}}</button>
  <button id="btn-reset" type="button" class="ghost">{{btn_reset}}</button>
  <span class="legend"><span class="dot gd"></span> {{label_gd}} &nbsp;<span class="dot bfgs"></span> {{label_bfgs}}</span>
</div>
<canvas id="cv" width="420" height="300" title="{{canvas_title}}"></canvas>
<div id="info" class="info">{{status_ready}}</div>

/* {{c_css_layout}} */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; color: #222; background: transparent; }
.controls { display: flex; align-items: center; gap: .5rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 13px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .45; cursor: default; }
.legend { font-size: .82rem; display: flex; align-items: center; gap: .3rem; margin-left: auto; }
.dot { display: inline-block; width: 10px; height: 10px; border-radius: 50%; }
.dot.gd   { background: #3a86ff; }
.dot.bfgs { background: #fb8500; }
canvas { display: block; border: 1px solid #dde; border-radius: 8px; max-width: 100%; }
.info { font-size: .88rem; font-weight: 600; margin-top: .4rem; min-height: 1.3em; }
.info.done-gd   { color: #3a86ff; }
.info.done-bfgs { color: #e76f00; }
.info.done-both { color: #0a7d33; }

// Code not found

Press Step to advance one iteration at a time, or Run to watch the full race. Notice how BFGS builds a better model of the local shape after just a few steps and then takes a nearly direct path, while gradient descent keeps bouncing off the steep walls. The step count shown in the corner makes the difference vivid.

What exactly does BFGS do on each step, and why does it work?

The secant condition. After moving from $\mathbf{x}_{k}$ to $\mathbf{x}_{k+1}$, we have two gradient vectors. Their difference $\mathbf{y}_{k} = \nabla f(\mathbf{x}_{k+1}) - \nabla f(\mathbf{x}_{k})$ and the step $\mathbf{s}_{k} = \mathbf{x}_{k+1} - \mathbf{x}_{k}$ must satisfy the secant condition: the new Hessian approximation $B_{k+1}$ should predict the observed gradient change, i.e. $B_{k+1} \mathbf{s}_{k} = \mathbf{y}_{k}$.

The minimal update. There are infinitely many matrices satisfying the secant condition. BFGS picks the one closest to $B_{k}$ in a natural matrix norm — a symmetric rank-2 correction:

$$B_{k+1} = B_{k} - \frac{B_{k}\mathbf{s}_{k}\mathbf{s}_{k}^{\top}B_{k}}{\mathbf{s}_{k}^{\top}B_{k}\mathbf{s}_{k}} + \frac{\mathbf{y}_{k}\mathbf{y}_{k}^{\top}}{\mathbf{y}_{k}^{\top}\mathbf{s}_{k}}$$

If $B_{k}$ is positive-definite and the curvature condition $\mathbf{y}_{k}^{\top}\mathbf{s}_{k} > 0$ holds (guaranteed by a Wolfe-condition line search), then $B_{k+1}$ stays positive-definite — meaning the search direction always points downhill.

Convergence. On smooth, strongly convex functions, BFGS achieves superlinear convergence: the error shrinks faster than any geometric rate, approaching quadratic convergence near the minimum. Gradient descent, by contrast, converges only linearly, at a rate governed by the condition number $\kappa$ of the Hessian.

L-BFGS. For large-scale problems, storing the full $n \times n$ matrix is too costly. Limited-memory BFGS stores only the last $m$ pairs $(\mathbf{s}_{k}, \mathbf{y}_{k})$ — typically $m \in \{5, \dots, 20\}$ — and reconstructs the product $B_{k}^{-1}\mathbf{g}$ implicitly in $O(mn)$ time. This is the optimizer used to train many large language models and neural networks, and it is also related to algorithms used in non-convex optimization.

Quasi-Newton methods are the workhorse of smooth, medium-to-large-scale optimization:

• Machine learning: L-BFGS trains logistic regression, SVMs, and small-to-medium neural networks. PyTorch's torch.optim.LBFGS is a standard option for full-batch settings. Some large language model fine-tuning recipes use L-BFGS variants for their fast convergence.
• Statistical fitting: maximum-likelihood estimation in mixed models, survival analysis, and structural equation models almost always uses a quasi-Newton solver internally — the user types fit() and BFGS is doing the work.
• Scientific computing: scipy.optimize.minimize defaults to L-BFGS-B (a bound-constrained variant) for unconstrained smooth problems. Molecular dynamics force-field minimization, geophysical inversion, and optimal control all rely on it.
• Engineering design: aerodynamic shape optimization, antenna design, and PDE-constrained optimization use adjoint methods to compute gradients cheaply, then feed those gradients to L-BFGS.

The key insight connecting all of these is the same: whenever you can compute gradients but not second derivatives, and the problem is smooth enough for a local quadratic model to be useful, quasi-Newton is the right tool. See also gradient descent's failure modes for the cases where it breaks down.

The core insight of BFGS is almost disarmingly simple: every time you take a step, the change in gradient tells you something about curvature. Collect those clues, keep the cheapest update that respects them, and you get a Hessian approximation that improves with every iteration — without ever computing a single second derivative.

The result is an algorithm that achieves superlinear convergence on smooth problems, costs only $O(n^{2})$ per step in its full form, and shrinks to $O(mn)$ per step in L-BFGS. It is the default optimizer in scipy, the go-to for statistical fitting, and the backbone of neural network training when full-batch optimization is viable.

The next time a model trains in seconds instead of hours, there is a good chance a rank-2 matrix update is quietly doing the heavy lifting behind the scenes.