R Hessian Matrix

Question

I need to create the Hessian matrix of a function given as:

func <- expression(sin(x+y)+cos(x-y))
vars <- c("x", "y")

I need the second order derivatives as expressions too, and I need to evaluate them lot of times, so I made a list of first order derivatives, and a list of list of second order derivatives.

funcD <- lapply(vars, function(v) D(func, v))
funcDD <- list(); for (i in 1:length(vars)) funcDD[[i]] <- lapply(vars, function(v) D(funcD[[i]], v))

So far, it works.

> funcDD
[[1]]
[[1]][[1]]
-(sin(x + y) + cos(x - y))

[[1]][[2]]
-(sin(x + y) - cos(x - y))


[[2]]
[[2]][[1]]
cos(x - y) - sin(x + y)

[[2]][[2]]
-(cos(x - y) + sin(x + y))

Now the questions: How can I create a matrix containing values of the evaluated expressions? Tried outer, didn't work.

> h <- outer(c(1:length(vars)), c(1:length(vars)), function(r, c) eval(funcDD[[r]][[c]], envir = list(x = 1, y = 2)))
Error in funcDD[[r]] : subscript out of bounds

Other question: Is there a more elegant way to store the second order derivative expressions? For instance is it possible to store expressions in a matrix instead of lists of lists?

Third question: Is it possible to get a vector of variables of an expression? Above I used vars <- c("x", "y") which I entered as input manually, is it necessary or is there a "get_variables"-like method?

Answer 1

The answer to question two is ,"mostly yes" and it offers an almost immediate answer to your question:

funcD <- sapply(vars, function(v) D(func, v))
funcDD <- matrix(list(), 2,2)
for (i in 1:length(vars)) 
        funcDD[,i] <- sapply(vars, function(v) D(funcD[[i]], v))
funcDD
#---------
     [,1]       [,2]      
[1,] Expression Expression
[2,] Expression Expression
> funcDD[1,1]
[[1]]
-(sin(x + y) + cos(x - y))

The "mostly" qualification is that one needs to use "list" rather than "expression" as the object type that the matrix is holding. Expressions don't really qualify as atomic objects, and you could easily extract the value and use it as a call, which might even be more convenient than having it as an expression:

> is.expression(funcDD[1,1])
[1] FALSE
> funcDD[1,1][[1]]
-(sin(x + y) + cos(x - y))
> class(funcDD[1,1][[1]])
[1] "call"

Turns out what was wanted was the same structure, so this calls each matrix element with the same specific vector as the evaluation environment and returns them all as a matrix.:

matrix(sapply(funcDD, eval, env=list(x=0, y=pi)), length(vars))
#---------
     [,1] [,2]
[1,]    1   -1
[2,]   -1    1