spectre.numdiff_scipy

Routines for numerical differentiation.

Module Contents

spectre.numdiff_scipy.group_dense(m, n, A)

spectre.numdiff_scipy.group_columns(A, order=0)

Group columns of a 2-D matrix for sparse finite differencing[1].

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

Aarray_like or sparse matrix, shape (m, n)

Matrix of which to group columns.

orderint, 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

groupsndarray 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

[1] (1,2)

A. Curtis, M. J. D. Powell, and J. Reid. On the estimation of sparse jacobian matrices. Journal of the Institute of Mathematics and its Applications, 13:117–120, 1974.

spectre.numdiff_scipy.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

funcallable

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.

x0array_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_stepNone 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_steparray_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.

f0None 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.

boundstuple 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, kwargstuple 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 Press et al.[2]. 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 Curtis et al.[1] 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

[2]

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes: The Art of Scientific Computing. Cambridge University Press, 3 edition, 2007. See Sec. 5.7.

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.])

spectre.numdiff_scipy.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

funcallable

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.

jaccallable

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.

x0array_like of shape (n,) or float

Point at which to estimate the derivatives. Float will be converted to 1-D array.

bounds2-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, kwargstuple 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

accuracyfloat

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