scipy.differentiate.

derivative#

scipy.differentiate.derivative(f, x, *, args=(), tolerances=None, maxiter=10, order=8, initial_step=0.5, step_factor=2.0, step_direction=0, preserve_shape=False, callback=None)[源代码]#

数值评估逐元素实标量函数的导数。

对于 f 的每个输出元素，derivative 使用有限差分微分逼近 f 在 x 对应元素处的一阶导数。

当 x、step_direction 和 args 包含（可广播的）数组时，此函数逐元素工作。

参数:

f可调用对象

需要求导的函数。签名必须是

f(xi: ndarray, *argsi) -> ndarray

f(xi, *argsi)，其中 xi 的每个元素都是有限实数，argsi 是一个元组，可能包含任意数量的与 xi 可广播的数组。f 必须是一个逐元素函数：对于有效的索引 j，每个标量元素 f(xi)[j] 必须等于 f(xi[j])。它不得改变数组 xi 或 argsi 中的数组。

x浮点数 array_like

求导的横坐标。必须与 args 和 step_direction 可广播。

argsarray_like 元组，可选

要传递给 f 的其他位置数组参数。数组必须彼此之间以及与 init 的数组可广播。如果所需根的可调用对象需要与 x 不可广播的参数，请使用 f 包裹该可调用对象，使得 f 只接受 x 和可广播的 *args。

tolerances浮点数字典，可选

绝对和相对容差。字典的有效键为

atol - 导数的绝对容差
rtol - 导数的相对容差

当 res.error < atol + rtol * abs(res.df) 时，迭代将停止。默认的 atol 是相应 dtype 的最小正规数，默认的 rtol 是相应 dtype 精度平方根。

orderint, 默认值：8

要使用的有限差分公式的（正整数）阶数。奇数将向上舍入为下一个偶数。

initial_step浮点数 array_like, 默认值：0.5

有限差分导数逼近的（绝对）初始步长。

step_factor浮点数，默认值：2.0

每次迭代中步长减小的因子；即，迭代 1 中的步长为 initial_step/step_factor。如果 step_factor < 1，则后续步骤将大于初始步长；如果不需要小于某个阈值的步长（例如，由于减法抵消错误），这可能很有用。

maxiterint, 默认值：10

要执行的算法的最大迭代次数。请参见注释。

step_direction整数 array_like

表示有限差分步进方向的数组（当 x 接近函数域的边界时使用。）必须与 x 和所有 args 可广播。如果为 0（默认值），则使用中心差分；如果为负值（例如 -1），则步长为非正值；如果为正值（例如 1），则所有步长为非负值。

preserve_shapebool, 默认值：False

在以下内容中，“f 的参数”指的是数组 xi 和 argsi 中的任何数组。令 shape 为 x 和 args 的所有元素的广播形状（这在概念上不同于传递给 f 的 xi` 和 ``argsi）。

当 preserve_shape=False（默认值）时，f 必须接受任何可广播形状的参数。
当 preserve_shape=True 时，f 必须接受形状为 shape 或 shape + (n,) 的参数，其中 (n,) 是正在评估函数的横坐标数。

在这两种情况下，对于 xi 中的每个标量元素 xi[j]，f 返回的数组必须在同一索引处包含标量 f(xi[j])。因此，输出的形状始终为输入 xi 的形状。

请参见示例。

callback可调用对象，可选

在第一次迭代之前和每次迭代之后调用的可选用户提供的函数。调用方式为 callback(res)，其中 res 是一个类似于 derivative 返回的 _RichResult （但包含所有变量的当前迭代值）。如果 callback 引发 StopIteration，则算法将立即终止，并且 derivative 将返回结果。callback 不得改变 res 或其属性。

返回:

res_RichResult

类似于 scipy.optimize.OptimizeResult 实例的对象，具有以下属性。这些描述的编写方式就好像这些值将是标量一样；但是，如果 f 返回数组，则输出将是具有相同形状的数组。

successbool 数组

算法成功终止的位置为 True（状态为 0）；否则为 False。

statusint 数组

表示算法退出状态的整数。

0 : 算法收敛到指定的容差。
-1 : 误差估计值增大，因此迭代终止。
-2 : 达到最大迭代次数。
-3 : 遇到非有限值。
-4 : 迭代被 callback 终止。
1 : 算法正常进行（仅在 callback 中）。

df浮点数数组

如果算法成功终止，则 f 在 x 处的导数。

error浮点数数组

误差估计值：当前导数估计值与前一次迭代估计值之差的绝对值。

nitint 数组

执行的算法迭代次数。

nfevint 数组

评估 f 的点的数量。

x浮点数数组

评估 f 的导数的值（在与 args 和 step_direction 广播后）。

另请参见

jacobian, hessian

说明

该实现受到了 jacobi [1]、numdifftools [2] 和 DERIVEST [3] 的启发，但其实现更直接地（并且可以说比较幼稚地）遵循了泰勒级数的理论。在第一次迭代中，导数使用阶数为 order 且最大步长为 initial_step 的有限差分公式进行估计。在后续的每次迭代中，最大步长减小 step_factor，并再次估计导数，直到达到终止条件。误差估计是当前导数近似值与前一次迭代的导数近似值之差的绝对值。

有限差分公式的模板设计为横坐标是“嵌套的”：在第一次迭代中，对 order + 1 个点计算 f 后，在后续的每次迭代中，f 仅在两个新点进行计算；有限差分公式所需的 order - 1 个先前计算的函数值被重用，并且两个函数值（在离 x 最远的点处计算的值）未被使用。

步长是绝对的。当步长相对于 x 的大小较小时，精度会损失；例如，如果 x 为 1e20，则默认的初始步长 0.5 无法被解析。因此，请考虑为较大的 x 值使用更大的初始步长。

默认容差在真导数恰好为零的点处难以满足。如果导数可能恰好为零，请考虑指定绝对容差（例如 atol=1e-12）以提高收敛性。

参考文献

[1]

Hans Dembinski (@HDembinski). jacobi. HDembinski/jacobi

[2]

Per A. Brodtkorb and John D’Errico. numdifftools. https://numdifftools.readthedocs.io/en/latest/

[3]

John D’Errico. DERIVEST: Adaptive Robust Numerical Differentiation. https://www.mathworks.com/matlabcentral/fileexchange/13490-adaptive-robust-numerical-differentiation

[4]

Numerical Differentition. Wikipedia. https://en.wikipedia.org/wiki/Numerical_differentiation

示例

计算几个点 x 处 np.exp 的导数。

>>> import numpy as np
>>> from scipy.differentiate import derivative
>>> f = np.exp
>>> df = np.exp  # true derivative
>>> x = np.linspace(1, 2, 5)
>>> res = derivative(f, x)
>>> res.df  # approximation of the derivative
array([2.71828183, 3.49034296, 4.48168907, 5.75460268, 7.3890561 ])
>>> res.error  # estimate of the error
array([7.13740178e-12, 9.16600129e-12, 1.17594823e-11, 1.51061386e-11,
       1.94262384e-11])
>>> abs(res.df - df(x))  # true error
array([2.53130850e-14, 3.55271368e-14, 5.77315973e-14, 5.59552404e-14,
       6.92779167e-14])

展示当步长减小时近似值的收敛性。每次迭代，步长减小 step_factor，因此，对于足够小的初始步长，每次迭代都会将误差减小 1/step_factor**order 倍，直到有限精度算术阻止进一步改进。

>>> import matplotlib.pyplot as plt
>>> iter = list(range(1, 12))  # maximum iterations
>>> hfac = 2  # step size reduction per iteration
>>> hdir = [-1, 0, 1]  # compare left-, central-, and right- steps
>>> order = 4  # order of differentiation formula
>>> x = 1
>>> ref = df(x)
>>> errors = []  # true error
>>> for i in iter:
...     res = derivative(f, x, maxiter=i, step_factor=hfac,
...                      step_direction=hdir, order=order,
...                      # prevent early termination
...                      tolerances=dict(atol=0, rtol=0))
...     errors.append(abs(res.df - ref))
>>> errors = np.array(errors)
>>> plt.semilogy(iter, errors[:, 0], label='left differences')
>>> plt.semilogy(iter, errors[:, 1], label='central differences')
>>> plt.semilogy(iter, errors[:, 2], label='right differences')
>>> plt.xlabel('iteration')
>>> plt.ylabel('error')
>>> plt.legend()
>>> plt.show()

../../_images/scipy-differentiate-derivative-1_00_00.png

>>> (errors[1, 1] / errors[0, 1], 1 / hfac**order)
(0.06215223140159822, 0.0625)

该实现在 x、step_direction 和 args 上进行了向量化处理。该函数在第一次迭代之前被调用一次以执行输入验证和标准化，然后在每次迭代中调用一次。

>>> def f(x, p):
...     f.nit += 1
...     return x**p
>>> f.nit = 0
>>> def df(x, p):
...     return p*x**(p-1)
>>> x = np.arange(1, 5)
>>> p = np.arange(1, 6).reshape((-1, 1))
>>> hdir = np.arange(-1, 2).reshape((-1, 1, 1))
>>> res = derivative(f, x, args=(p,), step_direction=hdir, maxiter=1)
>>> np.allclose(res.df, df(x, p))
True
>>> res.df.shape
(3, 5, 4)
>>> f.nit
2

默认情况下，preserve_shape 为 False，因此可调用对象 f 可以使用任何可广播的形状的数组进行调用。例如

>>> shapes = []
>>> def f(x, c):
...    shape = np.broadcast_shapes(x.shape, c.shape)
...    shapes.append(shape)
...    return np.sin(c*x)
>>>
>>> c = [1, 5, 10, 20]
>>> res = derivative(f, 0, args=(c,))
>>> shapes
[(4,), (4, 8), (4, 2), (3, 2), (2, 2), (1, 2)]

要理解这些形状的来源，并更好地理解 derivative 如何计算准确的结果，请注意，较高的 c 值对应于较高频率的正弦曲线。较高频率的正弦曲线会使函数的导数变化更快，因此需要更多的函数计算才能达到目标精度

>>> res.nfev
array([11, 13, 15, 17], dtype=int32)

初始的 shape，(4,)，对应于在单个横坐标和所有四个频率处计算函数值；这用于输入验证并确定存储结果的数组的大小和 dtype。下一个形状对应于在横坐标的初始网格和所有四个频率处计算函数值。连续调用函数会在两个更多横坐标处计算函数值，从而将近似值的有效阶数增加两个。但是，在以后的函数计算中，函数会在较少的频率处计算，因为相应的导数已经收敛到所需的容差。这节省了函数计算以提高性能，但这需要函数接受任何形状的参数。

“向量值”函数不太可能满足此要求。例如，考虑

>>> def f(x):
...    return [x, np.sin(3*x), x+np.sin(10*x), np.sin(20*x)*(x-1)**2]

这个被积函数与编写的 derivative 不兼容；例如，输出的形状将与 x 的形状不同。通过引入额外的参数，可以将此类函数转换为兼容的形式，但这会很不方便。在这种情况下，一个更简单的解决方案是使用 preserve_shape。

>>> shapes = []
>>> def f(x):
...     shapes.append(x.shape)
...     x0, x1, x2, x3 = x
...     return [x0, np.sin(3*x1), x2+np.sin(10*x2), np.sin(20*x3)*(x3-1)**2]
>>>
>>> x = np.zeros(4)
>>> res = derivative(f, x, preserve_shape=True)
>>> shapes
[(4,), (4, 8), (4, 2), (4, 2), (4, 2), (4, 2)]

这里，x 的形状为 (4,)。使用 preserve_shape=True，可以使用形状为 (4,) 或 (4, n) 的参数 x 调用该函数，这正是我们观察到的。