Skip to content

vectorized uni-root?

3 messages · ivo welch, R. Michael Weylandt

ivo welch
# 
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?

/iaw
R. Michael Weylandt
#
On Thu, Nov 8, 2012 at 3:05 PM, ivo welch <ivo.welch at gmail.com> wrote:
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
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: