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
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
Hi Ivo, The only problem is that uniroot() actually uses Brent's algorithm, and is based on the C code from netlib (there is also Fortran code available - zeroin.f). Brent's algorithm is more sophisticated than a simple bisection approach that you have vectorized. It combines bisection and secant methods along with some quadratic interpolation. The idea behind this hybrid approach (which by the way is a fundamental theme in all of numerical analysis) is to get faster convergence while not sacrificing the global convergence of bisection. So, the existing uniroot() will not be deprecated unless you can vectorize Brent's hybrid root-finding approach! Best, Ravi Ravi Varadhan, Ph.D. Assistant Professor The Center on Aging and Health Division of Geriatric Medicine & Gerontology Johns Hopkins University rvaradhan@jhmi.edu<mailto:rvaradhan@jhmi.edu> 410-502-2619 [[alternative HTML version deleted]]
The package rootoned on http://r-forge.r-project.org/R/?group_id=395 has an all-R version of zeroin (the algorithm of uniroot). This should also be in Rmpfr by Martin M., as it was set up for that use. I suspect it can be vectorized fairly easily. However, it may be simpler to write, or else abstract from that code, a vectorized false position or secant code (same formula, FP must have end points with opposite sign function values, so is `safer'). This avoids the bisection step and simplifies things, but may be a bit less efficient sometimes. JN Date: Fri, 9 Nov 2012 14:58:55 +0000 From: Ravi Varadhan <ravi.varadhan at jhu.edu> To: "'ivo.welch at gmail.com'" <ivo.welch at gmail.com> Cc: "r-help at r-project.org" <r-help at r-project.org> Subject: Re: [R] vectorized uni-root? Message-ID: <2F9EA67EF9AE1C48A147CB41BE2E15C32E8D67 at DOM-EB-MAIL1.win.ad.jhu.edu> Content-Type: text/plain Hi Ivo, The only problem is that uniroot() actually uses Brent's algorithm, and is based on the C code from netlib (there is also Fortran code available - zeroin.f). Brent's algorithm is more sophisticated than a simple bisection approach that you have vectorized. It combines bisection and secant methods along with some quadratic interpolation. The idea behind this hybrid approach (which by the way is a fundamental theme in all of numerical analysis) is to get faster convergence while not sacrificing the global convergence of bisection. So, the existing uniroot() will not be deprecated unless you can vectorize Brent's hybrid root-finding approach! Best, Ravi