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

From: Dimitris Rizopoulos <dimitris.rizopoulos_at_med.kuleuven.be>
Date: Thu, 27 Mar 2008 17:54:50 +0100

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


> 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. Received on Thu 27 Mar 2008 - 20:50:26 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 - 21:30:27 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.

list of date sections of archive