scipy.differentiate.

jacobian#

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

数值计算函数的雅可比矩阵。

参数:
f可调用对象

需要计算雅可比矩阵的函数。其签名必须为

f(xi: ndarray) -> ndarray

其中 xi 的每个元素都是有限实数。如果待微分函数接受额外的参数，请将其封装（例如使用 functools.partial 或 lambda）并将封装后的可调用对象传递给 jacobian。 f 不得修改数组 xi。有关向量化以及输入和输出维度，请参阅“备注”部分。

x浮点型类数组对象

计算雅可比矩阵的点。必须至少有一个维度。有关维度和向量化，请参阅“备注”部分。

tolerances浮点数字典，可选

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

  • atol - 导数的绝对容差

  • rtol - 导数的相对容差

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

maxiter整数，默认值：10

算法执行的最大迭代次数。请参阅“备注”部分。

order整数，默认值：8

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

initial_step浮点型类数组对象，默认值：0.5

有限差分导数近似的（绝对）初始步长。必须能与 x 和 step_direction 广播。

step_factor浮点数，默认值：2.0

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

step_direction整数类数组对象

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

返回:
res_RichResult

一个类似于 scipy.optimize.OptimizeResult 实例的对象，具有以下属性。描述假定值为标量；但是，如果 f 返回一个数组，则输出也将是相同形状的数组。

success布尔数组

当算法成功终止时（状态 0）为 True；否则为 False。

status整数数组

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

  • 0 : 算法已收敛到指定容差。

  • -1 : 误差估计增加，因此迭代终止。

  • -2 : 达到最大迭代次数。

  • -3 : 遇到非有限值。

df浮点数组

如果算法成功终止，则为函数 f 在点 x 处的雅可比矩阵。

error浮点数组

误差估计：当前雅可比矩阵估计值与上次迭代估计值之间差的幅度。

nit整数数组

算法执行的迭代次数。

nfev整数数组

函数 f 被评估的点的数量。

属性的每个元素都与 df 的相应元素相关联。例如，nfev 的元素 i 是为了计算 df 的元素 i 而评估函数 f 的点数。

另请参阅

derivative, hessian

备注

假设我们希望计算函数 \(f: \mathbf{R}^m \rightarrow \mathbf{R}^n\) 的雅可比矩阵。将变量 m 和 n 分别赋值为 \(m\) 和 \(n\) 的正整数值，并让 ... 表示任意的整数元组。如果我们要在一个单点处评估雅可比矩阵，那么

  • 参数 x 必须是形状为 (m,) 的数组

  • 参数 f 必须是向量化的，以接受形状为 (m, ...) 的数组。第一个轴表示 \(f\) 的 \(m\) 个输入；其余用于在一次调用中在多个点评估函数。

  • 参数 f 必须返回形状为 (n, ...) 的数组。第一个轴表示 \(f\) 的 \(n\) 个输出；其余用于在多个点评估函数的结果。

  • 结果对象的属性 df 将是形状为 (n, m) 的数组，即雅可比矩阵。

此函数也是向量化的，即可以在一次调用中在 k 个点评估雅可比矩阵。在这种情况下，x 将是形状为 (m, k) 的数组，f 将接受形状为 (m, k, ...) 的数组并返回形状为 (n, k, ...) 的数组，并且结果的 df 属性将具有形状 (n, m, k)。

假设所需的可调用对象 f_not_vectorized 未向量化；它只能接受形状为 (m,) 的数组。满足所需接口的一个简单解决方案是将 f_not_vectorized 封装如下：

def f(x):
    return np.apply_along_axis(f_not_vectorized, axis=0, arr=x)

或者，假设所需的可调用对象 f_vec_q 已向量化，但仅适用于形状为 (m, q) 的二维数组。为了满足所需接口，请考虑

def f(x):
    m, batch = x.shape[0], x.shape[1:]  # x.shape is (m, ...)
    x = np.reshape(x, (m, -1))  # `-1` is short for q = prod(batch)
    res = f_vec_q(x)  # pass shape (m, q) to function
    n = res.shape[0]
    return np.reshape(res, (n,) + batch)  # return shape (n, ...)

然后将封装后的可调用对象 f 作为 jacobian 的第一个参数传递。

参考文献

[1]

雅可比矩阵与行列式，维基百科，https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant

示例

Rosenbrock 函数将 \(\mathbf{R}^m \rightarrow \mathbf{R}\)；SciPy 实现 scipy.optimize.rosen 已向量化，可接受形状为 (m, p) 的数组并返回形状为 p 的数组。假设我们希望在 [0.5, 0.5, 0.5] 处评估雅可比矩阵（也称为梯度，因为函数返回标量）。

>>> import numpy as np
>>> from scipy.differentiate import jacobian
>>> from scipy.optimize import rosen, rosen_der
>>> m = 3
>>> x = np.full(m, 0.5)
>>> res = jacobian(rosen, x)
>>> ref = rosen_der(x)  # reference value of the gradient
>>> res.df, ref
(array([-51.,  -1.,  50.]), array([-51.,  -1.,  50.]))

作为具有多个输出的函数的示例，请考虑 [1] 中的示例 4。

>>> def f(x):
...     x1, x2, x3 = x
...     return [x1, 5*x3, 4*x2**2 - 2*x3, x3*np.sin(x1)]

真实的雅可比矩阵由以下给出

>>> def df(x):
...         x1, x2, x3 = x
...         one = np.ones_like(x1)
...         return [[one, 0*one, 0*one],
...                 [0*one, 0*one, 5*one],
...                 [0*one, 8*x2, -2*one],
...                 [x3*np.cos(x1), 0*one, np.sin(x1)]]

在任意点评估雅可比矩阵。

>>> rng = np.random.default_rng()
>>> x = rng.random(size=3)
>>> res = jacobian(f, x)
>>> ref = df(x)
>>> res.df.shape == (4, 3)
True
>>> np.allclose(res.df, ref)
True

在一次调用中在10个任意点评估雅可比矩阵。

>>> x = rng.random(size=(3, 10))
>>> res = jacobian(f, x)
>>> ref = df(x)
>>> res.df.shape == (4, 3, 10)
True
>>> np.allclose(res.df, ref)
True
本页目录