fzero initial bracketing [WAS: bug #31070]

[Carlo de Falco]

Dear all,

This discussion started as bug #31070 in the tracker,
but Jaroslav convinced me that this is not really a bug
but rather a suggestion for improvement so I am moving
the discussion here.

The algorithm used by fzero [1] for finding x s.t. f(x) = 0,
requires that the user provide an interval [a, b] such that f(a) * f(b) < 0

for compatibility with Matlab, though, fzero allows the user to just
provide an initial guess x0 for x and attempts at internally
constructing the initial
bracketing by looking for a and b in some 'neighborhood' of x0

To my knowledge there is no algorithm which is guaranteed to find such
initial bracketing,
so it is possible that the bracketing is not found in which case fzero
will throw an error and die.

Nonetheless, the use of fzero with an initial guess rather than an
initial bracketing as input is quite
common and there is a lot of Matlab code out there that relies on it,
even though other functions
would be probably more appropriate for this purpose.

Indeed, Matlab reportedly seems to be able to construct the initial
bracketing more often than Octave, so code
that relies on this form of fzero will not work out-of-the-box in Octave.
So I believe it makes sense to try and improve this aspect of fsolve
to avoid giving the (false) impression that
Matlab fzero is better than the Octave one.

The current approach in fzero is the following:
- set one end of the interval to x0
- set the other end to one of [0.9*x0, 1.1*x0, x0-1, x0+1, 0.5*x0
1.5*x0, -x0, 2*x0, -10*x0, 10*x0]
- if none of those work, gripe and die

One thing which is counter-intuitive about this approach is that one
guess that is closer
to the the actual solution does not work better than one that is
further away, e.g.:

>> fun = @(x) ((x-20)*40).^2-.1;
>> fzero (fun, 21)
ans =  20.0079056801977
>> fzero (fun, 20+.01)
error: fzero: not a valid initial bracketing
error: called from:
error:   /Users/carlo/Desktop/OF/octave/octave/scripts/optimization/fzero.m
at line 171, column 5


this is because fun has two roots that are closer to each other than
the minimum
resolution allowed by the bracketing heuristics.

What I don't like about this is that the user has no chance to control
that resolution (other than providing
the initial bracketing himself of course).

In my opinion it would be more consistent, for example, to link the
resolution of the initial
bracketing heuristics to the user provided tolerance parameter 'tolx'.

for example the approach in the attached patch tries, in addition to
the current guesses,
the current initial bracketing guesses:

[x0, x0+x0*tolx*2^k], tolx < tolx*2^k < 1
[x0-tolx*2^k, x0],  tolx < tolx*2^k < 1

i.e. by bisecting the interval [0, 2*x0] down to a minimum size of x0*tolx.

Applied to the example above this gives:

>> fzero (fun, 20+.01)
ans =  20.0079057081480

Of course the patch is very ugly, it's just meant to clarify what I
have in mind,
if you agree that I try to improve this behavior I'll make a proper changeset.

Sorry for writing such a long email about such a minor issue,
all the best,
Carlo


[1] Alefeld, G. E. and Potra, F. A. and Shi, Yixun
     Algorithm 748: enclosing zeros of continuous functions
     ACM TOMS 21 (3), 1995
     doi: http://doi.acm.org/10.1145/210089.210111


---------
diff -r 9f45b76c16e3 scripts/optimization/fzero.m
--- a/scripts/optimization/fzero.m      Sun Sep 19 18:50:30 2010 -0400
+++ b/scripts/optimization/fzero.m      Mon Sep 20 17:09:27 2010 +0200
@@ -148,12 +148,24 @@
     else
       aa = a;
     endif
+    found = false;
     for b = [0.9*aa, 1.1*aa, aa-1, aa+1, 0.5*aa 1.5*aa, -aa, 2*aa,
-10*aa, 10*aa]
       fb = fun (b); nfev += 1;
       if (sign (fa) * sign (fb) <= 0)
+        found = true;
         break;
       endif
     endfor
+    if (! found)
+      dx = aa*.5.^(floor(-log2 (tolx)):-1:0);
+      for b = [(aa+dx);(aa-dx)](:)'
+        fb = fun (b); nfev += 1;
+        if (sign (fa) * sign (fb) <= 0)
+          found = true;
+          break;
+        endif
+      endfor
+    endif
   endif

   if (b < a)
-----------------