From: Petr Savicky <savicky_at_praha1.ff.cuni.cz>

Date: Sat, 02 Apr 2011 23:25:08 +0200

}

print(range(of(x1, a)))

print(range(of(x2, a)))

R-help_at_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. Received on Sat 02 Apr 2011 - 21:32:16 GMT

Date: Sat, 02 Apr 2011 23:25:08 +0200

On Sat, Apr 02, 2011 at 08:24:07AM -0400, ivo welch wrote:

> curiosity---given that vector operations are so much faster than

*> scalar operations, would it make sense to make uniroot vectorized? if
**> I read the uniroot docs correctly, uniroot() calls an external C
**> routine which seems to be a scalar function. that must be slow. I am
**> thinking a vectorized version would be useful for an example such as
**>
**> of <- function(x,a) ( log(x)+x+a )
**> uniroot( of, c( 1e-7, 100 ), a=rnorm(1000000) )
*

Hi.

The slowest part of the solution using uniroot() is the repeated evaluation of the R level input function. If this can be vectorized, then a faster algorithm could be possible. The following is an example of a vectorized bisection, which is simpler and less efficient than "zeroin" used in uniroot().

of <- function(x,a) { log(x)+x+a }

a <- rnorm(1000)

x1 <- rep(1e-7, times=length(a))

x2 <- rep(100, times=length(a))

stopifnot(of(x1, a) < 0)

stopifnot(of(x2, a) > 0)

for (i in 1:60) {

x3 <- (x1 + x2)/2 pos <- of(x3, a) > 0 y1 <- ifelse(pos, x1, x3) y2 <- ifelse(pos, x3, x2) x1 <- y1 x2 <- y2

}

print(range(of(x1, a)))

print(range(of(x2, a)))

It can be more efficient to find approximations of the roots using a few iterations of the above approach and then switch to the Newton method, which can be vectorized easily.

Hope this helps for a start.

Petr Savicky.

R-help_at_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. Received on Sat 02 Apr 2011 - 21:32:16 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Sun 03 Apr 2011 - 01:20:27 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*