[R] vectorized uni-root?
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 + From: Ravi Varadhan ravi.varad...@jhu.edu To: 'ivo.we...@gmail.com' ivo.we...@gmail.com Cc: r-help@r-project.org r-help@r-project.org Subject: Re: [R] vectorized uni-root? Message-ID: 2f9ea67ef9ae1c48a147cb41be2e15c32e8...@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 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] vectorized uni-root?
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 rvarad...@jhmi.edumailto:rvarad...@jhmi.edu 410-502-2619 [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] vectorized uni-root?
On Thu, Nov 8, 2012 at 3:05 PM, ivo welch ivo.we...@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 __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] vectorized uni-root?
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 - ((fmid0)(flower0))|((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.weyla...@gmail.com wrote:
On Thu, Nov 8, 2012 at 3:05 PM, ivo welch ivo.we...@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
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
