[R] A faster way to compute finite-difference gradient of ascalarfunction of a large number of variables
Prof Brian Ripley
ripley at stats.ox.ac.uk
Thu Mar 27 18:23:46 CET 2008
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
More information about the R-help
mailing list