Quasi-Newton Methods

There are many variants of quasi-Newton methods. In all of them, the idea is to base the matrix (Subscript[B, k]) in the quadratic model

on an approximation of the Hessian matrix built up from the function and gradient values from some or all steps previously taken.

This loads a package that contains some utility functions:

<<Optimization`UnconstrainedProblems`

This shows a plot of the steps taken by the quasi-Newton method. The path is much less direct than for Newton’s method. The quasi-Newton method is used by default by (RefLink[FindMinimum, paclet:ref/FindMinimum]) for problems that are not sums of squares:

FindMinimumPlot[Cos[x^2 - 3 y] + Sin[x^2 + y^2],{{x,1},{y,1}}]

The first thing to notice about the path taken in this example is that it starts in the wrong direction. This direction is chosen because at the first step all the method has to go by is the gradient, and so it takes the direction of steepest descent. However, in subsequent steps, it incorporates information from the values of the function and gradient at the steps taken to build up an approximate model of the Hessian.

The methods used by the Wolfram Language are the Broyden–Fletcher–Goldfarb–Shanno (BFGS) updates and, for large systems, the limited-memory BFGS (L-BFGS) methods, in which the model (Subscript[B, k]) is not stored explicitly, but rather ((Subscript[B, k])^(-1))∇f((Subscript[x, k])) is calculated by gradients and step directions stored from past steps.

The BFGS method is implemented such that instead of forming the model Hessian (Subscript[B, k]) at each step, Cholesky factors (Subscript[L, k]) such that (Subscript[L, k]).((Subscript[L, k])^T)==(Subscript[B, k]) are computed so that only O((n^2)) operations are needed to solve the system (Subscript[B, k])(Subscript[s, k])==-∇f((Subscript[x, k])) [DS96] for a problem with n variables.

For large-scale sparse problems, the BFGS method can be problematic because, in general, the Cholesky factors (or the Hessian approximation (Subscript[B, k]) or its inverse) are dense, so the O((n^2)) memory and operations requirements become prohibitive compared to algorithms that take advantage of sparseness. The L-BFGS algorithm [NW99] forms an approximation to the inverse Hessian based on the last m past steps, which are stored. The Hessian approximation may not be as complete, but the memory and order of operations are limited to O(nm) for a problem with n variables. In the Wolfram Language, for problems over 250 variables, the algorithm is switched automatically to L-BFGS. You can control this with the method option "StepMemory"->m. With m==∞, the full BFGS method will always be used. Choosing an appropriate value of m is a tradeoff between speed of convergence and the work done per step. With m<3, you are most likely better off using a "conjugate gradient" algorithm.

This shows the same example function with the minimum computed using L-BFGS with m==5:

FindMinimumPlot[Cos[x^2 - 3 y] + Sin[x^2 + y^2],{{x,1},{y,1}}, Method->{"QuasiNewton", "StepMemory"->5}]

Quasi-Newton methods are chosen as the default in the Wolfram Language because they are typically quite fast and do not require computation of the Hessian matrix, which can be quite expensive both in terms of the symbolic computation and numerical evaluation. With an adequate line search, they can be shown to converge superlinearly [NW99] to a local minimum where the Hessian is positive definite. This means that

or, in other words, the steps keep getting smaller. However, for very high precision, this does not compare to the q-quadratic convergence rate of Newton's method.

This shows the number of steps and function evaluations required to find the minimum to high precision for the problem shown:

First[Timing[Block[{e = 0, s = 0},{{min, minpoint},{points}}  = Reap[FindMinimum[Cos[x^2 - 3 y] + Sin[x^2 + y^2],{{x,1.},{y,1.}}, Method->"QuasiNewton", WorkingPrecision->10000,StepMonitor:>(s++;Sow[{x,y}]), EvaluationMonitor:>e++]];Print[s," steps and ",e," evaluations"]]]]

Newton's method is able to find ten times as many digits with far fewer steps because of its quadratic convergence rate. However, the convergence with the quasi-Newton method is still superlinear since the ratio of the errors is clearly going to zero.

This makes a plot showing the ratios of the errors in the computation. The ratios of the errors are shown on a logarithmic scale so that the trend can clearly be seen over a large range of magnitudes:

exact = {x, y} /. Last[Solve[{x^2 + y^2 ==3  Pi/2, x^2 - 3 y == -Pi},{x,y}]];errs= Map[Norm[N[exact - #]]&, points];
ListPlot[Log[10, Drop[errs, 1]/Drop[errs, -1]]]

The following table summarizes the options you can use with quasi-Newton methods.

Method options for (RefLink[Method, paclet:ref/Method])->"QuasiNewton".