R help - Mar 2008 - A faster way to compute finite-difference gradient of a scalar function 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@jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 

----------------------------------------------------------------------------
--------

 


	[[alternative HTML version deleted]]

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