1"""Routines for numerical differentiation.""" 2import functools 3import numpy as np 4from numpy.linalg import norm 5 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 9 10 11def _adjust_scheme_to_bounds(x0, h, num_steps, scheme, lb, ub): 12 """Adjust final difference scheme to the presence of bounds. 13 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. 32 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'.") 49 50 if np.all((lb == -np.inf) & (ub == np.inf)): 51 return h, use_one_sided 52 53 h_total = h * num_steps 54 h_adjusted = h.copy() 55 56 lower_dist = x0 - lb 57 upper_dist = ub - x0 58 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 64 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) 71 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 76 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 81 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 86 87 return h_adjusted, use_one_sided 88 89 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. 95 96 Progressively smaller steps are used for larger floating point types. 97 98 Parameters 99 ---------- 100 f0_dtype: np.dtype 101 dtype of function evaluation 102 103 x0_dtype: np.dtype 104 dtype of parameter vector 105 106 method: {'2-point', '3-point', 'cs'} 107 108 Returns 109 ------- 110 EPS: float 111 relative step size. May be np.float16, np.float32, np.float64 112 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 121 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 128 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 134 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'}") 142 143 144def _compute_absolute_step(rel_step, x0, f0, method): 145 """ 146 Computes an absolute step from a relative step for finite difference 147 calculation. 148 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'} 157 158 Returns 159 ------- 160 h : float 161 The absolute step size 162 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 if rel_step is None: 170 rel_step = _eps_for_method(x0.dtype, f0.dtype, method) 171 sign_x0 = (x0 >= 0).astype(float) * 2 - 1 172 return rel_step * sign_x0 * np.maximum(1.0, np.abs(x0)) 173 174 175def _prepare_bounds(bounds, x0): 176 """ 177 Prepares new-style bounds from a two-tuple specifying the lower and upper 178 limits for values in x0. If a value is not bound then the lower/upper bound 179 will be expected to be -np.inf/np.inf. 180 181 Examples 182 -------- 183 >>> _prepare_bounds([(0, 1, 2), (1, 2, np.inf)], [0.5, 1.5, 2.5]) 184 (array([0., 1., 2.]), array([ 1., 2., inf])) 185 """ 186 lb, ub = [np.asarray(b, dtype=float) for b in bounds] 187 if lb.ndim == 0: 188 lb = np.resize(lb, x0.shape) 189 190 if ub.ndim == 0: 191 ub = np.resize(ub, x0.shape) 192 193 return lb, ub 194 195 196def group_columns(A, order=0): 197 """Group columns of a 2-D matrix for sparse finite differencing [1]_. 198 199 Two columns are in the same group if in each row at least one of them 200 has zero. A greedy sequential algorithm is used to construct groups. 201 202 Parameters 203 ---------- 204 A : array_like or sparse matrix, shape (m, n) 205 Matrix of which to group columns. 206 order : int, iterable of int with shape (n,) or None 207 Permutation array which defines the order of columns enumeration. 208 If int or None, a random permutation is used with `order` used as 209 a random seed. Default is 0, that is use a random permutation but 210 guarantee repeatability. 211 212 Returns 213 ------- 214 groups : ndarray of int, shape (n,) 215 Contains values from 0 to n_groups-1, where n_groups is the number 216 of found groups. Each value ``groups[i]`` is an index of a group to 217 which ith column assigned. The procedure was helpful only if 218 n_groups is significantly less than n. 219 220 References 221 ---------- 222 .. [1] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of 223 sparse Jacobian matrices", Journal of the Institute of Mathematics 224 and its Applications, 13 (1974), pp. 117-120. 225 """ 226 if issparse(A): 227 A = csc_matrix(A) 228 else: 229 A = np.atleast_2d(A) 230 A = (A != 0).astype(np.int32) 231 232 if A.ndim != 2: 233 raise ValueError("`A` must be 2-dimensional.") 234 235 m, n = A.shape 236 237 if order is None or np.isscalar(order): 238 rng = np.random.RandomState(order) 239 order = rng.permutation(n) 240 else: 241 order = np.asarray(order) 242 if order.shape != (n,): 243 raise ValueError("`order` has incorrect shape.") 244 245 A = A[:, order] 246 247 if issparse(A): 248 groups = group_sparse(m, n, A.indices, A.indptr) 249 else: 250 groups = group_dense(m, n, A) 251 252 groups[order] = groups.copy() 253 254 return groups 255 256 257def approx_derivative(fun, x0, method='3-point', rel_step=None, abs_step=None, 258 f0=None, bounds=(-np.inf, np.inf), sparsity=None, 259 as_linear_operator=False, args=(), kwargs={}): 260 """Compute finite difference approximation of the derivatives of a 261 vector-valued function. 262 263 If a function maps from R^n to R^m, its derivatives form m-by-n matrix 264 called the Jacobian, where an element (i, j) is a partial derivative of 265 f[i] with respect to x[j]. 266 267 Parameters 268 ---------- 269 fun : callable 270 Function of which to estimate the derivatives. The argument x 271 passed to this function is ndarray of shape (n,) (never a scalar 272 even if n=1). It must return 1-D array_like of shape (m,) or a scalar. 273 x0 : array_like of shape (n,) or float 274 Point at which to estimate the derivatives. Float will be converted 275 to a 1-D array. 276 method : {'3-point', '2-point', 'cs'}, optional 277 Finite difference method to use: 278 - '2-point' - use the first order accuracy forward or backward 279 difference. 280 - '3-point' - use central difference in interior points and the 281 second order accuracy forward or backward difference 282 near the boundary. 283 - 'cs' - use a complex-step finite difference scheme. This assumes 284 that the user function is real-valued and can be 285 analytically continued to the complex plane. Otherwise, 286 produces bogus results. 287 rel_step : None or array_like, optional 288 Relative step size to use. The absolute step size is computed as 289 ``h = rel_step * sign(x0) * max(1, abs(x0))``, possibly adjusted to 290 fit into the bounds. For ``method='3-point'`` the sign of `h` is 291 ignored. If None (default) then step is selected automatically, 292 see Notes. 293 abs_step : array_like, optional 294 Absolute step size to use, possibly adjusted to fit into the bounds. 295 For ``method='3-point'`` the sign of `abs_step` is ignored. By default 296 relative steps are used, only if ``abs_step is not None`` are absolute 297 steps used. 298 f0 : None or array_like, optional 299 If not None it is assumed to be equal to ``fun(x0)``, in this case 300 the ``fun(x0)`` is not called. Default is None. 301 bounds : tuple of array_like, optional 302 Lower and upper bounds on independent variables. Defaults to no bounds. 303 Each bound must match the size of `x0` or be a scalar, in the latter 304 case the bound will be the same for all variables. Use it to limit the 305 range of function evaluation. Bounds checking is not implemented 306 when `as_linear_operator` is True. 307 sparsity : {None, array_like, sparse matrix, 2-tuple}, optional 308 Defines a sparsity structure of the Jacobian matrix. If the Jacobian 309 matrix is known to have only few non-zero elements in each row, then 310 it's possible to estimate its several columns by a single function 311 evaluation [3]_. To perform such economic computations two ingredients 312 are required: 313 314 * structure : array_like or sparse matrix of shape (m, n). A zero 315 element means that a corresponding element of the Jacobian 316 identically equals to zero. 317 * groups : array_like of shape (n,). A column grouping for a given 318 sparsity structure, use `group_columns` to obtain it. 319 320 A single array or a sparse matrix is interpreted as a sparsity 321 structure, and groups are computed inside the function. A tuple is 322 interpreted as (structure, groups). If None (default), a standard 323 dense differencing will be used. 324 325 Note, that sparse differencing makes sense only for large Jacobian 326 matrices where each row contains few non-zero elements. 327 as_linear_operator : bool, optional 328 When True the function returns an `scipy.sparse.linalg.LinearOperator`. 329 Otherwise it returns a dense array or a sparse matrix depending on 330 `sparsity`. The linear operator provides an efficient way of computing 331 ``J.dot(p)`` for any vector ``p`` of shape (n,), but does not allow 332 direct access to individual elements of the matrix. By default 333 `as_linear_operator` is False. 334 args, kwargs : tuple and dict, optional 335 Additional arguments passed to `fun`. Both empty by default. 336 The calling signature is ``fun(x, *args, **kwargs)``. 337 338 Returns 339 ------- 340 J : {ndarray, sparse matrix, LinearOperator} 341 Finite difference approximation of the Jacobian matrix. 342 If `as_linear_operator` is True returns a LinearOperator 343 with shape (m, n). Otherwise it returns a dense array or sparse 344 matrix depending on how `sparsity` is defined. If `sparsity` 345 is None then a ndarray with shape (m, n) is returned. If 346 `sparsity` is not None returns a csr_matrix with shape (m, n). 347 For sparse matrices and linear operators it is always returned as 348 a 2-D structure, for ndarrays, if m=1 it is returned 349 as a 1-D gradient array with shape (n,). 350 351 See Also 352 -------- 353 check_derivative : Check correctness of a function computing derivatives. 354 355 Notes 356 ----- 357 If `rel_step` is not provided, it assigned as ``EPS**(1/s)``, where EPS is 358 determined from the smallest floating point dtype of `x0` or `fun(x0)`, 359 ``np.finfo(x0.dtype).eps``, s=2 for '2-point' method and 360 s=3 for '3-point' method. Such relative step approximately minimizes a sum 361 of truncation and round-off errors, see [1]_. Relative steps are used by 362 default. However, absolute steps are used when ``abs_step is not None``. 363 If any of the absolute steps produces an indistinguishable difference from 364 the original `x0`, ``(x0 + abs_step) - x0 == 0``, then a relative step is 365 substituted for that particular entry. 366 367 A finite difference scheme for '3-point' method is selected automatically. 368 The well-known central difference scheme is used for points sufficiently 369 far from the boundary, and 3-point forward or backward scheme is used for 370 points near the boundary. Both schemes have the second-order accuracy in 371 terms of Taylor expansion. Refer to [2]_ for the formulas of 3-point 372 forward and backward difference schemes. 373 374 For dense differencing when m=1 Jacobian is returned with a shape (n,), 375 on the other hand when n=1 Jacobian is returned with a shape (m, 1). 376 Our motivation is the following: a) It handles a case of gradient 377 computation (m=1) in a conventional way. b) It clearly separates these two 378 different cases. b) In all cases np.atleast_2d can be called to get 2-D 379 Jacobian with correct dimensions. 380 381 References 382 ---------- 383 .. [1] W. H. Press et. al. "Numerical Recipes. The Art of Scientific 384 Computing. 3rd edition", sec. 5.7. 385 386 .. [2] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of 387 sparse Jacobian matrices", Journal of the Institute of Mathematics 388 and its Applications, 13 (1974), pp. 117-120. 389 390 .. [3] B. Fornberg, "Generation of Finite Difference Formulas on 391 Arbitrarily Spaced Grids", Mathematics of Computation 51, 1988. 392 393 Examples 394 -------- 395 >>> import numpy as np 396 >>> from scipy.optimize import approx_derivative 397 >>> 398 >>> def f(x, c1, c2): 399 ... return np.array([x[0] * np.sin(c1 * x[1]), 400 ... x[0] * np.cos(c2 * x[1])]) 401 ... 402 >>> x0 = np.array([1.0, 0.5 * np.pi]) 403 >>> approx_derivative(f, x0, args=(1, 2)) 404 array([[ 1., 0.], 405 [-1., 0.]]) 406 407 Bounds can be used to limit the region of function evaluation. 408 In the example below we compute left and right derivative at point 1.0. 409 410 >>> def g(x): 411 ... return x**2 if x >= 1 else x 412 ... 413 >>> x0 = 1.0 414 >>> approx_derivative(g, x0, bounds=(-np.inf, 1.0)) 415 array([ 1.]) 416 >>> approx_derivative(g, x0, bounds=(1.0, np.inf)) 417 array([ 2.]) 418 """ 419 if method not in ['2-point', '3-point', 'cs']: 420 raise ValueError("Unknown method '%s'. " % method) 421 422 x0 = np.atleast_1d(x0) 423 if x0.ndim > 1: 424 raise ValueError("`x0` must have at most 1 dimension.") 425 426 lb, ub = _prepare_bounds(bounds, x0) 427 428 if lb.shape != x0.shape or ub.shape != x0.shape: 429 raise ValueError("Inconsistent shapes between bounds and `x0`.") 430 431 if as_linear_operator and not (np.all(np.isinf(lb)) 432 and np.all(np.isinf(ub))): 433 raise ValueError("Bounds not supported when " 434 "`as_linear_operator` is True.") 435 436 def fun_wrapped(x): 437 f = np.atleast_1d(fun(x, *args, **kwargs)) 438 if f.ndim > 1: 439 raise RuntimeError("`fun` return value has " 440 "more than 1 dimension.") 441 return f 442 443 if f0 is None: 444 f0 = fun_wrapped(x0) 445 else: 446 f0 = np.atleast_1d(f0) 447 if f0.ndim > 1: 448 raise ValueError("`f0` passed has more than 1 dimension.") 449 450 if np.any((x0 < lb) | (x0 > ub)): 451 raise ValueError("`x0` violates bound constraints.") 452 453 if as_linear_operator: 454 if rel_step is None: 455 rel_step = _eps_for_method(x0.dtype, f0.dtype, method) 456 457 return _linear_operator_difference(fun_wrapped, x0, 458 f0, rel_step, method) 459 else: 460 # by default we use rel_step 461 if abs_step is None: 462 h = _compute_absolute_step(rel_step, x0, f0, method) 463 else: 464 # user specifies an absolute step 465 sign_x0 = (x0 >= 0).astype(float) * 2 - 1 466 h = abs_step 467 468 # cannot have a zero step. This might happen if x0 is very large 469 # or small. In which case fall back to relative step. 470 dx = ((x0 + h) - x0) 471 h = np.where(dx == 0, 472 _eps_for_method(x0.dtype, f0.dtype, method) * 473 sign_x0 * np.maximum(1.0, np.abs(x0)), 474 h) 475 476 if method == '2-point': 477 h, use_one_sided = _adjust_scheme_to_bounds( 478 x0, h, 1, '1-sided', lb, ub) 479 elif method == '3-point': 480 h, use_one_sided = _adjust_scheme_to_bounds( 481 x0, h, 1, '2-sided', lb, ub) 482 elif method == 'cs': 483 use_one_sided = False 484 485 if sparsity is None: 486 return _dense_difference(fun_wrapped, x0, f0, h, 487 use_one_sided, method) 488 else: 489 if not issparse(sparsity) and len(sparsity) == 2: 490 structure, groups = sparsity 491 else: 492 structure = sparsity 493 groups = group_columns(sparsity) 494 495 if issparse(structure): 496 structure = csc_matrix(structure) 497 else: 498 structure = np.atleast_2d(structure) 499 500 groups = np.atleast_1d(groups) 501 return _sparse_difference(fun_wrapped, x0, f0, h, 502 use_one_sided, structure, 503 groups, method) 504 505 506def _linear_operator_difference(fun, x0, f0, h, method): 507 m = f0.size 508 n = x0.size 509 510 if method == '2-point': 511 def matvec(p): 512 if np.array_equal(p, np.zeros_like(p)): 513 return np.zeros(m) 514 dx = h / norm(p) 515 x = x0 + dx*p 516 df = fun(x) - f0 517 return df / dx 518 519 elif method == '3-point': 520 def matvec(p): 521 if np.array_equal(p, np.zeros_like(p)): 522 return np.zeros(m) 523 dx = 2*h / norm(p) 524 x1 = x0 - (dx/2)*p 525 x2 = x0 + (dx/2)*p 526 f1 = fun(x1) 527 f2 = fun(x2) 528 df = f2 - f1 529 return df / dx 530 531 elif method == 'cs': 532 def matvec(p): 533 if np.array_equal(p, np.zeros_like(p)): 534 return np.zeros(m) 535 dx = h / norm(p) 536 x = x0 + dx*p*1.j 537 f1 = fun(x) 538 df = f1.imag 539 return df / dx 540 541 else: 542 raise RuntimeError("Never be here.") 543 544 return LinearOperator((m, n), matvec) 545 546 547def _dense_difference(fun, x0, f0, h, use_one_sided, method): 548 m = f0.size 549 n = x0.size 550 J_transposed = np.empty((n, m)) 551 h_vecs = np.diag(h) 552 553 for i in range(h.size): 554 if method == '2-point': 555 x = x0 + h_vecs[i] 556 dx = x[i] - x0[i] # Recompute dx as exactly representable number. 557 df = fun(x) - f0 558 elif method == '3-point' and use_one_sided[i]: 559 x1 = x0 + h_vecs[i] 560 x2 = x0 + 2 * h_vecs[i] 561 dx = x2[i] - x0[i] 562 f1 = fun(x1) 563 f2 = fun(x2) 564 df = -3.0 * f0 + 4 * f1 - f2 565 elif method == '3-point' and not use_one_sided[i]: 566 x1 = x0 - h_vecs[i] 567 x2 = x0 + h_vecs[i] 568 dx = x2[i] - x1[i] 569 f1 = fun(x1) 570 f2 = fun(x2) 571 df = f2 - f1 572 elif method == 'cs': 573 f1 = fun(x0 + h_vecs[i]*1.j) 574 df = f1.imag 575 dx = h_vecs[i, i] 576 else: 577 raise RuntimeError("Never be here.") 578 579 J_transposed[i] = df / dx 580 581 if m == 1: 582 J_transposed = np.ravel(J_transposed) 583 584 return J_transposed.T 585 586 587def _sparse_difference(fun, x0, f0, h, use_one_sided, 588 structure, groups, method): 589 m = f0.size 590 n = x0.size 591 row_indices = [] 592 col_indices = [] 593 fractions = [] 594 595 n_groups = np.max(groups) + 1 596 for group in range(n_groups): 597 # Perturb variables which are in the same group simultaneously. 598 e = np.equal(group, groups) 599 h_vec = h * e 600 if method == '2-point': 601 x = x0 + h_vec 602 dx = x - x0 603 df = fun(x) - f0 604 # The result is written to columns which correspond to perturbed 605 # variables. 606 cols, = np.nonzero(e) 607 # Find all non-zero elements in selected columns of Jacobian. 608 i, j, _ = find(structure[:, cols]) 609 # Restore column indices in the full array. 610 j = cols[j] 611 elif method == '3-point': 612 # Here we do conceptually the same but separate one-sided 613 # and two-sided schemes. 614 x1 = x0.copy() 615 x2 = x0.copy() 616 617 mask_1 = use_one_sided & e 618 x1[mask_1] += h_vec[mask_1] 619 x2[mask_1] += 2 * h_vec[mask_1] 620 621 mask_2 = ~use_one_sided & e 622 x1[mask_2] -= h_vec[mask_2] 623 x2[mask_2] += h_vec[mask_2] 624 625 dx = np.zeros(n) 626 dx[mask_1] = x2[mask_1] - x0[mask_1] 627 dx[mask_2] = x2[mask_2] - x1[mask_2] 628 629 f1 = fun(x1) 630 f2 = fun(x2) 631 632 cols, = np.nonzero(e) 633 i, j, _ = find(structure[:, cols]) 634 j = cols[j] 635 636 mask = use_one_sided[j] 637 df = np.empty(m) 638 639 rows = i[mask] 640 df[rows] = -3 * f0[rows] + 4 * f1[rows] - f2[rows] 641 642 rows = i[~mask] 643 df[rows] = f2[rows] - f1[rows] 644 elif method == 'cs': 645 f1 = fun(x0 + h_vec*1.j) 646 df = f1.imag 647 dx = h_vec 648 cols, = np.nonzero(e) 649 i, j, _ = find(structure[:, cols]) 650 j = cols[j] 651 else: 652 raise ValueError("Never be here.") 653 654 # All that's left is to compute the fraction. We store i, j and 655 # fractions as separate arrays and later construct coo_matrix. 656 row_indices.append(i) 657 col_indices.append(j) 658 fractions.append(df[i] / dx[j]) 659 660 row_indices = np.hstack(row_indices) 661 col_indices = np.hstack(col_indices) 662 fractions = np.hstack(fractions) 663 J = coo_matrix((fractions, (row_indices, col_indices)), shape=(m, n)) 664 return csr_matrix(J) 665 666 667def check_derivative(fun, jac, x0, bounds=(-np.inf, np.inf), args=(), 668 kwargs={}): 669 """Check correctness of a function computing derivatives (Jacobian or 670 gradient) by comparison with a finite difference approximation. 671 672 Parameters 673 ---------- 674 fun : callable 675 Function of which to estimate the derivatives. The argument x 676 passed to this function is ndarray of shape (n,) (never a scalar 677 even if n=1). It must return 1-D array_like of shape (m,) or a scalar. 678 jac : callable 679 Function which computes Jacobian matrix of `fun`. It must work with 680 argument x the same way as `fun`. The return value must be array_like 681 or sparse matrix with an appropriate shape. 682 x0 : array_like of shape (n,) or float 683 Point at which to estimate the derivatives. Float will be converted 684 to 1-D array. 685 bounds : 2-tuple of array_like, optional 686 Lower and upper bounds on independent variables. Defaults to no bounds. 687 Each bound must match the size of `x0` or be a scalar, in the latter 688 case the bound will be the same for all variables. Use it to limit the 689 range of function evaluation. 690 args, kwargs : tuple and dict, optional 691 Additional arguments passed to `fun` and `jac`. Both empty by default. 692 The calling signature is ``fun(x, *args, **kwargs)`` and the same 693 for `jac`. 694 695 Returns 696 ------- 697 accuracy : float 698 The maximum among all relative errors for elements with absolute values 699 higher than 1 and absolute errors for elements with absolute values 700 less or equal than 1. If `accuracy` is on the order of 1e-6 or lower, 701 then it is likely that your `jac` implementation is correct. 702 703 See Also 704 -------- 705 approx_derivative : Compute finite difference approximation of derivative. 706 707 Examples 708 -------- 709 >>> import numpy as np 710 >>> from scipy.optimize import check_derivative 711 >>> 712 >>> 713 >>> def f(x, c1, c2): 714 ... return np.array([x[0] * np.sin(c1 * x[1]), 715 ... x[0] * np.cos(c2 * x[1])]) 716 ... 717 >>> def jac(x, c1, c2): 718 ... return np.array([ 719 ... [np.sin(c1 * x[1]), c1 * x[0] * np.cos(c1 * x[1])], 720 ... [np.cos(c2 * x[1]), -c2 * x[0] * np.sin(c2 * x[1])] 721 ... ]) 722 ... 723 >>> 724 >>> x0 = np.array([1.0, 0.5 * np.pi]) 725 >>> check_derivative(f, jac, x0, args=(1, 2)) 726 2.4492935982947064e-16 727 """ 728 J_to_test = jac(x0, *args, **kwargs) 729 if issparse(J_to_test): 730 J_diff = approx_derivative(fun, x0, bounds=bounds, sparsity=J_to_test, 731 args=args, kwargs=kwargs) 732 J_to_test = csr_matrix(J_to_test) 733 abs_err = J_to_test - J_diff 734 i, j, abs_err_data = find(abs_err) 735 J_diff_data = np.asarray(J_diff[i, j]).ravel() 736 return np.max(np.abs(abs_err_data) / 737 np.maximum(1, np.abs(J_diff_data))) 738 else: 739 J_diff = approx_derivative(fun, x0, bounds=bounds, 740 args=args, kwargs=kwargs) 741 abs_err = np.abs(J_to_test - J_diff) 742 return np.max(abs_err / np.maximum(1, np.abs(J_diff))) 743