如何使用 Numpy 计算导数？

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

Y = x2 + 1

使用 numpy？

比方说，我想要 x = 5的导数值..。

530920 次浏览
小开

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)
小开

根据你所要求的精度水平，你可以自己解决这个问题，使用简单的差异证明:

>>> (((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
小开
最佳答案

你有四个选择

  1. 有限差异
  2. 自动衍生工具
  3. 符号分化
  4. 手工计算导数。

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

如果你的问题足够简单，符号区分是理想的。符号方法最近越来越强大了。SymPy是一个优秀的项目，与 NumPy 集成良好。查看自动包装或 lambdify 功能或检查 关于一个类似问题的博客文章。

自动导数非常酷，不容易出现数字错误，但需要一些额外的库(谷歌为此，有一些很好的选择)。这是最健壮的，但也是最复杂/最难设置的选择。如果您可以将自己限制在 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.]
小开

我能想到的最直接的方法是使用 Numpy 梯度函数 numpy 梯度函数:

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) size 向量。

小开

我再加一个方法。

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)
小开

假设您想使用 numpy，您可以使用 严格的定义在任何点数值计算函数的导数:

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

你也可以使用 对称导数获得更好的效果:

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
小开

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

“ 高效地计算麻点码的导数。”

安装:

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))

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

小开

你可以使用 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
小开

为了计算一个数值函数的导数，使用这个二阶有限差分格式，如下所示: 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)在“科学计算的起源”课程中的讲座。


提交答案