polykin.math

hessian2

hessian2(
    f: Callable[[tuple[float, float]], float],
    x: tuple[float, float],
    h: float = 0.0,
) -> Float2x2Matrix

Calculate the numerical Hessian of a scalar function \(f(x,y)\) using the centered finite-difference scheme.

\[ H(x,y)=\begin{bmatrix} \frac{\partial^2f}{\partial x^2} & \frac{\partial^2f}{\partial x \partial y} \\ \frac{\partial^2f}{\partial y \partial x} & \frac{\partial^2f}{\partial y^2} \end{bmatrix} \]

where the partial derivatives are computed using the centered finite-difference schemes:

\[\begin{aligned} \frac{\partial^2f(x,y)}{\partial x^2} &= \frac{f(x+2h,y)-f(x,y)+f(x-2h,y)}{4 h^2} + O(h^2) \\ \frac{\partial^2f(x,y)}{\partial x \partial y} &= \frac{f(x+h,y+h)-f(x+h,y-h)-f(x-h,y+h)+f(x-h,y-h)}{4 h^2} + O(h^2) \end{aligned}\]

Although the matrix only contains 4 elements and is symmetric, a total of 9 function evaluations are performed.

References

• J. Martins and A. Ning. Engineering Design Optimization. Cambridge University Press, 2021.

PARAMETER	DESCRIPTION
f	Function to be diferentiated. TYPE: Callable[[tuple[float, float]], float]
x	Differentiation point. TYPE: tuple[float, float]
h	Finite-difference step. If 0, it will be set to the theoretical optimum value \(h = \sqrt[3]{3\epsilon}\). TYPE: float DEFAULT: 0.0

RETURNS	DESCRIPTION
Float2x2Matrix	Hessian matrix.

Examples:

Evaluate the numerical Hessian of f(x,y)=(x**2)*(y**3) at (2., -2.).

>>> from polykin.math import hessian2
>>> def fnc(x): return x[0]**2 * x[1]**3
>>> hessian2(fnc, (2., -2.))
array([[-15.99999951,  47.99999983],
       [ 47.99999983, -47.99999983]])

Source code in src/polykin/math/derivatives.py

def hessian2(f: Callable[[tuple[float, float]], float],
             x: tuple[float, float],
             h: float = 0.
             ) -> Float2x2Matrix:
    r"""Calculate the numerical Hessian of a scalar function $f(x,y)$ using the
    centered finite-difference scheme.

    $$
    H(x,y)=\begin{bmatrix}
    \frac{\partial^2f}{\partial x^2} & \frac{\partial^2f}{\partial x \partial y} \\ 
    \frac{\partial^2f}{\partial y \partial x} & \frac{\partial^2f}{\partial y^2} 
    \end{bmatrix}
    $$

    where the partial derivatives are computed using the centered
    finite-difference schemes:

    \begin{aligned}
    \frac{\partial^2f(x,y)}{\partial x^2} &= 
            \frac{f(x+2h,y)-f(x,y)+f(x-2h,y)}{4 h^2} + O(h^2) \\
    \frac{\partial^2f(x,y)}{\partial x \partial y} &=
            \frac{f(x+h,y+h)-f(x+h,y-h)-f(x-h,y+h)+f(x-h,y-h)}{4 h^2} + O(h^2)
    \end{aligned}

    Although the matrix only contains 4 elements and is symmetric, a total of
    9 function evaluations are performed.

    **References**

    *   J. Martins and A. Ning. Engineering Design Optimization. Cambridge
    University Press, 2021.

    Parameters
    ----------
    f : Callable[[tuple[float, float]], float]
        Function to be diferentiated.
    x : tuple[float, float]
        Differentiation point.
    h : float
        Finite-difference step. If `0`, it will be set to the theoretical
        optimum value $h = \sqrt[3]{3\epsilon}$.

    Returns
    -------
    Float2x2Matrix
        Hessian matrix.

    Examples
    --------
    Evaluate the numerical Hessian of f(x,y)=(x**2)*(y**3) at (2., -2.).
    >>> from polykin.math import hessian2
    >>> def fnc(x): return x[0]**2 * x[1]**3
    >>> hessian2(fnc, (2., -2.))
    array([[-15.99999951,  47.99999983],
           [ 47.99999983, -47.99999983]])
    """

    x0, x1 = x

    if h == 0:
        h = cbrt(3*eps)  # ~ 1e-5

    h0 = h*(1 + abs(x0))
    h1 = h*(1 + abs(x1))

    H = np.empty((2, 2))
    f0 = f(x)
    H[0, 0] = (f((x0 + 2*h0, x1)) - 2*f0 + f((x0 - 2*h0, x1)))/(4*h0**2)
    H[1, 1] = (f((x0, x1 + 2*h1)) - 2*f0 + f((x0, x1 - 2*h1)))/(4*h1**2)
    H[0, 1] = (f((x0 + h0, x1 + h1)) - f((x0 + h0, x1 - h1)) - f((x0 - h0, x1 + h1))
               + f((x0 - h0, x1 - h1)))/(4*h0*h1)
    H[1, 0] = H[0, 1]

    return H