[R] Is there a bisection method in R?

[Matt Shotwell]

I was just reading about the merge sort algorithm last night (BTW, here
is a fun link http://www.youtube.com/watch?v=t8g-iYGHpEA). There are
some interesting similarities in this context. Here's a recursive method
for bisection:

bisectMatt <- function(fn, lo, hi, tol = 1e-7, ...) {

    flo <- fn(lo, ...)
    fhi <- fn(hi, ...)

    if(flo * fhi > 0)
        stop("root is not bracketed by lo and hi")

    mid <- (lo + hi) / 2
    fmid <- fn(mid, ...)
    if(abs(fmid) <= tol || abs(hi-lo) <= tol)
        return(mid)


    if(fmid * fhi > 0)
        return(bisectMatt(fn, lo, mid, tol, ...))

    return(bisectMatt(fn, mid, hi, tol, ...))
} 

# Adapted from Ravi's original
bisectRavi <- function(fn, lo, hi, tol = 1e-7, ...) {

    flo <- fn(lo, ...) 
    fhi <- fn(hi, ...) 

    if (flo * fhi > 0) 
        stop("root is not bracketed by lo and hi")

    chg <- hi - lo
    while (abs(chg) > tol) {
        mid <- (lo + hi) / 2
        fmid <- fn(mid, ...)
        if (abs(fmid) <= tol) break
        if (flo * fmid < 0) hi <- mid 
        if (fhi * fmid < 0) lo <- mid
        chg <- hi - lo
    }

    return(mid)
}

testFn <- function(x, a) exp(-x) - a*x 

> system.time(bM <- bisectMatt(testFn, 0, 2, a=1))
   user  system elapsed 
  0.000   0.000   0.001 

> system.time(bR <- bisectRavi(testFn, 0, 2, a=1))
   user  system elapsed 
  0.000   0.000   0.001 

> bM
[1] 0.5671433

> bR
[1] 0.5671433

Of course, Ravi's version is better for production (and most likely
faster, though not significantly so in this example) because recursion
is more expensive than looping.

-Matt

On Fri, 2010-09-17 at 17:44 -0400, Ravi Varadhan wrote:
> Here is something simple (does not have any checks for bad input), yet
> should be adequate:
> 
> bisect <- function(fn, lower, upper, tol=1.e-07, ...) {
> f.lo <- fn(lower, ...) 
> f.hi <- fn(upper, ...) 
> feval <- 2
> 
> if (f.lo * f.hi > 0) stop("Root is not bracketed in the specified interval
> \n")
> chg <- upper - lower
> 
> while (abs(chg) > tol) {
> 	x.new <- (lower + upper) / 2
> 	f.new <- fn(x.new, ...)
> 	if (abs(f.new) <= tol) break
> 	if (f.lo * f.new < 0) upper <- x.new 
> 	if (f.hi * f.new < 0) lower <- x.new 
> 	chg <- upper - lower
> 	feval <- feval + 1
> }
> list(x = x.new, value = f.new, fevals=feval)
> }
> 
> # An example
> fn1 <- function(x, a) {
> exp(-x) - a*x 
> }
> 
> bisect(fn1, 0, 2, a=1)
>  
> bisect(fn1, 0, 2, a=2)
> 
> 
> Ravi.
> 
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
> Behalf Of Peter Dalgaard
> Sent: Friday, September 17, 2010 4:16 PM
> To: Gregory Gentlemen
> Cc: r-help at r-project.org
> Subject: Re: [R] Is there a bisection method in R?
> 
> On 09/17/2010 09:28 PM, Gregory Gentlemen wrote:
> > If uniroot is not a bisection method, then what function in R does use
> bisection?
> > 
> 
> Why do you assume that there is one? uniroot contains a better algorithm
> for finding bracketed roots.
> 
> It shouldn't be too hard to roll your own if you need one for
> pedagogical purposes.
> 

-- 
Matthew S. Shotwell
Graduate Student 
Division of Biostatistics and Epidemiology
Medical University of South Carolina