Second-Order & Quasi-Newton Methods
Newton's method, the Hessian, IRLS, and BFGS / L-BFGS.
On this page
Beginner: plain gradient descent only ever looks at the slope under its feet — the first derivative — and takes a step downhill proportional to that slope, scaled by a learning rate you have to pick. It has no idea whether the ground ahead is about to flatten out or steepen up; it finds out only after it gets there. Newton's method additionally looks at the curvature — the second derivative, or the Hessian matrix in more than one dimension — which tells it how the slope itself is changing. With both slope and curvature in hand, Newton's method can fit a parabola (a quadratic bowl) that locally matches the real function at the current point, and then jump straight to the bottom of that parabola in a single step, instead of nibbling downhill one small step at a time. Near the true minimum, where most smooth functions really do look almost exactly like a bowl, this is often dramatically faster than gradient descent.
Intermediate: the reason this works so well close to the optimum is a direct consequence of calculus, not a lucky coincidence — a smooth function's second-order Taylor expansion around its minimum is a quadratic bowl, up to an error term that shrinks faster than the quadratic term itself as you get closer. So "solve the quadratic approximation exactly" and "solve the real problem exactly" converge to being almost the same task once you're near the answer. The price is steep, though. Forming the Hessian for n parameters means computing an n×n matrix of second derivatives, and inverting (or solving a linear system with) that matrix costs roughly O(n\u00b3) — trivial for n in the tens, ruinous for a model with millions of parameters. Far from the optimum it's also risky in a specific, structural way: the local quadratic approximation can simply be a bad model of the real function out there, and Newton's raw update has no built-in notion of "downhill" — it only solves for where the local quadratic model is flat. Near a saddle point (common in non-convex deep-learning losses, where the Hessian has both positive and negative curvature directions), that stationary point of the local model can be a maximum or a saddle rather than a minimum, and the raw update can march the parameters straight toward it.
Advanced: quasi-Newton methods — most importantly BFGS and its limited-memory cousin L-BFGS — are the practical compromise the field actually uses. Instead of computing the true Hessian at all, they build up an approximation to the (inverse) Hessian purely from the gradients already computed at recent steps — the same information gradient descent already has lying around, just used more cleverly. Each step's change in position and change in gradient gives one more piece of curvature information "for free," and BFGS folds that piece into a running approximation that gets more curvature-aware over time, recovering much of Newton's fast convergence at a small fraction of the cost. L-BFGS pushes this further: rather than storing a full dense n×n approximate-Hessian matrix (which is itself O(n\u00b2) in memory — still too much for very large models), it keeps only a small, fixed-size window of the last handful of position and gradient differences, and reconstructs the effect of the approximate inverse Hessian from just those on the fly. That is exactly what "limited-memory" refers to, and it's what makes L-BFGS usable on problems full BFGS, let alone full Newton, never could touch.
Gradient descent's update rule:
Newton's method's update rule:
Line these up and the relationship is exact: Newton's method is gradient descent with the scalar learning rate replaced by an entire matrix, , the inverse Hessian. A scalar α can only ever stretch or shrink the gradient step by the same amount in every direction at once. The inverse Hessian instead rescales and rotates the step — taking large steps along directions the function is nearly flat in, small steps along directions it's steeply curved in, and correcting for how those directions interact — using real curvature information instead of one number picked in advance. That single substitution, one matrix in place of one scalar, is the entire formal difference between a first-order and a second-order optimizer.
Start from the second-order Taylor expansion of f around the current point , for a small step :
This right-hand side is a plain quadratic function of δ — a constant, plus a linear term, plus a quadratic term — so it can be minimized over δ exactly, the same way any quadratic can. Differentiate it with respect to δ. The linear term's gradient is just ∇f(θ_k); the quadratic term's gradient uses the standard identity for a symmetric matrix H, (the Hessian ∇²f(θ_k) is symmetric because mixed partial derivatives commute for any reasonably smooth f). So:
Setting this to zero — the first-order condition for a minimum of the quadratic approximation — gives a linear equation in δ:
Solving for δ (left-multiplying by the inverse Hessian, assuming it exists):
and taking that step, θ_{k+1} = θ_k + δ, is exactly Newton's update from the Formula section above. Every step of Newton's method is, precisely, "minimize the local quadratic model exactly," repeated.
Why this converges quadratically near a well-behaved minimum. The Taylor expansion used above is only an approximation — it drops a remainder term that is O(\u2016\u03b4\u2016\u00b3) for the function itself, which means the gradient's own linear approximation, ∇f(θ_k) + ∇²f(θ_k)δ, differs from the true gradient ∇f(θ_k + δ) by only O(‖δ‖²) (one derivative lower than the function's own remainder). Newton's step δ is chosen precisely to make that linear approximation of the gradient exactly zero — so the true gradient at the new point, ∇f(θ_{k+1}), is whatever is left over after that approximation error: O(‖δ‖²), i.e. of order the square of the distance just moved. Near a well-behaved minimum, where the Hessian stays roughly constant and positive definite, being O(‖δ‖²) away in gradient translates into being O(‖δ‖²) away in position too. In plain terms: each Newton step's error is roughly the square of the previous step's error, so the number of correct digits roughly doubles every iteration — a handful of steps can go from one correct digit to effectively machine precision. Gradient descent, by contrast, only ever multiplies its error by some fixed factor less than one each step (linear convergence) — a fixed fraction of progress per step, not a squaring of the remaining error, which is exactly why it typically needs far more iterations to reach the same precision.
Connection to IRLS. Apply this exact derivation to logistic regression's log-likelihood, and the resulting Newton update — worked out fully in the later Regression module once generalized linear models (GLMs) are introduced — turns out to be exactly the classical Iteratively Reweighted Least Squares (IRLS) algorithm: at each iteration, solve a weighted least-squares problem where the weights come from the current model's predicted probabilities. IRLS isn't a different algorithm from Newton's method wearing a different name by coincidence — it is Newton's method, specialized to this one log-likelihood, with its particular Hessian recognized as reducible to a weighted normal-equations solve.
Where this is used: this is precisely why logistic regression is so often fit with a Newton-type method (IRLS, or Newton-Raphson directly on the log-likelihood) rather than plain gradient descent — the underlying optimization problem is convex (lesson 2.2.1), so there's no saddle-point or non-convexity risk to worry about, and Newton's quadratic convergence typically gets a tight fit in a handful of iterations rather than the hundreds or thousands gradient descent might need. It's also exactly why L-BFGS is a common default solver for the maximum-likelihood fitting step underneath many classical statistical and machine learning models — it gets most of that same fast convergence without ever forming a Hessian.
The race replays automatically once on load. Drag the grey marker to a new point on the curve and release it to re-run the race from there — Newton's method (violet) reaches the bottom in a handful of large, curvature-aware jumps, while gradient descent (blue) crawls there in many small, fixed-size steps.
All three implementations below minimize the same one-dimensional loss f(x) = 0.5x\u00b2 + 0.05x\u2074 used in the diagram above, starting from x = 3. The first two hand-roll both Newton's method and gradient descent from an exact, manually-derived gradient and Hessian, printing the error after each step so the quadratic-vs-linear convergence-rate contrast is visible directly in the numbers. The third swaps in a real, production-grade optimizer library and compares a genuine second-order method against a genuine quasi-Newton one.
# Minimizing f(x) = 0.5*x^2 + 0.05*x^4, a well-behaved (convex, non-quadratic)
# 1D loss, comparing Newton's method against plain gradient descent -- both
# implemented with an exact, hand-derived gradient and Hessian, no autodiff.
def f(x):
return 0.5 * x**2 + 0.05 * x**4
def grad(x):
return x + 0.2 * x**3
def hess(x):
return 1 + 0.6 * x**2
x_true_min = 0.0
def newtons_method(x0, steps=6):
x = x0
print("Newton's method")
print(f" step 0: x={x:.6f} error={abs(x - x_true_min):.6f}")
for k in range(1, steps + 1):
x = x - grad(x) / hess(x)
print(f" step {k}: x={x:.6f} error={abs(x - x_true_min):.6f}")
return x
def gradient_descent(x0, lr=0.15, steps=40, print_every=5):
x = x0
print("Gradient descent")
print(f" step 0: x={x:.6f} error={abs(x - x_true_min):.6f}")
for k in range(1, steps + 1):
x = x - lr * grad(x)
if k % print_every == 0 or k == steps:
print(f" step {k}: x={x:.6f} error={abs(x - x_true_min):.6f}")
return x
newtons_method(3.0)
print()
gradient_descent(3.0)
# Typical output shows Newton's error shrinking roughly like it is being squared
# each step once it gets close (e.g. an error of ~0.11 becomes an error under
# 0.001 on the very next step), reaching machine-precision-level accuracy in
# about 4-5 steps -- while gradient descent's error shrinks by only a roughly
# constant fraction each step, and is still not fully converged after 40 steps.- Logistic regression fit via Newton's method / IRLS. Classical statistics packages (R's
glm(), Python'sstatsmodels) default to Newton-Raphson or IRLS for fitting logistic regression and other GLMs, because the log-likelihood is convex (2.2.1) and smooth enough that a handful of Newton iterations reaches a tight fit — exactly the scenario the Derivation section above builds toward. - L-BFGS as scikit-learn's default
LogisticRegressionsolver.solver='lbfgs'is the library's default precisely because it needs only gradients (cheap, and easy to get right for a convex loss), yet still converges in far fewer iterations than plain gradient descent would — the practical sweet spot this whole lesson is about, and it shows up as the literal default argument in one of the most-used ML libraries in existence. - Newton's method is essentially never used to train deep neural networks directly. A network with tens of millions or billions of parameters has a Hessian with that many squared entries — computing it, let alone inverting it, is completely infeasible at that scale (
O(n\u00b3)for the solve, on top ofO(n\u00b2)just to store it). This is exactly why the adaptive, per-parameter first-order optimizers from the previous lesson (2.2.3) — Adam, RMSProp, and their relatives — dominate deep learning instead: they approximate a little bit of curvature information (a per-parameter scale) far more cheaply than any true second-order method could. - Gauss-Newton and Levenberg-Marquardt are second-order-flavored methods purpose-built for nonlinear least-squares curve fitting — fitting a nonlinear model to data by minimizing squared residuals. They approximate the Hessian using only first-derivative (Jacobian) information specific to a sum-of-squares objective, which is cheaper than a true Hessian and numerically well-behaved for this one common problem shape. This is what runs, for example, behind
scipy.optimize.curve_fitand countless calibration and system-identification routines in engineering. - Natural gradient and K-FAC-style optimizers are attempts to bring curvature information back into large-scale training cheaply, by approximating the Hessian's structure (e.g. per-layer, or via the Fisher information matrix) rather than forming it exactly. They're mentioned here only as a forward pointer — genuinely promising, but still squarely in the research-and-specialized-tooling category rather than default practice (see the Expert note below).
- L-BFGS shows up throughout the scientific Python stack well beyond logistic regression — as the default or a standard option for maximum-likelihood fits of Gaussian process kernel hyperparameters, conditional random fields, and many other
scipy.optimize.minimizecall sites where gradients are cheap to compute but an exact Hessian is not.
- Running plain Newton's method far from the minimum on a non-convex loss and being surprised when it diverges or marches toward a saddle point or a maximum. The raw update has no concept of "downhill" — it solves for wherever the local quadratic approximation is stationary, which is a genuine minimum only if the Hessian there is positive definite. Near a saddle point (common in non-convex deep-learning losses, per the Beginner/Intermediate discussion above), the Hessian has negative eigenvalues, and Newton's step can happily move toward that saddle rather than away from it.
- Assuming "more curvature information is always better" while ignoring the very real
O(n\u00b3)-per-step cost of forming and solving with an exactn×nHessian. For a model with a few dozen parameters that's free; for a model with millions, it's a categorically bad trade — the previous lesson's adaptive first-order methods exist precisely because paying for exact curvature at that scale isn't worth it. - Confusing a quasi-Newton method's approximate curvature — built up purely from recent gradient differences — with the true Hessian. BFGS and L-BFGS's approximation is often excellent in practice, but its theoretical guarantees are weaker and more empirical than full Newton's method's guarantees near a true minimum; treating an L-BFGS run as "basically exact second-order information" overstates what it's actually doing.
Going deeper
Trust region methods are a more robust refinement of the basic Newton idea, aimed directly at the far-from-the-minimum failure mode in the Pitfall above. Instead of always taking the full Newton step wherever it points, a trust-region method only trusts its local quadratic model within some radius, solves a constrained version of the Newton problem inside that radius, and then checks how well the model's predicted decrease in the objective matched the actual decrease once the step was taken. A good match grows the trust radius for next time; a bad match shrinks it and the step is retried more cautiously. This single feedback loop — grow trust when the local model is being honest, shrink it when it isn't — is what turns "solve the local quadratic exactly, always" into something that degrades gracefully far from the optimum instead of diverging outright.
In modern deep learning, bringing curvature information back into large-scale training cheaply is an active area of ongoing research, not settled practice — K-FAC (Kronecker-Factored Approximate Curvature) approximates each layer's Fisher information matrix with a much cheaper Kronecker-factored structure, and Shampoo uses a related factored-preconditioner idea, both trying to recover some of second-order convergence's benefits without ever forming a true Hessian. Both have shown real gains in specific large-scale training setups, but neither has displaced Adam-family optimizers as the default choice the way L-BFGS has for classical convex fitting — treat this as a promising research direction to watch, not a settled recommendation.
On a non-convex loss, why can plain Newton's method actually move toward a saddle point instead of away from it, when gradient descent with a small step size tends not to?
Gradient descent always moves in the direction of steepest local decrease -- it only ever asks 'which way is downhill right now,' so on any direction where the loss is currently decreasing it keeps decreasing along that direction until the slope there flattens out. Newton's method asks a different question: 'where is the local quadratic model of the loss stationary' -- that is, where does the linear approximation of the gradient hit zero. That stationary point is a true minimum of the real function only when the Hessian at the current point is positive definite in every direction. Near a saddle point the Hessian has at least one negative eigenvalue, so the quadratic model itself curves downward in that direction and its stationary point sits on the far side of the saddle rather than at a minimum -- and Newton's update, having no separate notion of 'downhill' beyond solving for that stationary point, will step toward it anyway. Newton's method is only safe to run unmodified when convexity (a positive-definite Hessian, as guaranteed in the convex problems from 2.2.1) rules this failure mode out.
Newton's method replaces gradient descent's one-size-fits-all scalar step with a full curvature-aware matrix step, derived by exactly minimizing a local quadratic model of the loss — which is why it converges quadratically near a well-behaved minimum, but is O(n\u00b3)-per-step expensive and can misbehave badly (moving toward a saddle or maximum) far from one or on non-convex problems. Quasi-Newton methods, especially L-BFGS, are the practical compromise the field actually reaches for: approximate curvature built cheaply from recent gradients, getting most of Newton's speed without ever forming a Hessian — which is exactly why Newton-type methods and IRLS dominate convex fits like logistic regression, L-BFGS is scikit-learn's default logistic-regression solver, and true second-order methods are essentially absent from deep-learning training, where the previous lesson's adaptive first-order optimizers take over instead.