用有限差分法计算导数
fdm的Python项目详细描述
FDM: Finite Difference Methods
fdm用有限差分估计导数。
另请参见FDM.jl。
请参阅docs。
- Installation
- Multivariate Derivatives
- Scalar Derivatives
- Testing Sensitivities in a Reverse-Mode Automatic Differentation Framework
安装
pip install fdm
多元导数
fromfdmimportgradient,jacobian,jvp,hvp
为了便于说明，让我们考虑一个二次函数：
>>>a=np.random.randn(3,3);a=a@a.T>>>aarray([[3.57224794,0.22646662,-1.80432262],[0.22646662,4.72596213,3.46435663],[-1.80432262,3.46435663,3.70938152]])>>>deff(x):...return0.5*x@a@x
考虑以下输入值：
>>>x=np.array([1.0,2.0,3.0])
梯度
>>>grad=gradient(f)>>>grad(x)array([-1.38778668,20.07146076,16.25253519])>>>a@xarray([-1.38778668,20.07146076,16.25253519])
雅可比
>>>jac=jacobian(f)>>>jac(x)array([[-1.38778668,20.07146076,16.25253519]])>>>a@xarray([-1.38778668,20.07146076,16.25253519])
但是jacobian也适用于多值函数。
>>>deff2(x):...returna@x>>>jac2=jacobian(f2)>>>jac2(x)array([[3.57224794,0.22646662,-1.80432262],[0.22646662,4.72596213,3.46435663],[-1.80432262,3.46435663,3.70938152]])>>>aarray([[3.57224794,0.22646662,-1.80432262],[0.22646662,4.72596213,3.46435663],[-1.80432262,3.46435663,3.70938152]])
雅可比向量积（方向导数）
在标量情况下，jvp计算方向导数：
>>>v=np.array([0.5,0.6,0.7])# A direction>>>dir_deriv=jvp(f,v)>>>dir_deriv(x)22.725757753354657>>>np.sum(grad(x)*v)22.72575775335481
在多元情况下，jvp推广到雅可比向量积：
>>>prod=jvp(f2,v)>>>prod(x)array([0.65897811,5.37386023,3.77301973])>>>a@varray([0.65897811,5.37386023,3.77301973])
hessian向量积
>>>prod=hvp(f,v)>>>prod(x)array([[0.6589781,5.37386023,3.77301973]])>>>0.5*(a+a.T)@varray([0.65897811,5.37386023,3.77301973])
标量导数
>>>fromfdmimportcentral_fdm
让我们尝试用
二阶方法，其中我们知道np.sin是条件良好的。
>>>central_fdm(order=2,deriv=1,condition=1)(np.sin,1)-np.cos(1)4.307577627926662e-10
让我们试着用1估计np.sin在1的二阶导数
三阶方法。
>>>central_fdm(order=3,deriv=2,condition=1)(np.sin,1)+np.sin(1)-1.263876664436836e-07
嗯。
让我们检查一下这个三阶方法的准确性。
方法的步长和精度是在调用
FDM.estimate()。
>>>central_fdm(order=3,deriv=2,condition=1).estimate().acc8.733476581980376e-06
我们可能需要更精确一点。让我们检查一下 五阶方法。
>>>central_fdm(order=5,deriv=2,condition=1).estimate().acc7.343652562575155e-10
在1估计np.sin的二阶导数
五阶方法。
>>>central_fdm(order=5,deriv=2,condition=1)(np.sin,1)+np.sin(1)-9.145184609593571e-11
万岁！
最后，让我们验证一下增加订单确实可靠地增加了 准确度。
>>>foriinrange(3,11):...print(central_fdm(order=i,deriv=2,condition=1)(np.sin,1)+np.sin(1))-1.263876664436836e-076.341286606925678e-09-9.145184609593571e-112.7335911312320604e-126.588063428125679e-132.142730437526552e-132.057243264630415e-138.570921750106208e-14
在反向模式自动差分框架中测试灵敏度
考虑函数
defmul(a,b):returna*b
以及它的灵敏度
defs_mul(s_y,y,a,b):returns_y*b,a*s_y
灵敏度s_mul接受输出y的灵敏度s_y，
输出y，函数的参数mul；并返回一个元组
包含与a和b有关的灵敏度。
然后可以使用函数check_sensitivity断言
s_mul的实现是正确的：
>>>fromfdmimportcheck_sensitivity>>>check_sensitivity(mul,s_mul,(2,3))# Test at arguments `2` and `3`.
假设实现是错误的，例如
defs_mul_wrong(s_y,y,a,b):returns_y*b,b*s_y# Used `b` instead of `a` for the second sensitivity!
那么check_sensitivity应该抛出一个AssertionError：
>>>check_sensitivity(mul,s_mul,(2,3))AssertionError:Sensitivityofargument2offunction"mul"didnotmatchnumericalestimate.
欢迎加入QQ群-->： 979659372
