c10-7 (Numerical Recipes in C)

10.7 Variable Metric Methods in Multidimensions425*fret=dbrent(ax,xx,bx,f1dim,df1dim,TOL,&xmin);for (j=1;j<=n;j++) {Construct the vector results to return.xi[j] *= xmin;p[j] += xi[j];}free_vector(xicom,1,n);free_vector(pcom,1,n);#include "nrutil.h"extern int ncom;Defined in dlinmin.extern float *pcom,*xicom,(*nrfunc)(float []);extern void (*nrdfun)(float [], float []);float df1dim(float x){int j;float df1=0.0;float *xt,*df;xt=vector(1,ncom);df=vector(1,ncom);for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];(*nrdfun)(xt,df);for (j=1;j<=ncom;j++) df1 += df[j]*xicom[j];free_vector(df,1,ncom);free_vector(xt,1,ncom);return df1;}CITED REFERENCES AND FURTHER READING:Polak, E.

1971, Computational Methods in Optimization (New York: Academic Press), §2.3. [1]Jacobs, D.A.H. (ed.) 1977, The State of the Art in Numerical Analysis (London: AcademicPress), Chapter III.1.7 (by K.W. Brodlie). [2]Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§8.7.10.7 Variable Metric Methods inMultidimensionsThe goal of variable metric methods, which are sometimes called quasi-Newtonmethods, is not different from the goal of conjugate gradient methods: to accumulateinformation from successive line minimizations so that N such line minimizationslead to the exact minimum of a quadratic form in N dimensions. In that case, themethod will also be quadratically convergent for more general smooth functions.Both variable metric and conjugate gradient methods require that you are able tocompute your function’s gradient, or first partial derivatives, at arbitrary points.

Thevariable metric approach differs from the conjugate gradient in the way that it storesSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.

To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).}426Chapter 10.Minimization or Maximization of Functionslim Hi = A−1i→∞(10.7.1)Even better if the limit is achieved after N iterations instead of ∞.The reason that variable metric methods are sometimes called quasi-Newtonmethods can now be explained. Consider finding a minimum by using Newton’smethod to search for a zero of the gradient of the function. Near the current pointxi , we have to second orderf(x) = f(xi ) + (x − xi ) · ∇f(x i ) + 12 (x − xi ) · A · (x − xi )(10.7.2)so∇f(x) = ∇f(xi ) + A · (x − xi )(10.7.3)In Newton’s method we set ∇f(x) = 0 to determine the next iteration point:x − xi = −A−1 · ∇f(x i )(10.7.4)The left-hand side is the finite step we need take to get to the exact minimum; theright-hand side is known once we have accumulated an accurate H ≈ A−1 .The “quasi” in quasi-Newton is because we don’t use the actual Hessian matrixof f, but instead use our current approximation of it.

This is often better thanSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).and updates the information that is accumulated.

Instead of requiring intermediatestorage on the order of N , the number of dimensions, it requires a matrix of sizeN × N . Generally, for any moderate N , this is an entirely trivial disadvantage.On the other hand, there is not, as far as we know, any overwhelming advantagethat the variable metric methods hold over the conjugate gradient techniques, exceptperhaps a historical one. Developed somewhat earlier, and more widely propagated,the variable metric methods have by now developed a wider constituency of satisfiedusers.

Likewise, some fancier implementations of variable metric methods (goingbeyond the scope of this book, see below) have been developed to a greater level ofsophistication on issues like the minimization of roundoff error, handling of specialconditions, and so on. We tend to use variable metric rather than conjugate gradient,but we have no reason to urge this habit on you.Variable metric methods come in two main flavors.

One is the Davidon-FletcherPowell (DFP) algorithm (sometimes referred to as simply Fletcher-Powell). Theother goes by the name Broyden-Fletcher-Goldfarb-Shanno (BFGS). The BFGS andDFP schemes differ only in details of their roundoff error, convergence tolerances,and similar “dirty” issues which are outside of our scope [1,2] . However, it hasbecome generally recognized that, empirically, the BFGS scheme is superior in thesedetails. We will implement BFGS in this section.As before, we imagine that our arbitrary function f(x) can be locally approximated by the quadratic form of equation (10.6.1). We don’t, however, have anyinformation about the values of the quadratic form’s parameters A and b, exceptinsofar as we can glean such information from our function evaluations and lineminimizations.The basic idea of the variable metric method is to build up, iteratively, a goodapproximation to the inverse Hessian matrix A−1 , that is, to construct a sequenceof matrices Hi with the property,10.7 Variable Metric Methods in Multidimensions427using the true Hessian.

We can understand this paradoxical result by considering thedescent directions of f at xi . These are the directions p along which f decreases:∇f · p < 0. For the Newton direction (10.7.4) to be a descent direction, we must have∇f(xi ) · (x − xi ) = −(x − xi ) · A · (x − xi ) < 0(10.7.5)xi+1 − xi = A−1 · (∇fi+1 − ∇fi )(10.7.6)where ∇fj ≡ ∇f(xj ). Having made the step from xi to xi+1 , we might reasonablywant to require that the new approximation Hi+1 satisfy (10.7.6) as if it wereactually A−1 , that is,xi+1 − xi = Hi+1 · (∇fi+1 − ∇fi )(10.7.7)We might also imagine that the updating formula should be of the form Hi+1 =Hi + correction.What “objects” are around out of which to construct a correction term? Mostnotable are the two vectors xi+1 − xi and ∇fi+1 − ∇fi ; and there is also Hi .There are not infinitely many natural ways of making a matrix out of these objects,especially if (10.7.7) must hold! One such way, the DFP updating formula, is(xi+1 − xi ) ⊗ (xi+1 − xi )(xi+1 − xi ) · (∇fi+1 − ∇fi )[Hi · (∇fi+1 − ∇fi )] ⊗ [Hi · (∇fi+1 − ∇fi )]−(∇fi+1 − ∇fi ) · Hi · (∇fi+1 − ∇fi )Hi+1 = Hi +(10.7.8)where ⊗ denotes the “outer” or “direct” product of two vectors, a matrix: The ijcomponent of u⊗v is ui vj .

(You might want to verify that 10.7.8 does satisfy 10.7.7.)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.

To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).that is, A must be positive definite. In general, far from a minimum, we have noguarantee that the Hessian is positive definite. Taking the actual Newton step withthe real Hessian can move us to points where the function is increasing in value.The idea behind quasi-Newton methods is to start with a positive definite, symmetricapproximation to A (usually the unit matrix) and build up the approximating Hi ’sin such a way that the matrix Hi remains positive definite and symmetric.

Far fromthe minimum, this guarantees that we always move in a downhill direction. Close tothe minimum, the updating formula approaches the true Hessian and we enjoy thequadratic convergence of Newton’s method.When we are not close enough to the minimum, taking the full Newton stepp even with a positive definite A need not decrease the function; we may movetoo far for the quadratic approximation to be valid. All we are guaranteed is thatinitially f decreases as we move in the Newton direction. Once again we can usethe backtracking strategy described in §9.7 to choose a step along the direction ofthe Newton step p, but not necessarily all the way.We won’t rigorously derive the DFP algorithm for taking Hi into Hi+1 ; youcan consult [3] for clear derivations.