Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_numdiff.py: 6%
278 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Routines for numerical differentiation."""
2import functools
3import numpy as np
4from numpy.linalg import norm
6from scipy.sparse.linalg import LinearOperator
7from ..sparse import issparse, csc_matrix, csr_matrix, coo_matrix, find
8from ._group_columns import group_dense, group_sparse
11def _adjust_scheme_to_bounds(x0, h, num_steps, scheme, lb, ub):
12 """Adjust final difference scheme to the presence of bounds.
14 Parameters
15 ----------
16 x0 : ndarray, shape (n,)
17 Point at which we wish to estimate derivative.
18 h : ndarray, shape (n,)
19 Desired absolute finite difference steps.
20 num_steps : int
21 Number of `h` steps in one direction required to implement finite
22 difference scheme. For example, 2 means that we need to evaluate
23 f(x0 + 2 * h) or f(x0 - 2 * h)
24 scheme : {'1-sided', '2-sided'}
25 Whether steps in one or both directions are required. In other
26 words '1-sided' applies to forward and backward schemes, '2-sided'
27 applies to center schemes.
28 lb : ndarray, shape (n,)
29 Lower bounds on independent variables.
30 ub : ndarray, shape (n,)
31 Upper bounds on independent variables.
33 Returns
34 -------
35 h_adjusted : ndarray, shape (n,)
36 Adjusted absolute step sizes. Step size decreases only if a sign flip
37 or switching to one-sided scheme doesn't allow to take a full step.
38 use_one_sided : ndarray of bool, shape (n,)
39 Whether to switch to one-sided scheme. Informative only for
40 ``scheme='2-sided'``.
41 """
42 if scheme == '1-sided':
43 use_one_sided = np.ones_like(h, dtype=bool)
44 elif scheme == '2-sided':
45 h = np.abs(h)
46 use_one_sided = np.zeros_like(h, dtype=bool)
47 else:
48 raise ValueError("`scheme` must be '1-sided' or '2-sided'.")
50 if np.all((lb == -np.inf) & (ub == np.inf)):
51 return h, use_one_sided
53 h_total = h * num_steps
54 h_adjusted = h.copy()
56 lower_dist = x0 - lb
57 upper_dist = ub - x0
59 if scheme == '1-sided':
60 x = x0 + h_total
61 violated = (x < lb) | (x > ub)
62 fitting = np.abs(h_total) <= np.maximum(lower_dist, upper_dist)
63 h_adjusted[violated & fitting] *= -1
65 forward = (upper_dist >= lower_dist) & ~fitting
66 h_adjusted[forward] = upper_dist[forward] / num_steps
67 backward = (upper_dist < lower_dist) & ~fitting
68 h_adjusted[backward] = -lower_dist[backward] / num_steps
69 elif scheme == '2-sided':
70 central = (lower_dist >= h_total) & (upper_dist >= h_total)
72 forward = (upper_dist >= lower_dist) & ~central
73 h_adjusted[forward] = np.minimum(
74 h[forward], 0.5 * upper_dist[forward] / num_steps)
75 use_one_sided[forward] = True
77 backward = (upper_dist < lower_dist) & ~central
78 h_adjusted[backward] = -np.minimum(
79 h[backward], 0.5 * lower_dist[backward] / num_steps)
80 use_one_sided[backward] = True
82 min_dist = np.minimum(upper_dist, lower_dist) / num_steps
83 adjusted_central = (~central & (np.abs(h_adjusted) <= min_dist))
84 h_adjusted[adjusted_central] = min_dist[adjusted_central]
85 use_one_sided[adjusted_central] = False
87 return h_adjusted, use_one_sided
90@functools.lru_cache()
91def _eps_for_method(x0_dtype, f0_dtype, method):
92 """
93 Calculates relative EPS step to use for a given data type
94 and numdiff step method.
96 Progressively smaller steps are used for larger floating point types.
98 Parameters
99 ----------
100 f0_dtype: np.dtype
101 dtype of function evaluation
103 x0_dtype: np.dtype
104 dtype of parameter vector
106 method: {'2-point', '3-point', 'cs'}
108 Returns
109 -------
110 EPS: float
111 relative step size. May be np.float16, np.float32, np.float64
113 Notes
114 -----
115 The default relative step will be np.float64. However, if x0 or f0 are
116 smaller floating point types (np.float16, np.float32), then the smallest
117 floating point type is chosen.
118 """
119 # the default EPS value
120 EPS = np.finfo(np.float64).eps
122 x0_is_fp = False
123 if np.issubdtype(x0_dtype, np.inexact):
124 # if you're a floating point type then over-ride the default EPS
125 EPS = np.finfo(x0_dtype).eps
126 x0_itemsize = np.dtype(x0_dtype).itemsize
127 x0_is_fp = True
129 if np.issubdtype(f0_dtype, np.inexact):
130 f0_itemsize = np.dtype(f0_dtype).itemsize
131 # choose the smallest itemsize between x0 and f0
132 if x0_is_fp and f0_itemsize < x0_itemsize:
133 EPS = np.finfo(f0_dtype).eps
135 if method in ["2-point", "cs"]:
136 return EPS**0.5
137 elif method in ["3-point"]:
138 return EPS**(1/3)
139 else:
140 raise RuntimeError("Unknown step method, should be one of "
141 "{'2-point', '3-point', 'cs'}")
144def _compute_absolute_step(rel_step, x0, f0, method):
145 """
146 Computes an absolute step from a relative step for finite difference
147 calculation.
149 Parameters
150 ----------
151 rel_step: None or array-like
152 Relative step for the finite difference calculation
153 x0 : np.ndarray
154 Parameter vector
155 f0 : np.ndarray or scalar
156 method : {'2-point', '3-point', 'cs'}
158 Returns
159 -------
160 h : float
161 The absolute step size
163 Notes
164 -----
165 `h` will always be np.float64. However, if `x0` or `f0` are
166 smaller floating point dtypes (e.g. np.float32), then the absolute
167 step size will be calculated from the smallest floating point size.
168 """
169 # this is used instead of np.sign(x0) because we need
170 # sign_x0 to be 1 when x0 == 0.
171 sign_x0 = (x0 >= 0).astype(float) * 2 - 1
173 rstep = _eps_for_method(x0.dtype, f0.dtype, method)
175 if rel_step is None:
176 abs_step = rstep * sign_x0 * np.maximum(1.0, np.abs(x0))
177 else:
178 # User has requested specific relative steps.
179 # Don't multiply by max(1, abs(x0) because if x0 < 1 then their
180 # requested step is not used.
181 abs_step = rel_step * sign_x0 * np.abs(x0)
183 # however we don't want an abs_step of 0, which can happen if
184 # rel_step is 0, or x0 is 0. Instead, substitute a realistic step
185 dx = ((x0 + abs_step) - x0)
186 abs_step = np.where(dx == 0,
187 rstep * sign_x0 * np.maximum(1.0, np.abs(x0)),
188 abs_step)
190 return abs_step
193def _prepare_bounds(bounds, x0):
194 """
195 Prepares new-style bounds from a two-tuple specifying the lower and upper
196 limits for values in x0. If a value is not bound then the lower/upper bound
197 will be expected to be -np.inf/np.inf.
199 Examples
200 --------
201 >>> _prepare_bounds([(0, 1, 2), (1, 2, np.inf)], [0.5, 1.5, 2.5])
202 (array([0., 1., 2.]), array([ 1., 2., inf]))
203 """
204 lb, ub = [np.asarray(b, dtype=float) for b in bounds]
205 if lb.ndim == 0:
206 lb = np.resize(lb, x0.shape)
208 if ub.ndim == 0:
209 ub = np.resize(ub, x0.shape)
211 return lb, ub
214def group_columns(A, order=0):
215 """Group columns of a 2-D matrix for sparse finite differencing [1]_.
217 Two columns are in the same group if in each row at least one of them
218 has zero. A greedy sequential algorithm is used to construct groups.
220 Parameters
221 ----------
222 A : array_like or sparse matrix, shape (m, n)
223 Matrix of which to group columns.
224 order : int, iterable of int with shape (n,) or None
225 Permutation array which defines the order of columns enumeration.
226 If int or None, a random permutation is used with `order` used as
227 a random seed. Default is 0, that is use a random permutation but
228 guarantee repeatability.
230 Returns
231 -------
232 groups : ndarray of int, shape (n,)
233 Contains values from 0 to n_groups-1, where n_groups is the number
234 of found groups. Each value ``groups[i]`` is an index of a group to
235 which ith column assigned. The procedure was helpful only if
236 n_groups is significantly less than n.
238 References
239 ----------
240 .. [1] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
241 sparse Jacobian matrices", Journal of the Institute of Mathematics
242 and its Applications, 13 (1974), pp. 117-120.
243 """
244 if issparse(A):
245 A = csc_matrix(A)
246 else:
247 A = np.atleast_2d(A)
248 A = (A != 0).astype(np.int32)
250 if A.ndim != 2:
251 raise ValueError("`A` must be 2-dimensional.")
253 m, n = A.shape
255 if order is None or np.isscalar(order):
256 rng = np.random.RandomState(order)
257 order = rng.permutation(n)
258 else:
259 order = np.asarray(order)
260 if order.shape != (n,):
261 raise ValueError("`order` has incorrect shape.")
263 A = A[:, order]
265 if issparse(A):
266 groups = group_sparse(m, n, A.indices, A.indptr)
267 else:
268 groups = group_dense(m, n, A)
270 groups[order] = groups.copy()
272 return groups
275def approx_derivative(fun, x0, method='3-point', rel_step=None, abs_step=None,
276 f0=None, bounds=(-np.inf, np.inf), sparsity=None,
277 as_linear_operator=False, args=(), kwargs={}):
278 """Compute finite difference approximation of the derivatives of a
279 vector-valued function.
281 If a function maps from R^n to R^m, its derivatives form m-by-n matrix
282 called the Jacobian, where an element (i, j) is a partial derivative of
283 f[i] with respect to x[j].
285 Parameters
286 ----------
287 fun : callable
288 Function of which to estimate the derivatives. The argument x
289 passed to this function is ndarray of shape (n,) (never a scalar
290 even if n=1). It must return 1-D array_like of shape (m,) or a scalar.
291 x0 : array_like of shape (n,) or float
292 Point at which to estimate the derivatives. Float will be converted
293 to a 1-D array.
294 method : {'3-point', '2-point', 'cs'}, optional
295 Finite difference method to use:
296 - '2-point' - use the first order accuracy forward or backward
297 difference.
298 - '3-point' - use central difference in interior points and the
299 second order accuracy forward or backward difference
300 near the boundary.
301 - 'cs' - use a complex-step finite difference scheme. This assumes
302 that the user function is real-valued and can be
303 analytically continued to the complex plane. Otherwise,
304 produces bogus results.
305 rel_step : None or array_like, optional
306 Relative step size to use. If None (default) the absolute step size is
307 computed as ``h = rel_step * sign(x0) * max(1, abs(x0))``, with
308 `rel_step` being selected automatically, see Notes. Otherwise
309 ``h = rel_step * sign(x0) * abs(x0)``. For ``method='3-point'`` the
310 sign of `h` is ignored. The calculated step size is possibly adjusted
311 to fit into the bounds.
312 abs_step : array_like, optional
313 Absolute step size to use, possibly adjusted to fit into the bounds.
314 For ``method='3-point'`` the sign of `abs_step` is ignored. By default
315 relative steps are used, only if ``abs_step is not None`` are absolute
316 steps used.
317 f0 : None or array_like, optional
318 If not None it is assumed to be equal to ``fun(x0)``, in this case
319 the ``fun(x0)`` is not called. Default is None.
320 bounds : tuple of array_like, optional
321 Lower and upper bounds on independent variables. Defaults to no bounds.
322 Each bound must match the size of `x0` or be a scalar, in the latter
323 case the bound will be the same for all variables. Use it to limit the
324 range of function evaluation. Bounds checking is not implemented
325 when `as_linear_operator` is True.
326 sparsity : {None, array_like, sparse matrix, 2-tuple}, optional
327 Defines a sparsity structure of the Jacobian matrix. If the Jacobian
328 matrix is known to have only few non-zero elements in each row, then
329 it's possible to estimate its several columns by a single function
330 evaluation [3]_. To perform such economic computations two ingredients
331 are required:
333 * structure : array_like or sparse matrix of shape (m, n). A zero
334 element means that a corresponding element of the Jacobian
335 identically equals to zero.
336 * groups : array_like of shape (n,). A column grouping for a given
337 sparsity structure, use `group_columns` to obtain it.
339 A single array or a sparse matrix is interpreted as a sparsity
340 structure, and groups are computed inside the function. A tuple is
341 interpreted as (structure, groups). If None (default), a standard
342 dense differencing will be used.
344 Note, that sparse differencing makes sense only for large Jacobian
345 matrices where each row contains few non-zero elements.
346 as_linear_operator : bool, optional
347 When True the function returns an `scipy.sparse.linalg.LinearOperator`.
348 Otherwise it returns a dense array or a sparse matrix depending on
349 `sparsity`. The linear operator provides an efficient way of computing
350 ``J.dot(p)`` for any vector ``p`` of shape (n,), but does not allow
351 direct access to individual elements of the matrix. By default
352 `as_linear_operator` is False.
353 args, kwargs : tuple and dict, optional
354 Additional arguments passed to `fun`. Both empty by default.
355 The calling signature is ``fun(x, *args, **kwargs)``.
357 Returns
358 -------
359 J : {ndarray, sparse matrix, LinearOperator}
360 Finite difference approximation of the Jacobian matrix.
361 If `as_linear_operator` is True returns a LinearOperator
362 with shape (m, n). Otherwise it returns a dense array or sparse
363 matrix depending on how `sparsity` is defined. If `sparsity`
364 is None then a ndarray with shape (m, n) is returned. If
365 `sparsity` is not None returns a csr_matrix with shape (m, n).
366 For sparse matrices and linear operators it is always returned as
367 a 2-D structure, for ndarrays, if m=1 it is returned
368 as a 1-D gradient array with shape (n,).
370 See Also
371 --------
372 check_derivative : Check correctness of a function computing derivatives.
374 Notes
375 -----
376 If `rel_step` is not provided, it assigned as ``EPS**(1/s)``, where EPS is
377 determined from the smallest floating point dtype of `x0` or `fun(x0)`,
378 ``np.finfo(x0.dtype).eps``, s=2 for '2-point' method and
379 s=3 for '3-point' method. Such relative step approximately minimizes a sum
380 of truncation and round-off errors, see [1]_. Relative steps are used by
381 default. However, absolute steps are used when ``abs_step is not None``.
382 If any of the absolute or relative steps produces an indistinguishable
383 difference from the original `x0`, ``(x0 + dx) - x0 == 0``, then a
384 automatic step size is substituted for that particular entry.
386 A finite difference scheme for '3-point' method is selected automatically.
387 The well-known central difference scheme is used for points sufficiently
388 far from the boundary, and 3-point forward or backward scheme is used for
389 points near the boundary. Both schemes have the second-order accuracy in
390 terms of Taylor expansion. Refer to [2]_ for the formulas of 3-point
391 forward and backward difference schemes.
393 For dense differencing when m=1 Jacobian is returned with a shape (n,),
394 on the other hand when n=1 Jacobian is returned with a shape (m, 1).
395 Our motivation is the following: a) It handles a case of gradient
396 computation (m=1) in a conventional way. b) It clearly separates these two
397 different cases. b) In all cases np.atleast_2d can be called to get 2-D
398 Jacobian with correct dimensions.
400 References
401 ----------
402 .. [1] W. H. Press et. al. "Numerical Recipes. The Art of Scientific
403 Computing. 3rd edition", sec. 5.7.
405 .. [2] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
406 sparse Jacobian matrices", Journal of the Institute of Mathematics
407 and its Applications, 13 (1974), pp. 117-120.
409 .. [3] B. Fornberg, "Generation of Finite Difference Formulas on
410 Arbitrarily Spaced Grids", Mathematics of Computation 51, 1988.
412 Examples
413 --------
414 >>> import numpy as np
415 >>> from scipy.optimize._numdiff import approx_derivative
416 >>>
417 >>> def f(x, c1, c2):
418 ... return np.array([x[0] * np.sin(c1 * x[1]),
419 ... x[0] * np.cos(c2 * x[1])])
420 ...
421 >>> x0 = np.array([1.0, 0.5 * np.pi])
422 >>> approx_derivative(f, x0, args=(1, 2))
423 array([[ 1., 0.],
424 [-1., 0.]])
426 Bounds can be used to limit the region of function evaluation.
427 In the example below we compute left and right derivative at point 1.0.
429 >>> def g(x):
430 ... return x**2 if x >= 1 else x
431 ...
432 >>> x0 = 1.0
433 >>> approx_derivative(g, x0, bounds=(-np.inf, 1.0))
434 array([ 1.])
435 >>> approx_derivative(g, x0, bounds=(1.0, np.inf))
436 array([ 2.])
437 """
438 if method not in ['2-point', '3-point', 'cs']:
439 raise ValueError("Unknown method '%s'. " % method)
441 x0 = np.atleast_1d(x0)
442 if x0.ndim > 1:
443 raise ValueError("`x0` must have at most 1 dimension.")
445 lb, ub = _prepare_bounds(bounds, x0)
447 if lb.shape != x0.shape or ub.shape != x0.shape:
448 raise ValueError("Inconsistent shapes between bounds and `x0`.")
450 if as_linear_operator and not (np.all(np.isinf(lb))
451 and np.all(np.isinf(ub))):
452 raise ValueError("Bounds not supported when "
453 "`as_linear_operator` is True.")
455 def fun_wrapped(x):
456 f = np.atleast_1d(fun(x, *args, **kwargs))
457 if f.ndim > 1:
458 raise RuntimeError("`fun` return value has "
459 "more than 1 dimension.")
460 return f
462 if f0 is None:
463 f0 = fun_wrapped(x0)
464 else:
465 f0 = np.atleast_1d(f0)
466 if f0.ndim > 1:
467 raise ValueError("`f0` passed has more than 1 dimension.")
469 if np.any((x0 < lb) | (x0 > ub)):
470 raise ValueError("`x0` violates bound constraints.")
472 if as_linear_operator:
473 if rel_step is None:
474 rel_step = _eps_for_method(x0.dtype, f0.dtype, method)
476 return _linear_operator_difference(fun_wrapped, x0,
477 f0, rel_step, method)
478 else:
479 # by default we use rel_step
480 if abs_step is None:
481 h = _compute_absolute_step(rel_step, x0, f0, method)
482 else:
483 # user specifies an absolute step
484 sign_x0 = (x0 >= 0).astype(float) * 2 - 1
485 h = abs_step
487 # cannot have a zero step. This might happen if x0 is very large
488 # or small. In which case fall back to relative step.
489 dx = ((x0 + h) - x0)
490 h = np.where(dx == 0,
491 _eps_for_method(x0.dtype, f0.dtype, method) *
492 sign_x0 * np.maximum(1.0, np.abs(x0)),
493 h)
495 if method == '2-point':
496 h, use_one_sided = _adjust_scheme_to_bounds(
497 x0, h, 1, '1-sided', lb, ub)
498 elif method == '3-point':
499 h, use_one_sided = _adjust_scheme_to_bounds(
500 x0, h, 1, '2-sided', lb, ub)
501 elif method == 'cs':
502 use_one_sided = False
504 if sparsity is None:
505 return _dense_difference(fun_wrapped, x0, f0, h,
506 use_one_sided, method)
507 else:
508 if not issparse(sparsity) and len(sparsity) == 2:
509 structure, groups = sparsity
510 else:
511 structure = sparsity
512 groups = group_columns(sparsity)
514 if issparse(structure):
515 structure = csc_matrix(structure)
516 else:
517 structure = np.atleast_2d(structure)
519 groups = np.atleast_1d(groups)
520 return _sparse_difference(fun_wrapped, x0, f0, h,
521 use_one_sided, structure,
522 groups, method)
525def _linear_operator_difference(fun, x0, f0, h, method):
526 m = f0.size
527 n = x0.size
529 if method == '2-point':
530 def matvec(p):
531 if np.array_equal(p, np.zeros_like(p)):
532 return np.zeros(m)
533 dx = h / norm(p)
534 x = x0 + dx*p
535 df = fun(x) - f0
536 return df / dx
538 elif method == '3-point':
539 def matvec(p):
540 if np.array_equal(p, np.zeros_like(p)):
541 return np.zeros(m)
542 dx = 2*h / norm(p)
543 x1 = x0 - (dx/2)*p
544 x2 = x0 + (dx/2)*p
545 f1 = fun(x1)
546 f2 = fun(x2)
547 df = f2 - f1
548 return df / dx
550 elif method == 'cs':
551 def matvec(p):
552 if np.array_equal(p, np.zeros_like(p)):
553 return np.zeros(m)
554 dx = h / norm(p)
555 x = x0 + dx*p*1.j
556 f1 = fun(x)
557 df = f1.imag
558 return df / dx
560 else:
561 raise RuntimeError("Never be here.")
563 return LinearOperator((m, n), matvec)
566def _dense_difference(fun, x0, f0, h, use_one_sided, method):
567 m = f0.size
568 n = x0.size
569 J_transposed = np.empty((n, m))
570 h_vecs = np.diag(h)
572 for i in range(h.size):
573 if method == '2-point':
574 x = x0 + h_vecs[i]
575 dx = x[i] - x0[i] # Recompute dx as exactly representable number.
576 df = fun(x) - f0
577 elif method == '3-point' and use_one_sided[i]:
578 x1 = x0 + h_vecs[i]
579 x2 = x0 + 2 * h_vecs[i]
580 dx = x2[i] - x0[i]
581 f1 = fun(x1)
582 f2 = fun(x2)
583 df = -3.0 * f0 + 4 * f1 - f2
584 elif method == '3-point' and not use_one_sided[i]:
585 x1 = x0 - h_vecs[i]
586 x2 = x0 + h_vecs[i]
587 dx = x2[i] - x1[i]
588 f1 = fun(x1)
589 f2 = fun(x2)
590 df = f2 - f1
591 elif method == 'cs':
592 f1 = fun(x0 + h_vecs[i]*1.j)
593 df = f1.imag
594 dx = h_vecs[i, i]
595 else:
596 raise RuntimeError("Never be here.")
598 J_transposed[i] = df / dx
600 if m == 1:
601 J_transposed = np.ravel(J_transposed)
603 return J_transposed.T
606def _sparse_difference(fun, x0, f0, h, use_one_sided,
607 structure, groups, method):
608 m = f0.size
609 n = x0.size
610 row_indices = []
611 col_indices = []
612 fractions = []
614 n_groups = np.max(groups) + 1
615 for group in range(n_groups):
616 # Perturb variables which are in the same group simultaneously.
617 e = np.equal(group, groups)
618 h_vec = h * e
619 if method == '2-point':
620 x = x0 + h_vec
621 dx = x - x0
622 df = fun(x) - f0
623 # The result is written to columns which correspond to perturbed
624 # variables.
625 cols, = np.nonzero(e)
626 # Find all non-zero elements in selected columns of Jacobian.
627 i, j, _ = find(structure[:, cols])
628 # Restore column indices in the full array.
629 j = cols[j]
630 elif method == '3-point':
631 # Here we do conceptually the same but separate one-sided
632 # and two-sided schemes.
633 x1 = x0.copy()
634 x2 = x0.copy()
636 mask_1 = use_one_sided & e
637 x1[mask_1] += h_vec[mask_1]
638 x2[mask_1] += 2 * h_vec[mask_1]
640 mask_2 = ~use_one_sided & e
641 x1[mask_2] -= h_vec[mask_2]
642 x2[mask_2] += h_vec[mask_2]
644 dx = np.zeros(n)
645 dx[mask_1] = x2[mask_1] - x0[mask_1]
646 dx[mask_2] = x2[mask_2] - x1[mask_2]
648 f1 = fun(x1)
649 f2 = fun(x2)
651 cols, = np.nonzero(e)
652 i, j, _ = find(structure[:, cols])
653 j = cols[j]
655 mask = use_one_sided[j]
656 df = np.empty(m)
658 rows = i[mask]
659 df[rows] = -3 * f0[rows] + 4 * f1[rows] - f2[rows]
661 rows = i[~mask]
662 df[rows] = f2[rows] - f1[rows]
663 elif method == 'cs':
664 f1 = fun(x0 + h_vec*1.j)
665 df = f1.imag
666 dx = h_vec
667 cols, = np.nonzero(e)
668 i, j, _ = find(structure[:, cols])
669 j = cols[j]
670 else:
671 raise ValueError("Never be here.")
673 # All that's left is to compute the fraction. We store i, j and
674 # fractions as separate arrays and later construct coo_matrix.
675 row_indices.append(i)
676 col_indices.append(j)
677 fractions.append(df[i] / dx[j])
679 row_indices = np.hstack(row_indices)
680 col_indices = np.hstack(col_indices)
681 fractions = np.hstack(fractions)
682 J = coo_matrix((fractions, (row_indices, col_indices)), shape=(m, n))
683 return csr_matrix(J)
686def check_derivative(fun, jac, x0, bounds=(-np.inf, np.inf), args=(),
687 kwargs={}):
688 """Check correctness of a function computing derivatives (Jacobian or
689 gradient) by comparison with a finite difference approximation.
691 Parameters
692 ----------
693 fun : callable
694 Function of which to estimate the derivatives. The argument x
695 passed to this function is ndarray of shape (n,) (never a scalar
696 even if n=1). It must return 1-D array_like of shape (m,) or a scalar.
697 jac : callable
698 Function which computes Jacobian matrix of `fun`. It must work with
699 argument x the same way as `fun`. The return value must be array_like
700 or sparse matrix with an appropriate shape.
701 x0 : array_like of shape (n,) or float
702 Point at which to estimate the derivatives. Float will be converted
703 to 1-D array.
704 bounds : 2-tuple of array_like, optional
705 Lower and upper bounds on independent variables. Defaults to no bounds.
706 Each bound must match the size of `x0` or be a scalar, in the latter
707 case the bound will be the same for all variables. Use it to limit the
708 range of function evaluation.
709 args, kwargs : tuple and dict, optional
710 Additional arguments passed to `fun` and `jac`. Both empty by default.
711 The calling signature is ``fun(x, *args, **kwargs)`` and the same
712 for `jac`.
714 Returns
715 -------
716 accuracy : float
717 The maximum among all relative errors for elements with absolute values
718 higher than 1 and absolute errors for elements with absolute values
719 less or equal than 1. If `accuracy` is on the order of 1e-6 or lower,
720 then it is likely that your `jac` implementation is correct.
722 See Also
723 --------
724 approx_derivative : Compute finite difference approximation of derivative.
726 Examples
727 --------
728 >>> import numpy as np
729 >>> from scipy.optimize._numdiff import check_derivative
730 >>>
731 >>>
732 >>> def f(x, c1, c2):
733 ... return np.array([x[0] * np.sin(c1 * x[1]),
734 ... x[0] * np.cos(c2 * x[1])])
735 ...
736 >>> def jac(x, c1, c2):
737 ... return np.array([
738 ... [np.sin(c1 * x[1]), c1 * x[0] * np.cos(c1 * x[1])],
739 ... [np.cos(c2 * x[1]), -c2 * x[0] * np.sin(c2 * x[1])]
740 ... ])
741 ...
742 >>>
743 >>> x0 = np.array([1.0, 0.5 * np.pi])
744 >>> check_derivative(f, jac, x0, args=(1, 2))
745 2.4492935982947064e-16
746 """
747 J_to_test = jac(x0, *args, **kwargs)
748 if issparse(J_to_test):
749 J_diff = approx_derivative(fun, x0, bounds=bounds, sparsity=J_to_test,
750 args=args, kwargs=kwargs)
751 J_to_test = csr_matrix(J_to_test)
752 abs_err = J_to_test - J_diff
753 i, j, abs_err_data = find(abs_err)
754 J_diff_data = np.asarray(J_diff[i, j]).ravel()
755 return np.max(np.abs(abs_err_data) /
756 np.maximum(1, np.abs(J_diff_data)))
757 else:
758 J_diff = approx_derivative(fun, x0, bounds=bounds,
759 args=args, kwargs=kwargs)
760 abs_err = np.abs(J_to_test - J_diff)
761 return np.max(abs_err / np.maximum(1, np.abs(J_diff)))