當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

UA MATH575B 数值分析下II 牛顿算法

發布時間：2025/4/14 编程问答 51 豆豆

生活随笔收集整理的這篇文章主要介紹了 UA MATH575B 数值分析下II 牛顿算法小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

UA MATH575B 數值分析下II 牛頓算法

Pure Newton算法
- Damped Newton算法
- Levenberg-Marquardt算法
Quasi-Newton算法
- 割線法
- Broyden–Fletcher–Goldfarb–Shanno算法
示例

Pure Newton算法

考慮優化問題
$min_{x} f(x)$
Pure Newton算法其實就是把Newton-Raphson算法用來接優化的一階條件 $?f(x)=0\nabla f(x)=0$ ，遞推公式為
$xk+1=xk?H?1f(xk)?f(xk)x_{k+1} = x_{k} - H^{-1}f(x_{k}) \nabla f(x_{k})$
其中 $H^{-1}$ 代表Hessian矩陣的逆。因為 $?f(x)=0\nabla f(x)=0$ 的解不一定是全局最優解，它還可能是局部最優解或者鞍點，這個是Pure Newton算法最大的問題。還有就是這個下降的步長完全是由 $H?1f(xk)?f(xk)H^{-1}f(x_{k}) \nabla f(x_{k})$ 計算決定的，沒有控制步長的操作，可能導致下降路徑步與步之間變動非常大。下面兩種算法是針對這個缺點的改進。

Damped Newton算法

阻尼牛頓算法在牛頓算法的基礎上給步長加了一個阻尼因子 $t$
$xk+1=xk?tH?1f(xk)?f(xk)=xk+tΔxkx_{k+1} = x_{k} - tH^{-1}f(x_{k}) \nabla f(x_{k}) = x_{k} + t \Delta x_{k}$
在寫算法的時候這個 $t$ 用線性回溯搜索來找
$\argmin_{t \ge 0} f(x_k + t \Delta x_k)$
具體實現的時候可以先給定一個不算小的初值，逐漸減少直到
$f(xk+tΔxk)<f(xk)+α?f(xk)?tΔxkf(x_{k} + t \Delta x_{k}) < f(x_{k}) + \alpha \nabla f(x_{k}) \cdot t \Delta x_{k}$
$α\alpha$ 在[0,1)上取值，其意義在于當下降方向比較合適的時候，允許這一步多下降一點。

Levenberg-Marquardt算法

Levenberg-Marquardt算法在阻尼牛頓算法的基礎上又給Hessian矩陣加了一個松弛項
$xk+1=xk?t(Hf(xk)+μI)?1?f(xk)x_{k+1} = x_{k} - t(Hf(x_{k})+\mu I)^{-1} \nabla f(x_{k})$
當 $μ\mu$ 等于0時，這個算法變回阻尼牛頓算法，當 $μ\mu$ 趨近于無窮時，這個算法就成了梯度下降。

Quasi-Newton算法

割線法

相信大家都有這個感覺，牛頓算法要寫代碼就要算兩階偏導，簡單的Matlab算一算就好，復雜一點的Matlab算出來可能是錯的也可能算不出來，手算也會很麻煩。考慮一元函數，割線法就可以避免繁瑣的求導，用差商代替求導，幾何上就相當于割線代替切線，所以叫割線法。它的遞推公式長這樣，假設 $g = f^{'}$
$xk+1=xk?g(xk)g′(xk)x_{k+1} = x_{k} - \frac{g(x_{k})}{g'(x_{k})}$
割線法有兩個問題，收斂比較慢以及只能處理一元函數的優化。BFGS算法把割線法的思想推廣到多元，即想辦法用數值計算步驟替換求導，盡管無法解決收斂比牛頓算法慢的缺陷，但至少避免了多元函數求導的計算。

Broyden–Fletcher–Goldfarb–Shanno算法

BFGS算法的核心是在更新x的同時也更新Hessian矩陣的逆，從而避免求導和求逆的操作。這里給一個BFGS算法的簡單推導。用 $Δxk\Delta x_{k}$ 表示第 $k$ 步的步長， $C^k\hat{C}_{k}$ 表示第 $k$ 步對 $H^{-1}f(x_{k})$ 的近似，則
$Δxk=?H?1f(xk)?f(xk)≈?C^k?f(xk)\Delta x_{k} = -H^{-1}f(x_{k}) \nabla f(x_{k}) \approx -\hat{C}_{k}\nabla f(x_{k})$
類似阻尼牛頓法，給步長加個阻尼 $t$ ，定義 $dk=tΔxkd_k = t\Delta x_{k}$ ，阻尼因子同樣由線性回溯搜索確定。遞推方程為
$x_{k+1} = x_{k} + d_{k}$
現在推導 $C^k\hat{C}_{k}$ 的遞推方程。先定義梯度的差分
$gk=?f(xk+1)??f(xk)g_k = \nabla f(x_{k+1}) - \nabla f(x_{k})$
考慮梯度的Taylor展開
$?f(xk+1)??f(xk)≈Hf(xk)(xk+1?xk)C^k+1gk≈dk\nabla f(x_{k+1}) - \nabla f(x_{k}) \approx Hf(x_{k}) (x_{k+1}-x_{k}) \\ \hat{C}_{k+1} g_k \approx d_k$
老師上課推了，但我覺得好像有點簡單以為我會的，結果現在寫博客發現我錯了，我不會，先空在這兒吧，等我會了再補上。

示例

最小化函數
$H_2 (x,y)=exp?(8x-13y+21)+exp?(21y-13x+34)+0.0001exp?(x+y)$
這個函數的性質非常詭異，看個示意圖就知道了，它遞增非常快，數值非常大，如果畫圖的(x,y)區域設置得再大一點看到的就是一個平面加一個非常尖的峰。

figure(1); x=-.01:0.0005:0.01; % change region here y=-0.01:0.0005:0.01; [xx,yy]=meshgrid(x,y); zz=exp(8*xx-13*yy+21)+exp(21*yy-13*xx+34)+0.0001*exp(xx+yy); surf(xx,yy,zz); xlabel('x'); ylabel('y') zlabel('H_2(x,y)');

這個函數是個凸函數，最小值存在。先用最簡單的牛頓法試試看怎么樣。先定義這個函數，

function zz = H2(xx,yy) zz=exp(8*xx-13*yy+21)+exp(21*yy-13*xx+34)+0.0001*exp(xx+yy); end

然后定義函數計算牛頓法每一步更新的位移，

function xstep = PureNewtonStep(xx,yy) global M1 M2 M3; H = M1*exp(8*xx-13*yy+21)+M2*exp(21*yy-13*xx+34)+M3*exp(xx+yy); N1 = 8*exp(8*xx-13*yy+21)-13*exp(21*yy-13*xx+34)+0.0001*exp(xx+yy); N2 = -13*exp(8*xx-13*yy+21)+21*exp(21*yy-13*xx+34)+0.0001*exp(xx+yy); xstep = H\[N1;N2]; end

接下來寫一下牛頓法的迭代過程，設置更新的位移范數小于1e-8時停止，初值取(0,0)，記錄下每一步更新之后的位置和迭代次數

K = 0; % # of iterations TOL = 1e-8; % stop criteria xstep = 1; % descent step x_pureNewton = [0;0]; while norm(xstep) > TOLxstep = PureNewtonStep(x_pureNewton(1,end),x_pureNewton(2,end));x1 = x_pureNewton(:,end) - xstep;x_pureNewton = [x_pureNewton x1];K = K + 1; end

運行結果顯示需要迭代30次，最小值為2.116459796742611e+11，此時的位置是(19.571447286191784, 11.691722203621611)。(x, y)的路徑可以畫出來看看

下降的軌跡畫出來也很詭異，從原點出發以后直奔東北，拐個急彎之后又掉頭向西南移動。但事實上這是對的，因為這個函數最小值的通道比較狹窄（事實上這個最小值是局部最小，并不是全局最小值，要找到全局最小值需要更換不同初值調試）。再用BFGS來驗證一下這個結果，另外兩個算法我這次寫作業沒用上，有機會再試試吧。

function [x] = BFGS_H2(x) F = H2(x(1),x(2)); dF = Gradient_H2(x(1),x(2)); C = eye(length(x)); N = 0; while (norm(dF) > 1.e-12 && N < 32)d = -C*dF;while (H2(x(1)+d(1),x(2)+d(2)) >= F)d = d / 2.;if (norm(d) < eps)C = eye(length(x)); d = -dF;endendx = x + d; F = H2(x(1),x(2)); new_dF = Gradient_H2(x(1),x(2)); g = new_dF - dF; dF = new_dF;rho = 1. / (g'*d); mu = rho*(1. + rho*(g'*C*g));C = C - rho*(d*(g'*C) + (C*g)*d') + mu*d*d';N = N+1; end

BFGS算出來最小值點是(19.571447252082269, 11.691722181991503) ，最優的x與牛頓法前9個數字相同，最優的y與牛頓法前8個數字相同，基本是符合設定的容許誤差的。

總結

以上是生活随笔為你收集整理的UA MATH575B 数值分析下II 牛顿算法的全部內容，希望文章能夠幫你解決所遇到的問題。

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

上一篇： UA MATH575B 数值分析下III
下一篇： UA MATH574M 统计学习II 二