[Rd] RFC: a "safe" uniroot() function for future R

[Martin Maechler]

>>>>> Berend Hasselman <bhh at xs4all.nl>
>>>>>     on Mon, 3 Jun 2013 09:01:25 +0200 writes:

    > On 31-05-2013, at 10:34, Martin Maechler
    > <maechler at stat.math.ethz.ch> wrote:

    > …..
    >> From the current feedbacks, I'd come to propose / further
    >> discuss the following issues:
    >> 
    >> 1) the goal is to remain with one function uniroot()
    >> 
    >> 2) Instead of the 'Sig' = "sign(f'(x_0))" {not quite, but
    >> typically} with 4 different value classes, namely NULL,
    >> -1, 0, 1, (where +/- 1 are equivalent to any positive or
    >> negative finite number respectively),
    >> 
    >> we should either use a string with 4 different possible
    >> values or a {logical or NULL}, say 'upcrossing' (which
    >> also gives 4 values, NULL, NA, FALSE, TRUE).
    >> 
    >> 
    >> 3) [I'm not sure about this:] The new default of the
    >> 'Sig' replacement would correspond to the current Sig =
    >> NULL which does extend the search interval when that does
    >> not constitute a sign change.
    >> 
    >> Alternatively, implicitly proposed by Ravi Varadhan, the
    >> default would correspond to Sig = 0, i.e. to the current
    >> uniroot() behavior which signals an error as soon as the
    >> initial interval is not large enough.
    >> 
    >> Further opinions and suggestions for '2)' and '3)' are
    >> still very welcome!

    > I agree with Ravi. Ravi's example shows what can go wrong.
    > Let the default behaviour of a new uniroot() be the
    > current behaviour.  No harm will be done to any existent
    > applications.

    > An option would be to append a suggestion to the current
    > error message: "Consider setting argument Sig to non-zero
    > to expand the range of x" or something similar when Sig =
    > 0 and uniroot can't find a solution.

    > In addition I think that it should be possible to specify
    > an absolute lower and upper limit for x in order to avoid
    > evaluating a function for invalid values of its argument
    > (e.g. log(x) for negative x) Silly example:

    > xstart <- c(100,110)
    > unirootS(f, xstart, Sig=10,trace=2)

    > Berend

A few days before leaving for a conference and vacation, 
I have finally committed to the very first part of Berend's and
other suggestions, and my own plans (mentioned in this thread, a
month ago):

In R-devel:
   unirootS() no longer exists

and the (rendered) help file for  "new uniroot" currently starts with

------------------------------------------------------------ 

One Dimensional Root (Zero) Finding

Description:

     The function ‘uniroot’ searches the interval from ‘lower’ to
     ‘upper’ for a root (i.e., zero) of the function ‘f’ with respect
     to its first argument.

     Setting ‘extendInt’ to a non-‘"no"’ string, means searching for
     the correct ‘interval = c(lower,upper)’ if ‘sign(f(x))’ does not
     satisfy the requirements at the interval end points; see the
     ‘Details’ section.

Usage:

     uniroot(f, interval, ...,
             lower = min(interval), upper = max(interval),
             f.lower = f(lower, ...), f.upper = f(upper, ...),
             extendInt = c("no", "yes", "downX", "upX"), check.conv = FALSE,
             tol = .Machine$double.eps^0.25, maxiter = 1000, trace = 0)
     
Arguments:

       f: the function for which the root is sought.

interval: a vector containing the end-points of the interval to be
          searched for the root.

     ...: additional named or unnamed arguments to be passed to ‘f’

lower, upper: the lower and upper end points of the interval to be
          searched.

f.lower, f.upper: the same as ‘f(upper)’ and ‘f(lower)’, respectively.
          Passing these values from the caller where they are often
          known is more economical as soon as ‘f()’ contains
          non-trivial computations.

extendInt: character string specifying if the interval ‘c(lower,upper)’
          should be extended or directly produce an error when ‘f()’
          does not have differing signs at the endpoints.  The default,
          ‘"no"’, keeps the search interval and hence produces an
          error.

check.conv: logical indicating whether a convergence warning of the
          underlying ‘uniroot’ should be caught as an error and if
          non-convergence in ‘maxiter’ iterations should be an error
          instead of a warning.

     tol: the desired accuracy (convergence tolerance).

 maxiter: the maximum number of iterations.

   trace: integer number; if positive, tracing information is produced.
          Higher values giving more details.

Details:

     Note that arguments after ‘...’ must be matched exactly.

     Either ‘interval’ or both ‘lower’ and ‘upper’ must be specified:
     the upper endpoint must be strictly larger than the lower
     endpoint.  The function values at the endpoints must be of
     opposite signs (or zero), for ‘extendInt="no"’, the default.
     Otherwise, if ‘extendInt="yes"’, the interval is extended on both
     sides, in search of a sign change, i.e., until the search interval
     [l,u] satisfies f(l) * f(u) <= 0.

     If it is _known how_ f changes sign at the root x0, that is, if
     the function is increasing or decreasing there, ‘extendInt’ can
     (and typically should) be specified as ‘"upX"’ (for “upward
     crossing”) or ‘"downX"’, respectively.  Equivalently, define S:=
     +/- 1, to require S = sign(f(x0 + eps)) at the solution.  In that
     case, the search interval [l,u] possibly is extended to be such
     that S * f(l) <= 0 and S * f(u) >= 0.

-------------

As said above, this does not contain all the changes I've
planned, notably the exact interval extension algorithm *will*
change, on one hand according to what Duncan proposed
(economizing), but also in how the internal 'Delta' is
initialized.

You may want to look at the current partly new examples
and also provide your own, on this thread.
I plan to take it up after about July 20.

Martin