如何使用 Numpy 计算导数？答案

【问题标题】：How do I compute derivative using Numpy?如何使用 Numpy 计算导数？
【发布时间】：2012-03-26 16:50:34
【问题描述】：

如何计算函数的导数，例如

y = x²+1

使用numpy?

假设，我想要 x = 5 处的导数...

【问题讨论】：

你需要使用 Sympy：sympy.org/en/index.html Numpy 是 Python 的数值计算库
或者，您想要一种估算导数数值的方法吗？为此，您可以使用有限差分法，但请记住，它们往往会非常嘈杂。

标签： python math numpy

【解决方案1】：

你有四个选择

有限差分不需要外部工具，但容易出现数值错误，如果您处于多变量情况，可能需要一段时间。

如果您的问题足够简单，则符号区分是理想的选择。如今，符号方法变得非常强大。 SymPy 是一个很好的与 NumPy 集成的优秀项目。查看 autowrap 或 lambdify 函数或查看Jensen's blogpost about a similar question。

自动导数非常酷，不容易出现数字错误，但确实需要一些额外的库（谷歌为此，有一些不错的选择）。这是最强大但也是最复杂/最难设置的选择。如果您可以限制自己使用numpy 语法，那么Theano 可能是一个不错的选择。

这是一个使用 SymPy 的示例

In [1]: from sympy import *
In [2]: import numpy as np
In [3]: x = Symbol('x')
In [4]: y = x**2 + 1
In [5]: yprime = y.diff(x)
In [6]: yprime
Out[6]: 2⋅x

In [7]: f = lambdify(x, yprime, 'numpy')
In [8]: f(np.ones(5))
Out[8]: [ 2.  2.  2.  2.  2.]

【讨论】：

对不起，如果这看起来很愚蠢，3.Symbolic Differences 和 4.by hand contrasting 有什么区别？？
当我说“符号分化”时，我的意思是暗示这个过程是由计算机处理的。原则 3 和 4 的区别仅在于工作人员、计算机或程序员。由于一致性、可扩展性和惰性，3 优于 4。如果 3 找不到解决方案，则需要 4。
在第 7 行中，我们制作了 f，一个计算 y 对 x 的导数的函数。在 8 中，我们将这个导数函数应用于全为 1 的向量，并得到全为 2 的向量。这是因为，如第 6 行所述，yprime = 2*x.
为了完整起见，您也可以通过积分来进行微分（参见柯西积分公式），例如实现在mpmath 中（但不确定他们到底在做什么）。
有没有一种简单的方法可以在 numpy 中进行有限差分而不自己实现它？例如我想在预定义点找到函数的梯度。

【解决方案2】：

我能想到的最直接的方法是使用numpy's gradient function：

x = numpy.linspace(0,10,1000)
dx = x[1]-x[0]
y = x**2 + 1
dydx = numpy.gradient(y, dx)

这样，dydx 将使用中心差来计算，并且与 y 具有相同的长度，这与 numpy.diff 不同，numpy.diff 使用前向差并返回 (n-1) 大小向量。

【讨论】：

如果 dx 不是常数怎么办？
@weberc2，在这种情况下，您应该将一个向量除以另一个向量，但手动使用正向和反向导数分别处理边。
或者你可以用一个常数 dx 对 y 进行插值，然后计算梯度。
@Sparkler 感谢您的建议。如果我可以问 2 个小问题，(i) 为什么我们将 dx 传递给 numpy.gradient 而不是 x？ (ii) 我们也可以将你的最后一行做如下：dydx = numpy.gradient(y, numpy.gradient(x))？
从 v1.13 开始，可以使用数组作为第二个参数来指定非均匀间距。请参阅this page 的示例部分。

【解决方案3】：

NumPy 不提供计算导数的通用功能。但是，它可以处理多项式的简单特殊情况：

>>> p = numpy.poly1d([1, 0, 1])
>>> print p
   2
1 x + 1
>>> q = p.deriv()
>>> print q
2 x
>>> q(5)
10

如果您想以数值方式计算导数，则可以在绝大多数应用程序中使用中心差商。对于单点导数，公式类似于

x = 5.0
eps = numpy.sqrt(numpy.finfo(float).eps) * (1.0 + x)
print (p(x + eps) - p(x - eps)) / (2.0 * eps * x)

如果你有一个横坐标数组x 和一个对应的函数值数组y，你可以计算导数的近似值

numpy.diff(y) / numpy.diff(x)

【讨论】：

'为更一般的情况计算数值导数很容易'——我不同意，为一般情况计算数值导数非常困难。您只是选择了表现良好的函数。
2 之后是什么意思 >>>print p ?? （在第 2 行）
@DrStrangeLove：这就是指数。它旨在模拟数学符号。
@SvenMarnach 是最大指数吗？要不然是啥？？为什么它认为指数是2？我们只输入了系数...
@DrStrangeLove：输出应该读作1 * x**2 + 1。他们将2 放在上面的行中，因为它是一个指数。从远处看。

【解决方案4】：

假设您想使用numpy，您可以使用Rigorous definition 在任意点数值计算函数的导数：

def d_fun(x):
    h = 1e-5 #in theory h is an infinitesimal
    return (fun(x+h)-fun(x))/h

您也可以使用Symmetric derivative 以获得更好的结果：

def d_fun(x):
    h = 1e-5
    return (fun(x+h)-fun(x-h))/(2*h)

使用您的示例，完整代码应如下所示：

def fun(x):
    return x**2 + 1

def d_fun(x):
    h = 1e-5
    return (fun(x+h)-fun(x-h))/(2*h)

现在，您可以以数字方式在x=5 处找到导数：

In [1]: d_fun(5)
Out[1]: 9.999999999621423

【讨论】：

这是对导数基本定义的一个很好的应用。我一直想知道为什么我需要学习这个（除了理解割线接近切线的想法）。好吧，现在我知道了。
我会将函数添加为参数：def d_func(func, x):
@rb3652 首先，它用于导出所有导数规则。其次，它也经常用于数学证明。仅举几个应用程序。

【解决方案5】：

我再扔一个方法……

scipy.interpolate 的许多插值样条能够提供导数。因此，使用线性样条 (k=1)，样条的导数（使用 derivative() 方法）应该等效于前向差分。我不完全确定，但我相信使用三次样条导数类似于中心差分导数，因为它使用之前和之后的值来构造三次样条。

from scipy.interpolate import InterpolatedUnivariateSpline

# Get a function that evaluates the linear spline at any x
f = InterpolatedUnivariateSpline(x, y, k=1)

# Get a function that evaluates the derivative of the linear spline at any x
dfdx = f.derivative()

# Evaluate the derivative dydx at each x location...
dydx = dfdx(x)

【讨论】：

刚试过这个，我一直从这个函数 AxisError 得到错误：轴 -1 超出了维度 0 的数组的范围，我在社区上也没有看到任何答案，有什么帮助吗？
将您的问题作为新问题发布并在此处链接。可能需要提供一个导致您的错误发生的示例。我在 interp 函数中遇到的错误通常是因为输入的数据格式不正确 - 例如重复值、错误的维数、其中一个数组意外为空、数据未针对 x 排序或排序时不是有效的函数等。 scipy 可能错误地调用 numpy，但不太可能。检查 x.shape 和 y.shape。看看 np.interp() 是否有效 - 如果没有，它可能会提供更有用的错误。

【解决方案6】：

为了计算梯度，机器学习社区使用 Autograd：

“Efficiently computes derivatives of numpy code.”

安装：

pip install autograd

这是一个例子：

import autograd.numpy as np
from autograd import grad

def fct(x):
    y = x**2+1
    return y

grad_fct = grad(fct)
print(grad_fct(1.0))

它还可以计算复杂函数的梯度，例如多元函数。

【讨论】：

嗨，这个函数可以用来通过提供步长来区分两列数据吗？谢谢
这已在 3 年多前得到回答，但 autograd 不再被开发（只是维护）。 autograd 的主要开发人员现在将您指向github.com/google/jax

【解决方案7】：

您可以使用scipy，这非常简单：

scipy.misc.derivative(func, x0, dx=1.0, n=1, args=(), order=3)

求函数在某点的n阶导数。

在你的情况下：

from scipy.misc import derivative

def f(x):
    return x**2 + 1

derivative(f, 5, dx=1e-6)
# 10.00000000139778

【讨论】：

【解决方案8】：

根据您需要的精度水平，您可以使用简单的微分证明自行计算：

>>> (((5 + 0.1) ** 2 + 1) - ((5) ** 2 + 1)) / 0.1
10.09999999999998
>>> (((5 + 0.01) ** 2 + 1) - ((5) ** 2 + 1)) / 0.01
10.009999999999764
>>> (((5 + 0.0000000001) ** 2 + 1) - ((5) ** 2 + 1)) / 0.0000000001
10.00000082740371

我们实际上不能接受渐变的限制，但它有点有趣。不过你要小心，因为

>>> (((5+0.0000000000000001)**2+1)-((5)**2+1))/0.0000000000000001
0.0

【讨论】：

【解决方案9】：

要计算数值函数的导数，请使用如下所示的二阶有限差分格式： https://youtu.be/5QnToSn_oxk?t=1804

dx = 0.01
x = np.arange(-4, 4+dx, dx)
y = np.sin(x)
n = np.size(x)

yp = np.zeros(n)
yp[0] = (-3*y[0] + 4*y[1] - y[2]) / (2*dx)
yp[n-1] = (3 * y[n-1] - 4*y[n-2] + y[n-3]) / (2*dx)
for j in range(1,n-1):
    yp[j] = (y[j+1] - y[j-1]) / (2*dx)

或者，如果您想使用更高的顺序，请使用： https://youtu.be/5QnToSn_oxk?t=1374

所有内容均来自 Nathan Kutz 在“开始科学计算”课程的讲座。

【讨论】：