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

From: Christos Hatzis <christos.hatzis_at_nuverabio.com>
Date: Thu, 27 Mar 2008 12:36:56 -0400

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



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. Received on Thu 27 Mar 2008 - 16:35:31 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 27 Mar 2008 - 21:30:25 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive