Is there an R function or package for finding the gradient and Hessian of a vectorized equation?

Question

For example if I want to find the gradient and Hessian with respect to x of:

f = function(x,y,alpha,A,b){
  return((1/n)*(y-alpha*x)%*%(y-alpha*x) + (A%*%x-b)%*%(A%*%x-b))
}

Answer 1

CRAN package Deriv can compute symbolic derivatives of R functions.
In the code below I have removed the call to return in the posted function f.
Function DDeriv is a copy&paste from this RPubs post changed to use Deriv instead of base D that does not accept functions as its first argument.

library(Deriv)

DDeriv <- function(expr, name, order = 1){
  if(order < 1) stop("Order must be >= 1")
  if(order == 1) Deriv(expr, name)
  else DDeriv(Deriv(expr, name), name, order - 1)
}


f <- function(x,y,alpha,A,b){
  (1/n)*(y-alpha*x)%*%(y-alpha*x) + (A%*%x-b)%*%(A%*%x-b)
}

Using function Deriv directly:

Deriv(f, "x")
#function (x, y, alpha, A, b) 
#{
#    .e1 <- -alpha
#    .e3 <- A %*% x - b
#    .e5 <- y - alpha * x
#    (.e1 %*% .e5 + .e5 %*% .e1)/n + .e3 %*% A + A %*% .e3
#}

Using function DDeriv:

DDeriv(f, "x", 1)
#function (x, y, alpha, A, b) 
#{
#    .e1 <- -alpha
#    .e3 <- A %*% x - b
#    .e5 <- y - alpha * x
#    (.e1 %*% .e5 + .e5 %*% .e1)/n + .e3 %*% A + A %*% .e3
#}

DDeriv(f, "x", 2)
#function (x, y, alpha, A, b) 
#{
#    .e1 <- -alpha
#    2 * (.e1 %*% .e1/n) + 2 * A %*% A
#}