使用Python實現(xiàn)牛頓法求極值
對于一個多元函數(shù) 用牛頓法求其極小值的迭代格式為
其中 為函數(shù)
的梯度向量,
為函數(shù)
的Hesse(Hessian)矩陣。
上述牛頓法不是全局收斂的。為此可以引入阻尼牛頓法(又稱帶步長的牛頓法)。
我們知道,求極值的一般迭代格式為
其中 為搜索步長,
為搜索方向(注意所有的迭代格式都是先計算搜索方向,再計算搜索步長,如同瞎子下山一樣,先找到哪個方向可行下降,再決定下幾步)。
取下降方向 即得阻尼牛頓法,只不過搜索步長
不確定,需要用線性搜索技術(shù)確定一個較優(yōu)的值,比如精確線性搜索或者Goldstein搜索、Wolfe搜索等。特別地,當(dāng)
一直取為常數(shù)1時,就是普通的牛頓法。
以Rosenbrock函數(shù)為例,即有
于是可得函數(shù)的梯度
函數(shù) 的Hesse矩陣為
編寫Python代碼如下(使用版本為Python3.3):
""" Newton法 Rosenbrock函數(shù) 函數(shù) f(x)=100*(x(2)-x(1).^2).^2+(1-x(1)).^2 梯度 g(x)=(-400*(x(2)-x(1)^2)*x(1)-2*(1-x(1)),200*(x(2)-x(1)^2))^(T) """ import numpy as np import matplotlib.pyplot as plt def jacobian(x): return np.array([-400*x[0]*(x[1]-x[0]**2)-2*(1-x[0]),200*(x[1]-x[0]**2)]) def hessian(x): return np.array([[-400*(x[1]-3*x[0]**2)+2,-400*x[0]],[-400*x[0],200]]) 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=100*(x2-x1**2)**2+(1-x1)**2; # 給定的函數(shù) plt.contour(x1,x2,f,20) # 畫出函數(shù)的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 = 1 while i<imax and delta>10**(-5): p = -np.dot(np.linalg.inv(hessian(x)),jacobian(x)) x0 = x x = x + alpha*p W[:,i] = x delta = sum((x-x0)**2) print('第',i,'次迭代結(jié)果:') print(x,'\n') i=i+1 W=W[:,0:i] # 記錄迭代點 return W x0 = np.array([-1.2,1]) W=newton(x0) plt.plot(W[0,:],W[1,:],'g*',W[0,:],W[1,:]) # 畫出迭代點收斂的軌跡 plt.show()
上述代碼中jacobian(x)返回函數(shù)的梯度,hessian(x)返回函數(shù)的Hesse矩陣,用W矩陣記錄迭代點的坐標(biāo),然后畫出點的搜索軌跡。
可得輸出結(jié)果為
初始點為: [-1.2 1. ] 第 1 次迭代結(jié)果: [-1.1752809 1.38067416] 第 2 次迭代結(jié)果: [ 0.76311487 -3.17503385] 第 3 次迭代結(jié)果: [ 0.76342968 0.58282478] 第 4 次迭代結(jié)果: [ 0.99999531 0.94402732] 第 5 次迭代結(jié)果: [ 0.9999957 0.99999139] 第 6 次迭代結(jié)果: [ 1. 1.]
即迭代了6次得到了最優(yōu)解,畫出的迭代點的軌跡如下:
由于主要使用了Python的Numpy模塊來進(jìn)行計算,可以看出,代碼和最終的圖與Matlab是很相像的。
以上這篇使用Python實現(xiàn)牛頓法求極值就是小編分享給大家的全部內(nèi)容了,希望能給大家一個參考,也希望大家多多支持腳本之家。
相關(guān)文章
python多進(jìn)程下實現(xiàn)日志記錄按時間分割
這篇文章主要為大家詳細(xì)介紹了python多進(jìn)程下實現(xiàn)日志記錄按時間分割,具有一定的參考價值,感興趣的小伙伴們可以參考一下2019-07-07使用Django2快速開發(fā)Web項目的詳細(xì)步驟
這篇文章主要介紹了使用Django2快速開發(fā)Web項目的詳細(xì)步驟,文中通過示例代碼介紹的非常詳細(xì),對大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價值,需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧2019-01-01Python操作Redis數(shù)據(jù)庫的超詳細(xì)教程
大家應(yīng)該都知道redis是一個基于內(nèi)存的高效的鍵值型非關(guān)系數(shù)據(jù)庫,下面這篇文章主要給大家介紹了關(guān)于Python操作Redis的相關(guān)資料,文中通過實例代碼介紹的非常詳細(xì),需要的朋友可以參考下2022-06-06python使用pandas庫導(dǎo)入并保存excel、csv格式文件數(shù)據(jù)
CSV格式文件很方便各種工具之間傳遞數(shù)據(jù),平時工作過程之中會將數(shù)據(jù)保存為CSV格式,這篇文章主要介紹了python使用pandas庫導(dǎo)入并保存excel、csv格式文件數(shù)據(jù)的相關(guān)資料,需要的朋友可以參考下2017-12-12