Broyden–Fletcher–Goldfarb–Shanno algorithm

The optimization problem is to minimize ${\displaystyle f(\mathbf {x} )}$ , where ${\displaystyle \mathbf {x} }$  is a vector in ${\displaystyle \mathbb {R} ^{n}}$ , and ${\displaystyle f}$  is a differentiable scalar function. There are no constraints on the values that ${\displaystyle \mathbf {x} }$  can take.

The algorithm begins at an initial estimate ${\displaystyle \mathbf {x} _{0}}$  for the optimal value and proceeds iteratively to get a better estimate at each stage.

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

${\displaystyle B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k}),}$ 

where ${\displaystyle B_{k}}$  is an approximation to the Hessian matrixat${\displaystyle \mathbf {x} _{k}}$ , which is updated iteratively at each stage, and ${\displaystyle \nabla f(\mathbf {x} _{k})}$  is the gradient of the function evaluated at x_k. A line search in the direction p_k is then used to find the next point x_k+1 by minimizing ${\displaystyle f(\mathbf {x} _{k}+\gamma \mathbf {p} _{k})}$  over the scalar ${\displaystyle \gamma >0.}$ 

The quasi-Newton condition imposed on the update of ${\displaystyle B_{k}}$ is

${\displaystyle B_{k+1}(\mathbf {x} _{k+1}-\mathbf {x} _{k})=\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k}).}$ 

Let ${\displaystyle \mathbf {y} _{k}=\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}$  and ${\displaystyle \mathbf {s} _{k}=\mathbf {x} _{k+1}-\mathbf {x} _{k}}$ , then ${\displaystyle B_{k+1}}$  satisfies

${\displaystyle B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}}$ ,

which is the secant equation.

The curvature condition ${\displaystyle \mathbf {s} _{k}^{\top }\mathbf {y} _{k}>0}$  should be satisfied for ${\displaystyle B_{k+1}}$  to be positive definite, which can be verified by pre-multiplying the secant equation with ${\displaystyle \mathbf {s} _{k}^{T}}$ . If the function is not strongly convex, then the condition has to be enforced explicitly e.g. by finding a point x_k+1 satisfying the Wolfe conditions, which entail the curvature condition, using line search.

Instead of requiring the full Hessian matrix at the point ${\displaystyle \mathbf {x} _{k+1}}$  to be computed as ${\displaystyle B_{k+1}}$ , the approximate Hessian at stage k is updated by the addition of two matrices:

${\displaystyle B_{k+1}=B_{k}+U_{k}+V_{k}.}$ 

Both ${\displaystyle U_{k}}$  and ${\displaystyle V_{k}}$  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 ${\displaystyle B_{k+1}}$ , the update form can be chosen as ${\displaystyle B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\top }+\beta \mathbf {v} \mathbf {v} ^{\top }}$ . Imposing the secant condition, ${\displaystyle B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}}$ . Choosing ${\displaystyle \mathbf {u} =\mathbf {y} _{k}}$  and ${\displaystyle \mathbf {v} =B_{k}\mathbf {s} _{k}}$ , we can obtain:^[8]

${\displaystyle \alpha ={\frac {1}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}},}$ 

${\displaystyle \beta =-{\frac {1}{\mathbf {s} _{k}^{T}B_{k}\mathbf {s} _{k}}}.}$ 

Finally, we substitute ${\displaystyle \alpha }$  and ${\displaystyle \beta }$  into ${\displaystyle B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\top }+\beta \mathbf {v} \mathbf {v} ^{\top }}$  and get the update equation of ${\displaystyle B_{k+1}}$ :

${\displaystyle B_{k+1}=B_{k}+{\frac {\mathbf {y} _{k}\mathbf {y} _{k}^{\mathrm {T} }}{\mathbf {y} _{k}^{\mathrm {T} }\mathbf {s} _{k}}}-{\frac {B_{k}\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} }B_{k}^{\mathrm {T} }}{\mathbf {s} _{k}^{\mathrm {T} }B_{k}\mathbf {s} _{k}}}.}$ 

Consider the following unconstrained optimization problem $${\displaystyle {\begin{aligned}{\underset {\mathbf {x} \in \mathbb {R} ^{n}}{\text{minimize}}}\quad &f(\mathbf {x} ),\end{aligned}}}$$  where ${\displaystyle f:\mathbb {R} ^{n}\to \mathbb {R} }$  is a nonlinear objective function.

From an initial guess ${\displaystyle \mathbf {x} _{0}\in \mathbb {R} ^{n}}$  and an initial guess of the Hessian matrix ${\displaystyle B_{0}\in \mathbb {R} ^{n\times n}}$  the following steps are repeated as ${\displaystyle \mathbf {x} _{k}}$  converges to the solution:

1. Obtain a direction ${\displaystyle \mathbf {p} _{k}}$  by solving ${\displaystyle B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k})}$ .
2. Perform a one-dimensional optimization (line search) to find an acceptable stepsize ${\displaystyle \alpha _{k}}$  in the direction found in the first step. If an exact line search is performed, then ${\displaystyle \alpha _{k}=\arg \min f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})}$  . In practice, an inexact line search usually suffices, with an acceptable ${\displaystyle \alpha _{k}}$  satisfying Wolfe conditions.
3. Set ${\displaystyle \mathbf {s} _{k}=\alpha _{k}\mathbf {p} _{k}}$  and update ${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\mathbf {s} _{k}}$ .
4. ${\displaystyle \mathbf {y} _{k}={\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}}$ .
5. ${\displaystyle B_{k+1}=B_{k}+{\frac {\mathbf {y} _{k}\mathbf {y} _{k}^{\mathrm {T} }}{\mathbf {y} _{k}^{\mathrm {T} }\mathbf {s} _{k}}}-{\frac {B_{k}\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} }B_{k}^{\mathrm {T} }}{\mathbf {s} _{k}^{\mathrm {T} }B_{k}\mathbf {s} _{k}}}}$ .

Convergence can be determined by observing the norm of the gradient; given some ${\displaystyle \epsilon >0}$ , one may stop the algothim when ${\displaystyle ||\nabla f(\mathbf {x} _{k})||\leq \epsilon .}$ If${\displaystyle B_{0}}$  is initialized with ${\displaystyle B_{0}=I}$ , the first step will be equivalent to a gradient descent, but further steps are more and more refined by ${\displaystyle B_{k}}$ , the approximation to the Hessian.

The first step of the algorithm is carried out using the inverse of the matrix ${\displaystyle B_{k}}$ , which can be obtained efficiently by applying the Sherman–Morrison formula to the step 5 of the algorithm, giving

${\displaystyle B_{k+1}^{-1}=\left(I-{\frac {\mathbf {s} _{k}\mathbf {y} _{k}^{T}}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}}\right)B_{k}^{-1}\left(I-{\frac {\mathbf {y} _{k}\mathbf {s} _{k}^{T}}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}}\right)+{\frac {\mathbf {s} _{k}\mathbf {s} _{k}^{T}}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}}.}$ 

This can be computed efficiently without temporary matrices, recognizing that ${\displaystyle B_{k}^{-1}}$  is symmetric, and that ${\displaystyle \mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}\mathbf {y} _{k}}$  and ${\displaystyle \mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}}$  are scalars, using an expansion such as

${\displaystyle B_{k+1}^{-1}=B_{k}^{-1}+{\frac {(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}+\mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}\mathbf {y} _{k})(\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} })}{(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k})^{2}}}-{\frac {B_{k}^{-1}\mathbf {y} _{k}\mathbf {s} _{k}^{\mathrm {T} }+\mathbf {s} _{k}\mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}}{\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}}}.}$ 

Therefore, in order to avoid any matrix inversion, the inverse of the Hessian can be approximated instead of the Hessian itself: ${\displaystyle H_{k}{\overset {\operatorname {def} }{=}}B_{k}^{-1}.}$ ^[9]

From an initial guess ${\displaystyle \mathbf {x} _{0}}$  and an approximate inverted Hessian matrix ${\displaystyle H_{0}}$  the following steps are repeated as ${\displaystyle \mathbf {x} _{k}}$  converges to the solution:

1. Obtain a direction ${\displaystyle \mathbf {p} _{k}}$  by solving ${\displaystyle \mathbf {p} _{k}=-H_{k}\nabla f(\mathbf {x} _{k})}$ .
2. Perform a one-dimensional optimization (line search) to find an acceptable stepsize ${\displaystyle \alpha _{k}}$  in the direction found in the first step. If an exact line search is performed, then ${\displaystyle \alpha _{k}=\arg \min f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})}$  . In practice, an inexact line search usually suffices, with an acceptable ${\displaystyle \alpha _{k}}$  satisfying Wolfe conditions.
3. Set ${\displaystyle \mathbf {s} _{k}=\alpha _{k}\mathbf {p} _{k}}$  and update ${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\mathbf {s} _{k}}$ .
4. ${\displaystyle \mathbf {y} _{k}={\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}}$ .
5. ${\displaystyle H_{k+1}=H_{k}+{\frac {(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}+\mathbf {y} _{k}^{\mathrm {T} }H_{k}\mathbf {y} _{k})(\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} })}{(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k})^{2}}}-{\frac {H_{k}\mathbf {y} _{k}\mathbf {s} _{k}^{\mathrm {T} }+\mathbf {s} _{k}\mathbf {y} _{k}^{\mathrm {T} }H_{k}}{\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}}}}$ .

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

The BFGS update formula heavily relies on the curvature ${\displaystyle \mathbf {s} _{k}^{\top }\mathbf {y} _{k}}$  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 ^[11]) which modify ${\displaystyle \mathbf {s} _{k}}$  and/or ${\displaystyle \mathbf {y} _{k}}$  in order to obtain a more robust update.

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.^[12]
• InR, the BFGS algorithm (and the L-BFGS-B version that allows box constraints) is implemented as an option of the base function optim().^[13]
• InSciPy, the scipy.optimize.fmin_bfgs function implements BFGS.^[14] It is also possible to run BFGS using any of the L-BFGS algorithms by setting the parameter L to a very large number.
• InJulia, the Optim.jl package implements BFGS and L-BFGS as a solver option to the optimize() function (among other options).^[15]

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.
• LS-DYNA also uses BFGS to solve implicit Problems.

• Avriel, Mordecai (2003), Nonlinear Programming: Analysis and Methods, Dover Publishing, ISBN 978-0-486-43227-4
• Bonnans, J. Frédéric; Gilbert, J. Charles; Lemaréchal, Claude; Sagastizábal, Claudia A. (2006), "Newtonian Methods", Numerical Optimization: Theoretical and Practical Aspects (Second ed.), Berlin: Springer, pp. 51–66, ISBN 3-540-35445-X
• Fletcher, Roger (1987), Practical Methods of Optimization (2nd ed.), New York: John Wiley & Sons, ISBN 978-0-471-91547-8
• Luenberger, David G.; Ye, Yinyu (2008), Linear and nonlinear programming, International Series in Operations Research & Management Science, vol. 116 (Third ed.), New York: Springer, pp. xiv+546, ISBN 978-0-387-74502-2, MR 2423726
• Kelley, C. T. (1999), Iterative Methods for Optimization, Philadelphia: Society for Industrial and Applied Mathematics, pp. 71–86, ISBN 0-89871-433-8