[R] getting variable length numerical gradient

Randall R Schulz rschulz at sonic.net
Mon Sep 26 15:53:50 CEST 2005


Dimitris,

I'm new to R programming, and I'm trying to learn the proper way to do 
certain things. E.g., I had a piece of code with explicit iteration to 
apply some computations to a vector. It was pretty slow. I found a way 
to utilize R's built-in vectorization and it was sped up considerably.

So I want to ask about the code you supplied. Please see below.

(By the way, this message is best viewed using a mono-spaced font.)


On Sunday 25 September 2005 04:07, Dimitris Rizopoulos wrote:
> maybe you can find the following function useful (any comments are
> greatly appreciated):
>
> fd <- function(x, f, scalar = TRUE, ..., eps =
> sqrt(.Machine$double.neg.eps)){
>     f <- match.fun(f)
>     out <- if(scalar){
>         ...
>     } else{
>         n <- length(x)
>         res <- array(0, c(n, n))
>         f0 <- f(x, ...)
>         ex <- pmax(abs(x), 1)
>         for(i in 1:n){

This (following) statement will create a copy of the entire "x" vector 
on each iteration. It doesn't look like that's what you would want to 
do:

>             x. <- x

The computation described by this statement could be vectorized outside 
the loop:

>             x.[i] <- x[i] + eps * ex[i]

>             res[, i] <- c(f(x., ...) - f0) / (x.[i] - x[i])
>         }
>         res
>     }
>     out
> }

Offhand, I cannot tell for sure if the last line of that loop is 
vectorizable, but I have a hunch it is.

So at a minimum, it seems this fragment of your code:

        for(i in 1:n){
            x. <- x
            x.[i] <- x[i] + eps * ex[i]
            res[, i] <- c(f(x., ...) - f0) / (x.[i] - x[i])
        }

Could be more efficiently and succinctly replaced with this:

        x. <- x + eps * ex
        for (in in 1:n)
            res[, i] <- c(f(x., ...) - f0) / (x.[i] - x[i])


Could your someone else with R programming experience comment?


Thanks.


Randall Schulz




More information about the R-help mailing list