Re: Exit codes for fzero

[Rik]

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