与公众号同步更新，详细内容及相关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∈Rnmin​f(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−1​gk​(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} Hk​pk​=−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 Hk​pk​=−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−1​gk​，要求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−1​yk​=δ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−1​gk​，由式(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−1​gk​(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))−λgkT​Hk−1​gk​(B.15)

因Hk−1H_k^{-1}Hk−1​正定，故有gkTHk−1gk>0g_k^T H_k^{-1} g_k > 0gkT​Hk−1​gk​>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+1​yk​=δ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+1​yk​=Gk​yk​+Pk​yk​+Qk​yk​(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} Pk​yk​=δk​(B.20)

Qkyk=−Gkyk(B.21)Q_k y_k = -G_k y_k \quad \tag{B.21} Qk​yk​=−Gk​yk​(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​=δkT​yk​δ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​=−ykT​Gk​yk​Gk​yk​ykT​Gk​​(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​+δkT​yk​δk​δkT​​−ykT​Gk​yk​Gk​yk​ykT​Gk​​(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​=−Gk​gk​。

(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)+λk​pk​)=λ≥0min​f(x(k)+λpk​)

(5)置x(k+1)=x(k)+λkpkx^{(k+1)} = x^{(k)} + \lambda_k p_kx(k+1)=x(k)+λk​pk​。

(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​δk​yk​ykT​​−δkT​Bk​δk​Bk​δk​δkT​Bk​​(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_kBk​pk​=−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)+λk​pk​)=λ≥0min​f(x(k)+λpk​)

(5)置x(k+1)=x(k)+λkpkx^{(k+1)} = x^{(k)} + \lambda_k p_kx(k+1)=x(k)+λk​pk​。

(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−δkT​yk​δk​ykT​​)Gk​(I−δkT​yk​δk​ykT​​)T+δkT​yk​δ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 类算法。

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