下面定義了三個Python語言編寫的函數：函數表達式fun，梯度向量gfun，和海森矩陣hess。這三個表達式在后面各個算法的實現中會用到。

# 函數表達式fun fun = lambda x:100*(x[0]**2-x[1])**2 + (x[0]-1)**2# 梯度向量 gfun gfun = lambda x:np.array([400*x[0]*(x[0]**2-x[1])+2*(x[0]-1), -200*(x[0]**2-x[1])])# 海森矩陣 hess hess = lambda x:np.array([[1200*x[0]**2-400*x[1]+2, -400*x[0]],[-400*x[0],200]])

最速下降法的Python實現。

def gradient(fun,gfun,x0):#最速下降法求解無約束問題# x0是初始點，fun和gfun分別是目標函數和梯度maxk = 5000rho = 0.5sigma = 0.4k = 0epsilon = 1e-5while k<maxk:gk = gfun(x0)dk = -gkif np.linalg.norm(dk) < epsilon:breakm = 0mk = 0while m< 20:if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):mk = mbreakm += 1x0 += rho**mk*dkk += 1return x0,fun(x0),k

4.牛頓法
牛頓法的基本思想是用迭代點xkxk處的一階導數（梯度gkgk）和二階倒數（海森矩陣GkGk）對目標函數進行二次函數近似，然后把二次模型的極小點作為新的迭代點。

牛頓法推導

函數f(x)f(x)在xkxk處的泰勒展開式的前三項得
T(f,xk,3)=fk+gTk(x?xk)+12(x?xk)TGk(x?xk)
T(f,xk,3)=fk+gkT(x?xk)+12(x?xk)TGk(x?xk)

則其穩定點
?T=gk+Gk(x?xk)=0
?T=gk+Gk(x?xk)=0

若 GkGk非奇異，那么解上面的線性方程（標記其解為 xk+1xk+1）即得到牛頓法的迭代公式
xk+1=xk?G?1kgk
xk+1=xk?Gk?1gk

即優化方向為 dk=?G?1kgkdk=?Gk?1gk
求穩定點是用到了以下公式：

考慮y=xTAxy=xTAx，根據矩陣求導法則，有
dydx=Ax+ATxd2ydx2=A+AT
dydx=Ax+ATxd2ydx2=A+AT
注意實際中dkdk的是通過求解線性方程組Gkd=?gkGkd=?gk獲得的。

阻尼牛頓法及其Python實現
阻尼牛頓法采用了牛頓法的思想，唯一的差別是阻尼牛頓法并不直接采用xk+1=xk+dkxk+1=xk+dk，而是用線搜索技術獲得合適的步長因子αkαk,用公式xk+1=xk+δmdkxk+1=xk+δmdk計算xk+1xk+1。

阻尼牛頓法的算法描述：

step 1: 給定終止誤差0≤??1,δ∈(0,1),σ∈(0,0.5)0≤??1,δ∈(0,1),σ∈(0,0.5). 初始點x0∈Rnx0∈Rn. 令k←0k←0
step 2: 計算gk=?f(xk)gk=?f(xk). 若||gk||≤?||gk||≤?，停止迭代，輸出x?≈xkx?≈xk
step 3: 計算Gk=?2f(xk)Gk=?2f(xk)，并求解線性方程組得到解dkdk,
Gkd=?gk
Gkd=?gk

step 4: 記 mkmk是滿足下列不等式的最小非負整數 mm.
f(xk+βmdk)≤f(xk)+δβmgTkdk
f(xk+βmdk)≤f(xk)+δβmgkTdk

step 5: 令 αk=δmk,xk+1=xk+αkdk,k←k+1αk=δmk,xk+1=xk+αkdk,k←k+1，轉 step 2
阻尼牛頓法的Python實現：

def dampnm(fun,gfun,hess,x0):# 用阻尼牛頓法求解無約束問題#x0是初始點，fun，gfun和hess分別是目標函數值，梯度，海森矩陣的函數maxk = 500rho = 0.55sigma = 0.4k = 0epsilon = 1e-5while k < maxk:gk = gfun(x0)Gk = hess(x0)dk = -1.0*np.linalg.solve(Gk,gk)if np.linalg.norm(dk) < epsilon:breakm = 0mk = 0while m < 20:if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):mk = mbreakm += 1x0 += rho**mk*dkk += 1return x0,fun(x0),k

性能展示

作圖代碼：

iters = [] for i in xrange(np.shape(data)[0]):rt = dampnm(fun,gfun,hess,data[i])if rt[2] <=200:iters.append(rt[2])if i%100 == 0:print i,rt[2]plt.hist(iters,bins=50) plt.title(u'阻尼牛頓法迭代次數分布',{'fontname':'STFangsong','fontsize':18}) plt.xlabel(u'迭代次數',{'fontname':'STFangsong','fontsize':18}) plt.ylabel(u'頻率分布',{'fontname':'STFangsong','fontsize':18})

修正牛頓法法及其Python實現
牛頓法要求目標函數的海森矩陣G(x)=?2f(x)G(x)=?2f(x)在每個迭代點xkxk處是正定的，否則難以保證牛頓方向dk=?G?1kgkdk=?Gk?1gk是ff在xkxk處的下降方向。為克服這一缺陷，需要對牛頓法進行修正。
下面介紹兩種修正方法，分別是“牛頓-最速下降混合算法”和“修正牛頓法”。
牛頓-最速下降混合算法
該方法將牛頓法和最速下降法結合起來，基本思想是：當海森矩陣正定時，采用牛頓法確定的優化方向作為搜索方向；否則，即海森矩陣為非正定矩陣時，則采用負梯度方向作為搜索方向。

修正牛頓法
引入阻尼因子μk≥0μk≥0，即在每一迭代步適當選取參數μkμk，使得矩陣Ak=Gk+μkIAk=Gk+μkI正定，用AkAk代替GkGk來求解dkdk。

算法描述：

step 1: 給定終止誤差0≤??1,δ∈(0,1),σ∈(0,0.5)0≤??1,δ∈(0,1),σ∈(0,0.5). 初始點x0∈Rnx0∈Rn. 令k←0k←0
step 2: 計算gk=?f(xk)gk=?f(xk)，μk=||gk||1+τμk=||gk||1+τ. 若||gk||≤?||gk||≤?，停止迭代，輸出x?≈xkx?≈xk
step 3: 計算Gk=?2f(xk)Gk=?2f(xk)，并求解線性方程組得到解dkdk,
(Gk+μkI)d=?gk
(Gk+μkI)d=?gk

step 4: 記 mkmk是滿足下列不等式的最小非負整數 mm.
f(xk+βmdk)≤f(xk)+δβmgTkdk
f(xk+βmdk)≤f(xk)+δβmgkTdk

step 5: 令 αk=δmk,xk+1=xk+αkdk,k←k+1αk=δmk,xk+1=xk+αkdk,k←k+1，轉 step 2
修正牛頓法的Python實現代碼：

def revisenm(fun,gfun,hess,x0):# 用修正牛頓法求解無約束問題#x0是初始點，fun，gfun和hess分別是目標函數值，梯度，海森矩陣的函數maxk = 1e5n = np.shape(x0)[0]rho = 0.55sigma = 0.4tau = 0.0k = 0epsilon = 1e-5while k < maxk:gk = gfun(x0) ? ? ?if ?np.linalg.norm(gk) < epsilon:breakmuk = np.power(np.linalg.norm(gk),1+tau)Gk = hess(x0)Ak = Gk + muk*np.eye(n)dk = -1.0*np.linalg.solve(Ak,gk)m = 0mk = 0while m < 20:if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):mk = mbreakm += 1x0 += rho**mk*dkk += 1return x0,fun(x0),k

5.擬牛頓法
牛頓法的優點是具有二階收斂速度，缺點是：

但當海森矩陣G(xk)=?2f(x)G(xk)=?2f(x) 不正定時，不能保證所產生的方向是目標函數在xkxk處的下降方向。
特別地，當G(xk)G(xk)奇異時，算法就無法繼續進行下去。盡管修正牛頓法可以克服這一缺陷，但修正參數的取值很難把握，過大或過小都會影響到收斂速度。
牛頓法的每一步迭代都需要目標函數的海森矩陣G(xk)G(xk)，對于大規模問題其計算量是驚人的。
擬牛頓法的基本思想是用海森矩陣GkGk的某個近似矩陣BkBk取代GkGk. BkBk通常具有下面三個特點：

在某種意義下有Bk≈GkBk≈Gk ,使得相應的算法產生的方向近似于牛頓方向，確保算法具有較快的收斂速度。
對所有的kk，BkBk是正定的，從而使得算法所產生的方向是函數ff在xkxk處下降方向。
矩陣BkBk更新規則比較簡單
設 f:Rn→Rf:Rn→R在開集D?RnD?Rn上二次連續可微，那么ff在xk+1xk+1處的二次近似模型為：
f(x)≈f(xk+1)+gTk+1(x?xk+1)+12(x?xk+1)TGk+1(x?xk+1)
f(x)≈f(xk+1)+gk+1T(x?xk+1)+12(x?xk+1)TGk+1(x?xk+1)
對上式求導得
g(x)≈gk+1+Gk+1(x?xk+1)
g(x)≈gk+1+Gk+1(x?xk+1)

令 x=xkx=xk，位移 sk=xk+1?xksk=xk+1?xk，梯度差 yk=gk+1?gkyk=gk+1?gk，則有
Gk+1sk≈yk
Gk+1sk≈yk
.
因此，擬牛頓法中近似矩陣BkBk要滿足關系式
Bk+1sk=yk
Bk+1sk=yk

令 Hk+1=B?1k+1Hk+1=Bk+1?1，得到擬牛頓法的另一形式
Hk+1yk=sk
Hk+1yk=sk

其中 Hk+1Hk+1為海森矩陣逆矩陣的近似。上面兩個公式稱為 擬牛頓方程。
搜索方向由 dk=?Hkgkdk=?Hkgk或 Bkdk=?gkBkdk=?gk確定。根據近似矩陣的第三個特點，有
Bk+1=Bk+Ek,Hk+1=Hk+Dk
Bk+1=Bk+Ek,Hk+1=Hk+Dk

其中 Ek,DkEk,Dk為秩1或秩2矩陣。該公式稱為 校正規則。
通常將上面的三個式子（擬牛頓方程和校正規則）所確立的方法稱為擬牛頓法。從下面的擬牛頓法的幾個變種可以看出，不同的擬牛頓法的主要區別在于更新公式的不同。
DFP算法及其Python實現
DFP算法的校正公式如下：
Hk+1=Hk?HkykyTkHkyTkHkyk+sksTksTkyk
Hk+1=Hk?HkykykTHkykTHkyk+skskTskTyk
注意公式中的兩個分數結構，分母yTkHkykykTHkyk和sTkykskTyk是標量，分子HkykyTkHkHkykykTHk和sksTkskskT是與HkHk同型的矩陣，而且都是對稱矩陣。若HkHk正定，且sTkyk>0skTyk>0，則Hk+1Hk+1也正定。

當采用精確線搜索時，矩陣序列HkHk的正定新條件sTkyk>0skTyk>0可以被滿足。但對于ArmijoArmijo搜索準則來說，不能滿足這一條件，需要做如下修正：
Hk+1=???HkHk?HkykyTkHkyTkHkyk+sksTksTkyksTkyk≤0sTkyk>0???
Hk+1={HkskTyk≤0Hk?HkykykTHkykTHkyk+skskTskTykskTyk>0}
DFP算法描述：

step 1: 給定參數δ∈(0,1),σ∈(0,0.5)δ∈(0,1),σ∈(0,0.5)，初始點x0∈Rx0∈R,終止誤差0≤??10≤??1.初始對稱正定矩陣H0H0(通常取為G(x0)?1G(x0)?1或單位矩陣InIn).令k←0k←0
step 2: 計算gk=?f(xk)gk=?f(xk). 若||gk||≤?||gk||≤?，停止迭代，輸出x?≈xkx?≈xk
step 3: 計算搜索方向
dk=?Hkgk
dk=?Hkgk

step 4: 記 mkmk是滿足下列不等式的最小非負整數 mm.
f(xk+βmdk)≤f(xk)+δβmgTkdk
f(xk+βmdk)≤f(xk)+δβmgkTdk
令 αk=δmk,xk+1=xk+αkdkαk=δmk,xk+1=xk+αkdk
step 5: 由校正公式確定 Hk+1Hk+1
step 6： k←k+1k←k+1,轉 step 2
DFP算法的代碼實現

def dfp(fun,gfun,hess,x0):#功能：用DFP族算法求解無約束問題：min fun(x)#輸入：x0是初始點，fun,gfun分別是目標函數和梯度#輸出：x,val分別是近似最優點和最優解,k是迭代次數maxk = 1e5rho = 0.55sigma = 0.4epsilon = 1e-5k = 0n = np.shape(x0)[0]#海森矩陣可以初始化為單位矩陣Hk = np.linalg.inv(hess(x0)) #或者單位矩陣np.eye(n)while k < maxk :gk = gfun(x0)if np.linalg.norm(gk) < epsilon:breakdk = -1.0*np.dot(Hk,gk)m=0;mk=0while m < 20: # 用Armijo搜索求步長if fun(x0+rho**m*dk) < fun(x0)+sigma*rho**m*np.dot(gk,dk):mk = mbreakm += 1#print mk#DFP校正x = x0 + rho**mk*dksk = x - x0yk = gfun(x) - gk ??if np.dot(sk,yk) > 0:Hy = np.dot(Hk,yk)print Hysy = np.dot(sk,yk) #向量的點積yHy = np.dot(np.dot(yk,Hk),yk) # yHy是標量#表達式Hy.reshape((n,1))*Hy 中Hy是向量，生成矩陣Hk = Hk - 1.0*Hy.reshape((n,1))*Hy/yHy + 1.0*sk.reshape((n,1))*sk/syk += 1x0 = xreturn x0,fun(x0),k#分別是最優點坐標，最優值，迭代次數

由以上代碼可知，擬牛頓法只需要在迭代開始前計算一次海森矩陣，更一般的，根本就不計算海森矩陣，而是初始化為單位矩陣，在每次迭代過程中是不需按傳統方法（二階導數）計算海森矩陣的。

實際性能

相關代碼：

n = 50 x = y = np.linspace(-10,10,n) xx,yy = np.meshgrid(x,y) data = [[xx[i][j],yy[i][j]] for i in range(n) for j in range(n)] iters = [] for i in xrange(np.shape(data)[0]):rt = dfp(fun,gfun,hess,np.array(data[i]))if rt[2] <=200: # 提出迭代次數過大的iters.append(rt[2])if i%100 == 0:print i,rt[2]plt.hist(iters,bins=50) plt.title(u'DFP迭代次數分布',{'fontname':'STFangsong','fontsize':18}) plt.xlabel(u'迭代次數',{'fontname':'STFangsong','fontsize':18}) plt.ylabel(u'頻率分布',{'fontname':'STFangsong','fontsize':18})

BFGS算法及其Python實現
BFGS算法與DFP步驟基本相同，區別在于更新公式的差異
Bk+1=???BkBk?BkykyTkBkyTkBkyk+sksTksTkyksTkyk≤0sTkyk>0???
Bk+1={BkskTyk≤0Bk?BkykykTBkykTBkyk+skskTskTykskTyk>0}
BFGS算法的Python實現

def bfgs(fun,gfun,hess,x0):#功能：用BFGS族算法求解無約束問題：min fun(x)#輸入：x0是初始點，fun,gfun分別是目標函數和梯度#輸出：x,val分別是近似最優點和最優解,k是迭代次數 ?maxk = 1e5rho = 0.55sigma = 0.4epsilon = 1e-5k = 0n = np.shape(x0)[0]#海森矩陣可以初始化為單位矩陣Bk = eye(n)#np.linalg.inv(hess(x0)) #或者單位矩陣np.eye(n)while k < maxk:gk = gfun(x0)if np.linalg.norm(gk) < epsilon:breakdk = -1.0*np.linalg.solve(Bk,gk)m = 0mk = 0while m < 20: # 用Armijo搜索求步長if fun(x0+rho**m*dk) < fun(x0)+sigma*rho**m*np.dot(gk,dk):mk = mbreakm += 1#BFGS校正x = x0 + rho**mk*dksk = x - x0yk = gfun(x) - gk ??if np.dot(sk,yk) > 0: ? ?Bs = np.dot(Bk,sk)ys = np.dot(yk,sk)sBs = np.dot(np.dot(sk,Bk),sk)?Bk = Bk - 1.0*Bs.reshape((n,1))*Bs/sBs + 1.0*yk.reshape((n,1))*yk/ysk += 1x0 = xreturn x0,fun(x0),k#分別是最優點坐標，最優值，迭代次數?

實際性能

相關代碼：

iters = [] for i in xrange(np.shape(data)[0]):rt = bfgs(fun,gfun,hess,np.array(data[i]))if rt[2] <=200:iters.append(rt[2])if i%100 == 0:print i,rt[2]plt.hist(iters,bins=50) plt.title(u'BFGS迭代次數分布',{'fontname':'STFangsong','fontsize':18}) plt.xlabel(u'迭代次數',{'fontname':'STFangsong','fontsize':18}) plt.ylabel(u'頻率分布',{'fontname':'STFangsong','fontsize':18})

Broyden族算法及其Python實現
之前的BFGS和DFP校正都是由ykyk和BkskBksk(或 sksk和HkykHkyk)組成的秩2矩陣。而Broyden族算法采用了BFGS和DFP校正公式的凸組合，如下：
H?k+1=?kHBFGSk+1+(1??k)HDFPk+1=Hk?HkykyTkHkyTkHkyk+sksTksTkyk+?kvkvTk
Hk+1?=?kHk+1BFGS+(1??k)Hk+1DFP=Hk?HkykykTHkykTHkyk+skskTskTyk+?kvkvkT

其中 ?k∈[0,1]?k∈[0,1]， vkvk由下式定義：
vk=yTkHkyk??????√(skyTksk?HkykyTkHkyk)
vk=ykTHkyk(skykTsk?HkykykTHkyk)
Broyden族算法Python實現

def broyden(fun,gfun,hess,x0):#功能：用Broyden族算法求解無約束問題：min fun(x)#輸入：x0是初始點，fun,gfun分別是目標函數和梯度#輸出：x,val分別是近似最優點和最優解,k是迭代次數x0 = np.array(x0)maxk = 1e5rho = 0.55;sigma = 0.4;epsilon = 1e-5phi = 0.5;k=0;n = np.shape(x0)[0]Hk = np.linalg.inv(hess(x0))while k<maxk :gk = gfun(x0)if np.linalg.norm(gk) < epsilon:breakdk = -1*np.dot(Hk,gk)m=0;mk=0while m < 20: # 用Armijo搜索求步長if fun(x0+rho**m*dk) < fun(x0)+sigma*rho**m*np.dot(gk,dk):mk = mbreakm += 1#Broyden族校正x = x0 + rho**mk*dksk = x - x0yk = gfun(x) - gkHy = np.dot(Hk,yk)sy = np.dot(sk,yk)yHy = np.dot(np.dot(yk,Hk),yk)if(sy < 0.2 *yHy):theta = 0.8*yHy/(yHy-sy)sk = theta*sk + (1-theta)*Hysy = 0.2*yHyvk = np.sqrt(yHy)*(sk/sy-Hy/yHy)Hk = Hk - Hy.reshape((n,1))*Hy/yHy +sk.reshape((n,1))*sk/sy + phi*vk.reshape((n,1))*vkk += 1x0 = xreturn x0,fun(x0),k #分別是最優點坐標，最優值，迭代次數

實際性能

相關代碼：

n = 50 x = y = np.linspace(-10,10,n) xx,yy = np.meshgrid(x,y) data = [[xx[i][j],yy[i][j]] for i in range(n) for j in range(n)]iters = [] for i in xrange(np.shape(data)[0]):rt = lbfgs(fun,gfun,data[i])if rt[2] <=1000:iters.append(rt[2])if i%100 == 0:print i,rt[2]plt.hist(iters,bins=100) plt.title(u'L-BFGS法迭代次數分布',{'fontname':'STFangsong','fontsize':18}) plt.xlabel(u'迭代次數',{'fontname':'STFangsong','fontsize':18}) plt.ylabel(u'頻率分布',{'fontname':'STFangsong','fontsize':18})

L-BFGS算法及其Python實現
擬牛頓法(如BFGS算法)需要計算和存儲海森矩陣，其空間復雜度是n2n2，當nn很大時，其需要的內存量是非常大的。為了解決該內存問題，有限內存BFGS（即傳說中的L-BFGS算法）橫空出世。

基本BFGS算法的校正公式可以改寫成
Hk+1=(I?skyTksTkyk)Hk(I?yksTksTkyk)+sksTksTkyk
Hk+1=(I?skykTskTyk)Hk(I?ykskTskTyk)+skskTskTyk

記ρk=1/sTkykρk=1/skTyk，以及Vk=(I?ρkyksTk)Vk=(I?ρkykskT)，則上式可以寫成
Hk+1=VTkHkVk+ρksksTk
Hk+1=VkTHkVk+ρkskskT
給定一個初始矩陣H0H0(其取值后面有提到)，則由上式的遞推關系可以得到

H1=VT0H0V0+ρ0s0sT0
H1=V0TH0V0+ρ0s0s0T
H2=VT1H1V1+ρ1s1sT1=VT1(VT0H0V0+ρ0s0sT0)V1+ρ1s1sT1=VT1VT0H0V0V1+VT1ρ0s0sT0V1+ρ1s1sT1
H2=V1TH1V1+ρ1s1s1T=V1T(V0TH0V0+ρ0s0s0T)V1+ρ1s1s1T=V1TV0TH0V0V1+V1Tρ0s0s0TV1+ρ1s1s1T
???
???

更一般的，我們有
Hm+1=(VTKVTm?1?VT1VT0)H0(V0V1?Vm?1Vk)+(VTmVTm?1?VT2VT1)ρ0s0sT0(V1V2?Vm?1Vk)+(VTmVTm?1?VT3VT2)ρ1s1sT1(V2V3?Vm?1Vk)+?+(VTmVTm?1)ρm?2sm?2sTm?2(Vm?1Vk)VTkρm?1sm?1sTm?1Vk+ρmsmsTm
Hm+1=(VKTVm?1T?V1TV0T)H0(V0V1?Vm?1Vk)+(VmTVm?1T?V2TV1T)ρ0s0s0T(V1V2?Vm?1Vk)+(VmTVm?1T?V3TV2T)ρ1s1s1T(V2V3?Vm?1Vk)+?+(VmTVm?1T)ρm?2sm?2sm?2T(Vm?1Vk)VkTρm?1sm?1sm?1TVk+ρmsmsmT

在上式的右邊中， H0H0是由我們設置的，其余變量可由保存的最近 mm次迭代所形成的向量序列，如位移向量序列{sk}{sk}和一階導數差所形成的向量序列 {yk}{yk}獲得。也就是說，可由最近 mm次迭代的信息計算當前的海森矩陣的近似矩陣。
補充幾點：

H0=I?sTmymyTmymH0=I?smTymymTym，第一次迭代時，序列{ sksk}和{ ykyk}為空，則 H0=IH0=I
最初的若干次迭代， 序列{sksk}和{ykyk}的元素個數較小，會有n<mn<m，則將Hm+1Hm+1計算公式右邊的mm換成nn即可。
隨著迭代次數增加， 序列{sksk}和{ykyk}的元素個數增大，會有n>mn>m。由于我們只需要最近mm次迭代產生的序列元素，所以序列{sksk}和{ykyk}只需要保存最新的mm個元素即可，如果再有新的元素進入，則同時扔掉最舊的元素，以保證序列元素個數恒為mm。
這就是L-BFGS算法的思想。聰明的同學會問：你上面的公式不還是要計算海森矩陣的近似矩陣嗎？那豈不是還是需要n2n2規模的內存？
其實在實際中是不計算該矩陣的, 而且不是采用上面的方法，而是直接得到HkgkHkgk，而搜索方向就是dk=?Hkgkdk=?Hkgk。下面給出了計算的函數twoloop(s,y,ρρ,gk)的偽代碼，可知該函數內部的計算僅限于標量和向量，未出現矩陣。

函數參數s,ys,y即向量序列{sksk}和{ykyk}，ρρ為元素序列{ρkρk}，其元素ρk=1/sTkykρk=1/skTyk，gkgk是向量，為當前的一階導數，輸出為向量z=Hkgkz=Hkgk，即搜索方向的反方向
Functiontwoloop(s,y,ρ,gk)q=gkFori=k?1,k?2,…,k?mαi=ρisTiqq=q?αiyiHk=yTk?1sk?1/yTk?1yk?1z=HkqFori=k?m,k?m+1,…,k?1βi=ρiyTizz=z+si(αi?βi)Returnz
Functiontwoloop(s,y,ρ,gk)q=gkFori=k?1,k?2,…,k?mαi=ρisiTqq=q?αiyiHk=yk?1Tsk?1/yk?1Tyk?1z=HkqFori=k?m,k?m+1,…,k?1βi=ρiyiTzz=z+si(αi?βi)Returnz
給出L-BFGS的算法

可以看到其搜索方向dkdk是根據函數twolooptwoloop計算得到的。

L-BFGS算法Python實現

def twoloop(s, y, rho,gk):n = len(s) #向量序列的長度if np.shape(s)[0] >= 1:#h0是標量，而非矩陣h0 = 1.0*np.dot(s[-1],y[-1])/np.dot(y[-1],y[-1])else:h0 = 1a = empty((n,))q = gk.copy()?for i in range(n - 1, -1, -1):?a[i] = rho[i] * dot(s[i], q)q -= a[i] * y[i]z = h0*qfor i in range(n):b = rho[i] * dot(y[i], z)z += s[i] * (a[i] - b)return z ??def lbfgs(fun,gfun,x0,m=5):# fun和gfun分別是目標函數及其一階導數,x0是初值,m為儲存的序列的大小maxk = 2000rou = 0.55sigma = 0.4epsilon = 1e-5k = 0n = np.shape(x0)[0] #自變量的維度s, y, rho = [], [], []while k < maxk :gk = gfun(x0)if np.linalg.norm(gk) < epsilon:breakdk = -1.0*twoloop(s, y, rho,gk)m0=0;mk=0while m0 < 20: # 用Armijo搜索求步長if fun(x0+rou**m0*dk) < fun(x0)+sigma*rou**m0*np.dot(gk,dk):mk = m0breakm0 += 1x = x0 + rou**mk*dksk = x - x0yk = gfun(x) - gk ??if np.dot(sk,yk) > 0: #增加新的向量rho.append(1.0/np.dot(sk,yk))s.append(sk)y.append(yk)if np.shape(rho)[0] > m: #棄掉最舊向量rho.pop(0)s.pop(0)y.pop(0)k += 1x0 = xreturn x0,fun(x0),k#分別是最優點坐標，最優值，迭代次數

相關代碼：

n = 50 x = y = np.linspace(-10,10,n) xx,yy = np.meshgrid(x,y) data = [[xx[i][j],yy[i][j]] for i in range(n) for j in range(n)]iters = [] for i in xrange(np.shape(data)[0]):rt = lbfgs(fun,gfun,data[i])if rt[2] <=1000:iters.append(rt[2])if i%100 == 0:print i,rt[2]plt.hist(iters,bins=100) plt.title(u'L-BFGS法迭代次數分布',{'fontname':'STFangsong','fontsize':18}) plt.xlabel(u'迭代次數',{'fontname':'STFangsong','fontsize':18}) plt.ylabel(u'頻率分布',{'fontname':'STFangsong','fontsize':18})

如果對代碼有疑問，歡迎添加微信 ： 18565453898

創作挑戰賽新人創作獎勵來咯，堅持創作打卡瓜分現金大獎

總結

以上是生活随笔為你收集整理的L-BFGS算法/Broyden族/BFGS算法/阻尼牛顿法的Python实现代码的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。