与公众号同步更新，详细内容及相关ipynb文件在公众号中，公众号：AI入门小白

文章目录

牛顿法
拟牛顿法的思路
- DFP (Davidon-Fletcher- Powell) 算法(DFP algorithm)
- BFGS (Broyden-Fletcher-Goldfarl-Shanno) 算法(BFGS algorithm)
- Broyden 类算法(Broyden's algorithm)

牛顿法(Newton method) 和拟牛顿法(quasi-Newton method) 也是求解无约束最优化问题的常用方法，有收敛速度快的优点。牛顿法是迭代算法，每一步需要求解目标函数的黑塞矩阵的逆矩阵，计算比较复杂。拟牛顿法通过正定矩阵近似黑塞矩阵的逆矩阵或黑塞矩阵，简化了这一计算过程。

黑塞矩阵是一个多元函数的二阶偏导数构成的方阵，描述了函数的局部曲率。

牛顿法

考虑无约束最优化问题
min⁡x∈Rnf(x)(B.1)\min_{x \in R^n} f(x) \quad \tag{B.1} x∈Rnminf(x)(B.1)
其中x∗x^*x∗为目标函数的极小点。

假设f(x)f(x)f(x)具有二阶连续偏导数，若第kkk次迭代值为x(k)x^{(k)}x(k)，则可将f(x)f(x)f(x)在x(k)x^{(k)}x(k)附近进行二阶泰勒展开：
f(x)=f(x(k))+gkT(x−x(k))+12(x−x(k))TH(x(k))(x−x(k))(B.2)f(x) = f(x^{(k)}) + g_k^T (x - x^{(k)}) + \frac{1}{2} (x - x^{(k)})^T H(x^{(k)})(x - x^{(k)}) \quad \tag{B.2} f(x)=f(x(k))+gkT(x−x(k))+21(x−x(k))TH(x(k))(x−x(k))(B.2)
这里，gk=g(x(k))=▽f(x(k))g_k = g(x^{(k)}) = \triangledown f(x^{(k)})gk=g(x(k))=▽f(x(k))是f(x)f(x)f(x)的梯度向量在点x(k)x^{(k)}x(k)的值，H(x(k))H(x^{(k)})H(x(k))是f(x)f(x)f(x)的黑塞矩阵
H(x)=[∂2f∂xi∂xj]n×n(B.3)H(x) = \Bigg[\frac{\partial^2 f}{\partial x_i \partial x_j} \Bigg]_{n\times n} \quad \tag{B.3} H(x)=[∂xi∂xj∂2f]n×n(B.3)
在点x(k)x^{(k)}x(k)的值。函数f(x)f(x)f(x)有极值的必要条件是在极值点处一阶导数为0 ，即梯度向量为0。特别是当H(x(k))H(x^{(k)})H(x(k))是正定矩阵时，函数f(x)f(x)f(x)的极值为极小值。

牛顿法利用极小点的必要条件
▽f(x)=0(B.4)\triangledown f(x) = 0 \quad \tag{B.4} ▽f(x)=0(B.4)
每次选代中从点x(k)x^{(k)}x(k)开始，求目标函数的极小点，作为第k+1k + 1k+1次迭代值x(k+1)x^{(k+1)}x(k+1)。具体地，假设x(k+1)x^{(k+1)}x(k+1)满足:
▽f(x(k+1))=0(B.5)\triangledown f(x^{(k+1)}) = 0 \quad \tag{B.5} ▽f(x(k+1))=0(B.5)
由式(B.2) 有
▽f(x)=gk+Hk(x−x(k))(B.6)\triangledown f(x) = g_k + H_k(x - x^{(k)}) \quad \tag{B.6} ▽f(x)=gk+Hk(x−x(k))(B.6)
其中Hk=H(x(k))H_k = H(x^{(k)})Hk=H(x(k))。这样，式(B.5)成为
gk+Hk(x(k+1)−x(k))=0(B.7)g_k + H_k (x^{(k+1)} - x^{(k)}) = 0 \quad \tag{B.7} gk+Hk(x(k+1)−x(k))=0(B.7)
因此，
x(k+1)=x(k)−Hk−1gk(B.8)x^{(k+1)} = x^{(k)} - H_k^{-1} g_k \quad \tag{B.8} x(k+1)=x(k)−Hk−1gk(B.8)
或者
x(k+1)=x(k)+pk(B.9)x^{(k+1)} = x^{(k)} + p_k \quad \tag{B.9} x(k+1)=x(k)+pk(B.9)
其中，
Hkpk=−gk(B.10)H_k p_k = -g_k \quad \tag{B.10} Hkpk=−gk(B.10)
用式(B.8)作为迭代公式的算法就是牛顿法。

算法B.1(牛顿法)
输入：目标函数f(x)f(x)f(x)，梯度g(x)=▽f(x)g(x) = \triangledown f(x)g(x)=▽f(x)，黑塞矩阵H(x)H(x)H(x)，精度要求ε\varepsilonε；
输出：f(x)f(x)f(x)的极小点x∗x^*x∗。
(1)取初始点x(0)x^{(0)}x(0)，置k=0k = 0k=0。
(2)计算gk=g(x(k))g_k = g(x^{(k)})gk=g(x(k))。
(3)若∥gk∥<ε\lVert g_k \rVert < \varepsilon∥gk∥<ε，则停止计算，得近似解x∗=x(k)x^* = x^{(k)}x∗=x(k)。
(4)计算Hk=H(x(k))H_k = H(x^{(k)})Hk=H(x(k))，并求pkp_kpk
Hkpk=−gkH_k p_k = -g_k Hkpk=−gk
(5)置x(k+1)=x(k)+pkx^{(k+1)} = x^{(k)} + p_kx(k+1)=x(k)+pk。
(6)置k=k+1k = k+1k=k+1，转(2)

步骤(4) 求pk,pk=−Hk−1gkp_k, p_k = -H_k^{-1} g_kpk,pk=−Hk−1gk，要求Hk−1H_k^{-1}Hk−1，计算比较复杂，所以有其它改进的方法。

拟牛顿法的思路

在牛顿法的迭代中，需要计算黑塞矩阵的逆矩阵H−1H^{-1}H−1，这一计算比较复杂，考虑用一个nnn阶矩阵Gk=G(x(k))G_k = G(x^{(k)})Gk=G(x(k))来近似代替Hk−1=H−1(x(k))H_k^{-1} = H^{-1}(x^{(k)})Hk−1=H−1(x(k))。这就是拟牛顿法的基本想法。

先看牛顿法迭代中黑塞矩阵HkH_kHk满足的条件。首先，HkH_kHk满足以下关系。在式(B.6)中取x=x(k+1)x = x^{(k+1)}x=x(k+1)，即得
gk+1−gk=Hk(x(k+1)−x(k))(B.11)g_{k+1} - g_k = H_k(x^{(k+1)} - x^{(k)}) \quad \tag{B.11} gk+1−gk=Hk(x(k+1)−x(k))(B.11)
记yk=gk+1−gk,δk=x(k+1)−x(k)y_k = g_{k+1} - g_k, \delta_k = x^{(k+1)} - x^{(k)}yk=gk+1−gk,δk=x(k+1)−x(k)，则
yk=Hkδk(B.12)y_k = H_k \delta_k \quad \tag{B.12} yk=Hkδk(B.12)
或
Hk−1yk=δk(B.13)H_k^{-1} y_k = \delta_k \quad \tag{B.13} Hk−1yk=δk(B.13)
式(B.12) 或式(B.13) 称为拟牛顿条件。

如果HkH_kHk是正定的(Hk−1H_k^{-1}Hk−1也是正定的)，那么可以保证牛顿法搜索方向pkp_kpk是下降方向。这是因为搜索方向是pk=−Hk−1gkp_k = -H_k^{-1} g_kpk=−Hk−1gk，由式(B.8) 有
x=x(k)+λpk=x(k)−λHk−1gk(B.14)x = x^{(k)} + \lambda p_k = x^{(k)} - \lambda H_k^{-1} g_k \quad \tag{B.14} x=x(k)+λpk=x(k)−λHk−1gk(B.14)
所以f(x)f(x)f(x)在x(k)x^{(k)}x(k)的泰勒展开式(B.2) 可以近似写成:
f(x)=f(x(k))−λgkTHk−1gk(B.15)f(x) = f(x^{(k)}) - \lambda g_k^T H_k^{-1} g_k \quad \tag{B.15} f(x)=f(x(k))−λgkTHk−1gk(B.15)
因Hk−1H_k^{-1}Hk−1正定，故有gkTHk−1gk>0g_k^T H_k^{-1} g_k > 0gkTHk−1gk>0。当λ\lambdaλ为一个充分小的正数时，总有f(x)<f(x(k))f(x) < f(x^{(k)})f(x)<f(x(k))，也就是说pkp_kpk是下降方向。

拟牛顿法将GkG_kGk作为Hk−1H_k^{-1}Hk−1的近似，要求矩阵GkG_kGk满足同样的条件。首先，每次选代矩阵GkG_kGk是正定的。同时，GkG_kGk满足下面的拟牛顿条件:
Gk+1yk=δk(B.16)G_{k+1} y_k = \delta_k \quad \tag{B.16} Gk+1yk=δk(B.16)
按照拟牛顿条件选择GkG_kGk作为Hk−1H_k^{-1}Hk−1的近似或选择BkB_kBk作为HkH_kHk的近似的算法称为拟牛顿法。

按照拟牛顿条件，在每次选代中可以选择更新矩阵Gk+1G_{k+1}Gk+1：
Gk+1=Gk+ΔGk(B.17)G_{k+1} = G_k + \Delta G_k \quad \tag{B.17} Gk+1=Gk+ΔGk(B.17)
这种选择有一定的灵活性，因此有多种具体实现方法。下面介绍Broyden 类拟牛顿法。

DFP (Davidon-Fletcher- Powell) 算法(DFP algorithm)

DFP算法选择Gk+1G_{k+1}Gk+1的方法是，假设每一步迭代中矩阵Gk+1G_{k+1}Gk+1是由GkG_kGk加上两个附加项构成的，即
Gk+1=Gk+Pk+Qk(B.18)G_{k+1} = G_k + P_k +Q_k \quad \tag{B.18} Gk+1=Gk+Pk+Qk(B.18)
其中Pk,QkP_k, Q_kPk,Qk是待定矩阵。这时，
Gk+1yk=Gkyk+Pkyk+Qkyk(B.19)G_{k+1}y_k = G_k y_k + P_k y_k + Q_k y_k \quad \tag{B.19} Gk+1yk=Gkyk+Pkyk+Qkyk(B.19)
为使Gk+1G_{k+1}Gk+1满足拟牛顿条件，可使PkP_kPk和QkQ_kQk满足:
Pkyk=δk(B.20)P_k y_k = \delta_k \quad \tag{B.20} Pkyk=δk(B.20)
Qkyk=−Gkyk(B.21)Q_k y_k = -G_k y_k \quad \tag{B.21} Qkyk=−Gkyk(B.21)
事实上，不难找出这样的PkP_kPk和QkQ_kQk，例如取:
Pk=δkδkTδkTyk(B.22)P_k = \frac{\delta_k \delta_k^T}{\delta_k^T y_k} \quad \tag{B.22} Pk=δkTykδkδkT(B.22)
Qk=−GkykykTGkykTGkyk(B.23)Q_k = -\frac{G_k y_k y_k^T G_k}{y_k^T G_k y_k} \quad \tag{B.23} Qk=−ykTGkykGkykykTGk(B.23)
这样就可得到矩阵Gk+1G_{k+1}Gk+1的迭代公式:
Gk+1=Gk+δkδkTδkTyk−GkykykTGkykTGkyk(B.24)G_{k+1} = G_k + \frac{\delta_k \delta_k^T}{\delta_k^T y_k} - \frac{G_k y_k y_k^T G_k}{y_k^T G_k y_k} \quad \tag{B.24} Gk+1=Gk+δkTykδkδkT−ykTGkykGkykykTGk(B.24)
称为DFP 算法。

如果初始矩阵G0G_0G0是正定的，则迭代过程中的每个矩阵GkG_kGk都是正定的。

算法B.2(DFP算法)
输入：目标函数f(x)f(x)f(x)，梯度g(x)=∇f(x)g(x) = \nabla f(x)g(x)=∇f(x)，精度要求ε\varepsilonε；
输出：f(x)f(x)f(x)的极小点x∗x^*x∗。
(1)选定初始点x(0)x^{(0)}x(0)，取G0G_0G0为正定对称矩阵，置k=0k=0k=0。
(2)计算gk=g(x(k))g_k = g(x^{(k)})gk=g(x(k))。若∥gk∥<ε\lVert g_k \rVert < \varepsilon∥gk∥<ε，则停止计算，得近似解x∗=x(k)x^* = x^{(k)}x∗=x(k)；否则转(3)。
(3)置pk=−Gkgkp_k = -G_k g_kpk=−Gkgk。
(4)一维搜索:求λk\lambda_kλk使得
f(x(k)+λkpk)=min⁡λ≥0f(x(k)+λpk)f(x^{(k)} + \lambda_k p_k) = \min_{\lambda \geq 0} f(x^{(k)} + \lambda p_k) f(x(k)+λkpk)=λ≥0minf(x(k)+λpk)
(5)置x(k+1)=x(k)+λkpkx^{(k+1)} = x^{(k)} + \lambda_k p_kx(k+1)=x(k)+λkpk。
(6)计算gk+1=g(x(k+1))g_{k+1} = g(x^{(k+1)})gk+1=g(x(k+1))，若∥gk+1∥<ε\lVert g_{k+1} \rVert < \varepsilon∥gk+1∥<ε，则停止计算，得近似解x∗=x(k+1)x^* = x^{(k+1)}x∗=x(k+1)；否则，按式(B.24)算出Gk+1G_{k+1}Gk+1。
(7)置k=k+1k = k+1k=k+1，转(3)。

BFGS (Broyden-Fletcher-Goldfarl-Shanno) 算法(BFGS algorithm)

BFGS 算法是最流行的拟牛顿算法。
可以考虑用GkG_kGk逼近黑塞矩阵的逆矩阵H−1H^{-1}H−1，也可以考虑用BkB_kBk逼近黑塞矩阵HHH。
这时，相应的拟牛顿条件是
Bk+1δk=yk(B.25)B_{k+1} \delta_k = y_k \quad \tag{B.25} Bk+1δk=yk(B.25)
可以用同样的方法得到另一迭代公式。首先令
Bk+1=Bk+Pk+Qk(B.26)B_{k+1} = B_k + P_k + Q_k \quad \tag{B.26} Bk+1=Bk+Pk+Qk(B.26)
Bk+1δk=Bkδk+Pkδk+Qkδk(B.27)B_{k+1} \delta_k = B_k \delta_k + P_k \delta_k + Q_k \delta_k \quad \tag{B.27} Bk+1δk=Bkδk+Pkδk+Qkδk(B.27)
考虑使PkP_kPk和QkQ_kQk满足：
Pkδk=yk(B.28)P_k \delta_k = y_k \quad \tag{B.28} Pkδk=yk(B.28)
Qkδk=−Bkδk(B.29)Q_k \delta_k = -B_k \delta_k \quad \tag{B.29} Qkδk=−Bkδk(B.29)
找出适合条件的PkP_kPk和QkQ_kQk，得到BFGS算法矩阵Bk+1B_{k+1}Bk+1的迭代公式:
Bk+1=Bk+ykykTykTδk−BkδkδkTBkδkTBkδk(B.30)B_{k+1} = B_k + \frac{y_k y_k^T}{y_k^T \delta_k} - \frac{B_k \delta_k \delta_k^T B_k}{\delta_k^T B_k \delta_k} \quad \tag{B.30} Bk+1=Bk+ykTδkykykT−δkTBkδkBkδkδkTBk(B.30)
如果初始矩阵B0B_0B0是正定的，则迭代过程中的每个矩阵BkB_kBk都是正定的。

算法B.3(BFGS算法)
输入：目标函数f(x),g(x)=∇f(x)f(x),g(x) = \nabla f(x)f(x),g(x)=∇f(x)，精度要求ε\varepsilonε；
输出：f(x)f(x)f(x)的极小点x∗x^*x∗。
(1)选定初始点x(0)x^{(0)}x(0)，取B0B_0B0为正定对称矩阵，置k=0k=0k=0。
(2)计算gk=g(x(k))g_k = g(x^{(k)})gk=g(x(k))。若∥gk∥<ε\lVert g_k \rVert < \varepsilon∥gk∥<ε，则停止计算，得近似解x∗=x(k)x^* = x^{(k)}x∗=x(k)；否则转(3)。
(3)由Bkpk=−gkB_k p_k = -g_kBkpk=−gk求出pkp_kpk。
(4)一维搜索：求λk\lambda_kλk使得
f(x(k)+λkpk)=min⁡λ≥0f(x(k)+λpk)f(x^{(k)} + \lambda_k p_k) = \min_{\lambda \geq 0} f(x^{(k)} + \lambda p_k) f(x(k)+λkpk)=λ≥0minf(x(k)+λpk)
(5)置x(k+1)=x(k)+λkpkx^{(k+1)} = x^{(k)} + \lambda_k p_kx(k+1)=x(k)+λkpk。
(6)计算gk+1=g(x(k+1))g_{k+1} = g(x^{(k+1)})gk+1=g(x(k+1))，若∥gk+1∥<ε\lVert g_{k+1} \rVert < \varepsilon∥gk+1∥<ε，则停止计算，得近似解x∗=x(k+1)x^* = x^{(k+1)}x∗=x(k+1)；否则，按式(B.30)算出Bk+1B_{k+1}Bk+1。
(7)置k=k+1k = k + 1k=k+1，转(3)。

Broyden 类算法(Broyden’s algorithm)

Sherman-Morrison 公式:假设AAA是nnn阶可逆矩阵，u,vu, vu,v是nnn维向量，且A+uvTA + uv^TA+uvT也是可逆矩阵，则
(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u(A + uv^T)^{-1} = A^{-1} - \frac{A^{-1}uv^TA^{-1}}{1 + v^TA^{-1}u} (A+uvT)−1=A−1−1+vTA−1uA−1uvTA−1
可以从BFGS 算法矩阵BkB_kBk的迭代公式(B.30) 得到BFGS 算法关于GkG_kGk的迭代公式。事实上，若记Gk=Bk−1,Gk+1=Bk+1−1G_k = B_k^{-1}, G_{k+1} = B_{k+1}^{-1}Gk=Bk−1,Gk+1=Bk+1−1，那么对式(B.30) 两次应用Sherman-Morrison 公式即得
Gk+1=(I−δkykTδkTyk)Gk(I−δkykTδkTyk)T+δkδkTδkTyk(B.31)G_{k+1} = \Bigg(I - \frac{\delta_ky_k^T}{\delta_k^Ty_k} \Bigg)G_k\Bigg(I - \frac{\delta_ky_k^T}{\delta_k^Ty_k} \Bigg)^T + \frac{\delta_k\delta_k^T}{\delta_k^Ty_k} \quad \tag{B.31} Gk+1=(I−δkTykδkykT)Gk(I−δkTykδkykT)T+δkTykδkδkT(B.31)
称为BFGS 算法关于GkG_kGk的迭代公式。

由DFP 算法GkG_kGk的迭代公式(B.23) 得到的Gk+1G_{k+1}Gk+1记作GDFPG^{DFP}GDFP，由BFGS算法GkG_kGk的迭代公式(B.31) 得到的Gk+1G_{k+1}Gk+1记作GBFGSG^{BFGS}GBFGS，它们都满足方程拟牛顿条件式，所以它们的线性组合
Gk+1=αGDFP+(1−α)GBFGS(B.32)G_{k+1} = \alpha G^{DFP} + (1 - \alpha)G^{BFGS} \quad \tag{B.32} Gk+1=αGDFP+(1−α)GBFGS(B.32)
也满足拟牛顿条件式，而且是正定的。其中0≤α≤10 \leq \alpha \leq 10≤α≤1。这样就得到了一类拟牛顿法，称为Broyden 类算法。

数据来源：统计学习方法(第二版)

统计学习方法--牛顿法和拟牛顿法相关推荐

《统计学习方法》啃书辅助：附录B 牛顿法和拟牛顿法
<统计学习方法>啃书辅助:附录B 牛顿法和拟牛顿法 B.1 牛顿法 [基础知识]梯度 [基础知识]浅谈「正定矩阵」和「半正定矩阵」 - Xinyu Chen 的文章 - 知乎二阶泰勒展开 ...
机器学习优化算法中梯度下降,牛顿法和拟牛顿法的优缺点详细介绍
1.梯度下降法梯度下降法实现简单,当目标函数是凸函数时,梯度下降法的解是全局解.一般情况下,其解不保证是全局最优解,梯度下降法的速度也未必是最快的. 梯度下降法的优化思想:用当前位置负梯度方向作为搜 ...
梯度类算法原理：最速下降法、牛顿法和拟牛顿法
文章目录算法结构最速下降法牛顿法拟牛顿法算法结构梯度类算法有很多,本文主要学习最常见的3个算法:最速下降法.牛顿法和拟牛顿法.算法名称虽多,但是他们的算法结构都是一样的,可以描述为 (1) ...
梯度下降法、牛顿法和拟牛顿法——机器学习面试
梯度下降.牛顿.拟牛顿原理梯度下降牛顿法为Hesse矩阵参数更新的方程: 每一次迭代的更新方向都是当前点的牛顿方向,步长固定为1.每一次都需要计算一阶导数以及Hesse矩阵的逆矩阵,对于 ...
二阶偏微分方程组龙格库塔法_牛顿法和拟牛顿法——（书中附录B）
牛顿法(Newton method)和拟牛顿法(quasi-Newton method)也是求解无约束最优化问题的常用方法,具有收敛速度快的优点. 牛顿法是迭代算法,每一步需要求解目标函数的海赛矩阵的 ...
牛顿法（Newton Methods）、阻尼牛顿法和拟牛顿法
令 X=(x1,x2,⋯,xN)T∈RNX=(x_1,x_2,\cdots,x_N)^T \in {\bf R}^NX=(x1,x2,⋯,xN)T∈RN,目标函数 f:RN→Rf:{\bf R} ...
李航-统计学习方法-笔记-1：概论
写在前面本系列笔记主要记录<统计学习方法>中7种常用的机器学习分类算法,包括感知机,KNN,朴素贝叶斯,决策树,逻辑斯谛回归与最大熵模型,SVM,boosting. 课本还涉及到3种算法 ...
机器学习初学者手抄本：数学基础、机器学习经典算法、统计学习方法等
机器学习怎么学?当然是系统地学习了.没有时间这么办呢?利用碎片时间学习!很多人一天要花 2 个小时通勤,通勤路上有很多时间看手机.于是我把一些机器学习的基础知识做成了在线的机器学习手册,只需打开微信收 ...
李航老师《统计学习方法》的代码实现、课件、作业等相关资源的最全汇总
编辑 | Will 出品 | 字节AI 李航:毕业于日本京都大学电气电子工程系,日本东京大学获得计算机科学博士学位.1990年至2001年就职于日本NEC 公司中央研究所,任研究员,2001年至201 ...

统计学习方法--牛顿法和拟牛顿法