牛顿法与拟牛顿法摘记

1 牛顿法

\qquad 除了可以采用最速下降法求解“无约束最优化问题”，另一种常用的方法就是牛顿法。设 $f(\mathbf{x})$ 为二次可微的实函数，在第 $k$ 次迭代点 ${\mathbf{x}}^{(k)}$ 附近用二阶泰勒级数展开式 $\varphi (\mathbf{x})$ 来近似，也就是

$f(\mathbf{x})\approx \varphi (\mathbf{x})=f({\mathbf{x}}^{(k)})+\nabla f({\mathbf{x}}^{(k)}{)}^T(\mathbf{x}-{\mathbf{x}}^{(k)})+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}^{(k)}{)}^T{\nabla }^{2}f({\mathbf{x}}^{(k)})(\mathbf{x}-{\mathbf{x}}^{(k)}{)}^2$

\qquad 根据一阶必要条件，令 ${\varphi }^{\prime }(\mathbf{x})=0$，可得到

${\varphi }^{\prime }(\mathbf{x})=\nabla f({\mathbf{x}}^{(k)})+{\nabla }^{2}f({\mathbf{x}}^{(k)})(\mathbf{x}-{\mathbf{x}}^{(k)})=0$

\qquad 将 $\varphi (\mathbf{x})$ 的驻点作为下一个迭代点 ${\mathbf{x}}^{(k+1)}$ 就得到牛顿法的迭代公式

${\mathbf{x}}^{(k+1)}={\mathbf{x}}^{(k)}-{\nabla }^{2}f({\mathbf{x}}^{(k)}{)}^{-1}\nabla f({\mathbf{x}}^{(k)})$

\qquad 通常，也将牛顿方向定义为：$-{\nabla }^{2}f({\mathbf{x}}^{(k)}{)}^{-1}\nabla f({\mathbf{x}}^{(k)})$
\qquad
$\bullet$　实现代码

import numpy as np
import matplotlib.pyplot as plt
import time
def diffMat(f,x,h=0.001):	
    dx = np.array([h,0])
    dy = np.array([0,h])
    gx = (f(x+dx)-f(x-dx))/(2*h)
    gy = (f(x+dy)-f(x-dy))/(2*h)
    grad = np.asmatrix(np.array([gx,gy])).T       
    g2x = (f(x+dx) + f(x-dx) - 2*f(x))/(h*h)
    g2y = (f(x+dy) + f(x-dy) - 2*f(x))/(h*h)    
    gx2 = (f(x+dy+dx)-f(x+dy-dx))/(2*h)
    gx1 = (f(x-dy+dx)-f(x-dy-dx))/(2*h)    
    gy2 = (f(x+dx+dy)-f(x+dx-dy))/(2*h)
    gy1 = (f(x-dx+dy)-f(x-dx-dy))/(2*h)   
    gxy1 = (gx2 - gx1)/(2*h)
    gxy2 = (gy2 - gy1)/(2*h)
    hesse = np.matrix([[g2x,gxy1],[gxy2,g2y]])
    return grad,hesse
def newton(f,x0):
    k=1
    xk = x0
    while True:
        print('#',k)
        k = k+1        
        grad, hesse = diffMat(f,xk)
        d = np.asarray(hesse.I*grad).flatten()
        print('grad:\n',np.round(grad,4),'\nhesse:\n',np.round(hesse,4))
        xk = xk - d
        print('d:',d/np.linalg.norm(d))
        print('x[{}]:{}\n'.format(k,np.round(xk,4)))
        if np.linalg.norm(grad)<0.00001:
            break
    return xk
if __name__ == "__main__":       
    f = lambda x: (x[0]-1)**4 + (x[1]-0)**2 
    x0 = np.array([0,1],dtype='float')
    time0 = time.process_time()
    minval = newton(f,x0)
    time1 = time.process_time()
    print('x:',np.round(minval,4),'minval:',np.round(f(minval),4))
    print('time: %fs'%(time1-time0))

运行结果：

# 1
grad:
 [[-4.]
 [ 2.]] 
hesse:
 [[12.  0.]
 [ 0.  2.]]
d: [-0.316228    0.94868322]
x[2]:[0.3333 0.    ]
# 2
grad:
 [[-1.1852]
 [ 0.    ]] 
hesse:
 [[5.3333 0.    ]
 [0.     2.    ]]
d: [-1.00000000e+00  6.28995937e-10]
x[3]:[0.5556 0.    ]
（略）
# 12
grad:
 [[-0.]
 [ 0.]] 
hesse:
 [[1.6e-03 0.0e+00]
 [0.0e+00 2.0e+00]]
d: [-1.  0.]
x[13]:[0.9923 0.    ]

x: [0.9923 0.    ] minval: 0.0
time: 0.000000s

阻尼牛顿法

\qquad 牛顿法只是通过在 ${\mathbf{x}}^{(k)}$ 附近用二阶泰勒级数来近似目标函数，牛顿方向 $-{\nabla }^{2}f({\mathbf{x}}^{(k)}{)}^{-1}\nabla f({\mathbf{x}}^{(k)})$ 并不一定是下降方向，可能会使目标函数值增大。阻尼牛顿法，也称修正牛顿法，为了使目标函数值下降，增加了在牛顿方向上的一维搜索过程，其算法步骤为：
\qquad 步骤 $(1):$ 给定初始点 ${\mathbf{x}}^{(1)}\in R^n$，设置允许误差 $\epsilon >0$，令 $k=1$
\qquad 步骤 $(2):$ 计算搜索方向 ${\mathbf{d}}^{(k)}=-{\nabla }^{2}f({\mathbf{x}}^{(k)}{)}^{-1}\nabla f({\mathbf{x}}^{(k)})$
\qquad 步骤 $(3):$ 若 $\Vert \nabla f({\mathbf{x}}^{(k)})\Vert \le \epsilon$，则停止计算；
\qquad\qquad\qquad 否则，沿着 ${\mathbf{d}}^{(k)}$ 进行一维搜索求出 ${\lambda }_k$，使$f({\mathbf{x}}^{(k)}+{\lambda }_k{\mathbf{d}}^{(k)})=\underset{\lambda >0}{\min }f({\mathbf{x}}^{(k)}+\lambda {\mathbf{d}}^{(k)})$
\qquad 步骤 $(4):$ 令 ${\mathbf{x}}^{(k+1)}={\mathbf{x}}^{(k)}+{\lambda }_k{\mathbf{d}}^{(k)}$，且 $k=k+1$，转步骤 $(2)$
\qquad
$\bullet$　实现代码

import numpy as np
import matplotlib.pyplot as plt
import time
def diffMat(f,x,h=0.001):
（略）
def newton1d(f,xk,d,h=0.001):
    dif1 = (f(xk+h*d) - f(xk-h*d))/(2*h)
    dif2 = (f(xk+h*d) + f(xk-h*d) - 2*f(xk))/(h*h)
    deltax = dif1/dif2
    xk1 = xk - deltax * d
    return xk1
def dampednewton(f,x0):
    k=0
    xk = x0
    while True:
        k = k+1
        print('#',k)                
        grad, hesse = diffMat(f,xk)
        d = -np.asarray(hesse.I*grad).flatten()
        if np.linalg.norm(grad)<0.00001:
            break
        print('grad:\n',np.round(grad,4),'\nhesse:\n',np.round(hesse,4))
        print('d:',d/np.linalg.norm(d))
        xk = newton1d(f,xk,d)
        print('x[{}]:{}\n'.format(k,np.round(xk,4)))
    return xk
if __name__ == "__main__":       
    f = lambda x: (x[0]-1)**4 + (x[1]-0)**2 
    x0 = np.array([0,1],dtype='float')
    time0 = time.process_time()
    minval = dampednewton(f,x0)
    time1 = time.process_time()
    print('x:',np.round(minval,4),'minval:',np.round(f(minval),4))
    print('time: %fs'%(time1-time0))

运行结果：

# 1
grad:
 [[-4.]
 [ 2.]] 
hesse:
 [[12.  0.]
 [ 0.  2.]]
d: [ 0.316228   -0.94868322]
x[2]:[0.3333 0.    ]

# 2
grad:
 [[-1.1852]
 [ 0.    ]] 
hesse:
 [[5.3333 0.    ]
 [0.     2.    ]]
d: [ 1.00000000e+00 -1.33281943e-06]
x[3]:[0.5556 0.    ]

# 3
grad:
 [[-0.3512]
 [ 0.    ]] 
hesse:
 [[2.3704 0.    ]
 [0.     2.    ]]
d: [ 1.00000e+00 -3.53622e-12]
x[4]:[0.7037 0.    ]

# 4
grad:
 [[-0.1041]
 [ 0.    ]] 
hesse:
 [[ 1.0535 -0.    ]
 [-0.      2.    ]]
d: [1.00000000e+00 1.06224728e-13]
x[5]:[0.8025 0.    ]
（略）
# 11
grad:
 [[-0.]
 [ 0.]] 
hesse:
 [[0.0036 0.    ]
 [0.     2.    ]]
d: [ 1. -0.]
x[12]:[0.9884 0.    ]

# 12
x: [0.9884 0.    ] minval: 0.0
time: 0.000000s

\qquad

2 拟牛顿法

\qquad 牛顿法的优点是收敛很快，但是牛顿法的每次迭代过程中都需要计算二阶偏导数、求 $\text{Hesse}$ 矩阵的逆矩阵，而目标函数的 $\text{Hesse}$ 矩阵可能是非正定的。

\qquad 拟牛顿法$\text{(Quasi-Newton Method)}$是为了克服牛顿法的缺点而提出了“拟牛顿条件”，其基本思想是用“不包含二阶导数的矩阵”来近似牛顿法中 $\text{Hesse}$ 矩阵的逆矩阵。
\qquad

2.1 拟牛顿条件

\qquad 拟牛顿法是构造近似矩阵 $H_k$ 替代 “$\text{Hesse}$ 矩阵的逆矩阵”件，即：$H_k\approx {\nabla }^{2}f({\mathbf{x}}^{(k)}{)}^{-1}$。拟牛顿条件是构造 $H_k$ 替代 ${\nabla }^{2}f({\mathbf{x}}^{(k)}{)}^{-1}$ 执行牛顿法迭代运算时所需要满足的条件。

\qquad 使用一维搜索时，牛顿法的迭代公式为：

${\mathbf{x}}^{(k+1)}={\mathbf{x}}^{(k)}-{\lambda }_k{\nabla }^{2}f({\mathbf{x}}^{(k)}{)}^{-1}\nabla f({\mathbf{x}}^{(k)})$

\qquad 在第 $k$ 次迭代后得到了点 ${\mathbf{x}}^{(k+1)}$，将目标函数在点 ${\mathbf{x}}^{(k+1)}$ 进行二阶泰勒级数展开：

$f(\mathbf{x})\approx f({\mathbf{x}}^{(k+1)})+\nabla f({\mathbf{x}}^{(k+1)}{)}^T(\mathbf{x}-{\mathbf{x}}^{(k+1)})+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}^{(k+1)}{)}^T{\nabla }^{2}f({\mathbf{x}}^{(k+1)})(\mathbf{x}-{\mathbf{x}}^{(k+1)})$

\qquad 对上式两端取梯度，就得到：

$\nabla f(\mathbf{x})\approx \nabla f({\mathbf{x}}^{(k+1)})+{\nabla }^{2}f({\mathbf{x}}^{(k+1)})(\mathbf{x}-{\mathbf{x}}^{(k+1)})$

\qquad 将 $\mathbf{x}={\mathbf{x}}^{(k)}$ 带入上式：

$\nabla f({\mathbf{x}}^{(k)})\approx \nabla f({\mathbf{x}}^{(k+1)})+{\nabla }^{2}f({\mathbf{x}}^{(k+1)})({\mathbf{x}}^{(k)}-{\mathbf{x}}^{(k+1)})$

\qquad
\qquad 在上式中，令 ${\mathbf{p}}^{(k)}={\mathbf{x}}^{(k+1)}-{\mathbf{x}}^{(k)}$，${\mathbf{q}}^{(k)}=\nabla f({\mathbf{x}}^{(k+1)})-\nabla f({\mathbf{x}}^{(k)})$，那么：

${\mathbf{q}}^{(k)}\approx {\nabla }^{2}f({\mathbf{x}}^{(k+1)}){\mathbf{p}}^{(k)}$

\qquad
\qquad 若矩阵 ${\nabla }^{2}f({\mathbf{x}}^{(k+1)})$ 可逆，则：${\mathbf{p}}^{(k)}\approx {\nabla }^{2}f({\mathbf{x}}^{(k+1)}{)}^{-1}{\mathbf{q}}^{(k)}$

\qquad 若记 $H_{k+1}={\nabla }^{2}f({\mathbf{x}}^{(k+1)}{)}^{-1}$，那么

${\mathbf{p}}^{(k)}=H_{k+1}{\mathbf{q}}^{(k)}$ （拟牛顿条件）

\qquad 拟牛顿法就是为了构造拟牛顿法中的矩阵 $H_{k+1}$ 而提出的一类算法。
\qquad

2.2 秩1校正

\qquad 当 ${\nabla }^{2}f({\mathbf{x}}^{(k+1)}{)}^{-1}$ 是 $n$ 阶对称正定矩阵时，满足拟牛顿条件的近似矩阵 $H_{k+1}$ 也是 $n$ 阶对称正定矩阵。构造这样的近似矩阵，可以通过不断修正 $H_k$ 来构造下一个 $H_{k+1}$，通常令 $H_1=\mathbf{I}$，那么 $H_{k+1}$ 的修正过程为

$H_{k+1}=H_k+\Delta H_k$，　其中 $\Delta H_k$ 为校正矩阵。

\qquad
\qquad 确定 $\Delta H_k$ 的方法之一是，令 $\Delta H_k={\alpha }_k{\mathbf{z}}^{(k)}({\mathbf{z}}^{(k)}{)}^T$，其中 ${\alpha }_k$ 是一个常数，${\mathbf{z}}^{(k)}\in R^n$ 是 $n$ 维列向量。显然，矩阵 $\Delta H_k$ 的秩为 $1$，且 ${\mathbf{z}}^{(k)}$ 应该满足拟牛顿条件 ${\mathbf{p}}^{(k)}=H_{k+1}{\mathbf{q}}^{(k)}$，那么

$\begin{array}{rl}{\mathbf{p}}^{(k)}& =H_{k+1}{\mathbf{q}}^{(k)}\\ & =H_k{\mathbf{q}}^{(k)}+\Delta H_k{\mathbf{q}}^{(k)}\\ & =H_k{\mathbf{q}}^{(k)}+{\alpha }_k{\mathbf{z}}^{(k)}({\mathbf{z}}^{(k)}{)}^T{\mathbf{q}}^{(k)}\end{array}$

\qquad 由于 $({\mathbf{z}}^{(k)}{)}^T{\mathbf{q}}^{(k)}$ 是两个 $n$ 维列向量的内积，其值为常数，可得到：　${\mathbf{z}}^{(k)}=\frac{{\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)}}{{\alpha }_k({\mathbf{z}}^{(k)}{)}^T{\mathbf{q}}^{(k)}}$

\qquad 两边都左乘 $({\mathbf{q}}^{(k)}{)}^T$，可得到：

$({\mathbf{q}}^{(k)}{)}^T{\mathbf{z}}^{(k)}=({\mathbf{z}}^{(k)}{)}^T{\mathbf{q}}^{(k)}=\frac{({\mathbf{q}}^{(k)}{)}^T({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)})}{{\alpha }_k({\mathbf{z}}^{(k)}{)}^T{\mathbf{q}}^{(k)}}$

$⟹({\mathbf{q}}^{(k)}{)}^T({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)})={\alpha }_k(({\mathbf{z}}^{(k)}{)}^T{\mathbf{q}}^{(k)}{)}^2$

${\mathbf{x}}^T\mathbf{y}={\mathbf{y}}^T\mathbf{x},\forall \mathbf{x},\mathbf{y}\in R^n⟹({\mathbf{z}}^{(k)}{)}^T{\mathbf{q}}^{(k)}=({\mathbf{q}}^{(k)}{)}^T{\mathbf{z}}^{(k)}$

\qquad 可得到：

$\begin{array}{rl}\Delta H_k& ={\alpha }_k{\mathbf{z}}^{(k)}({\mathbf{z}}^{(k)}{)}^T\\ & ={\alpha }_k\frac{{\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)}}{{\alpha }_k({\mathbf{z}}^{(k)}{)}^T{\mathbf{q}}^{(k)}}{\left(\frac{{\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)}}{{\alpha }_k({\mathbf{z}}^{(k)}{)}^T{\mathbf{q}}^{(k)}}\right)}^T,({\mathbf{z}}^{(k)}{)}^T{\mathbf{q}}^{(k)}为常量\\ & ={\alpha }_k\frac{\left({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)}\right){\left({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)}\right)}^T}{{\alpha }_k{\left(({\mathbf{z}}^{(k)}{)}^T{\mathbf{q}}^{(k)}\right)}^2}\\ & ={\alpha }_k\frac{\left({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)}\right){\left({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)}\right)}^T}{({\mathbf{q}}^{(k)}{)}^T({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)})}\end{array}$

\qquad
\qquad 因此，$H_{k+1}$ 的修正过程称为秩1校正公式，也就是：

$H_{k+1}=H_k+\Delta H_k=H_k+{\alpha }_k\frac{\left({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)}\right){\left({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)}\right)}^T}{({\mathbf{q}}^{(k)}{)}^T({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)})}$

\qquad
\qquad 基于秩1校正公式的拟牛顿法计算步骤为：

$(1)$ 搜索方向为 ${\mathbf{d}}^{(k)}=-H_k\nabla f({\mathbf{x}}^{(k)}),H_1=\mathbf{I}$

$(2)$ 沿 ${\mathbf{d}}^{(k)}$ 进行一维搜索求出 ${\mathbf{x}}^{(k+1)}$

$(3)$ 求出 ${\mathbf{p}}^{(k)}={\mathbf{x}}^{(k+1)}-{\mathbf{x}}^{(k)}$ 以及 ${\mathbf{q}}^{(k)}=\nabla f({\mathbf{x}}^{(k+1)})-\nabla f({\mathbf{x}}^{(k)})$，计算近似矩阵

$H_{k+1}=H_k+{\alpha }_k\frac{\left({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)}\right){\left({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)}\right)}^T}{({\mathbf{q}}^{(k)}{)}^T({\mathbf{p}}^{(k)}-H_k{\mathbf{q}}^{(k)})}$

\qquad\quad 重新回到第 $(1)$ 步开始下一轮的计算，直到 $\Vert \nabla f({\mathbf{x}}^{(k)})\Vert <\epsilon$。
\qquad
$\bullet$　实现代码

import numpy as np
import matplotlib.pyplot as plt
import time
def diff1(f,x,h=0.001):
    dx = np.array([h,0])
    dy = np.array([0,h])
    gx = (f(x+dx)-f(x-dx))/(2*h)
    gy = (f(x+dy)-f(x-dy))/(2*h)
    return np.array([gx,gy])
def newton1d(f,xk,d,h=0.001):
    dif1 = (f(xk+h*d) - f(xk-h*d))/(2*h)
    dif2 = (f(xk+h*d) + f(xk-h*d) - 2*f(xk))/(h*h)
    deltax = dif1/dif2
    xk1 = xk - deltax * d
    return xk1
def newton(f,x0):
    k=1
    xk = x0
    Hk = np.eye(x0.shape[0])
    while True:
        gk = diff1(f,xk)        
        dk = -np.dot(Hk,gk)
        xk1 = newton1d(f,xk,dk)	
        gk1 = diff1(f,xk1)
        if np.linalg.norm(gk)<0.0001:
            break
        pk = xk1 - xk
        qk = gk1 - gk
        t1 = (pk-Hk.dot(qk)).reshape(x0.shape[0],-1)
        deltHk = t1.dot(t1.T)/qk.dot(pk-Hk.dot(qk))
        Hk = Hk + deltHk
        xk = xk1
        print('#{} min:{}'.format(k,np.round(xk,4)))
        k = k + 1
    return xk
if __name__ == "__main__":       
    f = lambda x: (x[0]-1)**4 + (x[1]-0)**2 
    x0 = np.array([0,1],dtype='float')
    time0 = time.process_time()
    minval = newton(f,x0)
    time1 = time.process_time()
    print('x:',np.round(minval,4),'minval:',np.round(f(minval),4))
    print('time: %fs'%(time1-time0))

运行结果：

#1 min:[0.4 0.8]
#2 min:[ 0.4084 -0.0043]
#3 min:[6.056e-01 1.000e-04]
#4 min:[ 0.7371 -0.    ]
#5 min:[0.8247 0.    ]
#6 min:[ 0.8831 -0.    ]
#7 min:[0.9221 0.    ]
#8 min:[ 0.9481 -0.    ]
#9 min:[0.9654 0.    ]
#10 min:[ 0.9769 -0.    ]
x: [ 0.9769 -0.    ] minval: 0.0
time: 0.000000s

\quad

2.3 DFP算法（变尺度法）

\qquad 变尺度法，也就是$\text{DFP}$法，其校正矩阵为：

$\Delta H_k=\frac{{\mathbf{p}}^{(k)}({\mathbf{p}}^{(k)}{)}^T}{({\mathbf{p}}^{(k)}{)}^T{\mathbf{q}}^{(k)}}-\frac{H_k{\mathbf{q}}^{(k)}({\mathbf{q}}^{(k)}{)}^T{H}_k}{({\mathbf{q}}^{(k)}{)}^T{H}_k{\mathbf{q}}^{(k)}}$

\qquad 近似矩阵 $H_{k+1}$ 的修正过程为：

$H_{k+1}=H_k+\Delta H_k=H_k+\frac{{\mathbf{p}}^{(k)}({\mathbf{p}}^{(k)}{)}^T}{({\mathbf{p}}^{(k)}{)}^T{\mathbf{q}}^{(k)}}-\frac{H_k{\mathbf{q}}^{(k)}({\mathbf{q}}^{(k)}{)}^T{H}_k}{({\mathbf{q}}^{(k)}{)}^T{H}_k{\mathbf{q}}^{(k)}}$

$\text{DFP}$法具有二次终止性，对于正定二次函数的目标函数，经过有限次搜索可以达到极小点。

\qquad
$\text{DFP}$法的计算步骤为：

$(1)$ 搜索方向为 ${\mathbf{d}}^{(k)}=-H_k\nabla f({\mathbf{x}}^{(k)}),H_1=\mathbf{I}$

$(2)$ 沿 ${\mathbf{d}}^{(k)}$ 进行一维搜索求出 ${\mathbf{x}}^{(k+1)}$

$(3)$ 求出 ${\mathbf{p}}^{(k)}={\mathbf{x}}^{(k+1)}-{\mathbf{x}}^{(k)}$ 以及 ${\mathbf{q}}^{(k)}=\nabla f({\mathbf{x}}^{(k+1)})-\nabla f({\mathbf{x}}^{(k)})$，计算近似矩阵

$H_{k+1}=H_k+\frac{{\mathbf{p}}^{(k)}({\mathbf{p}}^{(k)}{)}^T}{({\mathbf{p}}^{(k)}{)}^T{\mathbf{q}}^{(k)}}-\frac{H_k{\mathbf{q}}^{(k)}({\mathbf{q}}^{(k)}{)}^T{H}_k}{({\mathbf{q}}^{(k)}{)}^T{H}_k{\mathbf{q}}^{(k)}}$

\qquad\qquad 重新回到第 $(1)$ 步开始下一轮的计算，直到 $\Vert \nabla f({\mathbf{x}}^{(k)})\Vert <\epsilon$。
\qquad
$\bullet$　实现代码

import numpy as np
import matplotlib.pyplot as plt
import time
def diff1(f,x,h=0.001):
    dx = np.array([h,0])
    dy = np.array([0,h])
    gx = (f(x+dx)-f(x-dx))/(2*h)
    gy = (f(x+dy)-f(x-dy))/(2*h)
    return np.array([gx,gy])
def newton1d(f,xk,d,h=0.001):
    dif1 = (f(xk+h*d) - f(xk-h*d))/(2*h)
    dif2 = (f(xk+h*d) + f(xk-h*d) - 2*f(xk))/(h*h)
    deltax = dif1/dif2
    xk1 = xk - deltax * d
    return xk1
def newton(f,x0):
    k=1
    xk = x0
    Hk = np.eye(x0.shape[0])
    while True:
        gk = diff1(f,xk)        
        dk = -np.dot(Hk,gk)
        xk1 = newton1d(f,xk,dk)	
        gk1 = diff1(f,xk1)
        if np.linalg.norm(gk)<0.0001:
            break
        pk = xk1 - xk
        qk = gk1 - gk
        pkm = (xk1 - xk).reshape(x0.shape[0],-1)
        qkm = (gk1 - gk).reshape(x0.shape[0],-1)
        deltHk = pkm.dot(pkm.T)/pk.dot(qk) - Hk.dot(qkm.dot(qkm.T)).dot(Hk)/qk.dot(Hk.dot(qk))
        Hk = Hk + deltHk
        xk = xk1
        print('#{} min:{}'.format(k,np.round(xk,4)))
        k = k + 1
    return xk
if __name__ == "__main__":       
    f = lambda x: (x[0]-1)**4 + (x[1]-0)**2 
    # f = lambda x: 2*(x[0]-0)**2 + (x[1]-0)**2 - 4*x[0] + 2
    x0 = np.array([0,1],dtype='float')
    time0 = time.process_time()
    minval = newton(f,x0)
    time1 = time.process_time()
    print('x:',np.round(minval,4),'minval:',np.round(f(minval),4))
    print('time: %fs'%(time1-time0))

运行结果：

#1 min:[0.4 0.8]
#2 min:[ 0.4064 -0.0033]
#3 min:[0.6043 0.    ]
#4 min:[ 0.7362 -0.    ]
#5 min:[0.8241 0.    ]
#6 min:[ 0.8828 -0.    ]
#7 min:[0.9218 0.    ]
#8 min:[ 0.9479 -0.    ]
#9 min:[0.9653 0.    ]
#10 min:[ 0.9768 -0.    ]
x: [ 0.9768 -0.    ] minval: 0.0
time: 0.000000s

\qquad

2.4 BFGS公式

\qquad 在讨论拟牛顿条件时，首先构造出 ${\mathbf{q}}^{(k)}\approx {\nabla }^{2}f({\mathbf{x}}^{(k+1)}){\mathbf{p}}^{(k)}$ 的关系，然后记 $H_{k+1}={\nabla }^{2}f({\mathbf{x}}^{(k+1)}{)}^{-1}$，得到了 ${\mathbf{p}}^{(k)}=H_{k+1}{\mathbf{q}}^{(k)}$ 的拟牛顿条件，前文中的秩$1$校正和$\text{DFP}$公式都是基于 ${\mathbf{p}}^{(k)}=H_{k+1}{\mathbf{q}}^{(k)}$ 的拟牛顿条件。
\qquad
\qquad 而$\text{BFGS}$公式直接基于 ${\mathbf{q}}^{(k)}\approx {\nabla }^{2}f({\mathbf{x}}^{(k+1)}){\mathbf{p}}^{(k)}$ 的关系，用不含二阶导数的矩阵 $B_{k+1}$ 近似 ${\nabla }^{2}f({\mathbf{x}}^{(k+1)})$，从而有 ${\mathbf{q}}^{(k)}=B_{k+1}{\mathbf{p}}^{(k)}$，关于 $B_k$ 的修正公式为：

$B_{k+1}=B_k+\frac{{\mathbf{q}}^{(k)}({\mathbf{q}}^{(k)}{)}^T}{({\mathbf{q}}^{(k)}{)}^T{\mathbf{p}}^{(k)}}-\frac{B_k{\mathbf{p}}^{(k)}({\mathbf{p}}^{(k)}{)}^T{B}_k}{({\mathbf{p}}^{(k)}{)}^T{B}_k{\mathbf{p}}^{(k)}}$

\qquad 上式即为关于矩阵 $\text{B}$ 的 $\text{BFGS}$ 修正公式。
　
\qquad 由于 ${\mathbf{q}}^{(k)}=B_{k+1}{\mathbf{p}}^{(k)}$，因此 ${\mathbf{p}}^{(k)}=B_{k+1}^{-1}{\mathbf{q}}^{(k)}$，也就相当于令 $H_{k+1}=B_{k+1}^{-1}$ 时的拟牛顿法。
　
\qquad 利用$\text{Sherman-Morrison}$公式：

$(\mathbf{M}+\mathbf{u}{\mathbf{v}}^T{)}^{-1}={\mathbf{M}}^{-1}-\frac{{\mathbf{M}}^{-1}\mathbf{u}{\mathbf{v}}^T{\mathbf{M}}^{-1}}{1+{\mathbf{v}}^T{\mathbf{M}}^{-1}\mathbf{u}}$

\qquad 可以求出 $B_{k+1}^{-1}$，从而得到关于 $\mathbf{H}$ 的$\text{BFGS}$公式：

$H_{k+1}^{\text{BFGS}}=H_k+\left(1+\frac{({\mathbf{q}}^{(k)}{)}^T{H}_k{\mathbf{q}}^{(k)}}{({\mathbf{p}}^{(k)}{)}^T{\mathbf{q}}^{(k)}}\right)\frac{{\mathbf{p}}^{(k)}({\mathbf{p}}^{(k)}{)}^T}{({\mathbf{p}}^{(k)}{)}^T{\mathbf{q}}^{(k)}}-\frac{{\mathbf{p}}^{(k)}({\mathbf{q}}^{(k)}{)}^T{H}_k+H_k{\mathbf{q}}^{(k)}({\mathbf{p}}^{(k)}{)}^T}{({\mathbf{p}}^{(k)}{)}^T{\mathbf{q}}^{(k)}}$

\qquad
$\text{DFP}$公式和$\text{BFGS}$公式都是由 ${\mathbf{p}}^{(k)}$ 和 $H_k{\mathbf{q}}^{(k)}$ 构成的对称秩$2$校正，因此两个公式的加权组合具有相同的形式，所有这类修正公式的集合可以表示为：

$H_{k+1}^{\varphi }=(1-\varphi )H_{k+1}^{\text{DFP}}+\varphi H_{k+1}^{\text{BFGS}}$

\qquad 这类修正公式的全体，称为$\text{Broyden}$族。当 $\varphi =0$ 时为 $\text{DFP}$公式，当 $\varphi =1$ 时为 $\text{BFGS}$公式。