Broyden–Fletcher–Goldfarb–Shanno algorithm

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

The algorithm begins at an initial estimate $\backslash 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:

:$B\_k\; \backslash mathbf\{p\}\_k\; =\; -\backslash nabla\; f(\backslash mathbf\{x\}\_k),$

where $B\_k$ is an approximation to the Hessian matrix at $\backslash mathbf\{x\}\_k$, which is updated iteratively at each stage, and $\backslash nabla\; f(\backslash 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 $f(\backslash mathbf\{x\}\_k\; +\; \backslash gamma\backslash mathbf\{p\}\_k)$ over the scalar $\backslash gamma\; >\; 0.$

The quasi-Newton condition imposed on the update of $B\_k$ is

:$B\_\{k+1\}\; (\backslash mathbf\{x\}\_\{k+1\}\; -\; \backslash mathbf\{x\}\_k)\; =\; \backslash nabla\; f(\backslash mathbf\{x\}\_\{k+1\})\; -\; \backslash nabla\; f(\backslash mathbf\{x\}\_k).$

Let $\backslash mathbf\{y\}\_k\; =\; \backslash nabla\; f(\backslash mathbf\{x\}\_\{k+1\})\; -\; \backslash nabla\; f(\backslash mathbf\{x\}\_k)$ and $\backslash mathbf\{s\}\_k\; =\; \backslash mathbf\{x\}\_\{k+1\}\; -\; \backslash mathbf\{x\}\_k$, then $B\_\{k+1\}$ satisfies

:$B\_\{k+1\}\; \backslash mathbf\{s\}\_k\; =\; \backslash mathbf\{y\}\_k$,

which is the secant equation.

The curvature condition $\backslash mathbf\{s\}\_k^\backslash top\; \backslash mathbf\{y\}\_k\; >\; 0$ should be satisfied for $B\_\{k+1\}$ to be positive definite, which can be verified by pre-multiplying the secant equation with $\backslash 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 $\backslash mathbf\{x\}\_\{k+1\}$ to be computed as $B\_\{k+1\}$, the approximate Hessian at stage k is updated by the addition of two matrices:

:$B\_\{k+1\}\; =\; B\_k\; +\; U\_k\; +\; V\_k.$

Both $U\_k$ and $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 $B\_\{k+1\}$, the update form can be chosen as $B\_\{k+1\}\; =\; B\_k\; +\; \backslash alpha\backslash mathbf\{u\}\backslash mathbf\{u\}^\backslash top\; +\; \backslash beta\backslash mathbf\{v\}\backslash mathbf\{v\}^\backslash top$. Imposing the secant condition, $B\_\{k+1\}\; \backslash mathbf\{s\}\_k\; =\; \backslash mathbf\{y\}\_k$. Choosing $\backslash mathbf\{u\}\; =\; \backslash mathbf\{y\}\_k$ and $\backslash mathbf\{v\}\; =\; B\_k\; \backslash mathbf\{s\}\_k$, we can obtain:{{Citation | last1=Fletcher | first1=Roger | title=Practical methods of optimization | publisher=John Wiley & Sons | location=New York | edition=2nd | isbn=978-0-471-91547-8 | year=1987 | url-access=registration | url=https://archive.org/details/practicalmethods0000flet }}

:$\backslash alpha\; =\; \backslash frac\{1\}\{\backslash mathbf\{y\}\_k^T\; \backslash mathbf\{s\}\_k\},$

:$\backslash beta\; =\; -\backslash frac\{1\}\{\backslash mathbf\{s\}\_k^T\; B\_k\; \backslash mathbf\{s\}\_k\}.$

Finally, we substitute $\backslash alpha$ and $\backslash beta$ into $B\_\{k+1\}\; =\; B\_k\; +\; \backslash alpha\backslash mathbf\{u\}\backslash mathbf\{u\}^\backslash top\; +\; \backslash beta\backslash mathbf\{v\}\backslash mathbf\{v\}^\backslash top$ and get the update equation of $B\_\{k+1\}$:

:$B\_\{k+1\}\; =\; B\_k\; +\; \backslash frac\{\backslash mathbf\{y\}\_k\; \backslash mathbf\{y\}\_k^\{\backslash mathrm\{T\}\}\}\{\backslash mathbf\{y\}\_k^\{\backslash mathrm\{T\}\}\; \backslash mathbf\{s\}\_k\}\; -\; \backslash frac\{B\_k\; \backslash mathbf\{s\}\_k\; \backslash mathbf\{s\}\_k^\{\backslash mathrm\{T\}\}\; B\_k^\{\backslash mathrm\{T\}\}\; \}\{\backslash mathbf\{s\}\_k^\{\backslash mathrm\{T\}\}\; B\_k\; \backslash mathbf\{s\}\_k\}.$

Consider the following unconstrained optimization problem

From an initial guess $\backslash mathbf\{x\}\_0$ and an approximate inverted Hessian matrix $H\_0$ the following steps are repeated as $\backslash mathbf\{x\}\_k$ converges to the solution:

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

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 {{Citation needed|reason=there is a need to see in which works this was applied, esp. given the next sentence|date=May 2021}}. However, these quantities are technically defined by the true Hessian matrix, and the BFGS approximation may not converge to the true Hessian matrix.{{Cite journal | last1=Ge | last2=Powell| first1=Ren-pu | first2=M. J. D. | title=The Convergence of Variable Metric Matrices in Unconstrained Optimization | journal=Mathematical Programming |volume=27| year=1983 | issue=2|at=123 |doi=10.1007/BF02591941 | s2cid=8113073}}

The BFGS update formula heavily relies on the curvature $\backslash mathbf\{s\}\_k^\backslash top\; \backslash 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

{{citation

| title = Numerical Optimization

| author = Jorge Nocedal

| author2 = Stephen J. Wright

| year = 2006

}}

) which modify $\backslash mathbf\{s\}\_k$ and/or $\backslash 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.{{Cite web|title=GNU Scientific Library — GSL 2.6 documentation|url=https://www.gnu.org/software/gsl/doc/html/index.html|access-date=2020-11-22|website=www.gnu.org}}
• 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().{{Cite web|title=R: General-purpose Optimization|url=https://stat.ethz.ch/R-manual/R-devel/library/stats/html/optim.html|access-date=2020-11-22|website=stat.ethz.ch}}
• In SciPy, the scipy.optimize.fmin_bfgs function implements BFGS.{{Cite web|title=scipy.optimize.fmin_bfgs — SciPy v1.5.4 Reference Guide|url=https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin_bfgs.html|access-date=2020-11-22|website=docs.scipy.org}} It is also possible to run BFGS using any of the L-BFGS algorithms by setting the parameter L to a very large number. It is also one of the default methods used when running scipy.optimize.minimize with no constraints.{{Cite web|title=scipy.optimize.minimize — SciPy v1.5.4 Reference Guide|url=https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html|access-date=2025-01-22|website=docs.scipy.org}}
• In Julia, the [https://julianlsolvers.github.io/Optim.jl/stable/ Optim.jl] package implements BFGS and L-BFGS as a solver option to the optimize() function (among other options).{{cite web |title=Optim.jl Configurable options |url=https://julianlsolvers.github.io/Optim.jl/stable/#user/config/#solver-options |website=julianlsolvers}}

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.

{{Div col|colwidth=25em}}

• 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
• Compact quasi-Newton representation

{{Div col end}}

{{reflist}}

• {{Citation | last1=Avriel |first1=Mordecai |year=2003|title=Nonlinear Programming: Analysis and Methods|publisher= Dover Publishing|isbn= 978-0-486-43227-4}}
• {{Citation|last1=Bonnans |first1=J. Frédéric |last2=Gilbert |first2=J. Charles |last3=Lemaréchal |first3=Claude |author-link3=Claude Lemaréchal |last4=Sagastizábal |first4=Claudia A. |author4-link=Claudia Sagastizábal |title=Numerical Optimization: Theoretical and Practical Aspects |edition=Second |publisher=Springer |location=Berlin|year=2006 |isbn=3-540-35445-X |chapter=Newtonian Methods |pages=51–66 }}
• {{Citation | last1=Fletcher | first1=Roger | title=Practical Methods of Optimization | publisher=John Wiley & Sons | location=New York | edition=2nd | isbn=978-0-471-91547-8 | year=1987 | url-access=registration | url=https://archive.org/details/practicalmethods0000flet }}
• {{Citation|last1=Luenberger|first1=David G.|author-link1=David G. Luenberger|last2=Ye|first2=Yinyu|author-link2=Yinyu Ye|title=Linear and nonlinear programming|edition=Third|series=International Series in Operations Research & Management Science|volume=116|publisher=Springer|location=New York|year=2008|pages=xiv+546|isbn=978-0-387-74502-2| mr = 2423726}}
• {{citation |last=Kelley |first=C. T. |title=Iterative Methods for Optimization |location=Philadelphia |publisher=Society for Industrial and Applied Mathematics |year=1999 |isbn=0-89871-433-8 |pages=71–86 }}

{{Optimization algorithms|unconstrained}}

{{DEFAULTSORT:Broyden-Fletcher-Goldfarb-Shanno algorithm}}

Category:Optimization algorithms and methods