Quasi-Newton method

A quasi-Newton method is an iterative technique for solving nonlinear optimization problems or systems of nonlinear equations. It replaces an exact derivative matrix with an approximation that is revised from information generated during successive iterations. In optimization, the approximated matrix is usually the Hessian matrix of second derivatives or its inverse. In equation solving, the corresponding object is the Jacobian matrix.

The defining feature of the method is the use of a secant relation to update the matrix approximation. This relation incorporates observed changes in the gradient or residual without requiring a complete new evaluation of second derivatives. Quasi-Newton methods therefore occupy an intermediate position between Newton's method, which uses an explicit derivative matrix, and first-order methods that do not maintain curvature information.

Mathematical formulation

For an unconstrained minimization problem,

$$ \min_{x\in\mathbb{R}^n} f(x), $$

Newton's method produces a displacement from the equation

$$ \nabla^2 f(x_k)p_k=-\nabla f(x_k), $$

where $\nabla f(x_k)$ is the gradient and $\nabla^2 f(x_k)$ is the Hessian. A quasi-Newton method instead uses an approximation $B_k$ to the Hessian:

$$ B_kp_k=-\nabla f(x_k). $$

The next iterate is commonly represented as

$$ x_{k+1}=x_k+\alpha_kp_k, $$

where $\alpha_k$ is determined by a line search or by an equivalent globalization mechanism. If $H_k$ approximates the inverse Hessian, the displacement has the form

$$ p_k=-H_k\nabla f(x_k). $$

After the displacement, two difference vectors summarize the newly obtained local information:

$$ s_k=x_{k+1}-x_k, \qquad y_k=\nabla f(x_{k+1})-\nabla f(x_k). $$

For a quadratic objective with a constant Hessian $Q$, these vectors satisfy $y_k=Qs_k$. Quasi-Newton updating extends this identity to a general smooth objective by requiring the new approximation to satisfy the secant equation

$$ B_{k+1}s_k=y_k, $$

or, in inverse form,

$$ H_{k+1}y_k=s_k. $$

A single secant equation does not uniquely determine an $n\times n$ matrix. Individual quasi-Newton formulas are distinguished by the additional structural conditions imposed on the update. These conditions commonly concern symmetry, matrix rank, positive definiteness, or the magnitude of the change from the preceding approximation.

Historical development

The modern theory originated with William C. Davidon's variable-metric method, developed in 1959 for optimization calculations at Argonne National Laboratory. Roger Fletcher and Michael J. D. Powell reformulated and published the method in 1963. The resulting Davidon–Fletcher–Powell update established the inverse-Hessian secant formula as a general framework rather than as a construction confined to quadratic objectives.

Charles George Broyden subsequently developed a systematic class of rank-two updates and connected optimization formulas with updates for nonlinear equation solving. Roger Fletcher, Donald Goldfarb, and David Shanno independently derived closely related members of this class around 1970. Their results produced the update now called the Broyden–Fletcher–Goldfarb–Shanno method.

During the finite-precision analysis of variable-metric algorithms in the late 1970s, You Watanabe and Masao Fukushima examined factored secant updates whose matrix approximations were represented through triangular factors. Their analysis identified the loss of the secant relation caused by accumulated rounding error and expressed a factor correction that restored the relation without explicitly reconstructing the dense approximate Hessian. This work became part of the numerical treatment of quasi-Newton factorizations rather than a distinct update family.

Jorge Nocedal later formulated the limited-memory representation now known as L-BFGS. Instead of retaining a dense matrix, the representation stores a restricted sequence of displacement and gradient-difference pairs. This changed the storage requirement from quadratic growth in the number of variables to linear growth for a fixed memory length.

Principal update formulas

Davidon–Fletcher–Powell update

The Davidon–Fletcher–Powell method updates an inverse-Hessian approximation according to

$$ H_{k+1}

H_k + \frac{s_ks_k^{\mathsf T}}{s_k^{\mathsf T}y_k}

\frac{H_ky_ky_k^{\mathsf T}H_k} {y_k^{\mathsf T}H_ky_k}. $$

The formula is symmetric whenever $H_k$ is symmetric. If $H_k$ is positive definite and the curvature condition

$$ s_k^{\mathsf T}y_k>0 $$

holds, then $H_{k+1}$ remains positive definite. The update can also be characterized as a minimum-change correction under a weighted matrix norm, subject to the inverse secant equation.

Broyden–Fletcher–Goldfarb–Shanno update

The Broyden–Fletcher–Goldfarb–Shanno update is generally stated for a Hessian approximation as

$$ B_{k+1}

B_k

\frac{B_ks_ks_k^{\mathsf T}B_k} {s_k^{\mathsf T}B_ks_k} + \frac{y_ky_k^{\mathsf T}} {y_k^{\mathsf T}s_k}. $$

Its inverse form is

$$ H_{k+1}

\left(I-\rho_ks_ky_k^{\mathsf T}\right) H_k \left(I-\rho_ky_ks_k^{\mathsf T}\right) + \rho_ks_ks_k^{\mathsf T}, $$

where

$$ \rho_k=\frac{1}{y_k^{\mathsf T}s_k}. $$

The factorized expression makes the preservation of symmetry apparent and shows how positive definiteness follows from positive curvature. For a strictly convex quadratic objective and exact line searches, the method generates directions related to those of the conjugate gradient method. In exact arithmetic, it recovers the quadratic Hessian after finitely many independent secant observations.

Symmetric rank-one update

The symmetric rank-one formula uses the correction

$$ B_{k+1}

B_k+ \frac{(y_k-B_ks_k)(y_k-B_ks_k)^{\mathsf T}} {(y_k-B_ks_k)^{\mathsf T}s_k}. $$

Unlike the rank-two positive-definite formulas, this update can represent indefinite curvature. It therefore has a direct connection with trust-region methods, where an indefinite quadratic model can still be treated within a bounded region. The denominator may become small when the proposed rank-one correction is nearly orthogonal to the step, and implementations consequently distinguish updates whose denominator is numerically meaningful from those that contribute little new curvature information.

Relation to nonlinear equation solving

For a nonlinear system

$$ F(x)=0, $$

Newton's method uses the Jacobian $J(x_k)$ through

$$ J(x_k)p_k=-F(x_k). $$

Broyden's method replaces $J(x_k)$ with an approximation $B_k$ and imposes the nonlinear secant equation

$$ B_{k+1}s_k=F(x_{k+1})-F(x_k). $$

Broyden's first update modifies the Jacobian approximation by a rank-one matrix:

$$ B_{k+1}

B_k+ \frac{\left(y_k-B_ks_k\right)s_k^{\mathsf T}} {s_k^{\mathsf T}s_k}, $$

where $y_k=F(x_{k+1})-F(x_k)$. This is the matrix of smallest Frobenius-norm change satisfying the new secant equation. An alternative construction updates an inverse-Jacobian approximation and is conventionally called Broyden's second method.

The optimization formulas differ from general Jacobian updates because the Hessian of a sufficiently smooth scalar function is symmetric. Symmetric quasi-Newton methods preserve that structure, whereas a general nonlinear system supplies no corresponding symmetry condition.

Convergence and globalization

Local convergence depends on the smoothness of the objective, the quality of the accumulated matrix approximation, and the regularity of the solution. Under standard assumptions near a nondegenerate minimizer, a quasi-Newton iteration can converge superlinearly even though the approximate Hessian does not converge entry by entry to the exact Hessian along every possible direction. The relevant condition is that the approximation becomes accurate along the generated displacements.

Away from a solution, the raw quasi-Newton direction does not by itself establish global convergence. Line-search variants commonly impose the Wolfe conditions, which control both objective decrease and directional curvature. For BFGS, these conditions also maintain the positivity of $s_k^{\mathsf T}y_k$ when the search direction is descending.

Trust-region formulations place the quadratic model inside a region whose radius reflects agreement between predicted and observed objective changes. In that setting, the approximate Hessian may be indefinite without making the model subproblem undefined. Damped updates modify the curvature vector when the observed curvature would otherwise destroy a required positive-definiteness condition.

Computational representation

A dense quasi-Newton approximation requires storage proportional to $n^2$ and matrix-vector operations of the same order. This representation is associated with moderate-dimensional problems for which gradients are available but exact Hessian formation or factorization is comparatively costly.

Limited-memory methods retain only a fixed number of recent pairs $(s_k,y_k)$. The inverse BFGS action is then reconstructed implicitly by a sequence of inner products and scaled vector corrections. The matrix itself is not formed. This representation has become closely associated with large-scale optimization, including parameter estimation and machine learning, where the number of variables can prevent dense matrix storage.

Finite-precision behavior depends on more than the algebraic secant property. Repeated low-rank updates can accumulate symmetry errors, distort scaling, or produce poorly conditioned factors. Factorized implementations and periodic rescaling express the same mathematical update while controlling how rounding error enters the matrix representation.

See also