Computes the numerical Hessian of `functions`

or the symbolic Hessian of `characters`

.

```
hessian(f, var, params = list(), accuracy = 4, stepsize = NULL, drop = TRUE)
f %hessian% var
```

- f
array of

`characters`

or a`function`

returning a`numeric`

array.- var
vector giving the variable names with respect to which the derivatives are to be computed and/or the point where the derivatives are to be evaluated. See

`derivative`

.- params
`list`

of additional parameters passed to`f`

.- accuracy
degree of accuracy for numerical derivatives.

- stepsize
finite differences stepsize for numerical derivatives. It is based on the precision of the machine by default.

- drop
if

`TRUE`

, return the Hessian as a matrix and not as an`array`

for scalar-valued functions.

Hessian matrix for scalar-valued functions when `drop=TRUE`

, `array`

otherwise.

In Cartesian coordinates, the Hessian of a scalar-valued function \(F\) is the square matrix of second-order partial derivatives:

$$(H(F))_{ij} = \partial_{ij}F$$

When the function \(F\) is a tensor-valued function \(F_{d_1,\dots,d_n}\),
the `hessian`

is computed for each scalar component.

$$(H(F))_{d_1\dots d_n,ij} = \partial_{ij}F_{d_1\dots d_n}$$

It might be tempting to apply the definition of the Hessian as the Jacobian of the gradient to write it in arbitrary orthogonal coordinate systems. However, this results in a Hessian matrix that is not symmetric and ignores the distinction between vector and covectors in tensor analysis. The generalization to arbitrary coordinate system is not currently supported.

`%hessian%`

: binary operator with default parameters.

Guidotti, E. (2020). "calculus: High dimensional numerical and symbolic calculus in R". https://arxiv.org/abs/2101.00086

Other differential operators:
`curl()`

,
`derivative()`

,
`divergence()`

,
`gradient()`

,
`jacobian()`

,
`laplacian()`

```
### symbolic Hessian
hessian("x*y*z", var = c("x", "y", "z"))
#> [,1] [,2] [,3]
#> [1,] "0" "z" "y"
#> [2,] "z" "0" "x"
#> [3,] "y" "x" "0"
### numerical Hessian in (x=1, y=2, z=3)
f <- function(x, y, z) x*y*z
hessian(f = f, var = c(x=1, y=2, z=3))
#> [,1] [,2] [,3]
#> [1,] 1.100056e-10 3.00000e+00 2.000000e+00
#> [2,] 3.000000e+00 2.75014e-11 1.000000e+00
#> [3,] 2.000000e+00 1.00000e+00 1.222284e-11
### vectorized interface
f <- function(x) x[1]*x[2]*x[3]
hessian(f = f, var = c(1, 2, 3))
#> [,1] [,2] [,3]
#> [1,] 1.100056e-10 3.00000e+00 2.000000e+00
#> [2,] 3.000000e+00 2.75014e-11 1.000000e+00
#> [3,] 2.000000e+00 1.00000e+00 1.222284e-11
### symbolic vector-valued functions
f <- c("y*sin(x)", "x*cos(y)")
hessian(f = f, var = c("x","y"))
#> , , 1
#>
#> [,1] [,2]
#> [1,] "-(y * sin(x))" "cos(x)"
#> [2,] "0" "-sin(y)"
#>
#> , , 2
#>
#> [,1] [,2]
#> [1,] "cos(x)" "0"
#> [2,] "-sin(y)" "-(x * cos(y))"
#>
### numerical vector-valued functions
f <- function(x) c(sum(x), prod(x))
hessian(f = f, var = c(0,0,0))
#> , , 1
#>
#> [,1] [,2] [,3]
#> [1,] 1.790455e-13 0 0
#> [2,] 0.000000e+00 0 0
#>
#> , , 2
#>
#> [,1] [,2] [,3]
#> [1,] 0 1.790455e-13 0
#> [2,] 0 0.000000e+00 0
#>
#> , , 3
#>
#> [,1] [,2] [,3]
#> [1,] 0 0 1.790455e-13
#> [2,] 0 0 0.000000e+00
#>
### binary operator
"x*y^2" %hessian% c(x=1, y=3)
#> [,1] [,2]
#> [1,] 0 6
#> [2,] 6 2
```