Embedded Systems September 2000 Vol13_10

Issue link:

Contents of this Issue


Page 47 of 229

PROGRAMMER'S TOOLBOX e The bottom line is this: we want a method that combines the nicest fea- tures of both bisection and quadratic fit. Brent's method is such a method, but certain details, such as the halt criterion and the possibility of using quadratics exclusively, are cause for concern. the golden section. The function uses a new constant, golden, which I've added to my standard constant. cpp: double golden sqrt(5.0)+1.0)/2.0; Picking nits To summarize, the things that bother me about Brent's method are: • It can and will finish exclusively with quadratic fits • We cannot guarantee that the bracketing interval is small In addition to these concerns, I've found two other things to complain about. In both the C and Fortran ver- sions of Brent's method, I see that a small epsilon distance is used-ZEPS in Recipes, EPS in fmin.f. (Mittlemann and Spelluci use a machine-language fun ction to find the machine epsilon. I just stuck a reasonable number in its place.) This parameter is later u ed to set a value for tol3. Brent's method needs this epsilon value becau e it's looking for a certain maximum relative (not absolute) error in x. See tl1e line: tol1=eps*dabs(x)+tol3 If the predicted minimum of x just happens to be exactly zero, and tol3 were zero, Brent's method would fail because of a zero tolerance. Now, I don't claim to be tl1e world's leading expert in nwnerical methods, but I know that this step is just plain wrong. You don't try to use a relative error when the parameter you' re searching for can be zero. The step is also wholly unnecessary. Going in, we are given a span of x defin ed by the two end poin ts of the search range. Therefore we know exactly how wide this range is. We also know to expect to narrow it down only so far: roughly tl1e input range times the square root of the machine precision. Our algorithm can set an absolute tolerance at the begin- ning, and send ZEPS to the showers. The bottom line is thi s: we wa nt a meth od that combines th e ni cest features of both bisection and qua- drati c fit. Bren t ' me th od is such a method , but certain de tail s, such as th e halt crite rion and th e possibili ty of using quadrati cs exclu ive ly, are cause for conce rn . Instead of trying to psychoanalyze Brent and his a lgo- rithm, let's simply write o ne th at uses the same no tio ns, but does things our way. If it turns o ut look- ing very much like Brent's me th od, we should n ot be too surprised . I suspect, however, th at it will also di s- play significan t diffe rences, one of which I just me ntioned . One las t thought before we dive in: it occurrs to me that one reason algorithms like brent.c, fmin.f, and in fact most of the algorithms in Recipes seem so inscrutabl e is that th e ir authors wan ted th em to be time-efficient. At this point, 1 don ' t care abou t th at. I'd like to have numerical effi ciency, but I don 't care about counting clock This is the famous golden ratio, <1>. Numerically, it's approximately 1.618. Subtracting it from 2 gives 0.382, which is the value Press et al. call CGOLD, and Mittlemann and Spelucci call 1 c 1 • Don't be thrown by the fac t that they describe it as "the squared inverse of tl1e golden rati o." <1> shows funny behavior when squared or inverted. T he fun cti on gold_shrinkO in Listing 2 takes one step of th e gold- en section search. By tha t, I mean it begins with three points p0, p 1, and p2, and subdivides the inte rval on the "left," p0 •• p 1• It then re places e ithe r p0 or p1 and p2, and se ts up fo r the next ite rati on. Part of this se tup is to exchange p0 and p2, so we sh ouldn ' t make a ny assumptio n tha t, say, x0 < ~- Wha t eems to be th e "left" side of th e interval can really be le ft or righ t, on alte rnate inte rvals. As long as we' re this close, we cycl es. Therefo re, we can afford to sacrifice clock cycles fo r modular and tran - pare nt code. We ' ll worry about speed anothe r day. With that long-winded preamble, let's start W1iting some code. Getting it done I know I'm going to need some sub- fun ctions that do the kinds of opera- tions we' ll be needing, so I might as well write them first. The code is shown in Listing 2. Function swapO simply swaps two floating-point argu- ments. Press et al. use a macro for tl1i s; we can do tl1e same later. For now, a fun ction is just fine. Function gold() spli ts an interval into two unequal parts, according to 46 SEPTEMBER 2000 Embedded Systems Programming might as well write a full-blown golden section search routine. That's the one shown in Listing 2 as gold_searchO. As you can ee, fun ction gold_shrinkO does almost all of the work, so gold_searchO only needs to set things up properly, and decide when to quit. Note my use of an absolute error as a halt criterion. Figure 1 shows the performance of the algoritl1m. Note that it took 28 fun ction evaluations for a 104 improve- ment in accuracy, just about exactly what I predicted (this count includes sen1p evaluations of fi x)). The slope is also about six iterations per decade of improvement, again exactly what we might expect. Though this result was obtained fo r a specific problem, you can take the slope value to the bank. fil e

Articles in this issue

Archives of this issue

view archives of EETimes - Embedded Systems September 2000 Vol13_10