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

[Duncan Murdoch]

On 13-05-30 6:49 AM, Martin Maechler wrote:
>>>>>> Duncan Murdoch <murdoch.duncan at gmail.com>
>>>>>>      on Thu, 30 May 2013 05:27:57 -0400 writes:
>
>      > On Thu, May 30, 2013 at 4:18 AM, Martin Maechler <maechler at stat.math.ethz.ch
>      >> wrote:
>
>      >> With main R releases only happening yearly in spring, now is
>      >> good time to consider *and* discuss new features for what we
>      >> often call "R-devel" and more officially is
>      >> R Under development (unstable) (.....) -- "Unsuffered Consequences"
>      >>
>      >> Here is one such example I hereby expose to public scrutiny:
>      >>
>      >> A few minutes ago, I've committed the following to R-devel
>      >> (the 'trunk' in the svn repository for R):
>      >>
>      >> ------------------------------------------------------------------------
>      >> r62834 | maechler | 2013-05-30 10:01:33 +0200 (Thu, 30 May 2013) | 1 line
>      >> Changed paths:
>      >> M doc/NEWS.Rd
>      >> M src/library/stats/NAMESPACE
>      >> M src/library/stats/R/nlm.R
>      >> M src/library/stats/man/uniroot.Rd
>      >> M tests/Examples/stats-Ex.Rout.save
>      >>
>      >> new "safe" uniroot() =: unirootS()  [may change; see R-devel e-mail]
>      >> ------------------------------------------------------------------------
>      >>
>      >> The help file says
>      >>
>      >> ‘unirootS()’ is a “safe” version of ‘uniroot()’, built on
>      >> ‘uniroot()’, also useful as drop-in replacement of ‘uniroot()’ in
>      >> some cases.  “Safe” 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.
>      >>
>      >> We've had this function, called  safeUroot() in our package
>      >> copula for a while now, where an earlier not-exported version
>      >> has been in my package nor1mix even longer.
>      >> When I was tempted to introduce it into yet another CRAN package
>      >> I maintain,  I decided that this may be a sign that such a
>      >> simple [ utility for / generalization of ] uniroot() should
>      >> probably rather be added to R itself.
>      >>
>      >> The function definition, also visible in R's devel.sources, at the bottom
>      >> of
>      >> https://svn.r-project.org/R/trunk/src/library/stats/R/nlm.R ,
>      >> shows that unirootS() is a wrapper for uniroot() and
>      >> is in principle 100% back compatible to uniroot() itself for all
>      >> the cases where f(lower) and f(upper) are of differing sign and
>      >> hence uniroot() does not give a quick error.
>      >> unirootS() just has three new optional arguments, all with their
>      >> defaults set such as to remain uniroot() compatible.
>      >>
>      >> So, one option instead of the currently commited one would be to
>      >> adopt  unirootS() as "the new uniroot()" and rename current
>      >> uniroot to .uniroot() {and still export both}.
>      >>
>
>      > I would probably prefer this.
>
> I did originally, too, but then became less sure about possible CRAN
> checking "fallout"...
>
> Merging the two functions into one, and not keeping the original
> might be even the best solution, also easiest to maintain,
> including documentation.
>
>      >>
>      >> The capital "S" in the function name and the 'Sig' name is of
>      >> course quite a matter of taste, and this case corresponds to my
>      >> taste, but even that is part of the RFC.
>      >>
>      >>
>      >> unirootS <- function(f, interval, ...,
>      >> lower = min(interval), upper = max(interval),
>      >> f.lower = f(lower, ...), f.upper = f(upper, ...),
>      >> Sig = NULL, check.conv = FALSE,
>      >> tol = .Machine$double.eps^0.25, maxiter = 1000, trace  = 0)
>      >>
>
>      > A few comments:
>
> Thanks a lot, Duncan!
>
>      > 1.  I don't think the name "Sig" conveys the meaning of that parameter
>      > well.  If specified, it determines whether the function is increasing
>      > or decreasing at the root, so maybe "Increasing" or "Upcrossing" (with a
>      > logical value, default NA) would be better?
>
> Good point.  Note that it's completely unused if f() changes sign.

I don't think that's true (but maybe you've changed it since I 
downloaded; I'm offline right now).  For example, try

unirootS(function(x) x, c(-1, 1), Sig=-1)

This fails after 1000 iterations, whereas Sig=+1 is fine.

Duncan Murdoch


> In that case and if  f() is "down crossing" at x0 ,
> it would currently *not* warn even if  Sig == 1.
> For me, using  Upcrossing = TRUE (or similar) would rather
> suggest differently, i.e., I'd expect a warning when that
> requirement is not fulfilled.
>
> As you propose a default of NA, we *could* consider to warn in
> such cases, where the user prescribes the slope-sign at the crossing...
>
>      > 2.  In case 2 where the interval is expanded, wouldn't we save a bit of
>      > time by saving the initial values?  E.g. if Sig == 1 so we want an
>      > upcrossing,  but f.lower is positive, shouldn't we set upper to lower as we
>      > expand lower downwards?
>
> good point.
>
>      > 3.  Sometimes a user will want to force the solution to be between lower
>      > and upper, and will want to signal an error if they are not acceptable.
>      > If you do decide to merge this into uniroot that should be an option.
>
> yes, I agree.
> And that's actually the point for discussion if we
> should allow ourself to make this into uniroot():
> The back compatibility is only guaranteed in those case where
> old-uniroot() did *not* fail.
>
>      > 4.  It should count the search for the interval among the iterations, and
>      > quit if it can't find an interval in that time.  For example,
>
>      >   unirootS( function(x) 1, c(0,1) )
>
>      > never terminates.
>
> duh... of course!
> It shows I wrote the function for my own use, where example as the above
> would not happen.
>
> I've committed a simple change {and test} to catch such
> examples.
> To implement your proposal of *counting* the search
> iterations among the total is trivial if we decide to change
> the two functions into one uniroot(),
> whereas with the current  unirootS() -> uniroot(), this would
> need a bit of ugly code bloat... so I'd rather wait how this
> RFC develops.
>
> Martin
>
>
>
>      > Duncan Murdoch
>
>   [...........]
>
>      >> -----------
>      >>
>      >> As said, your comments are very welcome!
>      >> Note that I'm less interested in variations which gain 10-20% in
>      >> speed benchmarks, rather I'd appreciate proposals for changes
>      >> that give a "better" (in your sense) user interface.
>      >>
>      >> Martin Maechler, ETH Zurich
>