In the "powell" routine, p0 is the original point, pN is the point reached after minimizing in N directions, and pE is the extrapolated point (2 pN - p0). f0, fN and fE are the corresponding function values at these points.

When fE < f0, you may still keep the old set of directions if a certain condition is satisfied. In the NR implementation of powell, when this condition is satisfied, in addition to retaining the old direction set, you also continue the search from pN, not from pE, even if fE < fN. Why don't you continue the search from fE? It seems like this would speed things up.

Originally posted by paulp
In the "powell" routine, p0 is the original point, pN is the point reached after minimizing in N directions, and pE is the extrapolated point (2 pN - p0). f0, fN and fE are the corresponding function values at these points.

When fE < f0, you may still keep the old set of directions if a certain condition is satisfied. In the NR implementation of powell, when this condition is satisfied, in addition to retaining the old direction set, you also continue the search from pN, not from pE, even if fE < fN. Why don't you continue the search from fE? It seems like this would speed things up.

This should have been "Why don't you continue the search from pE?" (I'm still wondering why.)

I'm confused about your question. In the code, the call to linmin moves to the minimum of the new direction starting at the point p, and updates the point p to be the location of the new minimum. So you always start the search at the most recently found minimum.

I'm looking at the code from pages 336-337 of Numerical Recipes in Pascal (1989). If this has been changed in a more recent edition of Numerical Recipes in C, that would of course be helpful. The Pascal code inside the main loop goes something like this:

for i := 1 to n do
begin
...
linmin(p, xit^, n, fret);
...
end;

... {termination criterion}

{construct the extrapolated point and the average direction moved.
Then save the old starting point:}

for j := 1 to n do
begin
ptt^[j] := 2.0 * p[j] - pt^[j];
xit^[j] := p[j] - pt^[j];
pt^[j] := p[j];
end;

fptt := fnc(ptt^); {function value at extrapolated point}

if fptt < fp then
begin
t := 2.0 * (fp - 2.0 * fret + fptt) * sqr(fp - fret - del) - del * sqr(fp - fptt);

if t < 0 then
begin
linmin(p, xit^, n, fret);

for j := 1 to n do
xi[j, ibig] := xit^[j];
end;
end;


So in this code, even if fptt < fp, if t < 0, you don't start the next loop from
ptt, but instead start at p, even though fp > fptt.

OK, I understand what you're getting at. It turns out that Powell's method has a certain amount of heuristics, as described in the text. I tried replacing p by ptt every time the extrapolation is used, i.e. before the last call to linmin, as you suggested. I ran the code on a suite of test minimization problems. There were only small differences from the original results, about equally divided between more or fewer function evaluations. So it looks like it doesn't make much of a difference. Good question!

Thank you very much for investigating this! I will try the same experiment with the problems I've been working on. (I probably should have already tried this, but I wasn't confident enough.)