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: [EMAIL PROTECTED] > [mailto:[EMAIL PROTECTED] On Behalf Of Ravi Varadhan > Sent: Thursday, March 27, 2008 12:00 PM > To: [EMAIL PROTECTED] > 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: [EMAIL PROTECTED] > > Webpage: > http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html > > > > -------------------------------------------------------------- > -------------- > -------- > > > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@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@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.