博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
python实现牛顿法求解求解最小值(包括拟牛顿法)【最优化课程笔记】
阅读量:6175 次
发布时间:2019-06-21

本文共 3111 字,大约阅读时间需要 10 分钟。

什么是牛顿法

在第9章中介绍了一维搜索的牛顿法,什么是一维搜索的牛顿法?

首先介绍一下一维搜索

一维搜索

一维搜索其实也很简单,在许多迭代下降算法中,具有一个共同的特点,就是得到点x(k)后,需要按照某种规则确定一个方向d(k),再从x(k)出发,沿着d(k)的方向上求目标函数的极小点。从而得到x(k+1),重复以上做法,知道求得问题的解。这就是一维搜索。上面提到的d可以称作为步长因子。

一维搜索的方法有很多,归纳起来可以大体分成两类。

  • 试探法 : 这种方法需要按照某种方式找试探点,通过一系列试探点来确定极小点。

  • 函数逼近法,或是插值法 : 用某种简单的曲线逼近本来的函数曲线,通过求逼近函数的极小点来估计目标函数的极小点。

里面的种种门道也是颇为复杂,这里就讲讲一维搜索的牛顿法,这是函数逼近法的一种。

一维搜索的牛顿法

牛顿法的基本思想是,在极小点附近用二阶Taylor多项式近似目标函数f(x),进而求出极小点的估计值。

假设,考虑问题

$$min f(x), x\in\Re.$$

$$\phi(x) = f(x^{(k)}) + f^{'}(x^{(k)})(x-x^{(k)}) + \frac{1}{2}f^{''}(x^{(k)})(x-x^{(k)})^2$$

又令

$$\phi^{'}(x) = f^{'}(x^{(k)}) + f^{''}(x^{(k)})(x-x^{(k)}) = 0$$

得到 $ phi(x) $ 的驻点,记作x(k+1),则

$$x^{(k+1)} = x^{(k)} - \frac{f^{'}(x^{(k)})}{f^{''}(x^{(k)})}$$

有了以上公式的推导,相信已经能够理解,下面为了加深对这个算法的印象,毕竟大家都是程序员,需要用伪代码才能说服对方。

牛顿法的计算步骤:

(1)给定初始点 x(0),允许误差 delta > 0 ,置k = 0(2)若 

$$ f^{'}(x^{(k)}) < delta $$

则迭代停止,,得到结果。(3)计算点

$$x^{(k+1)} = x^{(k)} - \frac{f^{'}(x^{(k)})}{f^{''}(x^{(k)})}$$

置 k:=k+1,转步骤(2)

运用牛顿法是,初始点选择十分重要,如果初始点靠近极小点,则可能很快收敛,如果初始点远离极小点,迭代产生的点列可能不收敛于极小点。

牛顿法

一般来讲的牛顿法是使用导数的最优化方法,这里的牛顿法和一维搜索的牛顿法差别在于更新公式的不同,这里的迭代公式如下。

$$x^{(k)} = x^{(k)} - \nabla^{2}f(x^{(k)})^{-1} \nabla f(x^{(k)})$$

计算的步骤也是一样的。

其实误差delta可以设置是1-范式,也可以是2-范式。

coding time

"""Newton法Rosenbrock函数函数 f(x)梯度 g(x)hessen 矩阵"""import numpy as npimport matplotlib.pyplot as plt# 一阶导def jacobian(x):    return np.array([2*x[0]+3,2*x[1]+4])# 二阶导def hessian(x):    return np.array([[2,0],[0,2]])X1=np.arange(-1.5,1.5+0.05,0.05)X2=np.arange(-3.5,2+0.05,0.05)[x1,x2]=np.meshgrid(X1,X2)f=x1**2+x2**2+3*x1+4*x2-26; # 给定的函数plt.contour(x1,x2,f,20) # 画出函数的20条轮廓线def newton(x0):    print('初始点为:')    print(x0,'\n')    W=np.zeros((2,10**3))    i = 1    imax = 1000    W[:,0] = x0     x = x0    delta = 1    while i
0.1: p = -np.dot(np.linalg.inv(hessian(x)),jacobian(x)) print(jacobian(x)) print(hessian(x)) x0 = x x = x + p W[:,i] = x delta = sum((x-x0)) print('第'+str(i)+'次迭代结果:') print(x,'\n') i=i+1 W=W[:,0:i] # 记录迭代点 return Wx0 = np.array([1,1])W=newton(x0)plt.plot(W[0,:],W[1,:],'g*',W[0,:],W[1,:]) # 画出迭代点收敛的轨迹plt.show()

代码注释的很清楚了,就不解释了,下面再提一句,如果更新时乘上了alpha值,就成了阻尼牛顿法。

"""Newton法Rosenbrock函数函数 f(x)梯度 g(x)hessen 矩阵"""import numpy as npimport matplotlib.pyplot as plt# 一阶导def jacobian(x):    return np.array([2*x[0]+3,2*x[1]+4])# 二阶导def hessian(x):    return np.array([[2,0],[0,2]])X1=np.arange(-1.5,1.5+0.05,0.05)X2=np.arange(-3.5,2+0.05,0.05)[x1,x2]=np.meshgrid(X1,X2)f=x1**2+x2**2+3*x1+4*x2-26; # 给定的函数plt.contour(x1,x2,f,20) # 画出函数的20条轮廓线def newton(x0):    print('初始点为:')    print(x0,'\n')    W=np.zeros((2,10**3))    i = 1    imax = 1000    W[:,0] = x0     x = x0    delta = 1    alpha = 0.1    while i
0.1: p = -np.dot(np.linalg.inv(hessian(x)),jacobian(x)) print(jacobian(x)) print(hessian(x)) x0 = x x = x + alpha*p W[:,i] = x delta = sum((x-x0)) print('第'+str(i)+'次迭代结果:') print(x,'\n') i=i+1 W=W[:,0:i] # 记录迭代点 return Wx0 = np.array([1,1])W=newton(x0)plt.plot(W[0,:],W[1,:],'g*',W[0,:],W[1,:]) # 画出迭代点收敛的轨迹plt.show()

转载地址:http://klqba.baihongyu.com/

你可能感兴趣的文章
图片自适应
查看>>
amd cmd
查看>>
Linux下的uml画图工具
查看>>
xml返回数组数据
查看>>
约瑟夫问题总结
查看>>
spring mybatis 批量插入返回主键
查看>>
指针函数小用
查看>>
开源力量公开课第二十三期-从SVN到Git,次时代代码管理
查看>>
输入挂
查看>>
升级迁移前,存储过程统计各个用户下表的数据量,和迁移后的比对
查看>>
sql注入分类
查看>>
初识CSS选择器版本4
查看>>
[Hadoop in China 2011] 朱会灿:探析腾讯Typhoon云计算平台
查看>>
JavaScript之数组学习
查看>>
PHP 设置响应头来解决跨域问题
查看>>
CAS实现SSO单点登录原理
查看>>
博客园美化专用图片链接
查看>>
HDU_1969_二分
查看>>
高等代数葵花宝典—白皮书
查看>>
一种简单的图像修复方法
查看>>