按秩1法(详见博文《最优化方法Python计算:秩1拟牛顿法》)计算的修正矩阵 Q k + 1 = Q k + E k \boldsymbol{Q}_{k+1}=\boldsymbol{Q}_k+\boldsymbol{E}_k Qk+1=Qk+Ek无法保证其正定性。这时, d k + 1 = − Q k + 1 g k + 1 \boldsymbol{d}_{k+1}=-\boldsymbol{Q}_{k+1}\boldsymbol{g}_{k+1} dk+1=−Qk+1gk+1可能不是 f ( x ) f(\boldsymbol{x}) f(x)在 x k + 1 \boldsymbol{x}_{k+1} xk+1处的下降方向,致使算法失败。要摆脱秩1法的这一窘境,需另辟蹊径。很自然的想法是“秩2”修正:设 Q k \boldsymbol{Q}_k Qk对称正定,令修正矩阵
E k = Δ x k Δ x k ⊤ Δ x k ⊤ Δ g k − Q k Δ g k Δ g k ⊤ Q k Δ g k ⊤ Q k Δ g k \boldsymbol{E}_k=\frac{\Delta\boldsymbol{x}_k\Delta\boldsymbol{x}_k^\top}{\Delta\boldsymbol{x}_k^\top\Delta\boldsymbol{g}_k}-\frac{\boldsymbol{Q}_k\Delta\boldsymbol{g}_k\Delta\boldsymbol{g}_k^\top\boldsymbol{Q}_k}{\Delta\boldsymbol{g}_k^\top\boldsymbol{Q}_k\Delta\boldsymbol{g}_k} Ek=Δxk⊤ΔgkΔxkΔxk⊤−Δgk⊤QkΔgkQkΔgkΔgk⊤Qk
于是得到 Q k + 1 \boldsymbol{Q}_{k+1} Qk+1的秩2修正公式
Q k + 1 = Q k + E k = Q k + Δ x k Δ x k ⊤ Δ x k ⊤ Δ g k − Q k Δ g k Δ g k ⊤ Q k Δ g k ⊤ Q k Δ g k . ( 1 ) \boldsymbol{Q}_{k+1}=\boldsymbol{Q}_k+\boldsymbol{E}_k=\boldsymbol{Q}_k+\frac{\Delta\boldsymbol{x}_k\Delta\boldsymbol{x}_k^\top}{\Delta\boldsymbol{x}_k^\top\Delta\boldsymbol{g}_k}-\frac{\boldsymbol{Q}_k\Delta\boldsymbol{g}_k\Delta\boldsymbol{g}_k^\top\boldsymbol{Q}_k}{\Delta\boldsymbol{g}_k^\top\boldsymbol{Q}_k\Delta\boldsymbol{g}_k}.\quad\quad(1) Qk+1=Qk+Ek=Qk+Δxk⊤ΔgkΔxkΔxk⊤−Δgk⊤QkΔgkQkΔgkΔgk⊤Qk.(1)
利用式(1)作为正定阵 Q k + 1 \boldsymbol{Q}_{k+1} Qk+1修正公式的拟牛顿法是由Broyden,Fletcher,Goldfarb和Shanno在20世纪70年代各自独立提出来的,故常称为BFGS算法。可以证明:
定理1 (1)设 Q k \boldsymbol{Q}_k Qk对称正定,由式(1)确定的 Q k + 1 \boldsymbol{Q}_{k+1} Qk+1正定,当且仅当 Δ x k ⊤ Δ g k > 0 \Delta\boldsymbol{x}_k^\top\Delta\boldsymbol{g}_k>0 Δxk⊤Δgk>0。
(2)设目标函数 f ( x ) f(\boldsymbol{x}) f(x), x ∈ R n \boldsymbol{x}\in\text{R}^n x∈Rn一阶连续可微,且有极小值点 x 0 \boldsymbol{x}_0 x0。则BFGS算法每次迭代均有 Δ x k ⊤ Δ g k > 0 \Delta\boldsymbol{x}_k^\top\Delta\boldsymbol{g}_k>0 Δxk⊤Δgk>0。
BFGS算法是一个改进了的拟牛顿算法,读者可作为练习用Python实现BFGS算法。Python的scipy.optimize为用户提供了BFGS方法,只需要在调用minimize时将’BFGS’传递给method参数即可用BFGS方法计算目标函数的最优解。
例1 用scipy.optimize提供的BFGS方法计算Rosenbrock函数的最优解,给定初始点 x 1 = ( 100 100 ) \boldsymbol{x}_1=\begin{pmatrix}100\\100\end{pmatrix} x1=(100100)。
解:下列代码完成本例计算。
import numpy as np #导入numpy
from scipy.optimize import rosen, minimize #导入rosen, minimize
x=np.array([100,100]) #设置初始点
res=minimize(rosen,x,method='BFGS') #计算最优解
print(res)
运行程序,输出
fun: 1.8831204186846363e-11hess_inv: array([[0.49113161, 0.98272927],[0.98272927, 1.97132641]])jac: array([ 2.16488885e-06, -9.42470479e-07])message: 'Optimization terminated successfully.'nfev: 1488nit: 385njev: 496status: 0success: Truex: array([0.99999566, 0.99999131])
这意味着BFGS方法从初始点 x 1 = ( 100 100 ) \boldsymbol{x}_1=\begin{pmatrix}100\\100\end{pmatrix} x1=(100100)起,迭代385次算得Rosenbrock函数的最优解 x 0 \boldsymbol{x}_0 x0的近似值为 ( 0.99999566 0.99999131 ) \begin{pmatrix}0.99999566\\0.99999131\end{pmatrix} (0.999995660.99999131)。虽然运行效率未必优于秩1算法(见博文《最优化方法Python计算:秩1拟牛顿法》中例),但根据定理1,算法运行的可靠性得到了保证。