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

1 2 3 4 5 |

`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 |

`method.args` |
arguments passed to method. See |

`...` |
an additional arguments passed to |

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.

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

.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | ```
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)
# 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)
all.equal(hessc, jacc, tolerance = .Machine$double.eps)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.