polykin.math.derivatives

jacobian_forward

jacobian_forward(
    f: Callable[[FloatVector], FloatVector],
    x: FloatVector,
    *,
    fx: FloatVector | None = None,
    sclx: FloatVector | None = None,
    epsf: float | None = None
) -> FloatMatrix

Calculate the numerical Jacobian of a vector function \(\mathbf{f}(\mathbf{x})\) using the forward finite-difference scheme.

\[ J_{ij} = \frac{\partial f_i}{\partial x_j} = \frac{f_i(\mathbf{x} + \mathbf{e}_j h_j) - f_i(\mathbf{x})}{h_j} \]

The step size \(h_j\) is optimally determined according to the machine precision of the function values. Typically, the Jacobian is accurate to about half the number of reliable digits returned by the function.

If the function value at \(\mathbf{x}\) is provided, \(N\) function evaluations are required to compute the Jacobian, where \(N\) is the dimension of \(\mathbf{x}\).

References

• J.E. Dennis Jr., R.B. Schnabel, "Numerical Methods for Unconstrained Optimization and Nonlinear Equations", SIAM, 1996, p. 314.

PARAMETER	DESCRIPTION
f	Function to be differentiated. TYPE: Callable[[FloatVector], FloatVector]
x	Differentiation point. TYPE: FloatVector
fx	Function values at x, if available. TYPE: FloatVector | None DEFAULT: None
sclx	Scaling factors for x. Ideally, x[i]*sclx[i] is close to 1. By default, the factors are set internally based on the magnitudes of x. TYPE: FloatVector | None DEFAULT: None
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
FloatMatrix	Jacobian matrix.

Examples:

Evaluate the numerical jacobian of f(x) = x1**2 * x2**3 at (2, -2).

>>> from polykin.math import jacobian_forward
>>> import numpy as np
>>> def f(x): return x[0]**2 * x[1]**3
>>> jacobian_forward(f, np.array([2.0, -2.0]))
array([[-32.00000024,  47.99999928]])

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

def jacobian_forward(
    f: Callable[[FloatVector], FloatVector],
    x: FloatVector,
    *,
    fx: FloatVector | None = None,
    sclx: FloatVector | None = None,
    epsf: float | None = None,
) -> FloatMatrix:
    r"""Calculate the numerical Jacobian of a vector function
    $\mathbf{f}(\mathbf{x})$ using the forward finite-difference scheme.

    $$ J_{ij} = \frac{\partial f_i}{\partial x_j}
              = \frac{f_i(\mathbf{x} + \mathbf{e}_j h_j) - f_i(\mathbf{x})}{h_j} $$

    The step size $h_j$ is optimally determined according to the machine precision of the
    function values. Typically, the Jacobian is accurate to about half the number of
    reliable digits returned by the function.

    If the function value at $\mathbf{x}$ is provided, $N$ function evaluations are
    required to compute the Jacobian, where $N$ is the dimension of $\mathbf{x}$.

    **References**

    *   J.E. Dennis Jr., R.B. Schnabel, "Numerical Methods for Unconstrained Optimization
        and Nonlinear Equations", SIAM, 1996, p. 314.

    Parameters
    ----------
    f : Callable[[FloatVector], FloatVector]
        Function to be differentiated.
    x : FloatVector
        Differentiation point.
    fx : FloatVector | None
        Function values at `x`, if available.
    sclx : FloatVector | None
        Scaling factors for `x`. Ideally, `x[i]*sclx[i]` is close to 1. By
        default, the factors are set internally based on the magnitudes of `x`.
    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
    -------
    FloatMatrix
        Jacobian matrix.

    Examples
    --------
    Evaluate the numerical jacobian of f(x) = x1**2 * x2**3 at (2, -2).
    >>> from polykin.math import jacobian_forward
    >>> import numpy as np
    >>> def f(x): return x[0]**2 * x[1]**3
    >>> jacobian_forward(f, np.array([2.0, -2.0]))
    array([[-32.00000024,  47.99999928]])
    """
    fx = fx if fx is not None else f(x)
    sclx = np.abs(sclx) if sclx is not None else scalex(x)

    epsf = max(epsf, eps) if epsf is not None else eps
    h0 = math.sqrt(epsf)

    J = np.empty((fx.size, x.size))
    xh = x.copy()

    for i in range(x.size):
        h = h0 * max(abs(x[i]), 1 / sclx[i])
        xtemp = xh[i]
        xh[i] += h
        h = xh[i] - xtemp
        J[:, i] = (f(xh) - fx) / h
        xh[i] = xtemp

    return J