Calculate Hessian Matrix

Calculate a numerical approximation to the Hessian matrix of a function at a parameter value.

hessian(func, x, method="Richardson", method.args=list(), ...)

    # Default S3 method
hessian(func, x, method="Richardson",
        method.args=list(), ...)

Arguments

func

a function for which the first (vector) argument is used as a parameter vector.

x

the parameter vector first argument to func.

method

one of "Richardson" or "complex" indicating the method to use for the approximation.

method.args

arguments passed to method. See grad. (Arguments not specified remain with their default values.)

...

an additional arguments passed to func. WARNING: None of these should have names matching other arguments of this function.

Value

An n by n matrix of the Hessian of the function calculated at the point x.

Details

The function hessian calculates an numerical approximation to the n x n second derivative of a scalar real valued function with n-vector argument.

The argument method can be "Richardson" or "complex". Method "simple" is not supported.

For method "complex" the Hessian matrix is calculated as the Jacobian of the gradient. The function grad with method "complex" is used, and method.args is ignored for this (an eps of .Machine$double.eps is used). However, jacobian is used in the second step, with method "Richardson" and argument method.args is used for this. The default is method.args=list(eps=1e-4, d=0.1, zero.tol=sqrt(.Machine$double.eps/7e-7), r=4, v=2, show.details=FALSE). (These are the defaults for hessian with method "Richardson", which are slightly different from the defaults for jacobian with method "Richardson".) See addition comments in grad before choosing method "complex".

Methods "Richardson" uses genD and extracts the second derivative. For this method method.args=list(eps=1e-4, d=0.1, zero.tol=sqrt(.Machine$double.eps/7e-7), r=4, v=2, show.details=FALSE) is set as the default. hessian does one evaluation of func in order to do some error checking before calling genD, so the number of function evaluations will be one more than indicated for genD.

The argument side is not supported for second derivatives and since ... are passed to func there may be no error message if it is specified.

See also

jacobian, grad, genD

Examples

  sc2.f <- function(x){
    n <- length(x)
    sum((1:n) * (exp(x) - x)) / n
    }

  sc2.g <- function(x){
    n <- length(x)
    (1:n) * (exp(x) - 1) / n
    }

  x0 <- rnorm(5)
  hess <- hessian(func=sc2.f, x=x0)
  hessc <- hessian(func=sc2.f, x=x0, "complex")
  all.equal(hess, hessc, tolerance = .Machine$double.eps)
#> [1] "Mean relative difference: 7.495812e-11"
  
#  Hessian = Jacobian of the gradient
  jac  <- jacobian(func=sc2.g, x=x0)
  jacc <- jacobian(func=sc2.g, x=x0, "complex")
  all.equal(hess, jac, tolerance = .Machine$double.eps)
#> [1] "Mean relative difference: 8.21872e-11"
  all.equal(hessc, jacc, tolerance = .Machine$double.eps)
#> [1] "Mean relative difference: 1.390129e-14"