[R] A faster way to compute finite-difference gradient of ascalarfunction of a large number of variables
Ravi Varadhan
rvaradhan at jhmi.edu
Thu Mar 27 18:09:30 CET 2008
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!
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
More information about the R-help
mailing list