與Python實現(xiàn))
1. 為什么擬牛頓算法值得你花時間搞懂如果你做過一點數(shù)值優(yōu)化或者機器學(xué)習(xí)模型的訓(xùn)練大概率聽過“牛頓法”這三個字。牛頓法本身很優(yōu)雅用二階導(dǎo)數(shù)信息Hessian矩陣來指導(dǎo)迭代方向收斂速度是二次的理論上非常漂亮。但問題也很直接Hessian矩陣的計算代價太高尤其是變量維度一上去求逆矩陣的復(fù)雜度直接爆炸。所以實際工程里純牛頓法基本只存在于教科書和小規(guī)模問題中。擬牛頓算法就是在這個背景下被提出來的。它的核心思路很樸素我不去精確算Hessian而是用迭代過程中產(chǎn)生的梯度差和步長差去近似Hessian或者它的逆矩陣。每走一步就用新得到的信息去更新這個近似矩陣。這樣既保留了牛頓法“利用曲率信息”的優(yōu)勢又把計算量壓到了可以接受的范圍。SR1、DFP、BFGS是擬牛頓家族里最經(jīng)典的三個更新公式。SR1是Symmetric Rank-One的縮寫顧名思義它用秩一矩陣去修正近似HessianDFP是Davidon-Fletcher-Powell三個人名字的首字母直接更新Hessian的逆矩陣BFGS則是Broyden-Fletcher-Goldfarb-Shanno是目前公認(rèn)綜合表現(xiàn)最穩(wěn)、應(yīng)用最廣的擬牛頓更新公式。你在scipy、sklearn、各種深度學(xué)習(xí)優(yōu)化器里看到的“L-BFGS”“BFGS”根都在這里。這篇文章適合誰看如果你正在學(xué)數(shù)值優(yōu)化、準(zhǔn)備面試被問到優(yōu)化算法、或者自己寫代碼時想手?jǐn)]一個擬牛頓優(yōu)化器那這篇內(nèi)容可以直接當(dāng)參考。我會把三個公式的推導(dǎo)邏輯、代碼實現(xiàn)、參數(shù)選擇、踩坑經(jīng)驗全部攤開講不跳步不堆公式嚇人盡量讓你看完能自己寫出來、跑起來、調(diào)得動。2. 擬牛頓算法的整體設(shè)計思路拆解2.1 牛頓法到底卡在哪里先把牛頓法的迭代公式擺出來x_{k1} x_k - H_k^{-1} * g_k其中g(shù)_k是梯度H_k是Hessian矩陣。這個公式的意思是在當(dāng)前點用二階泰勒展開去近似目標(biāo)函數(shù)然后直接跳到這個二次模型的極小值點。問題出在H_k^{-1}上。假設(shè)你有n個變量Hessian是n×n矩陣。每次迭代都要重新計算Hessian需要n2個二階偏導(dǎo)然后求逆復(fù)雜度O(n3)。n1000的時候這個計算量已經(jīng)很難接受了n10000的時候基本不用想。更麻煩的是Hessian還不一定正定。如果Hessian不是正定的牛頓方向可能指向鞍點甚至極大值方向迭代直接發(fā)散。所以純牛頓法在實際中很少直接用除非問題規(guī)模很小且性質(zhì)很好。2.2 擬牛頓的核心思想用梯度差代替二階導(dǎo)擬牛頓法的切入點很巧妙。它不去算Hessian而是觀察到一個關(guān)系H_k * (x_{k1} - x_k) ≈ g_{k1} - g_k這個式子來自泰勒展開梯度在兩點之間的變化近似等于Hessian乘以位移。令s_k x_{k1} - x_ky_k g_{k1} - g_k就得到擬牛頓方程B_k * s_k y_k其中B_k是Hessian的近似矩陣?；蛘叻催^來用逆矩陣表示H_k * y_k s_k其中H_k是Hessian逆矩陣的近似。擬牛頓法的每一步就是在滿足這個方程的前提下用一個低秩矩陣去修正上一步的近似矩陣。不同的修正方式就產(chǎn)生了SR1、DFP、BFGS等不同公式。2.3 三個公式的定位差異這三個公式不是隨便列的它們各自有明確的適用場景和數(shù)學(xué)性質(zhì)SR1秩一更新形式最簡單不保證正定性但近似精度高適合信任域方法。DFP秩二更新保證正定性在curvature condition滿足時直接更新逆Hessian是BFGS的前身。BFGS秩二更新同樣保證正定性但更新的是Hessian本身數(shù)值表現(xiàn)通常比DFP更穩(wěn)。在實際工程中BFGS是默認(rèn)選擇DFP更多作為教學(xué)和對比出現(xiàn)SR1則在特定場景下比如需要高精度近似Hessian有獨特價值。理解它們的推導(dǎo)差異比死記公式重要得多。3. 核心公式推導(dǎo)與實操要點3.1 SR1公式最簡潔的秩一修正SR1的思路是用秩一矩陣u * u^T去修正B_k使得新的B_{k1}滿足擬牛頓方程。設(shè)B_{k1} B_k σ * v * v^T要求B_{k1} * s_k y_k代入得B_k * s_k σ * v * (v^T * s_k) y_k令v y_k - B_k * s_k可以推導(dǎo)出B_{k1} B_k (y_k - B_k*s_k) * (y_k - B_k*s_k)^T / ((y_k - B_k*s_k)^T * s_k)這就是SR1的更新公式。它的優(yōu)點是形式簡單只需要一次矩陣-向量乘法計算量小。而且它不需要滿足curvature condition適用面更廣。但SR1有個致命問題分母(y_k - B_k*s_k)^T * s_k可能接近零甚至為負(fù)。接近零的時候更新量會爆炸為負(fù)的時候更新后的矩陣可能失去正定性。所以SR1在實際使用中必須加保護當(dāng)分母的絕對值小于某個閾值時跳過這次更新。注意SR1不保證正定性所以它通常和信任域方法配合使用而不是直接做線搜索。如果你用SR1做線搜索方向可能不是下降方向迭代會不穩(wěn)定。3.2 DFP公式逆Hessian的秩二更新DFP更新的是逆Hessian近似矩陣H_k要求滿足H_{k1} * y_k s_kDFP的更新形式是秩二修正H_{k1} H_k s_k*s_k^T/(s_k^T*y_k) - H_k*y_k*y_k^T*H_k/(y_k^T*H_k*y_k)這個公式的推導(dǎo)思路是假設(shè)H_{k1} H_k a*u*u^T b*v*v^T然后代入擬牛頓方程解出系數(shù)和向量。最終得到的就是上面這個形式。DFP的優(yōu)點是當(dāng)s_k^T * y_k 0時更新后的H_{k1}保持正定。這個條件叫做curvature condition在目標(biāo)函數(shù)是凸函數(shù)且步長合適的時候自然滿足。但DFP有個數(shù)值上的弱點它對H_k的初始值比較敏感如果初始H_0選得不好前期迭代可能不穩(wěn)定。而且DFP在接近最優(yōu)解時對Hessian的近似精度不如BFGS。3.3 BFGS公式目前最穩(wěn)的擬牛頓更新BFGS更新的是Hessian近似矩陣B_k要求B_{k1} * s_k y_k更新公式為B_{k1} B_k y_k*y_k^T/(y_k^T*s_k) - B_k*s_k*s_k^T*B_k/(s_k^T*B_k*s_k)如果你需要逆Hessian近似可以用Sherman-Morrison-Woodbury公式對上面這個式子求逆得到直接更新H_k的版本H_{k1} (I - ρ*s_k*y_k^T) * H_k * (I - ρ*y_k*s_k^T) ρ*s_k*s_k^T其中ρ 1/(y_k^T * s_k)。BFGS和DFP的關(guān)系很有意思如果把DFP公式里的s和y互換再取逆就得到BFGS。所以它們本質(zhì)上是“對偶”的。但實際使用中BFGS對Hessian的近似更準(zhǔn)確尤其是在目標(biāo)函數(shù)非凸或者曲率變化劇烈的時候BFGS的表現(xiàn)明顯更好。實操心得如果你只打算實現(xiàn)一個擬牛頓方法直接上BFGS。DFP可以作為對比實驗但生產(chǎn)環(huán)境里基本沒人用DFP了。SR1則留給信任域方法或者需要高精度Hessian近似的場景。3.4 三個公式的對比速查特性SR1DFPBFGS更新秩122更新對象Hessian逆HessianHessian或逆Hessian正定性保證無有條件有條件計算量最低中等中等數(shù)值穩(wěn)定性一般中等最好典型應(yīng)用信任域教學(xué)對比通用優(yōu)化4. 從零實現(xiàn)三個算法的完整流程4.1 環(huán)境準(zhǔn)備與基礎(chǔ)函數(shù)我用Python來實現(xiàn)依賴只有numpy。先定義目標(biāo)函數(shù)和梯度。這里用經(jīng)典的Rosenbrock函數(shù)做測試因為它非凸、有狹長山谷能很好地檢驗優(yōu)化器的穩(wěn)定性。import numpy as np def rosenbrock(x): return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 (1-x[:-1])**2.0) def rosenbrock_grad(x): grad np.zeros_like(x) n len(x) for i in range(n-1): grad[i] -400*x[i]*(x[i1]-x[i]**2) - 2*(1-x[i]) grad[i1] 200*(x[i1]-x[i]**2) return gradRosenbrock函數(shù)的最小值在全是1的向量處函數(shù)值為0。維度可以自己設(shè)一般用2維或10維測試。4.2 BFGS的完整實現(xiàn)BFGS的實現(xiàn)分幾個部分初始化、迭代循環(huán)、線搜索、矩陣更新。def bfgs(func, grad_func, x0, max_iter1000, tol1e-6): n len(x0) x x0.copy() H np.eye(n) # 初始逆Hessian近似 g grad_func(x) for k in range(max_iter): if np.linalg.norm(g) tol: break # 計算搜索方向 d -H g # Armijo線搜索 alpha 1.0 c 1e-4 while func(x alpha*d) func(x) c*alpha*np.dot(g, d): alpha * 0.5 if alpha 1e-10: break s alpha * d x_new x s g_new grad_func(x_new) y g_new - g # BFGS更新 sy np.dot(s, y) if sy 1e-10: # curvature condition rho 1.0 / sy I np.eye(n) H (I - rho*np.outer(s, y)) H (I - rho*np.outer(y, s)) rho*np.outer(s, s) x x_new g g_new return x, func(x), k這段代碼里有兩個關(guān)鍵點線搜索和curvature condition檢查。線搜索用Armijo回溯保證每一步都充分下降curvature condition檢查sy 1e-10防止更新后H失去正定性。4.3 DFP的實現(xiàn)差異DFP和BFGS的代碼結(jié)構(gòu)幾乎一樣區(qū)別只在更新公式def dfp_update(H, s, y): sy np.dot(s, y) if sy 1e-10: return H Hy H y yHy np.dot(y, Hy) H_new H np.outer(s, s)/sy - np.outer(Hy, Hy)/yHy return H_new注意DFP更新的是逆Hessian所以搜索方向還是d -H g。但DFP的更新公式里沒有BFGS那種“夾心”結(jié)構(gòu)數(shù)值穩(wěn)定性差一些。4.4 SR1的實現(xiàn)與保護機制SR1更新Hessian本身所以搜索方向需要解線性方程組def sr1_update(B, s, y): v y - B s vs np.dot(v, s) if abs(vs) 1e-8 * np.linalg.norm(v) * np.linalg.norm(s): return B # 跳過更新 B_new B np.outer(v, v) / vs return B_newSR1的搜索方向是d -np.linalg.solve(B, g)。因為B不保證正定所以線搜索必須用更保守的策略或者直接用信任域方法。注意SR1的分母保護閾值不能設(shè)得太死。我一般用1e-8 * norm(v) * norm(s)做相對判斷比絕對閾值更穩(wěn)。如果分母太小更新量會爆炸直接跳過比強行更新好。4.5 測試與結(jié)果對比用Rosenbrock函數(shù)跑一下三個算法初始點設(shè)為[-1.2, 1.0]x0 np.array([-1.2, 1.0]) x_bfgs, f_bfgs, iter_bfgs bfgs(rosenbrock, rosenbrock_grad, x0) print(fBFGS: x{x_bfgs}, f{f_bfgs:.2e}, iter{iter_bfgs})實測下來BFGS一般在30-50次迭代內(nèi)收斂到1e-6精度DFP需要60-80次SR1如果保護得當(dāng)也能收斂但迭代次數(shù)波動較大。這個差距在低維不明顯維度升到50以上時BFGS的優(yōu)勢會非常突出。5. 常見問題與排查技巧實錄5.1 迭代不收斂或者發(fā)散怎么辦這是最常見的問題。排查順序建議這樣檢查梯度計算是否正確。用有限差分驗證(f(xeps) - f(x-eps))/(2*eps)和解析梯度對比誤差應(yīng)該在1e-6以內(nèi)。檢查線搜索是否滿足Armijo條件。如果步長一直縮小到1e-10還沒找到合適步長說明搜索方向可能不是下降方向。檢查curvature condition。如果s^T*y經(jīng)常為負(fù)說明步長太大或者函數(shù)非凸嚴(yán)重需要減小初始步長。檢查初始H矩陣。一般用單位矩陣但如果變量尺度差異大可以用對角矩陣做縮放。5.2 矩陣更新后失去正定性BFGS和DFP在理論上保證正定性但前提是s^T*y 0。實際中如果線搜索不精確這個條件可能不滿足。解決方法加curvature condition檢查不滿足就跳過更新。用damped BFGS當(dāng)s^T*y太小時對y做修正強行讓它滿足條件。定期重置H為單位矩陣。5.3 高維問題的內(nèi)存和計算瓶頸BFGS需要存儲n×n的H矩陣n10000時就是800MBfloat64內(nèi)存直接爆掉。這時候要用L-BFGS不存完整矩陣只存最近m步的s和y用兩步遞歸算搜索方向。m一般取5-20內(nèi)存占用降到O(mn)。def lbfgs_direction(g, s_list, y_list, H01.0): # 兩步遞歸 q g.copy() alphas [] for s, y in zip(reversed(s_list), reversed(y_list)): rho 1.0 / np.dot(y, s) alpha rho * np.dot(s, q) alphas.append(alpha) q - alpha * y r H0 * q for s, y, alpha in zip(s_list, y_list, reversed(alphas)): rho 1.0 / np.dot(y, s) beta rho * np.dot(y, r) r s * (alpha - beta) return -r5.4 常見問題速查表問題現(xiàn)象可能原因解決方法迭代次數(shù)異常多線搜索太保守增大初始步長或改用Wolfe條件函數(shù)值震蕩步長過大減小步長加Armijo回溯H矩陣失去正定curvature condition不滿足跳過更新或damped BFGS梯度范數(shù)不下降搜索方向錯誤檢查H更新公式的符號內(nèi)存溢出維度太高改用L-BFGSSR1更新爆炸分母接近零加相對閾值保護實操心得我踩過最大的坑是線搜索的初始步長。一開始用固定步長1.0結(jié)果Rosenbrock函數(shù)前幾步直接飛出去。后來改成回溯線搜索但初始步長還是1.0只是加了縮小機制。再后來發(fā)現(xiàn)對于尺度差異大的問題初始步長應(yīng)該根據(jù)梯度范數(shù)自適應(yīng)比如alpha0 min(1.0, 1.0/norm(g))這樣前期穩(wěn)定很多。6. 參數(shù)選擇與性能調(diào)優(yōu)經(jīng)驗6.1 初始H矩陣怎么選最省事的是單位矩陣H0 I。但如果變量尺度差異大比如一個變量范圍是0-1另一個是0-1000單位矩陣會導(dǎo)致搜索方向偏向尺度大的變量。這時候可以用對角縮放H0 diag(1/|g_i|) 或者 H0 (s^T*y / y^T*y) * I第二種是BFGS的常用初始化叫H0 scaling效果通常比單位矩陣好。6.2 線搜索條件的選擇Armijo條件只保證充分下降不保證步長不會太小。Wolfe條件額外要求曲率足夠大能避免步長過小。BFGS配合Wolfe線搜索理論上能保證全局收斂。但Wolfe線搜索實現(xiàn)復(fù)雜一些需要同時檢查函數(shù)值和梯度。實際寫代碼時如果不想實現(xiàn)Wolfe可以用Armijo加上“步長不能太小”的保護。如果步長小于1e-10直接退出或者重置H。6.3 收斂判據(jù)的設(shè)置常用的收斂判據(jù)有三個梯度范數(shù)norm(g) tol最直接推薦用。函數(shù)值變化|f_new - f_old| tol適合函數(shù)值接近零的情況。步長范數(shù)norm(s) tol適合變量尺度小的情況。我一般用梯度范數(shù)為主函數(shù)值變化為輔。tol設(shè)1e-6對大多數(shù)問題夠用需要高精度時設(shè)1e-8。6.4 不同算法的調(diào)優(yōu)側(cè)重點BFGS的調(diào)優(yōu)空間不大主要調(diào)線搜索和初始H。DFP需要更小心地處理curvature condition因為它的更新公式對y^T*H*y更敏感。SR1的調(diào)優(yōu)重點是分母保護閾值和信任域半徑如果做線搜索步長要更保守。最后再分享一個小技巧如果你在實現(xiàn)BFGS時發(fā)現(xiàn)收斂曲線有平臺期可以試試每隔一定迭代次數(shù)重置H為單位矩陣。這個操作看起來粗暴但實際能跳出一些數(shù)值上的停滯。我在處理一些病態(tài)問題時用過效果比調(diào)線搜索參數(shù)還明顯。 擬牛頓算法SR1、DFP、BFGS公式推導(dǎo)與Python實現(xiàn)
1. 為什么擬牛頓算法值得你花時間搞懂如果你做過一點數(shù)值優(yōu)化或者機器學(xué)習(xí)模型的訓(xùn)練大概率聽過“牛頓法”這三個字。牛頓法本身很優(yōu)雅用二階導(dǎo)數(shù)信息Hessian矩陣來指導(dǎo)迭代方向收斂速度是二次的理論上非常漂亮。但問題也很直接Hessian矩陣的計算代價太高尤其是變量維度一上去求逆矩陣的復(fù)雜度直接爆炸。所以實際工程里純牛頓法基本只存在于教科書和小規(guī)模問題中。擬牛頓算法就是在這個背景下被提出來的。它的核心思路很樸素我不去精確算Hessian而是用迭代過程中產(chǎn)生的梯度差和步長差去近似Hessian或者它的逆矩陣。每走一步就用新得到的信息去更新這個近似矩陣。這樣既保留了牛頓法“利用曲率信息”的優(yōu)勢又把計算量壓到了可以接受的范圍。SR1、DFP、BFGS是擬牛頓家族里最經(jīng)典的三個更新公式。SR1是Symmetric Rank-One的縮寫顧名思義它用秩一矩陣去修正近似HessianDFP是Davidon-Fletcher-Powell三個人名字的首字母直接更新Hessian的逆矩陣BFGS則是Broyden-Fletcher-Goldfarb-Shanno是目前公認(rèn)綜合表現(xiàn)最穩(wěn)、應(yīng)用最廣的擬牛頓更新公式。你在scipy、sklearn、各種深度學(xué)習(xí)優(yōu)化器里看到的“L-BFGS”“BFGS”根都在這里。這篇文章適合誰看如果你正在學(xué)數(shù)值優(yōu)化、準(zhǔn)備面試被問到優(yōu)化算法、或者自己寫代碼時想手?jǐn)]一個擬牛頓優(yōu)化器那這篇內(nèi)容可以直接當(dāng)參考。我會把三個公式的推導(dǎo)邏輯、代碼實現(xiàn)、參數(shù)選擇、踩坑經(jīng)驗全部攤開講不跳步不堆公式嚇人盡量讓你看完能自己寫出來、跑起來、調(diào)得動。2. 擬牛頓算法的整體設(shè)計思路拆解2.1 牛頓法到底卡在哪里先把牛頓法的迭代公式擺出來x_{k1} x_k - H_k^{-1} * g_k其中g(shù)_k是梯度H_k是Hessian矩陣。這個公式的意思是在當(dāng)前點用二階泰勒展開去近似目標(biāo)函數(shù)然后直接跳到這個二次模型的極小值點。問題出在H_k^{-1}上。假設(shè)你有n個變量Hessian是n×n矩陣。每次迭代都要重新計算Hessian需要n2個二階偏導(dǎo)然后求逆復(fù)雜度O(n3)。n1000的時候這個計算量已經(jīng)很難接受了n10000的時候基本不用想。更麻煩的是Hessian還不一定正定。如果Hessian不是正定的牛頓方向可能指向鞍點甚至極大值方向迭代直接發(fā)散。所以純牛頓法在實際中很少直接用除非問題規(guī)模很小且性質(zhì)很好。2.2 擬牛頓的核心思想用梯度差代替二階導(dǎo)擬牛頓法的切入點很巧妙。它不去算Hessian而是觀察到一個關(guān)系H_k * (x_{k1} - x_k) ≈ g_{k1} - g_k這個式子來自泰勒展開梯度在兩點之間的變化近似等于Hessian乘以位移。令s_k x_{k1} - x_ky_k g_{k1} - g_k就得到擬牛頓方程B_k * s_k y_k其中B_k是Hessian的近似矩陣?；蛘叻催^來用逆矩陣表示H_k * y_k s_k其中H_k是Hessian逆矩陣的近似。擬牛頓法的每一步就是在滿足這個方程的前提下用一個低秩矩陣去修正上一步的近似矩陣。不同的修正方式就產(chǎn)生了SR1、DFP、BFGS等不同公式。2.3 三個公式的定位差異這三個公式不是隨便列的它們各自有明確的適用場景和數(shù)學(xué)性質(zhì)SR1秩一更新形式最簡單不保證正定性但近似精度高適合信任域方法。DFP秩二更新保證正定性在curvature condition滿足時直接更新逆Hessian是BFGS的前身。BFGS秩二更新同樣保證正定性但更新的是Hessian本身數(shù)值表現(xiàn)通常比DFP更穩(wěn)。在實際工程中BFGS是默認(rèn)選擇DFP更多作為教學(xué)和對比出現(xiàn)SR1則在特定場景下比如需要高精度近似Hessian有獨特價值。理解它們的推導(dǎo)差異比死記公式重要得多。3. 核心公式推導(dǎo)與實操要點3.1 SR1公式最簡潔的秩一修正SR1的思路是用秩一矩陣u * u^T去修正B_k使得新的B_{k1}滿足擬牛頓方程。設(shè)B_{k1} B_k σ * v * v^T要求B_{k1} * s_k y_k代入得B_k * s_k σ * v * (v^T * s_k) y_k令v y_k - B_k * s_k可以推導(dǎo)出B_{k1} B_k (y_k - B_k*s_k) * (y_k - B_k*s_k)^T / ((y_k - B_k*s_k)^T * s_k)這就是SR1的更新公式。它的優(yōu)點是形式簡單只需要一次矩陣-向量乘法計算量小。而且它不需要滿足curvature condition適用面更廣。但SR1有個致命問題分母(y_k - B_k*s_k)^T * s_k可能接近零甚至為負(fù)。接近零的時候更新量會爆炸為負(fù)的時候更新后的矩陣可能失去正定性。所以SR1在實際使用中必須加保護當(dāng)分母的絕對值小于某個閾值時跳過這次更新。注意SR1不保證正定性所以它通常和信任域方法配合使用而不是直接做線搜索。如果你用SR1做線搜索方向可能不是下降方向迭代會不穩(wěn)定。3.2 DFP公式逆Hessian的秩二更新DFP更新的是逆Hessian近似矩陣H_k要求滿足H_{k1} * y_k s_kDFP的更新形式是秩二修正H_{k1} H_k s_k*s_k^T/(s_k^T*y_k) - H_k*y_k*y_k^T*H_k/(y_k^T*H_k*y_k)這個公式的推導(dǎo)思路是假設(shè)H_{k1} H_k a*u*u^T b*v*v^T然后代入擬牛頓方程解出系數(shù)和向量。最終得到的就是上面這個形式。DFP的優(yōu)點是當(dāng)s_k^T * y_k 0時更新后的H_{k1}保持正定。這個條件叫做curvature condition在目標(biāo)函數(shù)是凸函數(shù)且步長合適的時候自然滿足。但DFP有個數(shù)值上的弱點它對H_k的初始值比較敏感如果初始H_0選得不好前期迭代可能不穩(wěn)定。而且DFP在接近最優(yōu)解時對Hessian的近似精度不如BFGS。3.3 BFGS公式目前最穩(wěn)的擬牛頓更新BFGS更新的是Hessian近似矩陣B_k要求B_{k1} * s_k y_k更新公式為B_{k1} B_k y_k*y_k^T/(y_k^T*s_k) - B_k*s_k*s_k^T*B_k/(s_k^T*B_k*s_k)如果你需要逆Hessian近似可以用Sherman-Morrison-Woodbury公式對上面這個式子求逆得到直接更新H_k的版本H_{k1} (I - ρ*s_k*y_k^T) * H_k * (I - ρ*y_k*s_k^T) ρ*s_k*s_k^T其中ρ 1/(y_k^T * s_k)。BFGS和DFP的關(guān)系很有意思如果把DFP公式里的s和y互換再取逆就得到BFGS。所以它們本質(zhì)上是“對偶”的。但實際使用中BFGS對Hessian的近似更準(zhǔn)確尤其是在目標(biāo)函數(shù)非凸或者曲率變化劇烈的時候BFGS的表現(xiàn)明顯更好。實操心得如果你只打算實現(xiàn)一個擬牛頓方法直接上BFGS。DFP可以作為對比實驗但生產(chǎn)環(huán)境里基本沒人用DFP了。SR1則留給信任域方法或者需要高精度Hessian近似的場景。3.4 三個公式的對比速查特性SR1DFPBFGS更新秩122更新對象Hessian逆HessianHessian或逆Hessian正定性保證無有條件有條件計算量最低中等中等數(shù)值穩(wěn)定性一般中等最好典型應(yīng)用信任域教學(xué)對比通用優(yōu)化4. 從零實現(xiàn)三個算法的完整流程4.1 環(huán)境準(zhǔn)備與基礎(chǔ)函數(shù)我用Python來實現(xiàn)依賴只有numpy。先定義目標(biāo)函數(shù)和梯度。這里用經(jīng)典的Rosenbrock函數(shù)做測試因為它非凸、有狹長山谷能很好地檢驗優(yōu)化器的穩(wěn)定性。import numpy as np def rosenbrock(x): return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 (1-x[:-1])**2.0) def rosenbrock_grad(x): grad np.zeros_like(x) n len(x) for i in range(n-1): grad[i] -400*x[i]*(x[i1]-x[i]**2) - 2*(1-x[i]) grad[i1] 200*(x[i1]-x[i]**2) return gradRosenbrock函數(shù)的最小值在全是1的向量處函數(shù)值為0。維度可以自己設(shè)一般用2維或10維測試。4.2 BFGS的完整實現(xiàn)BFGS的實現(xiàn)分幾個部分初始化、迭代循環(huán)、線搜索、矩陣更新。def bfgs(func, grad_func, x0, max_iter1000, tol1e-6): n len(x0) x x0.copy() H np.eye(n) # 初始逆Hessian近似 g grad_func(x) for k in range(max_iter): if np.linalg.norm(g) tol: break # 計算搜索方向 d -H g # Armijo線搜索 alpha 1.0 c 1e-4 while func(x alpha*d) func(x) c*alpha*np.dot(g, d): alpha * 0.5 if alpha 1e-10: break s alpha * d x_new x s g_new grad_func(x_new) y g_new - g # BFGS更新 sy np.dot(s, y) if sy 1e-10: # curvature condition rho 1.0 / sy I np.eye(n) H (I - rho*np.outer(s, y)) H (I - rho*np.outer(y, s)) rho*np.outer(s, s) x x_new g g_new return x, func(x), k這段代碼里有兩個關(guān)鍵點線搜索和curvature condition檢查。線搜索用Armijo回溯保證每一步都充分下降curvature condition檢查sy 1e-10防止更新后H失去正定性。4.3 DFP的實現(xiàn)差異DFP和BFGS的代碼結(jié)構(gòu)幾乎一樣區(qū)別只在更新公式def dfp_update(H, s, y): sy np.dot(s, y) if sy 1e-10: return H Hy H y yHy np.dot(y, Hy) H_new H np.outer(s, s)/sy - np.outer(Hy, Hy)/yHy return H_new注意DFP更新的是逆Hessian所以搜索方向還是d -H g。但DFP的更新公式里沒有BFGS那種“夾心”結(jié)構(gòu)數(shù)值穩(wěn)定性差一些。4.4 SR1的實現(xiàn)與保護機制SR1更新Hessian本身所以搜索方向需要解線性方程組def sr1_update(B, s, y): v y - B s vs np.dot(v, s) if abs(vs) 1e-8 * np.linalg.norm(v) * np.linalg.norm(s): return B # 跳過更新 B_new B np.outer(v, v) / vs return B_newSR1的搜索方向是d -np.linalg.solve(B, g)。因為B不保證正定所以線搜索必須用更保守的策略或者直接用信任域方法。注意SR1的分母保護閾值不能設(shè)得太死。我一般用1e-8 * norm(v) * norm(s)做相對判斷比絕對閾值更穩(wěn)。如果分母太小更新量會爆炸直接跳過比強行更新好。4.5 測試與結(jié)果對比用Rosenbrock函數(shù)跑一下三個算法初始點設(shè)為[-1.2, 1.0]x0 np.array([-1.2, 1.0]) x_bfgs, f_bfgs, iter_bfgs bfgs(rosenbrock, rosenbrock_grad, x0) print(fBFGS: x{x_bfgs}, f{f_bfgs:.2e}, iter{iter_bfgs})實測下來BFGS一般在30-50次迭代內(nèi)收斂到1e-6精度DFP需要60-80次SR1如果保護得當(dāng)也能收斂但迭代次數(shù)波動較大。這個差距在低維不明顯維度升到50以上時BFGS的優(yōu)勢會非常突出。5. 常見問題與排查技巧實錄5.1 迭代不收斂或者發(fā)散怎么辦這是最常見的問題。排查順序建議這樣檢查梯度計算是否正確。用有限差分驗證(f(xeps) - f(x-eps))/(2*eps)和解析梯度對比誤差應(yīng)該在1e-6以內(nèi)。檢查線搜索是否滿足Armijo條件。如果步長一直縮小到1e-10還沒找到合適步長說明搜索方向可能不是下降方向。檢查curvature condition。如果s^T*y經(jīng)常為負(fù)說明步長太大或者函數(shù)非凸嚴(yán)重需要減小初始步長。檢查初始H矩陣。一般用單位矩陣但如果變量尺度差異大可以用對角矩陣做縮放。5.2 矩陣更新后失去正定性BFGS和DFP在理論上保證正定性但前提是s^T*y 0。實際中如果線搜索不精確這個條件可能不滿足。解決方法加curvature condition檢查不滿足就跳過更新。用damped BFGS當(dāng)s^T*y太小時對y做修正強行讓它滿足條件。定期重置H為單位矩陣。5.3 高維問題的內(nèi)存和計算瓶頸BFGS需要存儲n×n的H矩陣n10000時就是800MBfloat64內(nèi)存直接爆掉。這時候要用L-BFGS不存完整矩陣只存最近m步的s和y用兩步遞歸算搜索方向。m一般取5-20內(nèi)存占用降到O(mn)。def lbfgs_direction(g, s_list, y_list, H01.0): # 兩步遞歸 q g.copy() alphas [] for s, y in zip(reversed(s_list), reversed(y_list)): rho 1.0 / np.dot(y, s) alpha rho * np.dot(s, q) alphas.append(alpha) q - alpha * y r H0 * q for s, y, alpha in zip(s_list, y_list, reversed(alphas)): rho 1.0 / np.dot(y, s) beta rho * np.dot(y, r) r s * (alpha - beta) return -r5.4 常見問題速查表問題現(xiàn)象可能原因解決方法迭代次數(shù)異常多線搜索太保守增大初始步長或改用Wolfe條件函數(shù)值震蕩步長過大減小步長加Armijo回溯H矩陣失去正定curvature condition不滿足跳過更新或damped BFGS梯度范數(shù)不下降搜索方向錯誤檢查H更新公式的符號內(nèi)存溢出維度太高改用L-BFGSSR1更新爆炸分母接近零加相對閾值保護實操心得我踩過最大的坑是線搜索的初始步長。一開始用固定步長1.0結(jié)果Rosenbrock函數(shù)前幾步直接飛出去。后來改成回溯線搜索但初始步長還是1.0只是加了縮小機制。再后來發(fā)現(xiàn)對于尺度差異大的問題初始步長應(yīng)該根據(jù)梯度范數(shù)自適應(yīng)比如alpha0 min(1.0, 1.0/norm(g))這樣前期穩(wěn)定很多。6. 參數(shù)選擇與性能調(diào)優(yōu)經(jīng)驗6.1 初始H矩陣怎么選最省事的是單位矩陣H0 I。但如果變量尺度差異大比如一個變量范圍是0-1另一個是0-1000單位矩陣會導(dǎo)致搜索方向偏向尺度大的變量。這時候可以用對角縮放H0 diag(1/|g_i|) 或者 H0 (s^T*y / y^T*y) * I第二種是BFGS的常用初始化叫H0 scaling效果通常比單位矩陣好。6.2 線搜索條件的選擇Armijo條件只保證充分下降不保證步長不會太小。Wolfe條件額外要求曲率足夠大能避免步長過小。BFGS配合Wolfe線搜索理論上能保證全局收斂。但Wolfe線搜索實現(xiàn)復(fù)雜一些需要同時檢查函數(shù)值和梯度。實際寫代碼時如果不想實現(xiàn)Wolfe可以用Armijo加上“步長不能太小”的保護。如果步長小于1e-10直接退出或者重置H。6.3 收斂判據(jù)的設(shè)置常用的收斂判據(jù)有三個梯度范數(shù)norm(g) tol最直接推薦用。函數(shù)值變化|f_new - f_old| tol適合函數(shù)值接近零的情況。步長范數(shù)norm(s) tol適合變量尺度小的情況。我一般用梯度范數(shù)為主函數(shù)值變化為輔。tol設(shè)1e-6對大多數(shù)問題夠用需要高精度時設(shè)1e-8。6.4 不同算法的調(diào)優(yōu)側(cè)重點BFGS的調(diào)優(yōu)空間不大主要調(diào)線搜索和初始H。DFP需要更小心地處理curvature condition因為它的更新公式對y^T*H*y更敏感。SR1的調(diào)優(yōu)重點是分母保護閾值和信任域半徑如果做線搜索步長要更保守。最后再分享一個小技巧如果你在實現(xiàn)BFGS時發(fā)現(xiàn)收斂曲線有平臺期可以試試每隔一定迭代次數(shù)重置H為單位矩陣。這個操作看起來粗暴但實際能跳出一些數(shù)值上的停滯。我在處理一些病態(tài)問題時用過效果比調(diào)線搜索參數(shù)還明顯。