ed
        successfully.
    error : float array
        An estimate of the error: the magnitude of the difference between
        the current estimate of the Hessian and the estimate in the
        previous iteration.
    nfev : int array
        The number of points at which `f` was evaluated.

    Each element of an attribute is associated with the corresponding
    element of `ddf`. For instance, element ``[i, j]`` of `nfev` is the
    number of points at which `f` was evaluated for the sake of
    computing element ``[i, j]`` of `ddf`.

See Also
--------
derivative, jacobian

Notes
-----
Suppose we wish to evaluate the Hessian of a function
:math:`f: \mathbf{R}^m \rightarrow \mathbf{R}`, and we assign to variable
``m`` the positive integer value of :math:`m`. If we wish to evaluate
the Hessian at a single point, then:

- argument `x` must be an array of shape ``(m,)``
- argument `f` must be vectorized to accept an array of shape
  ``(m, ...)``. The first axis represents the :math:`m` inputs of
  :math:`f`; the remaining axes indicated by ellipses are for evaluating
  the function at several abscissae in a single call.
- argument `f` must return an array of shape ``(...)``.
- attribute ``dff`` of the result object will be an array of shape ``(m, m)``,
  the Hessian.

This function is also vectorized in the sense that the Hessian can be
evaluated at ``k`` points in a single call. In this case, `x` would be an
array of shape ``(m, k)``, `f` would accept an array of shape
``(m, ...)`` and return an array of shape ``(...)``, and the ``ddf``
attribute of the result would have shape ``(m, m, k)``. Note that the
axis associated with the ``k`` points is included within the axes
denoted by ``(...)``.

Currently, `hessian` is implemented by nesting calls to `jacobian`.
All options passed to `hessian` are used for both the inner and outer
calls with one exception: the `rtol` used in the inner `jacobian` call
is tightened by a factor of 100 with the expectation that the inner
error can be ignored. A consequence is that `rtol` should not be set
less than 100 times the precision of the dtype of `x`; a warning is
emitted otherwise.

References
----------
.. [1] Hessian matrix, *Wikipedia*,
       https://en.wikipedia.org/wiki/Hessian_matrix

Examples
--------
The Rosenbrock function maps from :math:`\mathbf{R}^m \rightarrow \mathbf{R}`;
the SciPy implementation `scipy.optimize.rosen` is vectorized to accept an
array of shape ``(m, ...)`` and return an array of shape ``...``. Suppose we
wish to evaluate the Hessian at ``[0.5, 0.5, 0.5]``.

>>> import numpy as np
>>> from scipy.differentiate import hessian
>>> from scipy.optimize import rosen, rosen_hess
>>> m = 3
>>> x = np.full(m, 0.5)
>>> res = hessian(rosen, x)
>>> ref = rosen_hess(x)  # reference value of the Hessian
>>> np.allclose(res.ddf, ref)
True

`hessian` is vectorized to evaluate the Hessian at multiple points
in a single call.

>>> rng = np.random.default_rng(4589245925010)
>>> x = rng.random((m, 10))
>>> res = hessian(rosen, x)
>>> ref = [rosen_hess(xi) for xi in x.T]
>>> ref = np.moveaxis(ref, 0, -1)
>>> np.allclose(res.ddf, ref)
True

)