Prev
Unit 11 · Numerical Optimisation
Newton and Quasi-Newton Methods
60m100 XP
Next
⛰️Numerical Optimisation

Newton and Quasi-Newton Methods

Newton's method for minimisation, BFGS, scipy.optimize.minimize and derivative-free Nelder–Mead

Argonne, 1959
William Davidon is a physicist at the Argonne National Laboratory near Chicago, using the laboratory's computer to fit theoretical models to experimental data — minimising functions of several variables. His runs with the available methods are so slow that, by his own account, the computer would often stop for maintenance before they finished. He finds a way to make each step smarter. Newton's method would use the Hessian matrix of second derivatives, which he cannot afford to compute. But every step taken reveals something about curvature: the change in gradient divided by the change in position. Davidon accumulates that information into a matrix that approximates the inverse Hessian, updating it after every step — a variable metric. His technical report, ANL-5990, is rejected by the journal he submits it to. Roger Fletcher and Michael Powell read it, analyse it, and publish an improved version in 1963 — the DFP method — and it transforms practical optimisation. In 1970 four researchers independently publish the refinement still in use, BFGS. Davidon's original paper finally appears in 1991, as the opening article of the first issue of the SIAM Journal on Optimization.
💡Newton's method for minimisation
Approximate ff near xk\mathbf x_k by its second-order Taylor model and minimise the model: ∇2f(xk) d=−∇f(xk),xk+1=xk+t d.\nabla^2f(\mathbf x_k)\,\mathbf d = -\nabla f(\mathbf x_k), \qquad \mathbf x_{k+1} = \mathbf x_k + t\,\mathbf d. This is Newton's method (lesson 2-4) applied to ∇f=0\nabla f = \mathbf 0. Near a minimum with positive definite Hessian it converges quadratically and is unaffected by bad scaling. Far away the Hessian may not be positive definite, so the step may not point downhill: practical codes fall back to the gradient or modify the Hessian, and use a line search. Damped Newton solves Rosenbrock from (−1.2,1)(-1.2, 1) in 22 iterations — against steepest descent's 10,916.
Derivation
The Newton step minimises a quadratic model
Keep the Taylor expansion of lesson 11-1 to second order, as a model of ff near xk\mathbf x_k:
m(d)=f(xk)+∇fTd+12dT∇2f d.m(\mathbf d) = f(\mathbf x_k) + \nabla f^{\mathsf T}\mathbf d + \tfrac12\mathbf d^{\mathsf T}\nabla^2f\,\mathbf d.
Its gradient with respect to d\mathbf d is ∇f+∇2f d\nabla f + \nabla^2f\,\mathbf d (as in lesson 4-3 for the energy). Setting it to zero gives the Newton step:
∇2f d=−∇f.\nabla^2f\,\mathbf d = -\nabla f.
If the Hessian is positive definite, the model is a bowl and d\mathbf d is its bottom; then ∇fTd=−dT∇2f d<0\nabla f^{\mathsf T}\mathbf d = -\mathbf d^{\mathsf T}\nabla^2f\,\mathbf d < 0, so d\mathbf d points downhill. If the Hessian is indefinite the model is a saddle, its 'bottom' may lie uphill — the failure the lab guards against by falling back to −∇f-\nabla f. For a quadratic ff the model is exact and one step lands on the minimum, whatever the condition number: Newton is immune to bad scaling.
💡BFGS
Quasi-Newton methods replace the Hessian by an approximation updated from gradients. With step s=xk+1−xk\mathbf s = \mathbf x_{k+1} - \mathbf x_k and gradient change y=∇fk+1−∇fk\mathbf y = \nabla f_{k+1} - \nabla f_k, require the new inverse-Hessian approximation to satisfy the secant equation Hk+1y=sH_{k+1}\mathbf y = \mathbf s — the same idea as Broyden's method (lesson 2-4) — while staying symmetric, positive definite and close to HkH_k. The BFGS update (Broyden, Fletcher, Goldfarb and Shanno, 1970): Hk+1=(I−ρsyT)Hk(I−ρysT)+ρ ssT,ρ=1/yTs.H_{k+1} = (I - \rho\mathbf s\mathbf y^{\mathsf T})H_k(I - \rho\mathbf y\mathbf s^{\mathsf T}) + \rho\,\mathbf s\mathbf s^{\mathsf T}, \qquad \rho = 1/\mathbf y^{\mathsf T}\mathbf s. Only gradients are needed; convergence is superlinear; each step costs O(n2)O(n^2). L-BFGS stores only the last few (s,y)(\mathbf s, \mathbf y) pairs and scales to millions of variables. BFGS is the default of scipy.optimize.minimize for unconstrained smooth problems.
Derivation
Checking the BFGS update
Where the secant equation comes from. For a quadratic ff the gradient changes linearly, ∇f(x+s)−∇f(x)=∇2f s\nabla f(\mathbf x + \mathbf s) - \nabla f(\mathbf x) = \nabla^2f\,\mathbf s, so y=∇2f s\mathbf y = \nabla^2f\,\mathbf s and the inverse Hessian maps y\mathbf y back to s\mathbf s. Any good approximation HH of the inverse Hessian should do the same along the step just taken: Hk+1y=sH_{k+1}\mathbf y = \mathbf s.
The update satisfies it. With ρ=1/yTs\rho = 1/\mathbf y^{\mathsf T}\mathbf s, the right-hand factor annihilates y\mathbf y:
(I−ρ ysT) y=y−ρ y (sTy)=y−y=0.(I - \rho\,\mathbf y\mathbf s^{\mathsf T})\,\mathbf y = \mathbf y - \rho\,\mathbf y\,(\mathbf s^{\mathsf T}\mathbf y) = \mathbf y - \mathbf y = \mathbf 0.
So the first term of Hk+1yH_{k+1}\mathbf y vanishes, and the second gives ρ s(sTy)=s\rho\,\mathbf s(\mathbf s^{\mathsf T}\mathbf y) = \mathbf s. The secant equation holds exactly.
It stays positive definite. The update is symmetric by construction, and it can be shown to keep HH positive definite whenever yTs>0\mathbf y^{\mathsf T}\mathbf s > 0 — the curvature along the step is positive, which a line search satisfying the Wolfe conditions guarantees. That is what keeps every BFGS direction −H∇f-H\nabla f pointing downhill.
📐Comparing methods on Rosenbrock
From (−1.2,1)(-1.2, 1) to tolerance 10−1010^{-10} with scipy.optimize.minimize:
  • Nelder–Mead (no derivatives): 132 iterations, 249 function evaluations;
  • conjugate gradient (nonlinear CG, gradients): 37 iterations, 80 evaluations;
  • BFGS: 34 iterations, 41 function and 41 gradient evaluations;
  • Newton-CG (gradient and Hessian): 85 iterations.
In ten dimensions the gap opens: Nelder–Mead needs 6,817 function evaluations, BFGS 112. Derivative information becomes more valuable the more variables there are.
If you cannot solve the proposed problem, try to solve first a simpler related problem.
– George Pólya, Mathematical Methods in Science (1977)
💡Nelder–Mead: no derivatives at all
Sometimes the objective is the output of a simulation, measured data or a black box with no usable derivative. John Nelder and Roger Mead's 1965 simplex method (not Dantzig's) keeps n+1n + 1 points — a triangle in 2-D — and replaces the worst by reflecting it through the others, expanding if that works well, contracting if not, shrinking as a last resort. It is robust and easy to use, but slow and without convergence guarantees in higher dimensions. Tuning a PI controller for a heater with dead time by minimising the integral of time-weighted absolute error over a simulated step response, Nelder–Mead takes 92 simulations to reach an ITAE of 11.49 — better than both the starting guess (79.96) and a Ziegler–Nichols-style rule (34.56).
💡Exact gradients: complex step and automatic differentiation
Finite differences (lesson 8-1) lose about half the digits: for f(x)=ex/sin⁡3x+cos⁡3xf(x) = e^x/\sqrt{\sin^3x + \cos^3x} at x=1.5x = 1.5 the best forward difference (h=10−8h = 10^{-8}) has a relative error of 2×10−82 \times 10^{-8}, and smaller hh makes it worse. Two ideas remove the subtraction altogether. The complex step (Lyness and Moler, 1967; revived by Squire and Trapp, 1998): for a real analytic ff, f′(x)≈Im⁡f(x+ih)/hf'(x) \approx \operatorname{Im} f(x + ih)/h. There is no difference of nearly equal numbers, so h=10−20h = 10^{-20} is fine — and the error is zero to machine precision. It needs code that accepts complex numbers. Automatic differentiation applies the chain rule to the program itself. Forward mode carries each intermediate value together with its derivative through every operation; it gives one directional derivative per pass. Reverse mode (Seppo Linnainmaa, 1970; known in machine learning as backpropagation) records the operations, then sweeps backwards through them, producing the whole gradient of a scalar function for a small multiple of the cost of evaluating it — whether it has 2 variables or 2 million. Finite differences would need n+1n + 1 evaluations. JAX, PyTorch and CasADi (widely used for trajectory optimisation) provide it; the adjoint methods of aerodynamic shape optimisation are reverse mode derived by hand. Exact gradients are what make BFGS and Newton work as advertised on large problems.
Derivation
Forward-mode automatic differentiation, by hand
Carry every quantity as a pair (value, derivative) and give each operation its calculus rule. For f(x)=xsin⁡xf(x) = x\sin x at x=2x = 2:
  1. x→(2, 1)x \to (2,\ 1) — the input, with dx/dx=1dx/dx = 1.
  2. s=sin⁡x→(sin⁡2, cos⁡2⋅1)=(0.909297, −0.416147)s = \sin x \to (\sin 2,\ \cos 2 \cdot 1) = (0.909297,\ -0.416147) — the chain rule.
  3. f=x⋅s→(2×0.909297, 1×0.909297+2×(−0.416147))=(1.818595, 0.077004)f = x \cdot s \to \big(2 \times 0.909297,\ 1 \times 0.909297 + 2 \times (-0.416147)\big) = (1.818595,\ 0.077004) — the product rule.
The second number is f′(2)=sin⁡2+2cos⁡2f'(2) = \sin 2 + 2\cos 2, exact to rounding: no step hh, no subtraction of nearly equal numbers. A program of any length is a sequence of such elementary operations, so the same bookkeeping differentiates it. Reverse mode records the operations on the way forward and then applies the chain rule backwards from the output, so one sweep gives the derivative with respect to every input at once — which is why it wins when there are millions of inputs and one output.
⚠Scaling and gradients
Optimisers work best when variables have similar magnitudes — a length in metres beside a thickness in metres spans three orders of magnitude; rescale them to order one. Supply analytic gradients when you can (jac=); finite-difference gradients (lesson 8-1) cost nn extra evaluations and lose half the digits. For simulation-based objectives, consider whether the simulation's own tolerance makes the objective noisy — gradient methods then stall, and derivative-free or surrogate methods may be better.
💡Where you will use this
Optimal Guidance: BFGS and SQP inside trajectory optimisers. Orbit determination and GNSS: nonlinear least squares (Unit 5) is a special Newton-type method. Control engineering: controller tuning on simulations. Structures and Aerodynamics: shape optimisation with adjoint gradients and L-BFGS. Data fitting: maximum-likelihood estimation.
📚Further reading
  • Jorge Nocedal and Stephen Wright, Numerical Optimization (2nd ed., Springer, 2006), chapter 6, 'Quasi-Newton Methods', and chapter 8, 'Calculating Derivatives' — BFGS and its positive definiteness, finite differences and automatic differentiation.
  • Andreas Griewank and Andrea Walther, Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation (2nd ed., SIAM, 2008).
  • William Davidon, 'Variable metric method for minimization', SIAM Journal on Optimization 1(1), 1–17 (1991) — the 1959 report, finally published.
↓But there is a problem
Every problem in this lesson was unconstrained: any x\mathbf x was allowed. Real design never is. A beam's stress must stay below the yield stress, a tank must hold a given volume, a budget is a budget, a trajectory must avoid the ground. Often the best design sits exactly on a limit. In 1947 George Dantzig, a mathematician working for the US Air Force, was asked to mechanise the planning of military logistics — what to move, where and when — under hundreds of limits on aircraft, crews and supplies. Every quantity in his problem was linear. That turned out to make it both enormously large and, remarkably, solvable.