polykin.math.derivatives

gradient_forward

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

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

\[ \frac{\partial f}{\partial x_i} = \frac{f(\mathbf{x} + \mathbf{e}_i h_i) - f(\mathbf{x})}{h_i} \]

The step size \(h_i\) is optimally determined according to the machine precision of the function values. Typically, the gradient 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 gradient, 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. 323.

PARAMETER	DESCRIPTION
f	Function to be differentiated. TYPE: Callable[[FloatVector], float]
x	Differentiation point. TYPE: FloatVector
fx	Function value at x, if available. TYPE: float | 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
FloatVector	Gradient vector.

Examples:

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

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

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

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

    $$ \frac{\partial f}{\partial x_i} =
       \frac{f(\mathbf{x} + \mathbf{e}_i h_i) - f(\mathbf{x})}{h_i} $$

    The step size $h_i$ is optimally determined according to the machine precision of the
    function values. Typically, the gradient 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 gradient, 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. 323.

    Parameters
    ----------
    f : Callable[[FloatVector], float]
        Function to be differentiated.
    x : FloatVector
        Differentiation point.
    fx : float | None
        Function value 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
    -------
    FloatVector
        Gradient vector.

    Examples
    --------
    Evaluate the numerical gradient of f(x) = x1**2 * x2**3 at (2, -2).
    >>> from polykin.math import gradient_forward
    >>> import numpy as np
    >>> def f(x): return x[0]**2 * x[1]**3
    >>> gradient_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)

    g = np.empty_like(x, dtype=float)
    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
        g[i] = (f(xh) - fx) / h
        xh[i] = xtemp

    return g