[R] Problems about Derivaties

[Gabor Grothendieck]

On 3/20/07, enrico.foscolo <enrico.foscolo at libero.it> wrote:
> Dear participants to the list,
>
> this is my problem: I want to obtain an expression that represents the second
> derivative of one function.
> With "deriv3" (package "stats") it is possible to evaluate the second
> derivative, but I do not know how I can get the (analytical) expression of this
> derivative.
>
> For example: Suppose that I have a function of this form:
> f(x,y)=x^3+y^3+(x^2)*(y^2). With "deriv3" I can evaluate the first derivative
> and the hessian matrix, as follows:
>
> > d<-deriv3(~x^3+y^3+(x^2)*(y^2),c("x","y"))
> > d[[1]]
> {
> .expr4 <- x^2
> .expr5 <- y^2
> .expr9 <- 2 * x
> .expr15 <- 2 * y
> .value <- x^3 + y^3 + .expr4 * .expr5
> .grad <- array(0, c(length(.value), 2), list(NULL, c("x",
> "y")))
> .hessian <- array(0, c(length(.value), 2, 2), list(NULL,
> c("x", "y"), c("x", "y")))
> .grad[, "x"] <- 3 * .expr4 + .expr9 * .expr5
> .hessian[, "x", "x"] <- 3 * .expr9 + 2 * .expr5
> .hessian[, "x", "y"] <- .hessian[, "y", "x"] <- .expr9 *
> .expr15
> .grad[, "y"] <- 3 * .expr5 + .expr4 * .expr15
> .hessian[, "y", "y"] <- 3 * .expr15 + .expr4 * 2
> attr(.value, "gradient") <- .grad
> attr(.value, "hessian") <- .hessian
> .value
> }
> > d[[1]][11]
> .hessian[, "x", "y"] <- .hessian[, "y", "x"] <- .expr9 * .expr15()
> > typeof(d[[1]][11])
> [1] "language"
> > is.call(d[[1]][11])
> [1] TRUE
>
> Is it possible to extract the formula of the second derivative of this function
> from the output of deriv3?
>
> I would like to use this expression as new function, which means that in this
> case, I would get the string ".hessian[, "x", "y"] <- .hessian[, "y", "x"] <-
> .expr9 * .expr15".
>

I am not really sure from your description exactly what the real problem
is but maybe this will help:

Hessian <- matrix(list(D(D(e, "x"), "x"), D(D(e, "x"), "y"), D(D(e,
"y"), "x"), D(D(e, "y"), "y")), 2)
Hessian[[1,1]]
Hessian[[1,2]]
Hessian[[2,1]]
Hessian[[2,2]]

# or

e <- expression(x^3+y^3+(x^2)*(y^2))
vars <- c("x", "y")
g <- expand.grid(vars, vars)
f <- function(x) D(D(e, x[1]), x[2])
Hessian <- matrix(apply(g, 1, f), 2)
Hessian[[1,1]]
Hessian[[1,2]]
Hessian[[2,1]]
Hessian[[2,2]]

# another possibility is the Ryacas package:

library(Ryacas)
yacas("HessianMatrix(x^3+y^3+(x^2)*(y^2), {x,y})")