Broyden–Fletcher–Goldfarb–Shanno algorithm

In numerical optimization, the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm is an iterative method for solving unconstrained nonlinear optimization problems. Like the related Davidon–Fletcher–Powell method, BFGS determines the descent direction by preconditioning the gradient with curvature information. It does so by gradually improving an approximation to the Hessian matrix of the loss function, obtained only from gradient evaluations (or approximate gradient evaluations) via a generalized secant method.

Since the updates of the BFGS curvature matrix do not require matrix inversion, its computational complexity is only , compared to  in Newton's method. Also in common use is L-BFGS, which is a limited-memory version of BFGS that is particularly suited to problems with very large numbers of variables (e.g., >1000). The BFGS-B variant handles simple box constraints.

The algorithm is named after Charles George Broyden, Roger Fletcher, Donald Goldfarb and David Shanno.

Rationale 
The optimization problem is to minimize , where  is a vector in , and  is a differentiable scalar function. There are no constraints on the values that  can take.

The algorithm begins at an initial estimate for the optimal value  and proceeds iteratively to get a better estimate at each stage.

The search direction pk at stage k is given by the solution of the analogue of the Newton equation:

where  is an approximation to the Hessian matrix at , which is updated iteratively at each stage, and  is the gradient of the function evaluated at xk.  A line search in the direction pk is then used to find the next point xk+1 by minimizing  over the scalar 

The quasi-Newton condition imposed on the update of  is

Let  and , then  satisfies 

, 

which is the secant equation. 

The curvature condition  should be satisfied for  to be positive definite, which can be verified by pre-multiplying the secant equation with . If the function is not strongly convex, then the condition has to be enforced explicitly e.g. by finding a point xk+1 satisfying the Wolfe conditions, which entail the curvature condition, using line search.

Instead of requiring the full Hessian matrix at the point  to be computed as , the approximate Hessian at stage k is updated by the addition of two matrices:

Both  and  are symmetric rank-one matrices, but their sum is a rank-two update matrix. BFGS and DFP updating matrix both differ from its predecessor by a rank-two matrix. Another simpler rank-one method is known as symmetric rank-one method, which does not guarantee the positive definiteness. In order to maintain the symmetry and positive definiteness of , the update form can be chosen as . Imposing the secant condition, . Choosing  and , we can obtain:

Finally, we substitute  and  into  and get the update equation of :

Algorithm 

From an initial guess  and an approximate Hessian matrix  the following steps are repeated as  converges to the solution:
 Obtain a direction  by solving .
 Perform a one-dimensional optimization (line search) to find an acceptable stepsize  in the direction found in the first step. If an exact line search is performed, then  . In practice, an inexact line search usually suffices, with an acceptable  satisfying Wolfe conditions.
 Set  and update .
 .
 .

 denotes the objective function to be minimized. Convergence can be checked by observing the norm of the gradient, .  If  is initialized with , the first step will be equivalent to a gradient descent, but further steps are more and more refined by , the approximation to the Hessian.

The first step of the algorithm is carried out using the inverse of the matrix , which can be obtained efficiently by applying the Sherman–Morrison formula to the step 5 of the algorithm, giving

 

This can be computed efficiently without temporary matrices, recognizing that  is symmetric,
and that  and  are scalars, using an expansion such as

 

Therefore, in order to avoid any matrix inversion, the inverse of the Hessian can be approximated instead of the Hessian itself: 

From an initial guess  and an approximate inverted Hessian matrix  the following steps are repeated as  converges to the solution:
 Obtain a direction  by solving .
 Perform a one-dimensional optimization (line search) to find an acceptable stepsize  in the direction found in the first step. If an exact line search is performed, then  . In practice, an inexact line search usually suffices, with an acceptable  satisfying Wolfe conditions.
 Set  and update .
 .
 .

In statistical estimation problems (such as maximum likelihood or Bayesian inference), credible intervals or confidence intervals for the solution can be estimated from the inverse of the final Hessian matrix . However, these quantities are technically defined by the true Hessian matrix, and the BFGS approximation may not converge to the true Hessian matrix.

Further developments 

The BFGS update formula heavily relies on the curvature  being strictly positive and bounded away from zero.
This condition is satisfied when we perform a line search with Wolfe conditions on a convex target.
However, some real-life applications (like Sequential Quadratic Programming methods) routinely produce negative or nearly-zero curvatures.
This can occur when optimizing a nonconvex target or when employing a trust-region approach instead of a line search.
It is also possible to produce spurious values due to noise in the target.

In such cases, one of the so-called damped BFGS updates can be used (see ) which modify  and/or  in order to obtain a more robust update.

Notable implementations 

Notable open source implementations are:

 ALGLIB implements BFGS and its limited-memory version in C++ and C#
 GNU Octave uses a form of BFGS in its fsolve function, with trust region extensions.
 The GSL implements BFGS as gsl_multimin_fdfminimizer_vector_bfgs2.
 In R, the BFGS algorithm (and the L-BFGS-B version that allows box constraints) is implemented as an option of the base function optim().
 In SciPy, the scipy.optimize.fmin_bfgs function implements BFGS. It is also possible to run BFGS using any of the L-BFGS algorithms by setting the parameter L to a very large number.
 In Julia, the Optim.jl package implements BFGS and L-BFGS as a solver option to the optimize() function (among other options).

Notable proprietary implementations include:

 The large scale nonlinear optimization software Artelys Knitro implements, among others, both BFGS and L-BFGS algorithms.
 In the MATLAB Optimization Toolbox, the fminunc function uses BFGS with cubic line search when the problem size is set to "medium scale."
 Mathematica includes BFGS.

See also 

 BHHH algorithm
 Davidon–Fletcher–Powell formula
 Gradient descent
 L-BFGS
 Levenberg–Marquardt algorithm
 Nelder–Mead method
 Pattern search (optimization)
 Quasi-Newton methods
 Symmetric rank-one

References

Further reading 
 
 
 
 
 

Optimization algorithms and methods