[R] A faster way to compute finite-difference gradient of a scalarfunction of a large number of variables

[Christos Hatzis]

Here is as solution that calculates derivatives using central differences by
appropriately embedding the vectors:

> grad.1
function(x, fn) {
x <- sort(x)
x.e <- head(embed(x, 2), -1)
y.e <- embed(fn(x), 3)
hh <- abs(diff(x.e[1, ]))
    y.e <- apply(y.e, 1, function(z) (z[1] - z[3])/(2 * hh))
cbind(x=x.e[, 1], grad=y.e)
}
> myfunc.1
function(x){
(exp(x) - x) / 10
}
> system.time(g.3 <- grad.1(p0, myfunc.1))
   user  system elapsed 
   0.06    0.00    0.07 
 
CAVEAT: this assumes that the x-values are equally spaced, i.e. same h
throughout, but this part can be easily modified to accommodate h1, h2 on
either side of x.

Also, your myfunc takes a vector input and returns a scalar.  If this is
what was intended, you will need to find a way to vectorize it.  

-Christos

> -----Original Message-----
> From: r-help-bounces at r-project.org 
> [mailto:r-help-bounces at r-project.org] On Behalf Of Ravi Varadhan
> Sent: Thursday, March 27, 2008 12:00 PM
> To: r-help at stat.math.ethz.ch
> Subject: [R] A faster way to compute finite-difference 
> gradient of a scalarfunction of a large number of variables
> 
> Hi All,
> 
>  
> 
> I would like to compute the simple finite-difference 
> approximation to the gradient of a scalar function of a large 
> number of variables (on the order of 1000).  Although a 
> one-time computation using the following function
> grad() is fast and simple enough, the overhead for repeated 
> evaluation of gradient in iterative schemes is quite 
> significant.  I was wondering whether there are better, more 
> efficient ways to approximate the gradient of a large scalar 
> function in R. 
> 
>  
> 
> Here is an example.  
> 
>  
> 
> grad <- function(x, fn=func, eps=1.e-07, ...){
> 
>       npar <- length(x)
> 
>       df <- rep(NA, npar)
> 
>       f <- fn(x, ...)
> 
>         for (i in 1:npar) {
> 
>             dx <- x
> 
>             dx[i] <- dx[i] + eps
> 
>             df[i] <- (fn(dx, ...) - f)/eps
> 
>         }
> 
>       df 
> 
> }
> 
>  
> 
>  
> 
> myfunc <- function(x){
> 
> nvec <- 1: length(x)
> 
> sum(nvec * (exp(x) - x)) / 10
> 
> }
> 
>  
> 
> myfunc.g <- function(x){
> 
> nvec <- 1: length(x)
> 
> nvec * (exp(x) - 1) / 10
> 
> }
> 
>  
> 
> p0 <- rexp(1000)
> 
> system.time(g.1 <- grad(x=p0, fn=myfunc))[1]
> 
> system.time(g.2 <- myfunc.g(x=p0))[1]
> 
> max(abs(g.2 - g.1))
> 
>  
> 
>  
> 
> Thanks in advance for any help or hints.
> 
>  
> 
> Ravi.
> 
> --------------------------------------------------------------
> --------------
> -------
> 
> Ravi Varadhan, Ph.D.
> 
> Assistant Professor, The Center on Aging and Health
> 
> Division of Geriatric Medicine and Gerontology 
> 
> Johns Hopkins University
> 
> Ph: (410) 502-2619
> 
> Fax: (410) 614-9625
> 
> Email: rvaradhan at jhmi.edu
> 
> Webpage:  
> http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
> 
>  
> 
> --------------------------------------------------------------
> --------------
> --------
> 
>  
> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> 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.
> 
>