[R] Computing High Order Derivatives (Numerically)

[Petr Savicky]

On Fri, Mar 23, 2012 at 12:35:57AM +0100, Gildas Mazo wrote:
> Dear R users,
> 
> Let f be a function over d variables x1,..,xd. I want to compute the k^th-order derivative with respect to x1,..,xk (k<=d). I have a by hand solution (see below) using an iterating code using D. However, I expect d to be high and f to be complicated. Then I want a vector x to be the input, instead of x1,..,xd. How to avoid the x1 <- x[1]; x2 <- x[2], etc steps in the code below? Moreover, D uses symbolic differentation and then eval evaluates the output to get a numerical result. But is there a way to compute the desired derivatives numerically directly (without using symbolic calculus at all)? Finally, what is the most efficient and fast way to get a numerical result for such derivatives?
> 
> Thank you very much in advance,
> Gildas
> 
> ### Code ###
> ### dif takes a function f, an order k, and a vector x as input. f must be a function of x1,..,xd with d >= k. The correspondance is done between xi and x[i]. The expression for f must be at the last row of the body function.
> dif <- function(f,k,x){
>   o <- list()
>   n <- length(body(f))
>   o[[1]] <- body(f)[[n]]
>   for (i in 1:k){
>     xi <- paste("x",i,sep="")
>     o[[i+1]] <- D(o[[i]],name=xi)
>   }
>   x1 <- x[1]
>   x2 <- x[2]
>   x3 <- x[3]
>   eval(o[[k+1]])
> }
> 
> ### Examples ###
> ## function to differentiate
> f <- function(x){
>   x1 <- x[1]
>   x2 <- x[2]
>   x3 <- x[3]
>   0.5*x1*x2*x3^2
> }
> ## derivative w.r.t. x1, x2 and x3 at the point (1,2,3).
> dif(f,3,c(1,2,3))
> 
> ### My Questions ###
> ## how to avoid to write by hand xi <- x[i] ??
> ## is there a way in R to compute such derivatives without using symbolic calculation but numerical compuation instead.

Hi.

For the first question, try the following

  dif <- function(f,k,x){
    o <- list()
    n <- length(body(f))
    o[[1]] <- body(f)[[n]]
    for (i in 1:k){
      xi <- paste("x",i,sep="")
      o[[i+1]] <- D(o[[i]],name=xi)
      assign(xi, x[i])
    }
    eval(o[[k+1]])
  }

For the second question, try the following.

  x <- c(1, 2, 3)
  k <- length(x)
  grid <- as.matrix(expand.grid(rep(list(c(0, 1)), times=k)))
  signs <- 1 - 2*(rowSums(1 - grid) %% 2)
  for (eps in 2^-(5:20)) {
      xeps <- eps*grid + rep(x, each=nrow(grid))
      print(sum(signs*apply(xeps, 1, FUN=f))/eps^k)
  }

  [1] 3.015625
  [1] 3.007812
  [1] 3.003906
  [1] 3.001953
  [1] 3.000977
  [1] 3.000488
  [1] 3.000244
  [1] 3.000122
  [1] 3
  [1] 3
  [1] 3
  [1] 3
  [1] 4
  [1] 0
  [1] 0
  [1] 0

If the above is computed in an exact arithmetic, then
with "eps" converging to zero, the result converges to
the required derivative. Since the numerical computations
are done with a rounding error, too small "eps" yields
a completely wrong result. The choice of a good "eps"
depends on the function and on "k". For a high "k", there
may even be no good "eps". See the considerations at

  http://en.wikipedia.org/wiki/Numerical_derivative

where the choice of "eps" is discussed in the simplest
case of a univariate function.

Hope this helps.

Petr Savicky.