f10-7

418 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).do12j=1,ncom df1dim=df1dim+df(j)*xicom(j) enddo 12 return END CITED REFERENCES AND FURTHER READING: Polak, E. 1971, Computational Methods inOptimization (NewYork: Academic Press), §2.3. [1] Jacobs,D.A.H.(ed.)1977, TheStateoftheArtinNumericalAnalysis (London:AcademicPress), Chapter III.1.7 (by K.W. Brodlie). [2] Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §8.7. 10.7 Variable Metric Methods in Multidimensions Thegoalof variablemetric methods,whicharesometimescalled quasi-Newton methods,isnotdifferentfromthegoalofconjugategradientmethods: toaccumulate information from successive line minimizations so that Nsuch line minimizations lead to the exact minimum of a quadratic form in Ndimensions. In that case, the methodwill also be quadraticallyconvergentformoregeneralsmoothfunctions. Bothvariablemetricandconjugategradientmethodsrequirethatyouareableto computeyourfunction’sgradient,orfirst partialderivatives,atarbitrarypoints. The variablemetricapproachdiffersfromtheconjugategradientinthewaythatit stores 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 size N×N. Generally,for any moderate N, this is an entirely trivial disadvantage. Ontheotherhand,thereisnot,asfarasweknow,anyoverwhelmingadvantage thatthevariablemetricmethodsholdovertheconjugategradienttechniques,except perhapsa historicalone. Developedsomewhatearlier,andmorewidelypropagated,thevariablemetricmethodshavebynowdevelopedawiderconstituencyofsatisfied users. Likewise, some fancier implementations of variable metric methods (going beyondthe scope of this book, see below) have been developedto a greater level ofsophistication on issues like the minimization of roundofferror, handling of special conditions,andso on. Wetendtousevariablemetricratherthanconjugategradient, but we have no reason to urge this habit on you. Variablemetricmethodscomeintwomainflavors. Oneisthe Davidon-Fletcher- Powell (DFP) algorithm (sometimes referred to as simply Fletcher-Powell ). The othergoesbythename Broyden-Fletcher-Goldfarb-Shanno(BFGS) .TheBFGSand DFP 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 has becomegenerallyrecognizedthat,empirically,theBFGSschemeissuperiorinthese details. We will implement BFGS in this section. As before, we imagine that our arbitrary function f(x)can be locally approx- imated by the quadratic form of equation (10.6.1). We don’t, however, have any 10.7VariableMetricMethodsinMultidimensions 419Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).information about the values of the quadratic form’s parameters Aandb, except insofar 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 good approximation to the inverse Hessian matrix A −1, that is, to construct a sequence of matrices Hiwith the property, lim i→∞Hi=A−1(10.7.1 ) Even better if the limit is achieved after Niterations instead of ∞. The reason that variable metric methods are sometimes called quasi-Newton methods can now be explained. Consider finding a minimum by using Newton’s method to search for a zero of the gradient of the function. Near the current pointx i, we have to second order f(x)= f(xi)+(x−xi)·∇f(xi)+1 2(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)=0to determine the next iteration point: x−xi=−A−1·∇f(xi)( 10.7.4 ) The left-hand side is the finite step we need take to get to the exact minimum; the right-handside is known once we have accumulated an accurate H≈A−1. The“quasi”inquasi-Newtonisbecausewedon’tusetheactualHessianmatrix off, but instead use our current approximation of it. This is often betterthan usingthetrueHessian. Wecanunderstandthisparadoxicalresultbyconsideringthedescent directions offatx i. These are the directions palong which fdecreases: ∇f·p<0. FortheNewtondirection(10.7.4)tobeadescentdirection,wemusthave ∇f(xi)·(x−xi)=−(x−xi)·A·(x−xi)<0( 10.7.5 ) which is true if Ais positive definite. In general, far from a minimum, we have no guarantee that the Hessian is positive definite. Taking the actual Newton step with the real Hessian can move us to points where the function is increasing in value. Theideabehindquasi-Newtonmethodsistostartwithapositivedefinite,symmetric approximation to A(usually the unit matrix) and build up the approximating H i’s in such a way that the matrix Hiremains positive definite and symmetric. Far from the minimum, this guarantees that we always move in a downhill direction. Close to the minimum, the updating formula approaches the true Hessian and we enjoythe quadratic convergence of Newton’s method. When we are not close enough to the minimum, taking the full Newton step peven with a positive definite Aneed not decrease the function; we may move too far for the quadratic approximation to be valid. All we are guaranteed is that initially fdecreases as we move in the Newton direction. Once again we can use the backtracking strategy described in §9.7 to choose a step along the direction of the Newton step p, but not necessarily all the way. 420 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).We won’t rigorously derive the DFP algorithm for taking HiintoHi+1; you can consult [3]for clear derivations. Following Brodlie (in [2]), we will give the following heuristic motivation of the procedure. Subtractingequation (10.7.4)at xi+1from that same equationat xigives xi+1−xi=A−1·(∇fi+1−∇ fi)( 10.7.6 ) where∇fj≡∇ f(xj). Having madethe step from xitoxi+1, we might reasonably want to require that the new approximation Hi+1satisfy (10.7.6) as if it were actuallyA−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 H i+1 = Hi+correction. What “objects” are around out of which to construct a correction term? Most notable are the two vectors xi+1−xiand∇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 updatingformula ,i s Hi+1 =Hi+(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)(10.7.8 ) where ⊗denotes the “outer” or “direct” product of two vectors, a matrix: The ij componentof u⊗visuivj. (Youmightwanttoverifythat10.7.8doessatisfy10.7.7.) TheBFGSupdatingformula is exactlythesame,but withoneadditionalterm, ··· +[ (∇fi+1−∇ fi)·Hi·(∇fi+1−∇ fi)]u⊗u (10.7.9 ) whereuis defined as the vector u≡(xi+1−xi) (xi+1−xi)·(∇fi+1−∇ fi) −Hi·(∇fi+1−∇ fi) (∇fi+1−∇ fi)·Hi·(∇fi+1−∇ fi)(10.7.10 ) (You might also verify that this satisfies 10.7.7.) You will have to take on faith — or else consult [3]for details of — the “deep” result that equation (10.7.8), with or without (10.7.9), does in fact convergeto A−1 inNsteps, if fis a quadratic form. Herenowistheroutine dfpminthatimplementsthequasi-Newtonmethod,and uses lnsrchfrom§9.7. As mentioned at the end of newtin§9.7, this algorithm can fail if your variables are badly scaled. 10.7VariableMetricMethodsinMultidimensions 421Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE dfpmin(p,n,gtol,iter,fret,func,dfunc) INTEGER iter,n,NMAX,ITMAX REAL fret,gtol,p(n),func,EPS,STPMX,TOLX PARAMETER (NMAX=50,ITMAX=200,STPMX=100.,EPS=3.e-8,TOLX=4.*EPS)EXTERNAL dfunc,func C USES dfunc,func,lnsrch Givenastarting point p(1:n)that isavectoroflength n,theBroyden-Fletcher-Goldfarb- Shanno variant of Davidon-Fletcher-Powell minimization is performed on a function func, usingitsgradientascalculatedbyaroutine dfunc. Theconvergencerequirementonzeroing the gradient is input as gtol. Returned quantities are p(1:n)(the location of the mini- mum), iter(thenumberofiterationsthatwereperformed),and fret(theminimumvalue of thefunction). Theroutine lnsrchiscalledtoperform approximate lineminimizations. Parameters: NMAXisthe maximumanticipated valueof n;ITMAXisthemaximumallowed number of iterations; STPMXis the scaled maximum step length allowed in line searches; TOLXis the convergence criterion on xvalues. INTEGER i,its,j LOGICAL check REAL den,fac,fad,fae,fp,stpmax,sum,sumdg,sumxi,temp,test, * dg(NMAX),g(NMAX),hdg(NMAX),hessin(NMAX,NMAX),* pnew(NMAX),xi(NMAX) fp=func(p) Calculate starting function valueand gradient, call dfunc(p,g)sum=0. do 12i=1,n and initializethe inverse Hessianto the unit matrix. do11j=1,n hessin(i,j)=0. enddo 11 hessin(i,i)=1.xi(i)=-g(i) Initial line direction. sum=sum+p(i)**2 enddo 12 stpmax=STPMX*max(sqrt(sum),float(n)) do27its=1,ITMAX Main loop over the iterations. iter=its call lnsrch(n,p,fp,g,xi,pnew,fret,stpmax,check,func) Thenewfunctionevaluationoccursin lnsrch;savethefunctionvaluein fpforthenext line search. It is usually safe to ignore the value of check. fp=fret do13i=1,n xi(i)=pnew(i)-p(i) Update the line direction, p(i)=pnew(i) and the current point. enddo 13 test=0. Test for convergence on ∆x. do14i=1,n temp=abs(xi(i))/max(abs(p(i)),1.) if(temp.gt.test)test=temp enddo 14 if(test.lt.TOLX)return do15i=1,n Save the old gradient, dg(i)=g(i) enddo 15 call dfunc(p,g) and get the new gradient. test=0. Test for convergence on zero gradient. den=max(fret,1.)do 16i=1,n temp=abs(g(i))*max(abs(p(i)),1.)/den if(temp.gt.test)test=temp enddo 16 if(test.lt.gtol)return do17i=1,n Compute difference of gradients, dg(i)=g(i)-dg(i) enddo 17 do19i=1,n and difference times current matrix. hdg(i)=0. 422 Chapter10. MinimizationorMaximizationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).do18j=1,n hdg(i)=hdg(i)+hessin(i,j)*dg(j) enddo 18 enddo 19 fac=0. Calculate dot products for the denominators. fae=0. sumdg=0.sumxi=0.do 21i=1,n fac=fac+dg(i)*xi(i) fae=fae+dg(i)*hdg(i)sumdg=sumdg+dg(i)**2sumxi=sumxi+xi(i)**2 enddo 21 if(fac.gt.sqrt(EPS*sumdg*sumxi))then Skip update if facnot sufficiently positive. fac=1./fac fad=1./fae do22i=1,n The vector that makes BFGS different from DFP: dg(i)=fac*xi(i)-fad*hdg(i) enddo 22 do24i=1,n The BFGS updating formula: do23j=i,n hessin(i,j)=hessin(i,j)+fac*xi(i)*xi(j) * -fad*hdg(i)*hdg(j)+fae*dg(i)*dg(j) hessin(j,i)=hessin(i,j) enddo 23 enddo 24 endifdo 26i=1,n Now calculate the next direction to go, xi(i)=0. do25j=1,n xi(i)=xi(i)-hessin(i,j)*g(j) enddo 25 enddo 26 enddo 27 and go back for another iteration. pause ’too many iterations in dfpmin’returnEND Quasi-Newton methods like dfpminwork well with the approximate line minimization done by lnsrch. The routines powell(§10.5) and frprmn(§10.6), however, need more accurate line minimization, which is carried out by the routine linmin. Advanced Implementationsof Variable Metric Methods Although rare, it can conceivably happen that roundoff errors cause the matrix Hito become nearly singular or non-positive-definite. This can be serious, because the supposedsearch directions might then not lead downhill, and because nearly singular H i’s tend to give subsequent Hi’s that are also nearly singular. There is a simple fix for this rare problem, the same as was mentioned in §10.4: In case of any doubt, you should restartthe algorithm at the claimed minimum point, and see if it goes anywhere. Simple, but not very elegant. Modern implementations of variable metricmethods deal with the problem in a more sophisticated way. Insteadofbuildingupanapproximationto A −1,itispossibletobuildupanapproximation ofAitself. Then, instead of calculating the left-hand side of (10.7.4) directly, one solves the set of linear equations A·(x−xi)=−∇ f(xi)( 10.7.11 ) At first glance this seems like a bad idea, since solving (10.7.11) is a process of order N3— and anyway, how does this help the roundoff problem? The trick is not to store Abut 10.8LinearProgrammingandtheSimplexMethod 423Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).rather a triangular decomposition of A, itsCholesky decomposition (cf.§2.9). The updating formula used for the Cholesky decomposition of Ais of order N2and can be arranged to guarantee that the matrix remains positive definite and nonsingular, even in the presence offinite roundoff. This method is due to Gill and Murray [1,2]. CITED REFERENCES AND FURTHER READING: Dennis,J.E., andSchnabel,R.B. 1983, NumericalMethods forUnconstrained Optimizationand Nonlinear Equations (Englewood Cliffs, NJ: Prentice-Hall). [1] Jacobs,D.A.H.(ed.)1977, TheStateoftheArtinNumericalAnalysis (London:AcademicPress), Chapter III.1, §§3–6 (by K. W. Brodlie). [2] Polak,E.1971, ComputationalMethodsinOptimization (NewYork:AcademicPress),pp.56ff.[3] Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), pp. 467–468. 10.8 Linear Programming and the Simplex Method The subject of linear programming , sometimes called linear optimization , concernsitselfwiththefollowingproblem: For Nindependentvariables x1,...,x N, maximize the function z=a01x1+a02x2+··· +a0NxN (10.8.1 ) subject to the primary constraints x1≥0,x 2≥0, ... x N≥0( 10.8.2 ) and simultaneously subject to M =m1+m2+m3additional constraints, m1of them of the form ai1x1+ai2x2+··· +aiNxN≤bi (bi≥0) i=1 ,...,m 1 (10.8.3 ) m2of them of the form aj1x1+aj2x2+··· +ajNxN≥bj≥0 j=m1+1 ,...,m 1+m2(10.8.4 ) and m3of them of the form ak1x1+ak2x2+··· +akNxN=bk≥0 k=m1+m2+1 ,...,m 1+m2+m3(10.8.5 ) The various aij’s can have either sign, or be zero. The fact that the b’s must all be nonnegative (as indicated by the final inequality in the above three equations) is a matter of convention only, since you can multiply any contrary inequality by −1. There is no particular significance in the number of constraints Mbeing less than, equal to, or greater than the number of unknowns N.