spectre.numdiff_scipy ===================== .. py:module:: spectre.numdiff_scipy .. autoapi-nested-parse:: Routines for numerical differentiation. Module Contents --------------- .. py:function:: group_dense(m, n, A) .. py:function:: group_columns(A, order=0) Group columns of a 2-D matrix for sparse finite differencing\ :footcite:p:`Curtis1974`. Two columns are in the same group if in each row at least one of them has zero. A greedy sequential algorithm is used to construct groups. Parameters ---------- A : array_like or sparse matrix, shape (m, n) Matrix of which to group columns. order : int, iterable of int with shape (n,) or None Permutation array which defines the order of columns enumeration. If int or None, a random permutation is used with `order` used as a random seed. Default is 0, that is use a random permutation but guarantee repeatability. Returns ------- groups : ndarray of int, shape (n,) Contains values from 0 to n_groups-1, where n_groups is the number of found groups. Each value ``groups[i]`` is an index of a group to which ith column assigned. The procedure was helpful only if n_groups is significantly less than n. References ---------- .. footbibliography:: .. py:function:: approx_derivative(fun, x0, dofs=None, method='3-point', *args, rel_step=None, abs_step=None, f0=None, bounds=(-np.inf, np.inf), kwargs={}) Compute finite difference approximation of the derivatives of a vector-valued function. If a function maps from R^n to R^m, its derivatives form m-by-n matrix called the Jacobian, where an element (i, j) is a partial derivative of f[i] with respect to x[j]. Parameters ---------- fun : callable Function of which to estimate the derivatives. The argument x passed to this function is ndarray of shape (n,) (never a scalar even if n=1). It must return 1-D array_like of shape (m,) or a scalar. x0 : array_like of shape (n,) or float Point at which to estimate the derivatives. Float will be converted to a 1-D array. method : {'3-point', '2-point', 'cs'}, optional Finite difference method to use: - '2-point' - use the first order accuracy forward or backward difference. - '3-point' - use central difference in interior points and the second order accuracy forward or backward difference near the boundary. - 'cs' - use a complex-step finite difference scheme. This assumes that the user function is real-valued and can be analytically continued to the complex plane. Otherwise, produces bogus results. rel_step : None or array_like, optional Relative step size to use. If None (default) the absolute step size is computed as ``h = rel_step * sign(x0) * max(1, abs(x0))``, with `rel_step` being selected automatically, see Notes. Otherwise ``h = rel_step * sign(x0) * abs(x0)``. For ``method='3-point'`` the sign of `h` is ignored. The calculated step size is possibly adjusted to fit into the bounds. abs_step : array_like, optional Absolute step size to use, possibly adjusted to fit into the bounds. For ``method='3-point'`` the sign of `abs_step` is ignored. By default relative steps are used, only if ``abs_step is not None`` are absolute steps used. f0 : None or array_like, optional If not None it is assumed to be equal to ``fun(x0)``, in this case the ``fun(x0)`` is not called. Default is None. bounds : tuple of array_like, optional Lower and upper bounds on independent variables. Defaults to no bounds. Each bound must match the size of `x0` or be a scalar, in the latter case the bound will be the same for all variables. Use it to limit the range of function evaluation. Bounds checking is not implemented when `as_linear_operator` is True. args, kwargs : tuple and dict, optional Additional arguments passed to `fun`. Both empty by default. The calling signature is ``fun(x, *args, **kwargs)``. Returns ------- J : {ndarray, sparse matrix, LinearOperator} Finite difference approximation of the Jacobian matrix. If `as_linear_operator` is True returns a LinearOperator with shape (m, n). Otherwise it returns a dense array or sparse matrix depending on how `sparsity` is defined. If `sparsity` is None then a ndarray with shape (m, n) is returned. If `sparsity` is not None returns a csr_matrix with shape (m, n). For sparse matrices and linear operators it is always returned as a 2-D structure, for ndarrays, if m=1 it is returned as a 1-D gradient array with shape (n,). See Also -------- check_derivative : Check correctness of a function computing derivatives. Notes ----- If `rel_step` is not provided, it assigned as ``EPS**(1/s)``, where EPS is determined from the smallest floating point dtype of `x0` or `fun(x0)`, ``np.finfo(x0.dtype).eps``, s=2 for '2-point' method and s=3 for '3-point' method. Such relative step approximately minimizes a sum of truncation and round-off errors, see :footcite:t:`Press2007`. Relative steps are used by default. However, absolute steps are used when ``abs_step is not None``. If any of the absolute or relative steps produces an indistinguishable difference from the original `x0`, ``(x0 + dx) - x0 == 0``, then a automatic step size is substituted for that particular entry. A finite difference scheme for '3-point' method is selected automatically. The well-known central difference scheme is used for points sufficiently far from the boundary, and 3-point forward or backward scheme is used for points near the boundary. Both schemes have the second-order accuracy in terms of Taylor expansion. Refer to :footcite:t:`Curtis1974` for the formulas of 3-point forward and backward difference schemes. For dense differencing when m=1 Jacobian is returned with a shape (n,), on the other hand when n=1 Jacobian is returned with a shape (m, 1). Our motivation is the following: a) It handles a case of gradient computation (m=1) in a conventional way. b) It clearly separates these two different cases. b) In all cases np.atleast_2d can be called to get 2-D Jacobian with correct dimensions. References ---------- .. footbibliography:: Examples -------- >>> import numpy as np >>> from scipy.optimize._numdiff import approx_derivative >>> >>> def f(x, c1, c2): ... return np.array([x[0] * np.sin(c1 * x[1]), x[0] * np.cos(c2 * x[1])]) >>> x0 = np.array([1.0, 0.5 * np.pi]) >>> approx_derivative(f, x0, args=(1, 2)) array([[ 1., 0.], [-1., 0.]]) Bounds can be used to limit the region of function evaluation. In the example below we compute left and right derivative at point 1.0. >>> def g(x): ... return x**2 if x >= 1 else x >>> x0 = 1.0 >>> approx_derivative(g, x0, bounds=(-np.inf, 1.0)) array([ 1.]) >>> approx_derivative(g, x0, bounds=(1.0, np.inf)) array([ 2.]) .. py:function:: check_derivative(fun, jac, x0, bounds=(-np.inf, np.inf), args=(), kwargs={}) Check correctness of a function computing derivatives (Jacobian or gradient) by comparison with a finite difference approximation. Parameters ---------- fun : callable Function of which to estimate the derivatives. The argument x passed to this function is ndarray of shape (n,) (never a scalar even if n=1). It must return 1-D array_like of shape (m,) or a scalar. jac : callable Function which computes Jacobian matrix of `fun`. It must work with argument x the same way as `fun`. The return value must be array_like or sparse matrix with an appropriate shape. x0 : array_like of shape (n,) or float Point at which to estimate the derivatives. Float will be converted to 1-D array. bounds : 2-tuple of array_like, optional Lower and upper bounds on independent variables. Defaults to no bounds. Each bound must match the size of `x0` or be a scalar, in the latter case the bound will be the same for all variables. Use it to limit the range of function evaluation. args, kwargs : tuple and dict, optional Additional arguments passed to `fun` and `jac`. Both empty by default. The calling signature is ``fun(x, *args, **kwargs)`` and the same for `jac`. Returns ------- accuracy : float The maximum among all relative errors for elements with absolute values higher than 1 and absolute errors for elements with absolute values less or equal than 1. If `accuracy` is on the order of 1e-6 or lower, then it is likely that your `jac` implementation is correct. See Also -------- approx_derivative : Compute finite difference approximation of derivative. Examples -------- >>> import numpy as np >>> from scipy.optimize._numdiff import check_derivative >>> >>> >>> def f(x, c1, c2): ... return np.array([x[0] * np.sin(c1 * x[1]), x[0] * np.cos(c2 * x[1])]) >>> def jac(x, c1, c2): ... return np.array( ... [ ... [np.sin(c1 * x[1]), c1 * x[0] * np.cos(c1 * x[1])], ... [np.cos(c2 * x[1]), -c2 * x[0] * np.sin(c2 * x[1])], ... ] ... ) >>> >>> x0 = np.array([1.0, 0.5 * np.pi]) >>> check_derivative(f, jac, x0, args=(1, 2)) 2.4492935982947064e-16