Thank you, Dimitris & Christos. Yes, "myfunc" is a scalar function that needs to be minimized over a high-dimensional parameter space. I was afraid that there might be no better way, apart from coding in C. Thanks, Dimitris, for confirming my fear!
Best regards, 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 ---------------------------------------------------------------------------- -------- -----Original Message----- From: Dimitris Rizopoulos [mailto:[EMAIL PROTECTED] Sent: Thursday, March 27, 2008 12:55 PM To: [EMAIL PROTECTED]; 'Ravi Varadhan' Cc: [EMAIL PROTECTED] Subject: Re: [R] A faster way to compute finite-difference gradient of ascalarfunction of a large number of variables If 'myfunc' is a vector function and can be vectorized in R, then it is even faster to use the following: grad.vec <- function(x, fn, ..., eps = sqrt(.Machine$double.neg.eps)){ x1 <- x + eps * pmax(abs(x), 1) x2 <- x - eps * pmax(abs(x), 1) (fn(x1, ...) - fn(x2, ...)) / (x1 - x2) } 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, ])) apply(y.e, 1, function(z) (z[1] - z[3])/(2 * hh)) } myfunc.1 <- function(x){ (exp(x) - x) / 10 } p0 <- rexp(1000) system.time(for(i in 1:500) out1 <- grad.1(p0, myfunc.1)) system.time(for(i in 1:500) out2 <- grad.vec(p0, myfunc.1)) but for scalar functions I don't think that there is any other way than the one Ravi already has used (maybe doing the whole computation in C??). Best, Dimitris ---- Dimitris Rizopoulos Biostatistical Centre School of Public Health Catholic University of Leuven Address: Kapucijnenvoer 35, Leuven, Belgium Tel: +32/(0)16/336899 Fax: +32/(0)16/337015 Web: http://med.kuleuven.be/biostat/ http://www.student.kuleuven.be/~m0390867/dimitris.htm ----- Original Message ----- From: "Christos Hatzis" <[EMAIL PROTECTED]> To: "'Ravi Varadhan'" <[EMAIL PROTECTED]>; <[EMAIL PROTECTED]> Sent: Thursday, March 27, 2008 5:36 PM Subject: Re: [R] A faster way to compute finite-difference gradient of ascalarfunction of a large number of variables > 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. > Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm ______________________________________________ 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.
