Nash - Compact Numerical Methods for Computers (523163), страница 38
Текст из файла (страница 38)
A Nelder-Mead minimisation procedure (cont.)for i := 1 to n do if 0.1*abs(Bvec[i])>step then step := 0.1*abs(Bvec[i]);writeln(‘Stepsize computed as ‘,step);for j := 2 to n1 do {main loop to build polytope} {STEP 5}begin {STEP 6}action := ‘BUILD‘;for i := 1 to n do P[i,j] := Bvec[i]; {set the parameters}{alternative strategy -- variable step size -- in the build phase}{ step := 0.1*abs(Bvec[j-l])+0.001; }{Note the scaling and avoidance of zero stepsize.}trystep := step; {trial step -- STEP 7}while P[j-1,j]=Bvecu[j-1] dobeginP[j-1,j] := Bvec[j-1]+trystep; trystep := trystep*10;end; (while]size := size+trystep; {to compute a size measure for polytope -- STEP 8}end; {loop on j for parameters}oldsize := size; {to save the size measure -- STEP 9}calcvert := true; {must calculate vertices when polytope is new}shrinkfail := false; {initialize shrink failure flag so we don’t havefalse convergence}repeat {main loop for Nelder-Mead operations -- STEP 10}if calcvert thenbeginfor j := 1 to nl do {compute the function at each vertex}beginif j<>L then {We already have function value for L(owest) vertex.}beginfor i := 1 to n do Bvec[i] := P[i,j]; {get the parameter values}f := fininfn(n,Bvec,Workdata,notcomp); {function calculation}if notcomp then f := big; funcount := funcount+1; P[nl,j] := f;end; {if j<>L clause}end; {loop on j to compute polytope vertices}calcvert := false; {remember to reset flag so we don’t calculatevertices every cycle of algorithm}end; {calculation of vertices}{STEP 11: find the highest and lowest vertices in current polytope}VL := P[nl,L]; {supposedly lowest value}VH := VL; {highest value must hopefully be higher}H := L; {pointer to highest vertex initialized to L}{Now perform the search}for j := 1 to nl dobeginif j<>L thenbeginf := P[nl,j]; {function value at vertex j}if f<VL thenbeginL := j; VL := f; {save new ‘low’}end;if f>VH thenbegin175176Compact numerical methods for computersAlgorithm 19.
A Nelder-Mead minimisation procedure (cont.)H := j; VH := f; {save new ‘high’}end;end; {if j<>L}end; {search for highest and lowest}{STEP 12: test and display current polytope information}if VH>VL+convtol thenbegin {major cycle of the method}str(funcount:5,tstr); {this is purely for orderly display of resultsin aligned columns}writeln(action,tstr,’‘,VH,’ ‘,VL);VN := beta*VL+(1.0-beta)*VH;{interpolate to get “next to highest” function value -- there aremany options here, we have chosen a fairly conservative one.}for i := 1 to n do {compute centroid of all but point H -- STEP 13}begintemp := -P[i,H]; {leave out point H by subtraction}for j := 1 to nl do temp := temp+P[i,j];P[i,C] := temp/n; {centroid parameter i}end, {loop on i for centroid}for i := 1 to n do {compute reflection in Bvec[] -- STEP 14}Bvec[i] := (1.0+alpha)*P[i,C]-alpha*P[i,H];f := fminfn(n,Bvec,Workdata,notcomp); {function value at refln point}if notcomp then f := big; {When function is not computable, a verylarge value is assigned.}funcount := funcount+l; {increment function count}action := ‘REFLECTION ’; {label the action taken}VR := f; {STEP 15: test if extension should be tried}if VR<VL thenbegin {STEP 16: try extension}P[nl,C] := f; {save the function value at reflection point}for i := 1 to n dobeginf := gamma*Bvec[i]+(1-gamma)*P[i,C];P[i,C] := Bvec[i]; {save the reflection point in case we need it}Bvec[i] := f;end; {loop on i for extension point}f := fminfn(n,Bvec,Workdata,notcomp); {function calculation}if notcomp then f := big; funcount := funcount+1;if f<VR then {STEP 17: test extension}begin {STEP 18: save extension point, replacing H}for i := 1 to n do P[i,H] := Bvec[i];P[n1,H] := f; {and its function value}action := ‘EXTENSION’; {change the action label}end {replace H}elsebegin {STEP 19: save reflection point}for i := 1 to n do P[i,H] := P[i,C];P[n1,H] := VR; {save reflection function value}end {save reflection point in H; note action is still reflection}end {try extension}else {reflection point not lower than current lowest point}begin {reduction and shrink -- STEP 20}Direct search methodsAlgorithm 19.
A Nelder-Mead minimisation procedure (cont.)action := ‘HI-REDUCTION‘; {default to hi-side reduction}if VR<VH then {save reflection -- then try reduction on lo-sideif function value not also < VN}begin {STEP 21: replace H with reflection point}for i := 1 to n do P[i,H] := Bvec[i];P[n1,H] := VR; {and save its function value}action := ‘LO-REDUCTION‘; {re-label action taken}end; {R replaces H so reduction on lo-side}{STEP 22: carry out the reduction step}for i := 1 to n do Bvec[i] := (l-beta)*P[i,H]+beta*P[i,C];f := fminfn(n,Bvec,Workdata,notcomp); {function calculation}if notcomp then f := big; funcount := funcount+1;{STEP 23: test reduction point}if f<P[nl,H] then {replace H -- may be old R in this case,so we do not use VH in this comparison}begin {STEP 24: save new point}for i := 1 to n do P[i,H] := Bvec[i];P[nl,H] := f; {and its function value, which may now not bethe highest in polytope}end {replace H}else {not a new point during reduction}{STEP 25: test for failure of all tactics so far to reduce thefunction value.
Note that this cannot be an ‘else’ statementfrom the ‘if VR<VH’ since we have used that statement in STEP21 to save the reflection point as a prelude to lo-reductiontactic, a step which is omitted when we try the hi-reduction,which has failed to reduce the function value.}if VR>=VH then {hi-reduction has failed to find a point lowerthan H, and reflection point was also higher}begin {STEP 26: shrink polytope toward point L}action := ‘SHRINK’;calcvert := true; {must recalculate the vertices after this}size := 0.0;for j := 1 to nl dobeginif j<>L then {ignore the low vertex}for i := l to n dobeginP[i,j] := beta*(P[i,j]-P[i,L])+P[i,L]; {note the form ofexpression used to avoid rounding errors}size := size+abs(P[i,j]-P[i,L]);end; {loop on i and if j<>L}end; {loop on j}if size<oldsize then {STEP 27 -- test if shrink reduced size}begin {the new polytope is ‘smaller*, so we can proceed}shrinkfail := false; {restart after shrink}oldsize := size;endelse {shrink failed -- polytope has not shrunk}begin {STEP 28 -- exit on failure}writeln(‘Polytope size measure not decreased in shrink’);shrinkfail := true;177178Compact numerical methods for computersAlgorithm 19.
A Nelder-Mead minimisation procedure (cont.)end;{else shrink failed}end; {if VR>=VH -- shrink}end; {reduction}end; {if VH>VL+...}{STEP 29 -- end of major cycle of method}until ((VH<=VL+convtol) or shrinkfail );{STEP 30: if progress made, or polytope shrunk successfully, tryanother major cycle from STEP 10}end; {begin minimisation}{STEP 31} {save best parameters and function value found}writeln(‘Exiting from Alg19.pas Nelder Mead polytope minimiser’);writeIn(‘‘,funcount,’ function evaluations used’);Fmin := P[nl,L]; {save best value found}for i := 1 to n do X[i] := P[i,L];if shrinkfail then fail := true;{STEP 32} {exit}end; {alg19.pas == nmmin}14.3.
AN AXIAL SEARCH PROCEDUREThe following procedure is a device designed primarily to generate informationconcerning the shape of the surface S(b) near the minimum, which will be labelledb*. In some instances, where a minimisation algorithm has converged prematurely, the axial search may reveal a point having a lower function value than thatat the supposed minimum. The step taken along each axis is somewhat arbitrary.O’Neill (1971), for instance, in applying the axial search after a Nelder-Meadprogram, uses ±0·1001 times the initial steps used to generate the simplex.However, these latter increments must be supplied by the program user.
I preferto adjust the parameter bi by an increments = e ( |b i | +e )(13.17)where e is some small number such as the square root of the machine precision.Section 18.2 gives a discussion of this choice. Its principal advantage is that theincrement always adjusts the parameter. Alternatively. I have employed theassignment(14.18)bi := bi (1±0·001)unless these fail to change the parameter, in which case I usebi := ±0·001.(14.19)The latter axial search was used in the set of tests used to compare algorithms forfunction minimisation which were included in the first edition and reported in §18.4of both editions.
What follows below reflects our current usage. including somemeasures of the curvature and symmetry of the functional surface near the presumedminimum.Direct search methodsAlgorithm 20. Axial searchprocedure axissrch(n: integer; {the number of parameters in thefunction to be minimised}var Bvec: rvector; {the parameter values oninput (Bvec) and output (X) from minmeth}var Fmin: real; {‘minimum’ function value}var lowerfn: boolean; {set true if lower valueis found during the axial search}Workdata: probdata); {user defined data area}{alg20.pas == axial search verification of function minimumfor a function of n parameters.Note: in this version, function evaluations are not counted.Copyright 1988 J.C.Nash}varcradius, eps, f, fplus, step, temp, tilt : real;i : integer;notcomp : boolean;beginwriteln(‘alg20.pas--axial search’);eps := calceps; {machine precision}eps := sqrt(eps); {take its square root for a step scaling}writeln(‘Axis’:6,’ Stepsize ‘: 14,‘function + ‘: 14,‘function - ‘: 14,’ rad.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
