Re: Exit codes for fzero

[Jaroslav Hajek]

On Wed, Feb 17, 2010 at 2:57 AM, Rik <address@hidden> wrote:
> Jaroslav Hajek wrote:
>
>>>
>>> For a patch, I considered a simple level test, such as fval <
>>> SOME_TOLERANCE or (fb - fa) < SOME_TOLERANCE.  These would catch one type
>>> of singularity where the right- or left-hand limit was undefined and going
>>> to infinity.  However, it wouldn't catch a simple jump discontinuity such
>>> as a step function from -1 to +1.  In that case both limits would exist and
>>> both would be well-bounded but the two limits do not converge on the same
>>> value.   Instead, my patch looks at the numerically calculated slope
>>> between the last two computed points of fzero.  If the slope is extremely
>>> high this *may* indicate a discontinuity of either the infinite or bounded
>>> sort.
>>>
>>> A draft of the patch is attached.  If it looks good and I get no other
>>> comments I will commit it in a few days.
>>>
>>
>> I have some comments. This looks good, but I think you should more
>> carefully consider what an "extremely high" slope is.
>> In particular, I don't like the absolute constant 1e6, because it
>> makes the test scaling-dependent. Yes, absolute scaling independence
>> is impossible due to limited FP range, but 1e6 is way too small. I
>> think it would be better if the slope of the final bracketing were
>> compared to the slope of the initial bracketing.
> There is some scale dependence in the patch.  I used a comparison to
> max(1e6, 0.5/tolX) for the determination of "extremely steep".  This breaks
> the comparison into a fixed zone of 1e6 and a scale dependent zone.
>
> My goal was to catch a jump discontinuity of 1 unit or a delta_y of 1.  The
> delta_x that fzero arrives at is basically the tolerance in x that the user
> has specified.  The slope for any jump discontinuity of this sort is
> therefore 1/delta_x which is why I set one of the bounds as 0.5/tolX.  This
> scales with the user's requirements.  For the default tolerance of 1e-8
> this is 5E7.  I don't know if this is as large as you want, but if it is
> increased it won't catch the jump discontinuity.  For comparison, a
> function diverging to +/- infinity such as tan(x) or 1/x has a slope of
> 2e16 at the default tolerance.
>
> This test can be fooled of course.  The following straight line with large
> enough slope does it.
> [x, fval, info, out] = fzero (@(x) 1e8*x-pi, [-1,1], struct("TolX",1e-8))
>
> which returns an info of -2 with the new patch.  This seems okay to me
> because the test was only meant to be a warning about singularities rather
> than a definitive judgment.  Moreover, there is a rather straightforward
> test to clarify the matter.  I can keep calling fzero with tighter
> tolerances and if fzero continues to report the solution as a singularity
> it very likely is.  Using that technique
>
> [x, fval, info, out] = fzero (@(x) 1e8*x-pi, [-1,1], struct("TolX",1e-9))
>
> properly returns info=1 indicating that it has converged on a true zero.
>
> The problem with the above setting is that low tolerances result in totally
> unacceptable results.  If I'm only looking for the rough location of a zero
> I might use a tolX of 1e-2 which would result in a maximum slope of just 50
> before classifying the result as a singularity.  I can think of lots of
> functions that have slopes above 50.  Thus, for sloppy tolerances I put a
> bound of 1e6 before classifying something as a singularity.  This bound
> won't catch tan(x) at a tolx of 1e-2, but just catches it at a tolx of
> 1e-3.  I'll happily change this if you a better scale-dependent algorithm.
>
>>
>> Btw., Matlab 2007 misclassifies for instance this simple example:
>>
>>>> [x, fval, info, out] = fzero (@(x) sign(x) + (x==0), [-5,4])
> Interesting, the new patched fzero catches this function.
>
>
> --Rik
>

What about the following? Just compare the final slope relatively to
the starting one. This should achieve scaling-independency in both x
and f.

# HG changeset patch
# User Jaroslav Hajek <address@hidden>
# Date 1266405408 -3600
# Node ID 11bdf1fec1aacbdbf5844aa9583888fdcee8e4df
# Parent  83fa590b8a09d9318b3c649901e156717d9a6c2a
[mq]: patch.fzero

diff --git a/scripts/optimization/fzero.m b/scripts/optimization/fzero.m
--- a/scripts/optimization/fzero.m
+++ b/scripts/optimization/fzero.m
@@ -38,6 +38,8 @@
 ## Maximum number of iterations or function evaluations has been exhausted.
 ## @item -1
 ## The algorithm has been terminated from user output function.
+## @item -2
+## The algorithm may have converged to a singular point.
 ## @end itemize
 ## @seealso{optimset, fsolve}
 ## @end deftypefn
@@ -133,6 +135,8 @@
     error ("fzero:bracket", "fzero: not a valid initial bracketing");
   endif

+  slope0 = (fb - fa) / (b - a);
+
   itype = 1;

   if (abs (fa) < abs (fb))
@@ -281,6 +285,13 @@
     endif
   endwhile

+  ## Check solution for a singularity by examining slope
+  if (info == 1)
+    if ((b - a) != 0 && abs ((fb - fa)/(b - a) / slope0) > max (1e6, 0.5/tolx))
+      info = -2;
+    endif
+  endif
+
   output.iterations = niter;
   output.funcCount = nfev;
   output.bracket = [a, b];


this classifies correctly your example as well as the one that fools Matlab.

what do you think?

-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz