Skip to content

getting variable length numerical gradient

6 messages · Antonio, Fabio Di Narzo, Dimitris Rizopoulos, Randall R Schulz

Original
Antonio, Fabio Di Narzo
# 
Hi all.
I have a numerical function f(x), with x being a vector of generic
size (say k=4), and I wanna take the numerically computed gradient,
using deriv or numericDeriv (or something else).

My difficulties here are that in deriv and numericDeric the function
is passed as an expression, and one have to pass the list of variables
involved as a char vector... So, it's a pure R programming question.


Have a nice sunday,
Antonio, Fabio Di Narzo.
Dimitris Rizopoulos
# 
maybe you can find the following function useful (any comments are 
greatly appreciated):

fd <- function(x, f, scalar = TRUE, ..., eps = 
sqrt(.Machine$double.neg.eps)){
    f <- match.fun(f)
    out <- if(scalar){
        if(length(f0 <- f(x, ...)) != length(x))
            stop("'f' must be vectorized")
        x. <- x + eps * pmax(abs(x), 1)
        c(f(x., ...) - f0) / (x. - x)
    } else{
        n <- length(x)
        res <- array(0, c(n, n))
        f0 <- f(x, ...)
        ex <- pmax(abs(x), 1)
        for(i in 1:n){
            x. <- x
            x.[i] <- x[i] + eps * ex[i]
            res[, i] <- c(f(x., ...) - f0) / (x.[i] - x[i])
        }
        res
    }
    out
}


## Examples

x <- seq(-3.3, 3.3, 0.1)
all.equal(fd(x, pnorm, mean = 0.5), dnorm(x, mean = 0.5))


# Approximate the Hessian matrix for a logistic regression

# the score vector function
gn <- function(b, y, X){
    p <- as.vector(plogis(X %*% b))
    -colSums(X * (y - p))
}

# We simulate some data and fit the logistic regression
n <- 800
x1 <- runif(n,-3, 3); x2 <- runif(n, -3, 3)
pr <- plogis(0.8 + 0.4 * x1 - 0.3 * x2)
y <- rbinom(n, 1, pr)
fm <- glm(y ~ x1 + x2, binomial)

## The Hessian using forward difference approximation
fd(fm$coef, gn, scalar = FALSE, y = y, X = cbind(1, x1, x2))

## The true Hessian
solve(summary(fm)$cov.unscaled)


I hope it helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
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://www.med.kuleuven.be/biostat/
     http://www.student.kuleuven.be/~m0390867/dimitris.htm


----- Original Message ----- 
From: "Antonio, Fabio Di Narzo" <antonio.fabio at gmail.com>
To: <R-help at stat.math.ethz.ch>
Sent: Sunday, September 25, 2005 11:37 AM
Subject: [R] getting variable length numerical gradient
Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
Antonio, Fabio Di Narzo
# 
Tnx very much Dimitris,
your code does what I need. I've just adapted it to my needs (e.g., I
don't deal with scalar functions), and so solved my problem.

Given this, is there a way to use the deriv function in the base
package, within this context (variable length vector of indipendent
variables)?

Best,
Antonio, Fabio Di Narzo.
On 9/25/05, Dimitris Rizopoulos <dimitris.rizopoulos at med.kuleuven.be> wrote:
Randall R Schulz
# 
Dimitris,

I'm new to R programming, and I'm trying to learn the proper way to do 
certain things. E.g., I had a piece of code with explicit iteration to 
apply some computations to a vector. It was pretty slow. I found a way 
to utilize R's built-in vectorization and it was sped up considerably.

So I want to ask about the code you supplied. Please see below.

(By the way, this message is best viewed using a mono-spaced font.)
On Sunday 25 September 2005 04:07, Dimitris Rizopoulos wrote:
This (following) statement will create a copy of the entire "x" vector 
on each iteration. It doesn't look like that's what you would want to 
do:
The computation described by this statement could be vectorized outside 
the loop:
Offhand, I cannot tell for sure if the last line of that loop is 
vectorizable, but I have a hunch it is.

So at a minimum, it seems this fragment of your code:

?? ?? ?? ?? for(i in 1:n){
?? ?? ?? ?? ?? ?? x. <- x
?? ?? ?? ?? ?? ?? x.[i] <- x[i] + eps * ex[i]
?? ?? ?? ?? ?? ?? res[, i] <- c(f(x., ...) - f0) / (x.[i] - x[i])
?? ?? ?? ?? }

Could be more efficiently and succinctly replaced with this:

        x. <- x + eps * ex
        for (in in 1:n)
            res[, i] <- c(f(x., ...) - f0) / (x.[i] - x[i])


Could your someone else with R programming experience comment?


Thanks.


Randall Schulz
Dimitris Rizopoulos
# 
Randall,

thanks for your comments; however, you have to take into account what 
is the purpose of the function here! The goal is to approximate 
*partial* derivatives numerically, using in fact the definition of the 
partial derivatives. If you recall this definition I hope that you can 
see why I change the ith element of the x vector and not the whole 
one. You could also test your approach with the original one in the 
logistic regression example and see the difference.

I hope it is more clear now.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
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://www.med.kuleuven.be/biostat/
     http://www.student.kuleuven.be/~m0390867/dimitris.htm


----- Original Message ----- 
From: "Randall R Schulz" <rschulz at sonic.net>
To: "R Help" <R-Help at stat.math.ethz.ch>
Sent: Monday, September 26, 2005 3:53 PM
Subject: Re: [R] getting variable length numerical gradient


Dimitris,

I'm new to R programming, and I'm trying to learn the proper way to do
certain things. E.g., I had a piece of code with explicit iteration to
apply some computations to a vector. It was pretty slow. I found a way
to utilize R's built-in vectorization and it was sped up considerably.

So I want to ask about the code you supplied. Please see below.

(By the way, this message is best viewed using a mono-spaced font.)
On Sunday 25 September 2005 04:07, Dimitris Rizopoulos wrote:
This (following) statement will create a copy of the entire "x" vector
on each iteration. It doesn't look like that's what you would want to
do:
The computation described by this statement could be vectorized 
outside
the loop:
Offhand, I cannot tell for sure if the last line of that loop is
vectorizable, but I have a hunch it is.

So at a minimum, it seems this fragment of your code:

for(i in 1:n){
x. <- x
x.[i] <- x[i] + eps * ex[i]
res[, i] <- c(f(x., ...) - f0) / (x.[i] - x[i])
}

Could be more efficiently and succinctly replaced with this:

        x. <- x + eps * ex
        for (in in 1:n)
            res[, i] <- c(f(x., ...) - f0) / (x.[i] - x[i])


Could your someone else with R programming experience comment?


Thanks.


Randall Schulz

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! 
http://www.R-project.org/posting-guide.html


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
Randall R Schulz
# 
Dimitris,
On Monday 26 September 2005 07:16, Dimitris Rizopoulos wrote:
Yes. Thanks for clearing up my misunderstanding.
Randall Schulz