polykin.math.derivatives

hessian2_centered

hessian2_centered(
    f: Callable[[tuple[float, float]], float],
    x: tuple[float, float],
    epsf: float | None = None,
) -> 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 differentiated. TYPE: Callable[[tuple[float, float]], float]
x	Differentiation point. TYPE: tuple[float, float]
epsf	Machine precision of the function values. If None, machine precision of 64-bit floating-point type is assumed. If the number of reliable base-10 digits in the results returned by the function is \(n\), then epsf is approximately \(10^{-n}\). TYPE: float | None DEFAULT: None

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_centered
>>> def f(x): return x[0]**2 * x[1]**3
>>> hessian2_centered(f, (2.0, -2.0))
array([[-15.99999951,  47.99999983],
       [ 47.99999983, -47.99999983]])

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

def hessian2_centered(
    f: Callable[[tuple[float, float]], float],
    x: tuple[float, float],
    epsf: float | None = None,
) -> 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 differentiated.
    x : tuple[float, float]
        Differentiation point.
    epsf : float | None
        Machine precision of the function values. If `None`, machine precision of 64-bit
        floating-point type is assumed. If the number of reliable base-10 digits in the
        results returned by the function is $n$, then `epsf` is approximately $10^{-n}$.

    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_centered
    >>> def f(x): return x[0]**2 * x[1]**3
    >>> hessian2_centered(f, (2.0, -2.0))
    array([[-15.99999951,  47.99999983],
           [ 47.99999983, -47.99999983]])
    """
    x0, x1 = x

    epsf = max(epsf, eps) if epsf is not None else eps
    h = math.cbrt(3 * epsf)
    h0 = h * max(1.0, abs(x0))
    h1 = h * max(1.0, abs(x1))

    H = np.empty((2, 2))
    fx = f(x)
    H[0, 0] = (f((x0 + 2 * h0, x1)) - 2 * fx + f((x0 - 2 * h0, x1))) / (4 * h0**2)
    H[1, 1] = (f((x0, x1 + 2 * h1)) - 2 * fx + 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