[R] A faster way to compute finite-difference gradient of ascalarfunction of a large number of variables

[Dimitris Rizopoulos]

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