深入浅出最优化(4) 拟牛顿法

1 拟牛顿法的数学基础

对于牛顿法，我们保留其快速收敛性，同时克服牛顿法黑森矩阵需要正定的问题以及避免计算黑森矩阵以减少计算量，我们提出了拟牛顿法。

假定当前点为 $x_{k+1}$ ，若我们用已得到的 $x_k,x_{k+1}$ 及其一阶导数信息 $g_k,g_{k+1}$ ，构造一个正定矩阵 $B_{k+1}$ 作为 $G_{k+1}$ 的近似。这样根据牛顿法下降方向的产生公式 $\nabla^2f(x_{k+1})d_{k+1}=-\nabla f(x_{k+1})$ ，近似可以得到 $B_{k+1}d_{k+1}=-g_{k+1}$ 。设 $H_{k+1}=B_{k+1}^{-1}$ ，则可以得到 $d_{k+1}=-H_{k+1}g_{k+1}$ 。

$B_{k+1}$ 的得出是很理想的，我们希望它能根据已得到的 $x_k,x_{k+1}$ 及其一阶导数信息 $g_k,g_{k+1}$ 用纯迭代的方法得到，而避免黑森矩阵的计算。但想要求得这样的 $B_{k+1}$ ，首先要知道什么样的 $B_{k+1}$ 可以作为 $G_{k+1}$ 的近似矩阵。这个限制条件是保证 $x_k$ 与 $x_{k+1}$ 满足梯度的割线方程。我们知道目标函数梯度函数上两点间的割线方程为 $\nabla f(x_k)-\nabla f(x_{k+1})=G_{k+1}(x_k-x_{k+1})$ ，那么作为近似的 $B_{k+1}$ 也必须使得这个方程成立，即 $\nabla f(x_k)-\nabla f(x_{k+1})=B_{k+1}(x_k-x_{k+1})$ 。

记 $y_k=\nabla f(x_k)-\nabla f(x_{k+1})$ ， $s_k=x_k-x_{k+1}$ ，则有拟牛顿方程 $B_{k+1}s_k=y_k$ 和第二拟牛顿方程 $s_k=H_{k+1}y_k$ 。

2 $B_{k+1}$ 的计算方法

2.1 BFGS修正公式

给定初始正定矩阵 $B_0$ ，之后有 $B_{k+1}=B_k+\Delta_k$ ，假设 $\Delta_k$ 是秩为2的矩阵，则根据矩阵分解有 $\Delta_k=a_ku_ku_k^T+b_kv_kv_k^T$ ，其中 $a_k,b_k$ 是任意常数， $u_k,v_k$ 是任意向量，则 $B_{k+1}=B_k+a_ku_ku_k^T+b_ku_ku_k^T$ ，代入拟牛顿方程得到 $a_ku_ku_k^Ts_k+b_ku_ku_k^Ts_k=y_k-B_ks_k$ 。不妨设 $u_k=y_k,v_k=B_ks_k$ ，则 $a_k=\frac{1}{u_k^Ts_k}=\frac{1}{y_k^Ts_k},b_k=\frac{1}{v_k^Ts_k}=-\frac{1}{s_k^TB_ks_k}$ 。得到BFGS修正公式：

$B_{k+1}=B_k-\frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}+\frac{y_ky_k^T}{y_k^Ts_k}$

这个公式可以保证每步产生的 $B_k$ 的正定。（证明见附录1）

2.2 DFP修正公式

$H_k=B_k^{-1}$ ，又根据第二拟牛顿方程 $s_k=H_{k+1}y_k$ ，和BFGS公式使用同样方法推导可得： $H_{k+1}=H_k-\frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\frac{s_ks_k^T}{s_k^Ty_k}$

3 拟牛顿法的步骤

给定初始点 $x_0\in R^n$ ， $B_0$ 或 $H_0$ ，精度 $\epsilon>0$ ，令 $k=0$
若满足终止准则，则得解 $x_k$ ，算法终止
解线性方程组 $B_k d_k+\nabla f(x_k)=0$ ，得解 $d_k$ ，或计算 $d_k=-H_kg_k$
由线性搜索计算步长 $\alpha_k$ ，令 $x_{k+1}=x_k+\alpha_kd_k$
用修正公式得到 $B_{k+1}$ 或 $H_{k+1}$ ，使其满足拟牛顿方程， $k=k+1$ ，转2

4 性能评估

收敛性：因为保证了 $B_k$ 正定，则根据古典牛顿法的条件，每一步均会产生下降方向，所以满足收敛性
收敛速度：若在最优解处梯度为0、黑森矩阵正定，则当且仅当 $\displaystyle\lim_{k→∞}\frac{||(B_k-\nabla^2f(x^*))||}{||d_k||}=0$ 时，算法超线性收敛。（证明见附录2）
二次终止性：显然满足，而且因为无需步长搜索，每步步长为1，对于二次函数可以一次取到最优解
计算量：小，只需计算两点梯度即可
储存空间：大，对于n元目标函数，需要储存一个大小为n*n的矩阵

5 实战测试

对于本文集的第一篇文章[深入浅出最优化(1) 最优化问题概念与基本知识]中提出的最小二乘问题， $x_1,x_2,x_3,x_4$ 的初值均在 $[-2,2]$ 的范围内随机生成，总共生成100组起点。统计迭代成功（在1000步内得到最优解且单次步长搜索迭代次数不超过1000次）的样本的平均迭代步数、平均迭代时间和得到的最优解及残差平方和最小值。

BFGS：

平均迭代步数	平均迭代时间	最优解	残差平方和最小值
18.1	0.28s	$x_1=0.1891~x_2=0.2820~x_3=0.1463~ x_4=0.1744$	$1.5939\times10^{-4}$

DFP：

平均迭代步数	平均迭代时间	最优解	残差平方和最小值
102.05	1.86s	$x_1=0.1689~x_2=1.0273~x_3=0.3685~ x_4=0.4585$	$2.6031\times10^{-4}$

代码实现

本博客所有代码在https://github.com/HarmoniaLeo/optimization-in-a-nutshell开源，如果帮助到你，请点个star，谢谢这对我真的很重要！

你可以在上面的GitHub链接或本文集的第一篇文章深入浅出最优化(1) 最优化问题概念与基本知识中找到Function.py和lagb.py

BFGS：

import numpy as np
from Function import Function   #定义法求导工具
from lagb import *  #线性代数工具库
from scipy import linalg

n=2 #x的长度

def myFunc(x):  #x是一个包含所有参数的列表
    return x[0]**2 + 2*x[1]**2 + 2*x[0] - 6*x[1] +1 #目标方程

k=0
x=np.zeros(n)   #初始值点
e=0.001
sigma=0.4
rho=0.55
beta=1
tar=Function(myFunc)
B=np.array(np.eye(n))
while tar.norm(x)>e:
    a=1
    d=-np.linalg.solve(B,tar.grad(x))
    if not (tar.value(x+a*d)<=tar.value(x)+rho*a*dot(turn(tar.grad(x)),d) and \
        dot(turn(tar.grad(x+a*d)),d)>=sigma*dot(turn(tar.grad(x)),d)):
        a=beta
        while tar.value(x+a*d)>tar.value(x)+rho*a*dot(turn(tar.grad(x)),d):
            a*=rho
        while dot(turn(tar.grad(x+a*d)),d)<sigma*dot(turn(tar.grad(x)),d):
            a1=a/rho
            da=a1-a
            while tar.value(x+(a+da)*d)>tar.value(x)+rho*(a+da)*dot(turn(tar.grad(x)),d):
                da*=rho
            a+=da
    x1=x+a*d
    if tar.norm(x1)<e:
        print(k+1)
        x=x1
        break
    s=x-x1
    y=tar.grad(x)-tar.grad(x1)
    B=B-muldot(B,s,turn(s),B)/muldot(turn(s),B,s)+dot(y,turn(y))/dot(turn(y),s)
    x=x1
    k+=1
    print(k)
print(x)

DFP：

import numpy as np
from Function import Function   #定义法求导工具
from lagb import *  #线性代数工具库
from scipy import linalg

n=2 #x的长度

def myFunc(x):  #x是一个包含所有参数的列表
    return x[0]**2 + 2*x[1]**2 + 2*x[0] - 6*x[1] +1 #目标方程

k=0
x=np.zeros(n)   #初始值点
e=0.001
sigma=0.4
rho=0.55
beta=1
tar=Function(myFunc)
H=np.array(np.eye(n))
while tar.norm(x)>e:
    a=1
    d=-dot(H,tar.grad(x))
    if not (tar.value(x+a*d)<=tar.value(x)+rho*a*dot(turn(tar.grad(x)),d) and \
        dot(turn(tar.grad(x+a*d)),d)>=sigma*dot(turn(tar.grad(x)),d)):
        a=beta
        while tar.value(x+a*d)>tar.value(x)+rho*a*dot(turn(tar.grad(x)),d):
            a*=rho
        while dot(turn(tar.grad(x+a*d)),d)<sigma*dot(turn(tar.grad(x)),d):
            a1=a/rho
            da=a1-a
            while tar.value(x+(a+da)*d)>tar.value(x)+rho*(a+da)*dot(turn(tar.grad(x)),d):
                da*=rho
            a+=da
    x1=x+a*d
    if tar.norm(x1)<e:
        print(k+1)
        x=x1
        break
    s=x-x1
    y=tar.grad(x)-tar.grad(x1)
    H=H-muldot(H,y,turn(y),H)/muldot(turn(y),H,y)+dot(s,turn(s))/dot(turn(s),y)
    x=x1
    k+=1
    print(k)
print(x)