On Thu, 27 Mar 2008, Ravi Varadhan wrote: > 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!
But there is C code available to do it: no one has remembered what is in 'Writing R Extensions' where it crops up in several places. See numericDeriv in stats, for one piece of C code with an R interface. > > 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. > -- Brian D. Ripley, [EMAIL PROTECTED] Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595 ______________________________________________ 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.
