R help - Apr 2011 - uniroot speed and vectorization?

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) )

I would have timed this, but I would have used a 'for' loop, which is
probably not the "R way" of doing this.  has someone already written a
package that does this?

/iaw

----
Ivo Welch (ivo.welch at brown.edu, ivo.welch at gmail.com)

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.