最优化方法Python计算：秩1拟牛顿法

对二次型函数$f(\mathbf{x})=\frac{1}{2}{\mathbf{x}}^{\top }\mathbf{Hx}-{\mathbf{x}}^{\top }\mathbf{b}$，其中$\mathbf{H}$对称正定。$f(\mathbf{x})$的梯度$\mathbf{g}(\mathbf{x})=\nabla f(\mathbf{x})=\mathbf{Hx}-\mathbf{b}$，Hesse阵${\nabla }^{2}f(\mathbf{x})=\mathbf{H}$。取定初始点${\mathbf{x}}_1$，牛顿法的迭代式为${\mathbf{x}}_{k+1}={\mathbf{x}}_k-{\mathbf{H}}^{-1}{\mathbf{g}}_k$，其中${\mathbf{g}}_k=\mathbf{g}({\mathbf{x}}_k)$。由${\mathbf{g}}_k={\mathbf{Hx}}_k-\mathbf{b}$，有
$${\mathbf{g}}_{k+1}-{\mathbf{g}}_k={\mathbf{Hx}}_{k+1}-{\mathbf{Hx}}_k=\mathbf{H}({\mathbf{x}}_{k+1}-{\mathbf{x}}_k)$$
或
$${\mathbf{H}}^{-1}({\mathbf{g}}_{k+1}-{\mathbf{g}}_k)={\mathbf{x}}_{k+1}-{\mathbf{x}}_k.$$
若记
$$\{\begin{array}{l}\Delta {\mathbf{x}}_k={\mathbf{x}}_{k+1}-{\mathbf{x}}_k\\ \Delta {\mathbf{g}}_k={\mathbf{g}}_{k+1}-{\mathbf{g}}_k\end{array}$$
则称
$$\Delta {\mathbf{g}}_k=\mathbf{H}\Delta \mathbf{x}$$
或
$${\mathbf{H}}^{-1}\Delta {\mathbf{g}}_k=\Delta {\mathbf{x}}_k$$
为牛顿方程。
牛顿法能正确工作的一个重要前提是初始点${\mathbf{x}}_1$需与最优解点${\mathbf{x}}_0$充分接近。因为这时迭代式${\mathbf{x}}_{k+1}={\mathbf{x}}_k-{\mathbf{H}}_k^{-1}{\mathbf{g}}_k$可能使得$f({\mathbf{x}}_{k+1})$未必小于$f({\mathbf{x}}_k)$，$k=1,2,\cdots$。这可以通过选取合理的步长${\alpha }_k$加以改善，即把迭代式改为
$${\mathbf{x}}_{k+1}={\mathbf{x}}_k-{\alpha }_k{\mathbf{H}}_k^{-1}{\mathbf{g}}_k$$
其中${\mathbf{g}}_k=\nabla f({\mathbf{x}}_k)$，${\mathbf{H}}_k={\nabla }^{2}f({\mathbf{x}}_k)$。选取合适的${\alpha }_k$，譬如，${\alpha }_k$为${\varphi }_k(\alpha )=f({\mathbf{x}}_k-\alpha {\mathbf{H}}_k^{-1}{\mathbf{g}}_k)$的最小值点，使得$f({\mathbf{x}}_{k+1})<f({\mathbf{x}}_k)$，$k=1,2,\cdots$。另一方面，有如下事实
定理1 设函数$f(\mathbf{x}),\mathbf{x}\in {\text{R}}^n$一阶连续可微，$\mathbf{Q}\in {\text{R}}^{n\times n}$，${\mathbf{Q}}^{\top }=\mathbf{Q}$且正定。对${\mathbf{x}}_k\in {\text{R}}^n$，若${\mathbf{g}}_k=\nabla f({\mathbf{x}}_k)\ne \mathbf{o}$，记${\alpha }_k=\arg \underset{\alpha >0}{\min }f({\mathbf{x}}_k-\alpha \mathbf{Q}{\mathbf{g}}_k)$，${\mathbf{x}}_{k+1}={\mathbf{x}}_k-{\alpha }_k\mathbf{Q}{\mathbf{g}}_k$。则
$$f({\mathbf{x}}_{k+1})<f({\mathbf{x}}_k).$$
对于非二次型函数$f(\mathbf{x})$而言，其Hesse矩阵$\mathbf{H}={\nabla }^{2}f(\mathbf{x})$，将依赖于$\mathbf{x}$。当$n$很大时，在牛顿算法的迭代计算中更新${\mathbf{H}}_k={\nabla }^{2}f({\mathbf{x}}_k)$会消耗大量资源。我们的目标是对迭代点${\mathbf{x}}_k$寻求${\mathbf{H}}_k$的一个近似对称正定矩阵${\mathbf{Q}}_k$，使得在迭代中更新${\mathbf{Q}}_k$比更新${\mathbf{H}}_k$更高效，且矩阵${\mathbf{Q}}_k$应满足牛顿方程${\mathbf{Q}}_k^{-1}\Delta {\mathbf{g}}_{k-1}\approx \Delta {\mathbf{x}}_{k-1}$或${\mathbf{g}}_{k-1}\approx {\mathbf{Q}}_k\Delta {\mathbf{x}}_{k-1}$。其中，$\Delta {\mathbf{x}}_{k-1}={\mathbf{x}}_k-{\mathbf{x}}_{k-1}$，$\Delta {\mathbf{g}}_{k-1}={\mathbf{g}}_k-{\mathbf{g}}_{k-1}$。这样改造牛顿法所得算法称为拟牛顿法。
实现拟牛顿法最简单的是秩1法。
定理2 设目标函数$f(\mathbf{x})$，$\mathbf{x}\in \in {\text{R}}^n$有最小值点${\mathbf{x}}_0$且一阶连续可微。$k=1$时，令对称正定阵${\mathbf{Q}}_1=\mathbf{I}$。对$k>1$，令${\mathbf{Q}}_k^{-1}={\mathbf{Q}}_{k-1}^{-1}+{\mathbf{E}}_{k-1}$，其中，
$${\mathbf{E}}_{k-1}=\frac{(\Delta {\mathbf{x}}_{k-1}-{\mathbf{Q}}_{k-1}^{-1}\Delta {\mathbf{g}}_{k-1})(\Delta {\mathbf{x}}_{k-1}-{\mathbf{Q}}_{k-1}^{-1}\Delta {\mathbf{g}}_{k-1}{)}^{\top }}{(\Delta {\mathbf{x}}_{k-1}-{\mathbf{Q}}_{k-1}^{-1}\Delta {\mathbf{g}}_{k-1}{)}^{\top }\Delta {\mathbf{g}}_{k-1}}$$
则迭代式
$${\mathbf{Q}}_k^{-1}={\mathbf{Q}}_{k-1}^{-1}+{\mathbf{E}}_{k-1}$$
满足牛顿方程${\mathbf{Q}}_k^{-1}\Delta {\mathbf{g}}_{k-1}\approx \Delta {\mathbf{x}}_{k-1}$。
利用定理1和定理2。得到秩1法搜索点序列的迭代式
$${\mathbf{x}}_{k+1}={\mathbf{x}}_k-{\alpha }_k{\mathbf{Q}}_k^{-1}{\mathbf{g}}_k$$
其中，${\alpha }_k=\arg \underset{\alpha >0}{\min }f({\mathbf{x}}_k-\alpha \mathbf{Q}{\mathbf{g}}_k)$。下列Python函数实现秩1算法。

import numpy as np
from scipy.optimize import minimize_scalar,OptimizeResult
def rank1(f, x1, gtol, **options):
    n=x1.size
    xk=x1
    gk=grad(f,xk)
    Qk=np.eye(n)
    dk=-np.matmul(Qk,gk)
    phi=lambda a: f(xk+a*dk)
    k=1
    while np.linalg.norm(gk)>gtol:
        ak=minimize_scalar(phi).x
        dx=ak*dk
        xk=xk+dx
        dg=grad(f,xk)-gk
        vk=dx-np.matmul(Qk,dg)
        Ek=np.matmul(vk.reshape(n,1),vk.reshape(1,n))/np.dot(vk,dg)
        Qk+=Ek
        gk+=dg
        dk=-np.matmul(Qk,gk)
        k+=1
    bestx=xk
    besty=f(bestx)
    return OptimizeResult(fun=besty, x=bestx, nit=k)

程序中第3~24行定义的rank1函数实现对称秩1算法。参数f表示目标函数$f(\mathbf{x})$，x1表示初始点${\mathbf{x}}_1$，eps表示容错误差$\epsilon$，options实现minimize与本函数的信息交换机制。
第4~10行执行初始化操作：第4行读取自变量维数n。第5行将表示迭代点的xk初始化为x1。第6行调用计算梯度的函数grad（见博文《最优化方法Python计算：n元函数梯度与Hesse阵的数值计算》）计算目标函数的梯度
$${\mathbf{g}}_k=\nabla f({\mathbf{x}}_k)$$
赋予变量gk。第7行调用numpy的eye函数将正定阵${\mathbf{Q}}_k^{-1}$初始化为单位阵${\mathbf{I}}_n$赋予Qk。第8行调用numpy的matmul函数计算
$${\mathbf{d}}_k=-{\mathbf{Q}}_k^{-1}{\mathbf{g}}_k$$
赋予表示搜索刚想的变量dk。第9行用{\bf{lambda}}运算符定义一元函数
$$\varphi (\alpha )=f({\mathbf{x}}_k+\alpha {\mathbf{d}}_k)$$
为phi。
第11~21行执行迭代的while循环中，第12行调用scipy.optimize的minimize_scalar函数（第2行导入），计算
$${\alpha }_k=\arg \underset{\alpha >0}{\min }\varphi (\alpha )$$
赋予ak。第13行计算
$$\Delta \mathbf{x}={\mathbf{x}}_{k+1}-{\mathbf{x}}_k={\mathbf{x}}_k+{\alpha }_k{\mathbf{d}}_k-{\mathbf{x}}_k={\alpha }_k{\mathbf{d}}_k$$
赋予dx。第14行计算
$${\mathbf{x}}_{k+1}={\mathbf{x}}_k+\alpha {\mathbf{d}}_k={\mathbf{x}}_k+\Delta {\mathbf{x}}_k$$
更新迭代点xk。第15行计算
$$\Delta {\mathbf{g}}_k={\mathbf{g}}_{k+1}-{\mathbf{g}}_k$$
dg。第16行调用numpy的matmul函数计算向量
$${\mathbf{v}}_k=\Delta {\mathbf{x}}_k-{\mathbf{Q}}_k\Delta {\mathbf{g}}_k$$
赋予vk。第17行计算矩阵
$${\mathbf{E}}_k=\frac{{\mathbf{v}}_k^{\top }{\mathbf{v}}_k}{{\mathbf{v}}_k\Delta {\mathbf{g}}^{\top }}$$
赋予Ek。第18行计算
$${\mathbf{Q}}_{k+1}^{-1}={\mathbf{Q}}_k^{-1}+{\mathbf{E}}_k$$
更新Qk。第19行计算
$${\mathbf{g}}_{k+1}=\Delta {\mathbf{g}}_k-{\mathbf{g}}_k$$
更新gk。第20行计算
$${\mathbf{Q}}_{k+1}^{-1}{\mathbf{g}}_{k+1}$$
更新dk。\par
第22~24用$f({\mathbf{x}}_k)$，${\mathbf{x}}_k$以及$k$构造OptimizeResult（第2行导入）对象并返回。以下例子说明函数rank1的应用。
例1 设初始点${\mathbf{x}}_1=\left(\begin{array}{c}100\\ 100\end{array}\right)$，用rank1函数计算Rosenbrock函数的最优解，给定容错误差$\epsilon =10^{-8}$。
解：下列代码完成本例计算。

import numpy as np                                      #导入numpy
from scipy.optimize import rosen, minimize              #导入rosen, minimize
x=np.array([100,100])                                   #设置初始点
res=minimize(rosen,x,method=rank1,options={'gtol':1e-8})#计算最优解
print(res)

第3行设置初始点${\mathbf{x}}_1=\left(\begin{array}{c}100\\ 100\end{array}\right)$。第4行调用scipy.optimize的minimize函数（第2行导入），传递rank1给参数method，计算rosen（第2行导入）表示的Rosenbrock函数的最优解。运行程序，输出

 fun: 4.889083789243004e-26
 nit: 49
   x: array([1., 1.])

即以${\mathbf{x}}_1=\left(\begin{array}{c}100\\ 100\end{array}\right)$为初始点，$\epsilon =10^{-8}$为容错误差，实现对称秩1算法的rank1函数迭代49次，算得Rosenbrock函数的最优解${\mathbf{x}}_0=\left(\begin{array}{c}1\\ 1\end{array}\right)$。需提请注意的是，对牛顿法而言，以${\mathbf{x}}_1=\left(\begin{array}{c}100\\ 100\end{array}\right)$为初始点是得不到最优解${\mathbf{x}}_0$的，读者可自行验证。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！