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.[1]
The BFGS method belongs to quasi-Newton methods, a class of hill-climbing optimization techniques that seek a stationary point of a (preferably twice continuously differentiable) function. For such problems, a necessary condition for optimality is that the gradient be zero. Newton's method and the BFGS methods are not guaranteed to converge unless the function has a quadratic Taylor expansion near an optimum. However, BFGS has proven to have good performance even for non-smooth optimizations.[2]
In quasi-Newton methods, the Hessian matrix of second derivatives is not computed. Instead, the Hessian matrix is approximated using updates specified by gradient evaluations (or approximate gradient evaluations). Quasi-Newton methods are generalizations of the secant method to find the root of the first derivative for multidimensional problems. In multi-dimensional problems, the secant equation does not specify a unique solution, and quasi-Newton methods differ in how they constrain the solution. The BFGS method is one of the most popular members of this class.[3] 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[4] variant handles simple box constraints.
The algorithm is named after Charles George Broyden, Roger Fletcher, Donald Goldfarb and David Shanno.
Contents
1 Rationale
2 Algorithm
3 Implementations
4 See also
5 Notes
6 Bibliography
7 External links
Rationale
The optimization problem is to minimize f(x){displaystyle f(mathbf {x} )}, where x{displaystyle mathbf {x} } is a vector in Rn{displaystyle mathbb {R} ^{n}}, and f{displaystyle f} is a differentiable scalar function. There are no constraints on the values that x{displaystyle mathbf {x} } can take.
The algorithm begins at an initial estimate for the optimal value x0{displaystyle mathbf {x} _{0}} 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:
- Bkpk=−∇f(xk),{displaystyle B_{k}mathbf {p} _{k}=-nabla f(mathbf {x} _{k}),}
where Bk{displaystyle B_{k}} is an approximation to the Hessian matrix, which is updated iteratively at each stage, and ∇f(xk){displaystyle nabla f(mathbf {x} _{k})} 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 f(xk+αpk){displaystyle f(mathbf {x} _{k}+alpha mathbf {p} _{k})} over the scalar α>0{displaystyle alpha >0}.
The quasi-Newton condition imposed on the update of Bk{displaystyle B_{k}} is
- Bk+1(xk+1−xk)=∇f(xk+1)−∇f(xk).{displaystyle B_{k+1}(mathbf {x} _{k+1}-mathbf {x} _{k})=nabla f(mathbf {x} _{k+1})-nabla f(mathbf {x} _{k}).}
Let yk=∇f(xk+1)−∇f(xk){displaystyle mathbf {y} _{k}=nabla f(mathbf {x} _{k+1})-nabla f(mathbf {x} _{k})} and sk=xk+1−xk{displaystyle mathbf {s} _{k}=mathbf {x} _{k+1}-mathbf {x} _{k}}, then Bk+1{displaystyle B_{k+1}} satisfies Bk+1sk=yk{displaystyle B_{k+1}mathbf {s} _{k}=mathbf {y} _{k}}, which is the secant equation. The curvature condition sk⊤yk>0{displaystyle mathbf {s} _{k}^{top }mathbf {y} _{k}>0} should be satisfied. If the function is not strongly convex, then the condition has to be enforced explicitly.
Instead of requiring the full Hessian matrix at the point xk+1 to be computed as Bk+1, the approximate Hessian at stage k is updated by the addition of two matrices:
- Bk+1=Bk+Uk+Vk.{displaystyle B_{k+1}=B_{k}+U_{k}+V_{k}.}
Both Uk and Vk 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 Bk+1{displaystyle B_{k+1}}, the update form can be chosen as Bk+1=Bk+αuu⊤+βvv⊤{displaystyle B_{k+1}=B_{k}+alpha mathbf {u} mathbf {u} ^{top }+beta mathbf {v} mathbf {v} ^{top }}. Imposing the secant condition, Bk+1sk=yk{displaystyle B_{k+1}mathbf {s} _{k}=mathbf {y} _{k}}. Choosing u=yk{displaystyle mathbf {u} =mathbf {y} _{k}} and v=Bksk{displaystyle mathbf {v} =B_{k}mathbf {s} _{k}}, we can obtain:[5]
- α=1ykTsk,{displaystyle alpha ={frac {1}{mathbf {y} _{k}^{T}mathbf {s} _{k}}},}
- β=−1skTBksk.{displaystyle beta =-{frac {1}{mathbf {s} _{k}^{T}B_{k}mathbf {s} _{k}}}.}
Finally, we substitute α{displaystyle alpha } and β{displaystyle beta } into Bk+1=Bk+αuu⊤+βvv⊤{displaystyle B_{k+1}=B_{k}+alpha mathbf {u} mathbf {u} ^{top }+beta mathbf {v} mathbf {v} ^{top }} and get the update equation of Bk+1{displaystyle B_{k+1}}:
- Bk+1=Bk+ykykTykTsk−BkskskTBkTskTBksk.{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}}}.}
Algorithm
From an initial guess x0{displaystyle mathbf {x} _{0}} and an approximate Hessian matrix B0{displaystyle B_{0}} the following steps are repeated as xk{displaystyle mathbf {x} _{k}} converges to the solution:
- Obtain a direction pk{displaystyle mathbf {p} _{k}} by solving Bkpk=−∇f(xk){displaystyle B_{k}mathbf {p} _{k}=-nabla f(mathbf {x} _{k})}.
- Perform a one-dimensional optimization (line search) to find an acceptable stepsize αk{displaystyle alpha _{k}} in the direction found in the first step, so αk=argminf(xk+αpk){displaystyle alpha _{k}=arg min f(mathbf {x} _{k}+alpha mathbf {p} _{k})}.
- Set sk=αkpk{displaystyle mathbf {s} _{k}=alpha _{k}mathbf {p} _{k}} and update xk+1=xk+sk{displaystyle mathbf {x} _{k+1}=mathbf {x} _{k}+mathbf {s} _{k}}.
yk=∇f(xk+1)−∇f(xk){displaystyle mathbf {y} _{k}={nabla f(mathbf {x} _{k+1})-nabla f(mathbf {x} _{k})}}.
Bk+1=Bk+ykykTykTsk−BkskskTBkTskTBksk{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}}}}.
f(x){displaystyle f(mathbf {x} )} denotes the objective function to be minimized. Convergence can be checked by observing the norm of the gradient, ||∇f(xk)||{displaystyle ||nabla f(mathbf {x} _{k})||}. If B0{displaystyle B_{0}}is initialized with B0=I{displaystyle B_{0}=I}, the first step will be equivalent to a gradient descent, but further steps are more and more refined by Bk{displaystyle B_{k}}, the approximation to the Hessian.
The first step of the algorithm is carried out using the inverse of the matrix Bk{displaystyle B_{k}}, which can be obtained efficiently by applying the Sherman–Morrison formula to the step 5 of the algorithm, giving
- Bk+1−1=(I−skykTykTsk)Bk−1(I−ykskTykTsk)+skskTykTsk.{displaystyle B_{k+1}^{-1}=left(I-{frac {s_{k}y_{k}^{T}}{y_{k}^{T}s_{k}}}right)B_{k}^{-1}left(I-{frac {y_{k}s_{k}^{T}}{y_{k}^{T}s_{k}}}right)+{frac {s_{k}s_{k}^{T}}{y_{k}^{T}s_{k}}}.}
This can be computed efficiently without temporary matrices, recognizing that Bk−1{displaystyle B_{k}^{-1}} is symmetric,
and that ykTBk−1yk{displaystyle mathbf {y} _{k}^{mathrm {T} }B_{k}^{-1}mathbf {y} _{k}} and skTyk{displaystyle mathbf {s} _{k}^{mathrm {T} }mathbf {y} _{k}} are scalars, using an expansion such as
- Bk+1−1=Bk−1+(skTyk+ykTBk−1yk)(skskT)(skTyk)2−Bk−1ykskT+skykTBk−1skTyk.{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}}}.}
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.
Implementations
The GSL implements BFGS as gsl_multimin_fdfminimizer_vector_bfgs2. Ceres Solver implements both BFGS and L-BFGS.
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.
Octave uses BFGS with a double-dogleg approximation to the cubic line search.
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 the MATLAB Optimization Toolbox, the fminunc function uses BFGS with cubic line search when the problem size is set to "medium scale."
A high-precision arithmetic version of BFGS (pBFGS), implemented in C++ and integrated with the high-precision arithmetic package ARPREC is robust against numerical instability (e.g. round-off errors).
Another C++ implementation of BFGS (along with L-BFGS, L-BFGS-B, CG, and Newton's method) using Eigen (C++ library) is available on GitHub under the MIT License here.
The ensmallen C++ library has implementation of several BFGS variants, provided under the BSD license.
BFGS and L-BFGS are also implemented in C as part of the open-source Gnu Regression, Econometrics and Time-series Library (gretl).
The large scale nonlinear optimization software Artelys Knitro implements, among others, both BFGS and L-BFGS algorithms.
See also
- Quasi-Newton methods
- Davidon–Fletcher–Powell formula
- L-BFGS
- Gradient descent
- Nelder–Mead method
- Levenberg–Marquardt algorithm
- Pattern search (optimization)
- BHHH algorithm
- Symmetric rank-one
Notes
^ Fletcher, Roger (1987), Practical methods of optimization (2nd ed.), New York: John Wiley & Sons, ISBN 978-0-471-91547-8.mw-parser-output cite.citation{font-style:inherit}.mw-parser-output .citation q{quotes:"""""""'""'"}.mw-parser-output .citation .cs1-lock-free a{background:url("//upload.wikimedia.org/wikipedia/commons/thumb/6/65/Lock-green.svg/9px-Lock-green.svg.png")no-repeat;background-position:right .1em center}.mw-parser-output .citation .cs1-lock-limited a,.mw-parser-output .citation .cs1-lock-registration a{background:url("//upload.wikimedia.org/wikipedia/commons/thumb/d/d6/Lock-gray-alt-2.svg/9px-Lock-gray-alt-2.svg.png")no-repeat;background-position:right .1em center}.mw-parser-output .citation .cs1-lock-subscription a{background:url("//upload.wikimedia.org/wikipedia/commons/thumb/a/aa/Lock-red-alt-2.svg/9px-Lock-red-alt-2.svg.png")no-repeat;background-position:right .1em center}.mw-parser-output .cs1-subscription,.mw-parser-output .cs1-registration{color:#555}.mw-parser-output .cs1-subscription span,.mw-parser-output .cs1-registration span{border-bottom:1px dotted;cursor:help}.mw-parser-output .cs1-ws-icon a{background:url("//upload.wikimedia.org/wikipedia/commons/thumb/4/4c/Wikisource-logo.svg/12px-Wikisource-logo.svg.png")no-repeat;background-position:right .1em center}.mw-parser-output code.cs1-code{color:inherit;background:inherit;border:inherit;padding:inherit}.mw-parser-output .cs1-hidden-error{display:none;font-size:100%}.mw-parser-output .cs1-visible-error{font-size:100%}.mw-parser-output .cs1-maint{display:none;color:#33aa33;margin-left:0.3em}.mw-parser-output .cs1-subscription,.mw-parser-output .cs1-registration,.mw-parser-output .cs1-format{font-size:95%}.mw-parser-output .cs1-kern-left,.mw-parser-output .cs1-kern-wl-left{padding-left:0.2em}.mw-parser-output .cs1-kern-right,.mw-parser-output .cs1-kern-wl-right{padding-right:0.2em}
^ Lewis, Adrian S.; Overton, Michael (2009), Nonsmooth optimization via BFGS (PDF)
^ Nocedal & Wright (2006), page 24
^ Byrd, Richard H.; Lu, Peihuang; Nocedal, Jorge; Zhu, Ciyou (1995), "A Limited Memory Algorithm for Bound Constrained Optimization", SIAM Journal on Scientific Computing, 16 (5): 1190–1208, CiteSeerX 10.1.1.645.5814, doi:10.1137/0916069
^ Fletcher, Roger (1987), Practical methods of optimization (2nd ed.), New York: John Wiley & Sons, ISBN 978-0-471-91547-8
Bibliography
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), Numerical optimization: Theoretical and practical aspects, Universitext (Second revised ed. of translation of 1997 French ed.), Berlin: Springer-Verlag, doi:10.1007/978-3-540-35447-5, ISBN 978-3-540-35445-1, MR 2265882
Broyden, C. G. (1970), "The convergence of a class of double-rank minimization algorithms", Journal of the Institute of Mathematics and Its Applications, 6: 76–90, doi:10.1093/imamat/6.1.76
Fletcher, R. (1970), "A New Approach to Variable Metric Algorithms", Computer Journal, 13 (3): 317–322, doi:10.1093/comjnl/13.3.317
Fletcher, Roger (1987), Practical methods of optimization (2nd ed.), New York: John Wiley & Sons, ISBN 978-0-471-91547-8
Goldfarb, D. (1970), "A Family of Variable Metric Updates Derived by Variational Means", Mathematics of Computation, 24 (109): 23–26, doi:10.1090/S0025-5718-1970-0258249-6
Luenberger, David G.; Ye, Yinyu (2008), Linear and nonlinear programming, International Series in Operations Research & Management Science, 116 (Third ed.), New York: Springer, pp. xiv+546, ISBN 978-0-387-74502-2, MR 2423726
Nocedal, Jorge; Wright, Stephen J. (2006), Numerical Optimization (2nd ed.), Berlin, New York: Springer-Verlag, ISBN 978-0-387-30303-1
Shanno, David F. (July 1970), "Conditioning of quasi-Newton methods for function minimization", Mathematics of Computation, 24 (111): 647–656, doi:10.1090/S0025-5718-1970-0274029-X, MR 0274029
Shanno, David F.; Kettler, Paul C. (July 1970), "Optimal conditioning of quasi-Newton methods", Mathematics of Computation, 24 (111): 657–664, doi:10.1090/S0025-5718-1970-0274030-6, MR 0274030
External links
Source code of high-precision BFGS A C++ source code of BFGS with high-precision arithmetic