[R] Computing High Order Derivatives (Numerically)

[Gildas Mazo]

Dear Petr Savicky,

this helped indeed. Thank you very much.

Gildas

----- Mail original -----
> De: "Petr Savicky" <savicky at cs.cas.cz>
> À: r-help at r-project.org
> Envoyé: Vendredi 23 Mars 2012 09:39:37
> Objet: Re: [R] Computing High Order Derivatives (Numerically)
> 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.
> 
> ______________________________________________
> 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.

-- 
Gildas Mazo
PhD student
MISTIS team at INRIA
Grenoble, France