Broyden–Fletcher–Goldfarb–Shanno Algorithm

Broyden–Fletcher–Goldfarb–Shanno algorithm (BFGS) is an iterative quasi-Newton method in numerical optimization and scientific computing for solving unconstrained nonlinear minimization problems min_x∈ℝⁿ f(x) where f is twice continuously differentiable. Discovered independently in 1970 by Charles George Broyden, Roger Fletcher, Donald Goldfarb, and David F. Shanno, the BFGS algorithm bridges the gap between first-order gradient descent—which requires only O(n) operations per step but converges slowly on ill-conditioned valleys—and second-order Newton's method, which achieves quadratic convergence but requires evaluating and factoring the n × n Hessian matrix ∇²f(x_k) at a cost of O(n³) operations per step.

Instead of computing second derivatives directly, BFGS accumulates curvature information across successive steps by enforcing the secant equation on the observed step displacement s_k = x_k+1 − x_k and gradient change y_k = ∇f(x_k+1) − ∇f(x_k). By applying a symmetric rank-two update directly to the inverse Hessian approximation H_k ≈ [∇²f(x_k)]⁻¹, BFGS computes each search direction in O(n²) operations without solving a linear system via Cholesky decomposition or LU decomposition, while achieving local superlinear convergence.

At iteration k, BFGS constructs a local quadratic model of the objective function around x_k:

m_k(p) = f(x_k) + ∇f(x_k)^Tp + (1/2) p^TB_kp,

where B_k ∈ ℝ^n×n is a symmetric positive-definite matrix approximating the true Hessian ∇²f(x_k). Minimizing m_k(p) yields the quasi-Newton search direction p_k = −B_k⁻¹∇f(x_k) = −H_k∇f(x_k), where H_k = B_k⁻¹ denotes the inverse Hessian approximation. After stepping to x_k+1 = x_k + α_kp_k along a step length α_k > 0, the updated Hessian approximation B_k+1 is required to mimic the action of the true Hessian along the most recent displacement s_k = x_k+1 − x_k:

B_k+1 s_k = y_k, or equivalently H_k+1 y_k = s_k,

where y_k = ∇f(x_k + s_k) − ∇f(x_k). Because the secant equation imposes only n linear constraints on the n(n + 1)/2 independent entries of the symmetric matrix H_k+1, infinitely many symmetric matrices satisfy H_k+1y_k = s_k.

Variational Derivation and the BFGS Update Formulas

Whereas the earlier Davidon–Fletcher–Powell (DFP) formula minimizes the weighted Frobenius distance ‖W^1/2(B − B_k)W^1/2‖_F on the forward Hessian B, BFGS solves the dual variational problem on the inverse Hessian H:

min_H ‖W^1/2(H − H_k)W^1/2‖_F subject to H = H^T and Hy_k = s_k,

where W is any symmetric positive-definite weight matrix satisfying Ws_k = y_k (such as the average Gauss–Newton Hessian). Solving this constrained matrix minimization yields the BFGS inverse Hessian update:

H_k+1 = (I − ρ_ks_ky_k^T) H_k (I − ρ_ky_ks_k^T) + ρ_ks_ks_k^T, where ρ_k = 1 / (y_k^Ts_k).

Applying the Sherman–Morrison–Woodbury matrix inversion identity produces the corresponding forward Hessian BFGS update, which modifies B_k by a rank-two correction:

B_k+1 = B_k − (B_ks_ks_k^TB_k) / (s_k^TB_ks_k) + (y_ky_k^T) / (y_k^Ts_k).

Hereditary Positive-Definiteness via the Wolfe Conditions

For p_k+1 = −H_k+1∇f(x_k+1) to remain a strict descent direction, H_k+1 must inherit positive-definiteness from H_k. For any non-zero vector z ∈ ℝⁿ, forming the quadratic form z^TH_k+1z shows that H_k+1 is positive-definite if and only if the curvature condition holds:

s_k^Ty_k > 0.

Geometrically, s_k^Ty_k = α_k [∇f(x_k+1)^Tp_k − ∇f(x_k)^Tp_k] > 0 means that the directional slope along p_k becomes less negative between x_k and x_k+1. This condition is automatically enforced whenever the step size α_k satisfies the Wolfe line-search conditions:

1. Armijo sufficient decrease: f(x_k + α_kp_k) ≤ f(x_k) + c₁α_k∇f(x_k)^Tp_k (typically c₁ = 10⁻⁴)
2. Curvature condition: ∇f(x_k + α_kp_k)^Tp_k ≥ c₂∇f(x_k)^Tp_k (typically c₂ = 0.9).

Because 0 < c₁ < c₂ < 1 and ∇f(x_k)^Tp_k < 0, any Wolfe step guarantees s_k^Ty_k ≥ (1 − c₂) |α_k∇f(x_k)^Tp_k| > 0. By the Dennis–Moré theorem (1974), BFGS with unit step length α_k = 1 accepted asymptotically achieves Q-superlinear convergence (lim_k→∞ ‖x_k+1 − x‖ / ‖x_k − x‖ = 0) and self-corrects poor initial Hessian scalings far more rapidly than DFP.

When the number of variables n exceeds 10⁴—as in neural network training, seismic waveform inversion, numerical weather variational data assimilation, or electronic-structure geometry relaxation—storing a dense n × n matrix H_k is impossible. Introduced by Jorge Nocedal in 1980, Limited-Memory BFGS (L-BFGS) never forms H_k explicitly. Instead, it stores only the m most recent vector pairs {(s_i, y_i)}_i=k−m^k−1 (with m typically between 5 and 20) and computes the matrix-vector product H_k∇f(x_k) directly using Nocedal's two-loop recursion in 4mn multiplications, starting from the scaled identity seed H_k⁽⁰⁾ = [(s_k−1^Ty_k−1)/(y_k−1^Ty_k−1)] I. Its bound-constrained extension L-BFGS-B combines gradient projection along the box constraints l ≤ x ≤ u with subspace minimization over the free variables.

Why is BFGS preferred over the DFP update in numerical software?

Although DFP and BFGS are mathematical duals (interchanging B ↔ H and s_k ↔ y_k), numerical analysis by Powell and Byrd–Nocedal showed that when H_k accidentally acquires a tiny eigenvalue, DFP struggles to recover and can stall for hundreds of iterations, whereas the rank-one term +ρ_ks_ks_k^T in the BFGS inverse update immediately restores step length along s_k.

How should the initial inverse Hessian matrix H0 be chosen?

Setting H₀ = I takes a first step identical to steepest descent, after which scaling H₀ ← [(s₀^Ty₀)/(y₀^Ty₀)] I before applying the first BFGS update matches the inverse Rayleigh quotient of the true Hessian along s₀, making the second step scale-invariant.

1. Broyden, C. G. (1970). The convergence of a class of double-rank minimization algorithms 1. General considerations. Journal of the Institute of Mathematics and Its Applications, 6(1), 76–90. https://doi.org/10.1093/imamat/6.1.76
2. Fletcher, R. (1970). A new approach to variable metric algorithms. The Computer Journal, 13(3), 317–322. https://doi.org/10.1093/comjnl/13.3.317
3. Goldfarb, D. (1970). A family of variable-metric methods derived by variational means. Mathematics of Computation, 24(109), 23–26. https://doi.org/10.1090/S0025-5718-1970-0258249-6
4. Shanno, D. F. (1970). Conditioning of quasi-Newton methods for function minimization. Mathematics of Computation, 24(111), 647–656. https://doi.org/10.1090/S0025-5718-1970-0274029-X
5. Liu, D. C., & Nocedal, J. (1989). On the limited memory BFGS method for large scale optimization. Mathematical Programming, 45(1), 503–528. https://doi.org/10.1007/BF01589116