[R] vectorized uni-root?

[ivo welch]

hi michael---this can be done better than with outer().  A vectorized
bisection search that works with an arbitrary number of arguments is

bisection <- function(f, lower, upper, ..., numiter =100, tolerance =
.Machine$double.eps^0.25 ) {
  stopifnot(length(lower) == length(upper))

  flower <- f(lower, ...); fupper <- f(upper, ...)
  for (n in 1:numiter) {
    mid <- (lower+upper)/2
    fmid <- f(mid, ...)
    if (all(abs(fmid) < tolerance)) break
    samesign <- ((fmid<0)&(flower<0))|((fmid>=0)&(flower>=0))
    lower <- ifelse(samesign, mid, lower )
    flower <- ifelse(samesign, fmid, flower )
    upper <- ifelse(!samesign, mid, upper )
    fupper <- ifelse(!samesign, fmid, fupper )
  }
  return(list( mid=mid, fmid=fmid, lower=lower, upper=upper,
flower=flower, fupper=fupper, n=n ))
}

test.function <- function(x, a) (x-a)*(x**2-a)
K <- 5
print( bisection( test.function, lower=rep(-100,K), upper=rep(100,K), a=1:K ))

this is almost a direct generalization of the simpler uniroot()
function that is currently in R.  if any of the R developers is
reading this, maybe they can add a version of this function and
deprecate the non-vectorized version.

/iaw



On Thu, Nov 8, 2012 at 9:53 AM, R. Michael Weylandt
<michael.weylandt at gmail.com> wrote:
> On Thu, Nov 8, 2012 at 3:05 PM, ivo welch <ivo.welch at gmail.com> wrote:
>> dear R experts--- I have (many) unidimensional root problems.  think
>>
>>   loc.of.root <- uniroot( f= function(x,a) log( exp(a) + a) + a,
>> c(.,9e10), a=rnorm(1) ) $root
>>
>> (for some coefficients a, there won't be a solution; for others, it
>> may exceed the domain.  implied volatilities in various Black-Scholes
>> formulas and variant formulas are like this, too.)
>>
>> except I don't need 1 root, but a few million.  to get so many roots,
>> I can use a for loop, or an lapply or mclapply.  alternatively,
>> because f is a vectorized function in both x and a, evaluations for a
>> few million values will be cheap.  so, one could probably write a
>> clever bisection search here.
>>
>> does a "vectorized" uniroot exist already?
>
> Not to my knowledge, but I think you're on track for a nice
> quick-and-dirty if your function is cheap like the above. Make a 2D
> grid of results with outer() and interpolate to find roots as needed.
> For smooth objective functions, it will likely be cheaper to increase
> precision than to loop over optimization routines written in R.
>
> Cheers,
> Michael