写在前面

SciPy的optimize模块提供了许多数值优化算法,下面对其中的一些记录。

非线性方程组求解

SciPy中对非线性方程组求解是fslove()函数,它的调用形式一般为fslove(fun, x0),fun是计算非线性方程组的误差函数,它需要一个参数x,fun依靠x来计算线性方程组的每个方程的值(或者叫误差),x0是x的一个初始值。

"""
计算非线性方程组:
  5x1+3 = 0
  4x0^2-2sin(x1x2)=0
  x1x2-1.5=0
"""
## 误差函数
def fun(x):
  x0,x1,x2 = x.tolist()
  return[5*x1+3,4x0^2-2sin(x1x2),x1x2-1.5]

result = optimize.fsolve(fun,[1,1,1])
## result
[-0.70622057  -0.6  -2.5]

在计算非线性方程中的解时,比如像坐标上升算法,其中需要用到未知数的导数,同样,scipy的fslove()也提供了fprime参数传递未知数的雅各比矩阵从而加速计算,传递的雅各比矩阵每一行时某一方程对各个未知数的导数。对于上面的例子,我们可以写下如下的雅各比矩阵传入。

def j(x):
  x0,x1,x2 = x.tolist()
  return[[0,5,0],[8*x0,-2*x2*cos(x1*x2],[0,x2,x1]]

result = optimize.fsolve(fun,[1,1,1],fprime=j)
#result
[-0.70622057  -0.6  -2.5]

scipy的内部在实现fslove时应该时应该是利用了坐标上升算法或者梯度相关优化算法,但本人没有考证,有兴趣的可以看看源码。

最小二乘拟合

关于最小二乘算法的理论这里并不想谈,网上解释的文章也挺多,在 optimize模块中,可以使用leastsq()对数据进行最小二乘拟合计算。 leastsq()的用法很简单,只需要将计箅误差的函数和待确定参数的初始值传递给它即可。

x = np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78])
y = np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])
def residual(p):
  k,b = p
  return y-(k*x+b)
r = optimize.leastsq(residual,[1,0])
k,b = r[0]
# print k
.613495349193
# print b
.79409254326
def func(x,p):
  """
    计算的正弦波 :A*sin(2*pi*k*x+theta)
  """
  A,k,theta = p
  return A*sin(2*np.pi*k*x+theta)

def redis(p,y,x):
  return y-func(x,p)

x = np.linspace(0,2*np.pi,100)
A,k,theta = 10,0.34,np.pi/6
y0 = func(x,[A,k,theta])
# 加入噪声
np.random.seed(0)
y1 = y0+2*np.random.randn(len(x))
p0 = [7,0.40,0]
# p0是A,k,theta的初始值,y1,x要拟合的数据
plsq = optimize.leastsq(redis, p0,args=(y1,x))
print [A,k,theta] #真是的参数值
print plsq[0]  #拟合后的参数值

对于像正弦波或者余弦波的曲线拟合,optimize提供curve_fit()函数,它的使用方式和leastq()稍有不同,它直接计算曲线的值,比如上面的拟合正弦波可以用cureve_fit()来写。

def func2(x,p):
  """
    计算的正弦波 :A*sin(2*pi*k*x+theta)
  """
  A,k,theta = p
  return A*sin(2*np.pi*k*x+theta)
ret,_=optimize.curve_fit(func2,x,y1,p0=p0)

该函数有一个缺点就是对于初始值敏感,如果初始频率和真实频率值差太多,会导致最后无法收敛到真是频率。

局部最小值

optimize模块还提供了常用的最小值算法如:Nelder-Mead、Powell、CG、BFGS、Newton-CG等,在这些最小值计算时,往往会传入一阶导数矩阵(雅各比矩阵)或者二阶导数矩阵(黑塞矩阵)从而加速收敛,这些最优化算法往往不能保证收敛到全局最小值,大部分会收敛到局部极小值。这些函数的调用方式为:

optimize.minimize(target_fun,init_val,method,jac,hess) 
target_fun:函数的表达式计算; 
init_val:初始值; 
method:最小化的算法; 
jac:雅各比矩阵 
hess:黑塞矩阵。

全局最小值算法

全局最小值使用optimize.basinhopping()来实现,这个函数首先要定义一个误差计算方式,比如平方误差函数,niter时迭代的次数,最后还需要一个局部极小值优化方法,minimizer_kwargs传入。比如上面的正弦函数拟合:

def func1(x,p):
  """
    计算的正弦波 :A*sin(2*pi*k*x+theta)
  """
  A,k,theta = p
  return A*sin(2*np.pi*k*x+theta)
def func_error(p,y,x):
  return np.sum((y-func1(x,p)**2)
result = optimize.basinhopping(func_error,[1,1,1],niter=10,
              minimizer_kwargs={"method":"L-BFGS-B",
                        "args":(y1,x1)})
## [1,1,1]是传入的初始值,args是需要拟合的数据

以上这篇python科学计算之scipy——optimize用法就是小编分享给大家的全部内容了,希望能给大家一个参考,也希望大家多多支持。

标签:
python,scipy,optimize

免责声明:本站文章均来自网站采集或用户投稿,网站不提供任何软件下载或自行开发的软件! 如有用户或公司发现本站内容信息存在侵权行为,请邮件告知! 858582#qq.com
白云城资源网 Copyright www.dyhadc.com

评论“python科学计算之scipy——optimize用法”

暂无“python科学计算之scipy——optimize用法”评论...

《魔兽世界》大逃杀!60人新游玩模式《强袭风暴》3月21日上线

暴雪近日发布了《魔兽世界》10.2.6 更新内容,新游玩模式《强袭风暴》即将于3月21 日在亚服上线,届时玩家将前往阿拉希高地展开一场 60 人大逃杀对战。

艾泽拉斯的冒险者已经征服了艾泽拉斯的大地及遥远的彼岸。他们在对抗世界上最致命的敌人时展现出过人的手腕,并且成功阻止终结宇宙等级的威胁。当他们在为即将于《魔兽世界》资料片《地心之战》中来袭的萨拉塔斯势力做战斗准备时,他们还需要在熟悉的阿拉希高地面对一个全新的敌人──那就是彼此。在《巨龙崛起》10.2.6 更新的《强袭风暴》中,玩家将会进入一个全新的海盗主题大逃杀式限时活动,其中包含极高的风险和史诗级的奖励。

《强袭风暴》不是普通的战场,作为一个独立于主游戏之外的活动,玩家可以用大逃杀的风格来体验《魔兽世界》,不分职业、不分装备(除了你在赛局中捡到的),光是技巧和战略的强弱之分就能决定出谁才是能坚持到最后的赢家。本次活动将会开放单人和双人模式,玩家在加入海盗主题的预赛大厅区域前,可以从强袭风暴角色画面新增好友。游玩游戏将可以累计名望轨迹,《巨龙崛起》和《魔兽世界:巫妖王之怒 经典版》的玩家都可以获得奖励。