From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>

Date: Thu, 27 Mar 2008 17:23:46 +0000 (GMT)

Date: Thu, 27 Mar 2008 17:23:46 +0000 (GMT)

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: rvaradhan_at_jhmi.edu
**>
**> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
**>
**>
**>
**> ----------------------------------------------------------------------------
**> --------
**>
**>
**> -----Original Message-----
**> From: Dimitris Rizopoulos [mailto:dimitris.rizopoulos_at_med.kuleuven.be]
**> Sent: Thursday, March 27, 2008 12:55 PM
**> To: christos.hatzis_at_nuverabio.com; 'Ravi Varadhan'
**> Cc: r-help_at_stat.math.ethz.ch
**> 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" <christos.hatzis_at_nuverabio.com>
**> To: "'Ravi Varadhan'" <rvaradhan_at_jhmi.edu>; <r-help_at_stat.math.ethz.ch>
**> 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: 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.
**>>
**>
**>
**> Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
**>
**> ______________________________________________
**> 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.
**>
*

-- Brian D. Ripley, ripley_at_stats.ox.ac.uk 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_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 - 17:32:07 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 - 18:30:24 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.
*