Re: Exit codes for fzero

[Jaroslav Hajek]

On Wed, Feb 10, 2010 at 11:42 PM, Rik <address@hidden> wrote:
> Jaroslav Hajek wrote:
>
>>> One thought I had would be a single sanity check at the end of fzero to
>>> verify that fval is somewhere near zero.  That would easily catch this case
>>> where the end value of 7x10^14 cannot be confused with zero.  Is there
>>> another better way to get an error exit code when the algorithm has failed?
>>>
>>> --Rik
>>>
>>>
>>
>> Hi Rik,
>> your explanation is correct - the algorithm falls back to bisection
>> when inverse interpolation fails. Basically, the function is assumed
>> to be continuous, so based on that assumption the algorithm
>> incorrectly deduces that a root is bracketed (as if the -Inf and Inf
>> were connected).
>> This is where the actual bracketing interval may be useful. After
>> fsolve "converges", you can check out the final bracketing, and if the
>> jump in values is excessive, consider it to be a discontinuity. I
>> don't think it's much useful to put such a check directly into fzero;
>> the best way to react is probably problem-dependent.
>>
>
> My hope was to avoid coding checks around every function call.  When I call
> an ordinary system function I can use the exit code alone to determine if
> something went wrong.  In the current case I would need to code something 
> like:
>
> if ((INFO == 1) && (abs(FEVAL) < SOME_TOLERANCE)) { result was good }
>

You can easily encapsulate this in your own function. You can even
replace the existing fzero, like this:

fzero = @my_fzero;

...

function [z, fval, info, out] = my_fzero (varargin)
[x, fval, info, out] = fzero (varargin {:});
if ((info == 1) && (abs(fval) > SOME_TOLERANCE))
  info = -2;
endif
endfunction

...
replace fzero in current scope
old_fzero = @fzero;
fzero = @(varargin) my_



> If this is the way that everyone should be verifying the results of fzero
> then it makes more sense to locate this code within fzero.m rather than
> duplicated  in hundreds of user scripts.
>

Why everyone? If you know (or assume) that your function is
continuous, you need no check.


> My second rationale is that I want to be able to use fzero as a black box
> without understanding all of its internals.  When investigating new
> functions I won't necessarily know in advance whether a function is
> pathological or not.  Or I might be exploring a parameter space and have
> the same function with random parameters for which only a few values
> produce bad functions.  In this case I would like to just try the function
> and have fzero itself report whether it was successful since I can't guess
> at all of the ways that the algorithm might fail.
>

The problem is that fzero itself can't really tell the truth. Consider
this modification:

octave:5> [x,fval,info,output] = fzero(@(x) ifelse (abs(x-pi) < 1e-15,
1e30*(x-pi), 1./(x-pi)), 3)
x =  3.1416
fval =  8.8818e+14
info =  1
output =
{
  iterations =  85
  funcCount =  87
  bracket =

     3.1416   3.1416

  bracketf =

    -8.8818e+14   8.8818e+14

}

the result looks similar, but in fact this function has a root at x =
pi, and is continuous everywhere. Yes, this is pretty artificial, but
you can construct more realistic examples.


> Maybe I'm wishing for something that just doesn't exist.  Can anyone with
> Matlab run 'fzero(@(x) 1./(x-pi), 3)' and see whether the exit code reports
> success or failure?
>
> --Rik
>

>> [x, f, info] = fzero(@(x) 1./(x-pi), 3)

x =

    3.1416


f =

  -7.5060e+14


info =

    -5

yes, apparently the exit code indicates a problem. Let's try the second example:

>>  [x,fval,info] = fzero(@(x) ifelse (abs(x-pi) < 1e-15, 1e30*(x-pi), 
>> 1./(x-pi)), 3)

x =

    3.1416


fval =

  -7.5060e+14


info =

    -5


It's wrong here. So it's not like there is something magical, probably
just a test like you proposed above.
But I'm not against it; if you want to, propose a patch.

regards

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