polykin.math.derivatives

jacobian_forward

jacobian_forward(
    f: Callable[[FloatVector], FloatVector],
    x: FloatVector,
    fx: FloatVector | None = None,
    sclx: FloatVector | None = None,
    ndigit: int | 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 number of reliable digits of \(\mathbf{f}\) and the magnitude and scale of each \(\mathbf{x}\) component.

Typically, the first ndigit/2 digits of the Jacobian are accurate.

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 diferentiated. 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
ndigit	Number of reliable base-10 digits in the values returned by f. This parameter is optional when the function is evaluated analytically, but is essential when the function involves numerical procedures (such as root finding or ODE integration). If None, machine precision is assumed. TYPE: int | 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,
    ndigit: int | 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 number of
    reliable digits of $\mathbf{f}$ and the magnitude and scale of each
    $\mathbf{x}$ component.

    Typically, the first `ndigit/2` digits of the Jacobian are accurate.

    **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 diferentiated.
    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`.
    ndigit : int | None
        Number of reliable base-10 digits in the values returned by `f`. This
        parameter is optional when the function is evaluated analytically,
        but is essential when the function involves numerical procedures (such
        as root finding or ODE integration). If `None`, machine precision is
        assumed.

    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 = sclx if sclx is not None else scalex(x)

    η = eps if ndigit is None else 10**(-ndigit)
    h0 = np.sqrt(η)

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

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

    return J