Computing the Jacobian matrix for x^2*y, 5*x+sin(y)

Question

If I have an R^n --> R^m function, how can I create its Jacobian matrix?

For instance:

expression( x^2*y, 5*x+sin(y) ) # f : R^2 --> R^2

I'd like to have a matrix of expressions like:

expression(2*x*y) expression(x^2)
expression(5) expression(cos(y))

Is there any solution?

Answer 1

At the beginning your question seem a bit vague (although we did offer more specifics later). There is a deriv function that allows simple symbolic differentiation. Since you did not offer any data at first, I assumed symbolic results were you wanted until the second comment.

ex <- expression( x^2*y, 5*x+sin(y) )
sapply(ex, deriv, c("x", "y") )
#-------result--------------------
expression({
    .expr1 <- x^2
    .value <- .expr1 * y
    .grad <- array(0, c(length(.value), 2L), list(NULL, c("x", 
        "y")))
    .grad[, "x"] <- 2 * x * y
    .grad[, "y"] <- .expr1
    attr(.value, "gradient") <- .grad
    .value
}, {
    .value <- 5 * x + sin(y)
    .grad <- array(0, c(length(.value), 2L), list(NULL, c("x", 
        "y")))
    .grad[, "x"] <- 5
    .grad[, "y"] <- cos(y)
    attr(.value, "gradient") <- .grad
    .value
})

Alternatively you could have use this, which would return a list of functions:

sapply(ex, deriv, c("x", "y") , func=TRUE)

[Edit 1: Answering first comment] If you want an object to allow individual matrix functions try this:

res <- sapply(ex, function(ex1) lapply(c("x","y") , 
                            function(arg) deriv(ex1, arg, func=TRUE) ) )

#--------------
> res[1,1]
[[1]]
function (x) 
{
    .value <- x^2 * y
    .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
    .grad[, "x"] <- 2 * x * y
    attr(.value, "gradient") <- .grad
    .value
}

Or for the expression version:

res <- sapply(ex, function(ex1) lapply(c("x","y") , 
                                   function(arg) deriv(ex1, arg) ) )
#----------------------
> res[1,1]
[[1]]
expression({
    .value <- x^2 * y
    .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
    .grad[, "x"] <- 2 * x * y
    attr(.value, "gradient") <- .grad
    .value
})

[Edit 2: answering second comment] To evaluate an expression, you pass it to eval with an appropriate environment:

> eval( res[1,1][[1]], envir=list(x=4,y=5) )
[1] 80
attr(,"gradient")
      x
[1,] 40

There was a niggling wringkle of needing to extract the first (and only) element of the list items in the matrix by first using "[" and then using "[[".