与公众号同步更新,详细内容及相关ipynb文件在公众号中,公众号:AI入门小白
文章目录
牛顿法拟牛顿法的思路DFP (Davidon-Fletcher- Powell) 算法(DFP algorithm)BFGS (Broyden-Fletcher-Goldfarl-Shanno) 算法(BFGS algorithm)Broyden 类算法(Broyden's algorithm)牛顿法(Newton method) 和拟牛顿法(quasi-Newton method) 也是求解无约束最优化问题的常用方法,有收敛速度快的优点。牛顿法是迭代算法,每一步需要求解目标函数的黑塞矩阵的逆矩阵,计算比较复杂。拟牛顿法通过正定矩阵近似黑塞矩阵的逆矩阵或黑塞矩阵,简化了这一计算过程。
黑塞矩阵是一个多元函数的二阶偏导数构成的方阵,描述了函数的局部曲率。
牛顿法
考虑无约束最优化问题
minx∈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 类算法。
数据来源:统计学习方法(第二版)