python计算函数梯度,梯度算子的计算函数

如何使用python计算常微分方程？

常用形式

odeint(func, y0, t,args,Dfun)

一般这种形式就够用了。

下面是官方的例子，求解的是

D(D(y1))-t*y1=0

为了方便，采取D=d/dt。如果我们令初值

y1(0) = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)

D(y1)(0) = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)

这个微分方程的解y1=airy(t)。

令D(y1)=y0，就有这个常微分方程组。

D(y0)=t*y1

D(y1)=y0

Python求解该微分方程。

from scipy.integrate import odeint

from scipy.special import gamma, airy

y1_0 = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)

y0_0 = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)

y0 = [y0_0, y1_0]

def func(y, t):

... return [t*y[1],y[0]]

def gradient(y,t):

... return [[0,t],[1,0]]

x = arange(0,4.0, 0.01)

t = x

ychk = airy(x)[0]

y = odeint(func, y0, t)

y2 = odeint(func, y0, t, Dfun=gradient)

print ychk[:36:6]

[ 0.355028 0.339511 0.324068 0.308763 0.293658 0.278806]

print y[:36:6,1]

[ 0.355028 0.339511 0.324067 0.308763 0.293658 0.278806]

print y2[:36:6,1]

[ 0.355028 0.339511 0.324067 0.308763 0.293658 0.278806]

得到的解与精确值相比，误差相当小。

=======================================================================================================

args是额外的参数。

用法请参看下面的例子。这是一个洛仑兹曲线的求解，并且用matplotlib绘出空间曲线图。（来自《python科学计算》）

from scipy.integrate import odeint

import numpy as np

def lorenz(w, t, p, r, b):

# 给出位置矢量w，和三个参数p, r, b 计算出

# dx/dt, dy/dt, dz/dt 的值

x, y, z = w

# 直接与lorenz 的计算公式对应

return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])

t = np.arange(0, 30, 0.01) # 创建时间点

# 调用ode 对lorenz 进行求解, 用两个不同的初始值

track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))

track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0))

# 绘图

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt

fig = plt.figure()

ax = Axes3D(fig)

ax.plot(track1[:,0], track1[:,1], track1[:,2])

ax.plot(track2[:,0], track2[:,1], track2[:,2])

plt.show()

===========================================================================

scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)

计算常微分方程（组）

使用 FORTRAN库odepack中的lsoda解常微分方程。这个函数一般求解初值问题。

参数：

func : callable(y, t0, ...) 计算y在t0 处的导数。

y0 : 数组 y的初值条件（可以是矢量）

t : 数组 为求出y，这是一个时间点的序列。初值点应该是这个序列的第一个元素。

args : 元组 func的额外参数

Dfun : callable(y, t0, ...) 函数的梯度(Jacobian)。即雅可比多项式。

col_deriv : boolean. True，Dfun定义列向导数（更快），否则Dfun会定义横排导数

full_output : boolean 可选输出，如果为True 则返回一个字典，作为第二输出。

printmessg : boolean 是否打印convergence 消息。

返回： y : array, shape (len(y0), len(t))

数组，包含y值，每一个对应于时间序列中的t。初值y0 在第一排。

infodict : 字典，只有full_output == True 时，才会返回。

字典包含额为的输出信息。

键值：

‘hu’ vector of step sizes successfully used for each time step.

‘tcur’ vector with the value of t reached for each time step. (will always be at least as large as the input times).

‘tolsf’ vector of tolerance scale factors, greater than 1.0, computed when a request for too much accuracy was detected.

‘tsw’ value of t at the time of the last method switch (given for each time step)

‘nst’ cumulative number of time steps

‘nfe’ cumulative number of function evaluations for each time step

‘nje’ cumulative number of jacobian evaluations for each time step

‘nqu’ a vector of method orders for each successful step.

‘imxer’index of the component of largest magnitude in the weighted local error vector (e / ewt) on an error return, -1 otherwise.

‘lenrw’ the length of the double work array required.

‘leniw’ the length of integer work array required.

‘mused’a vector of method indicators for each successful time step: 1: adams (nonstiff), 2: bdf (stiff)

其他参数，官方网站和文档都没有明确说明。相关的资料，暂时也找不到。

相关问答

Q1: Python气象数据处理与绘图(2)：常用数据计算方法

对于气象绘图来讲，第一步是对数据的处理，通过各类公式，或者统计方法将原始数据处理为目标数据。

按照气象统计课程的内容，我给出了一些常用到的统计方法的对应函数：

在计算气候态，区域平均时均要使用到求均值函数，对应NCL中的dim_average函数，在python中通常使用np.mean()函数

numpy.mean(a, axis, dtype)

假设a为[time,lat,lon]的数据，那么

需要特别注意的是，气象数据中常有缺测，在NCL中，使用求均值函数会自动略过，而在python中，当任意一数与缺测(np.nan)计算的结果均为np.nan，比如求[1,2,3,4，np.nan]的平均值，结果为np.nan

因此，当数据存在缺测数据时，通常使用np.nanmean()函数，用法同上，此时[1,2,3,4，np.nan]的平均值为(1+2+3+4)/4 = 2.5

同样的，求某数组最大最小值时也有np.nanmax(), np.nanmin()函数来补充np.max(), np.min()的不足。

其他很多np的计算函数也可以通过在前边加‘nan’来使用。

另外，

也可以直接将a中缺失值全部填充为0。

np.std(a, axis, dtype)

用法同np.mean()

在NCL中有直接求数据标准化的函数dim_standardize()

其实也就是一行的事，根据需要指定维度即可。

皮尔逊相关系数：

相关可以说是气象科研中最常用的方法之一了，numpy函数中的np.corrcoef(x, y)就可以实现相关计算。但是在这里我推荐scipy.stats中的函数来计算相关系数：

这个函数缺点和有点都很明显，优点是可以直接返回相关系数R及其P值，这避免了我们进一步计算置信度。而缺点则是该函数只支持两个一维数组的计算，也就是说当我们需要计算一个场和一个序列的相关时，我们需要循环来实现。

其中a[time,lat,lon]，b[time]

(NCL中为regcoef()函数)

同样推荐Scipy库中的stats.linregress(x,y)函数：

slop: 回归斜率

intercept：回归截距

r_value： 相关系数

p_value： P值

std_err： 估计标准误差

直接可以输出P值，同样省去了做置信度检验的过程，遗憾的是仍需同相关系数一样循环计算。

Q2: Python怎么做最优化

一、概观scipy中的optimize子包中提供了常用的最优化算法函数实现。我们可以直接调用这些函数完成我们的优化问题。optimize中函数最典型的特点就是能够从函数名称上看出是使用了什么算法。下面optimize包中函数的概览：1.非线性最优化fmin -- 简单Nelder-Mead算法fmin_powell -- 改进型Powell法fmin_bfgs -- 拟Newton法fmin_cg -- 非线性共轭梯度法fmin_ncg -- 线性搜索Newton共轭梯度法leastsq -- 最小二乘2.有约束的多元函数问题fmin_l_bfgs_b ---使用L-BFGS-B算法fmin_tnc ---梯度信息fmin_cobyla ---线性逼近fmin_slsqp ---序列最小二乘法nnls ---解|| Ax - b ||_2 for x=03.全局优化anneal ---模拟退火算法brute --强力法4.标量函数fminboundbrentgoldenbracket5.拟合curve_fit-- 使用非线性最小二乘法拟合6.标量函数求根brentq ---classic Brent (1973)brenth ---A variation on the classic Brent（1980）ridder ---Ridder是提出这个算法的人名bisect ---二分法newton ---牛顿法fixed_point7.多维函数求根fsolve ---通用broyden1 ---Broyden’s first Jacobian approximation.broyden2 ---Broyden’s second Jacobian approximationnewton_krylov ---Krylov approximation for inverse Jacobiananderson ---extended Anderson mixingexcitingmixing ---tuned diagonal Jacobian approximationlinearmixing ---scalar Jacobian approximationdiagbroyden ---diagonal Broyden Jacobian approximation8.实用函数line_search ---找到满足强Wolfe的alpha值check_grad ---通过和前向有限差分逼近比较检查梯度函数的正确性二、实战非线性最优化fmin完整的调用形式是：fmin(func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None)不过我们最常使用的就是前两个参数。一个描述优化问题的函数以及初值。后面的那些参数我们也很容易理解。如果您能用到，请自己研究。下面研究一个最简单的问题，来感受这个函数的使用方法：f(x)=x**2-4*x+8，我们知道，这个函数的最小值是4，在x=2的时候取到。from scipy.optimize import fmin #引入优化包def myfunc(x):return x**2-4*x+8 #定义函数x0 = [1.3] #猜一个初值xopt = fmin(myfunc, x0) #求解print xopt #打印结果运行之后，给出的结果是：Optimization terminated successfully.Current function value: 4.000000Iterations: 16Function evaluations: 32[ 2.00001953]程序准确的计算得出了最小值，不过最小值点并不是严格的2，这应该是由二进制机器编码误差造成的。除了fmin_ncg必须提供梯度信息外，其他几个函数的调用大同小异，完全类似。我们不妨做一个对比：from scipy.optimize import fmin,fmin_powell,fmin_bfgs,fmin_cgdef myfunc(x):return x**2-4*x+8x0 = [1.3]xopt1 = fmin(myfunc, x0)print xopt1printxopt2 = fmin_powell(myfunc, x0)print xopt2printxopt3 = fmin_bfgs(myfunc, x0)print xopt3printxopt4 = fmin_cg(myfunc,x0)print xopt4给出的结果是：Optimization terminated successfully.Current function value: 4.000000Iterations: 16Function evaluations: 32[ 2.00001953]Optimization terminated successfully.Current function value: 4.000000Iterations: 2Function evaluations: 531.99999999997Optimization terminated successfully.Current function value: 4.000000Iterations: 2Function evaluations: 12Gradient evaluations: 4[ 2.00000001]Optimization terminated successfully.Current function value: 4.000000Iterations: 2Function evaluations: 15Gradient evaluations: 5[ 2.]我们可以根据给出的消息直观的判断算法的执行情况。每一种算法数学上的问题，请自己看书学习。个人感觉，如果不是纯研究数学的工作，没必要搞清楚那些推导以及定理云云。不过，必须了解每一种算法的优劣以及能力所及。在使用的时候，不妨多种算法都使用一下，看看效果分别如何，同时，还可以互相印证算法失效的问题。在from scipy.optimize import fmin之后，就可以使用help(fmin)来查看fmin的帮助信息了。帮助信息中没有例子，但是给出了每一个参数的含义说明，这是调用函数时候的最有价值参考。有源码研究癖好的，或者当你需要改进这些已经实现的算法的时候，可能需要查看optimize中的每种算法的源代码。在这里：https:/ / github. com/scipy/scipy/blob/master/scipy/optimize/optimize.py聪明的你肯定发现了，顺着这个链接往上一级、再往上一级，你会找到scipy的几乎所有源码！

Q3: python如何实现求函数的在一个连续区间的最值？

如果函数是确定的，可以用导数的方法进行计算，但是如果函数是不确定的，就需要用优化的方法来处理了，比如常用的梯度上升法，模拟退火等，希望可以帮到你。

Q4: 如何通过Python进行深度学习？

作者 | Vihar Kurama

编译 | 荷叶

来源 | 云栖社区

摘要：深度学习背后的主要原因是人工智能应该从人脑中汲取灵感。本文就用一个小例子无死角的介绍一下深度学习！

人脑模拟

深度学习背后的主要原因是人工智能应该从人脑中汲取灵感。此观点引出了“神经网络”这一术语。人脑中包含数十亿个神经元，它们之间有数万个连接。很多情况下，深度学习算法和人脑相似，因为人脑和深度学习模型都拥有大量的编译单元（神经元），这些编译单元（神经元）在独立的情况下都不太智能，但是当他们相互作用时就会变得智能。

我认为人们需要了解到深度学习正在使得很多幕后的事物变得更好。深度学习已经应用于谷歌搜索和图像搜索，你可以通过它搜索像“拥抱”这样的词语以获得相应的图像。－杰弗里·辛顿

神经元

神经网络的基本构建模块是人工神经元，它模仿了人类大脑的神经元。这些神经元是简单、强大的计算单元，拥有加权输入信号并且使用激活函数产生输出信号。这些神经元分布在神经网络的几个层中。

inputs 输入 outputs 输出 weights 权值 activation 激活

人工神经网络的工作原理是什么？

深度学习由人工神经网络构成，该网络模拟了人脑中类似的网络。当数据穿过这个人工网络时，每一层都会处理这个数据的一方面，过滤掉异常值，辨认出熟悉的实体，并产生最终输出。

输入层：该层由神经元组成，这些神经元只接收输入信息并将它传递到其他层。输入层的图层数应等于数据集里的属性或要素的数量。输出层：输出层具有预测性，其主要取决于你所构建的模型类型。隐含层：隐含层处于输入层和输出层之间，以模型类型为基础。隐含层包含大量的神经元。处于隐含层的神经元会先转化输入信息，再将它们传递出去。随着网络受训练，权重得到更新，从而使其更具前瞻性。

神经元的权重

权重是指两个神经元之间的连接的强度或幅度。你如果熟悉线性回归的话，可以将输入的权重类比为我们在回归方程中用的系数。权重通常被初始化为小的随机数值，比如数值0-1。

前馈深度网络

前馈监督神经网络曾是第一个也是最成功的学习算法。该网络也可被称为深度网络、多层感知机（MLP）或简单神经网络，并且阐明了具有单一隐含层的原始架构。每个神经元通过某个权重和另一个神经元相关联。

该网络处理向前处理输入信息，激活神经元，最终产生输出值。在此网络中，这称为前向传递。

inputlayer 输入层 hidden layer 输出层 output layer 输出层

激活函数

激活函数就是求和加权的输入到神经元的输出的映射。之所以称之为激活函数或传递函数是因为它控制着激活神经元的初始值和输出信号的强度。

用数学表示为：

我们有许多激活函数，其中使用最多的是整流线性单元函数、双曲正切函数和solfPlus函数。

激活函数的速查表如下：

反向传播

在网络中，我们将预测值与预期输出值相比较，并使用函数计算其误差。然后，这个误差会传回这个网络，每次传回一个层，权重也会根绝其导致的误差值进行更新。这个聪明的数学法是反向传播算法。这个步骤会在训练数据的所有样本中反复进行，整个训练数据集的网络更新一轮称为一个时期。一个网络可受训练数十、数百或数千个时期。

prediction error 预测误差

代价函数和梯度下降

代价函数度量了神经网络对给定的训练输入和预期输出“有多好”。该函数可能取决于权重、偏差等属性。

代价函数是单值的，并不是一个向量，因为它从整体上评估神经网络的性能。在运用梯度下降最优算法时，权重在每个时期后都会得到增量式地更新。

兼容代价函数

用数学表述为差值平方和：

target 目标值 output 输出值

权重更新的大小和方向是由在代价梯度的反向上采取步骤计算出的。

其中η 是学习率

其中Δw是包含每个权重系数w的权重更新的向量，其计算方式如下：

target 目标值 output 输出值

图表中会考虑到单系数的代价函数

initial weight 初始权重 gradient 梯度 global cost minimum 代价极小值

在导数达到最小误差值之前，我们会一直计算梯度下降，并且每个步骤都会取决于斜率（梯度）的陡度。

多层感知器（前向传播）

这类网络由多层神经元组成，通常这些神经元以前馈方式（向前传播）相互连接。一层中的每个神经元可以直接连接后续层的神经元。在许多应用中，这些网络的单元会采用S型函数或整流线性单元（整流线性激活）函数作为激活函数。

现在想想看要找出处理次数这个问题，给定的账户和家庭成员作为输入

要解决这个问题，首先，我们需要先创建一个前向传播神经网络。我们的输入层将是家庭成员和账户的数量，隐含层数为1， 输出层将是处理次数。

将图中输入层到输出层的给定权重作为输入：家庭成员数为2、账户数为3。

现在将通过以下步骤使用前向传播来计算隐含层（i，j）和输出层（k）的值。

步骤：

1， 乘法－添加方法。

2， 点积（输入＊权重）。

3，一次一个数据点的前向传播。

4， 输出是该数据点的预测。

i的值将从相连接的神经元所对应的输入值和权重中计算出来。

i = (2 * 1) + (3* 1) → i = 5

同样地，j = (2 * -1) + (3 * 1) → j =1

K = (5 * 2) + (1* -1) → k = 9

Python中的多层感知器问题的解决

激活函数的使用

为了使神经网络达到其最大预测能力，我们需要在隐含层应用一个激活函数，以捕捉非线性。我们通过将值代入方程式的方式来在输入层和输出层应用激活函数。

这里我们使用整流线性激活（ReLU）:

用Keras开发第一个神经网络

关于Keras：

Keras是一个高级神经网络的应用程序编程接口，由Python编写，能够搭建在TensorFlow,CNTK,或Theano上。

使用PIP在设备上安装Keras,并且运行下列指令。

在keras执行深度学习程序的步骤

1，加载数据；

2，创建模型；

3，编译模型；

4，拟合模型；

5，评估模型。

开发Keras模型

全连接层用Dense表示。我们可以指定层中神经元的数量作为第一参数，指定初始化方法为第二参数，即初始化参数，并且用激活参数确定激活函数。既然模型已经创建，我们就可以编译它。我们在底层库（也称为后端）用高效数字库编译模型，底层库可以用Theano或TensorFlow。目前为止，我们已经完成了创建模型和编译模型，为进行有效计算做好了准备。现在可以在PIMA数据上运行模型了。我们可以在模型上调用拟合函数f（），以在数据上训练或拟合模型。

我们先从KERAS中的程序开始，

神经网络一直训练到150个时期，并返回精确值。

Q5: OpenCV-Python系列四：图像分割(2)--梯度

上一期提到的图像阈值处理，不仅可以实现获取你想要的目标区域(作为mask使用)，还可以帮你获取图像的边缘信息，那关于图像边缘，本期将从另外的角度来处理。

对边缘信息与背景差异较大的场景，你也可以使用threshold分割，不过若阈值不好选取，Laplacian梯度算子就不失为一直尝试方案，而且上网看看，关于Laplacian算子还可以用来判断图像的模糊程度，这个在相机的自动对焦当中，是否可以尝试判断下？

不过处理的效果并不理想，图像低灰阶部分边缘信息丢失严重。

对于sobel，laplacian算子我们可以使用cv2.filter2D()来实现，配置相应的核模板即可，如实现提取水平方向边缘信息：

你可以依据实际的应用需求来配置提取边缘的角度信息，这里以45度角(垂直向下逆时针旋转45度)为例：

对此，你可以采用下面的方式来解决：