拟牛顿法(变尺度算法)【python实现_numpy】

1.题目：$$\begin{gathered} minf(x_1,x_2)=2x_1^2+x_2^2-4x_1+2 \\ where{x}^{(0)}=(2,1{)}^T,H_1=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],ϵ=0.01 \end{gathered}$$
2. DFP算法
{\quad} step1：给定${x\text{x}x}^0\in {R\text{R}R}^n，$允许误差$ϵ>0。$

{\quad} step2：令${H\text{H}H}_0:={E\text{E}E}_n，$计算${g\text{g}g}^0:=\nabla f({x\text{x}x}^0)，$令$k:=0。$

{\quad} step3：令$p^k:=-{H\text{H}H}_k{g\text{g}g}^k。$

{\quad} step4：$求{\lambda }_k:f({x\text{x}x}^k+{\lambda }_k{p\text{p}p}^k)=\underset{\lambda >0}{\min }f({x\text{x}x}^k+\lambda {p\text{p}p}^k),令{x\text{x}x}^{x+1}:={x\text{x}x}^k+{\lambda }_k{p\text{p}p}^k,计算{g\text{g}g}^{k+1}:=\nabla f({x\text{x}x}^{k+1})。$

{\quad} step5：若$||{g\text{g}g}^{k+1}||<ϵ,则停止，打印{x\text{x}x}^{k+1};否则，转step6。$

{\quad} step6：若$k=n-1$,则令${x\text{x}x}^0:={x\text{x}x}^{k+1},转step2；否则，转step7。$

{\quad} step7：计算
$$\begin{gathered} \Delta {g\text{g}g}^k={g\text{g}g}^{k+1}-{g\text{g}g}^k,\Delta {x\text{x}x}^k={x\text{x}x}^{k+1}-{x\text{x}x}^k \\ {H\text{H}H}_{k+1}:={H\text{H}H}_k+\Delta {g\text{g}g}^k={H\text{H}H}_k+\frac{\Delta {x\text{x}x}^k(\Delta {x\text{x}x}^k{)}^T}{(\Delta {x\text{x}x}^k{)}^T\Delta {g\text{g}g}^k}-\frac{{H\text{H}H}_k\Delta {g\text{g}g}^k(\Delta {g\text{g}g}^k{)}^T{H\text{H}H}_k}{(\Delta {g\text{g}g}^k{)}^T{H\text{H}H}_k\Delta {g\text{g}g}^k} \end{gathered}$$ $令k:=k+1,转step3$
3. BFGS算法
{\quad} 仅需将DFP的step7中计算$H_{k+1}$的方式改变即可。
$$\begin{gathered} H_{k+1}=H_k+{\beta }_k\frac{\Delta {x\text{x}x}^k(\Delta {x\text{x}x}^k{)}^T}{(\Delta {x\text{x}x}^k{)}^T\Delta {g\text{g}g}^k}-\frac{\Delta {x\text{x}x}^k(\Delta {g\text{g}g}^k{)}^T{H}_k+H_k\Delta {g\text{g}g}^k(\Delta {x\text{x}x}^k{)}^T}{(\Delta {x\text{x}x}^k{)}^T\Delta {g\text{g}g}^k} \\ 其中, \\ {\beta }_k=1+\frac{(\Delta {g\text{g}g}^k{)}^T{H}_k\Delta {g\text{g}g}^k}{(\Delta {x\text{x}x}^k{)}^T\Delta {g\text{g}g}^k} \end{gathered}$$4.python实现
{\quad} 4.1代码

import numpy as np
import math

def get_g(x):
    a = [4*x[0,0]-4,2*x[1,0]]
    b = math.sqrt(a[0]**2+a[1]**2)
    c = np.mat([[a[0]],[a[1]]])
    return [c,b]

def get_x(x,p):
    a = p[0,0]*2-x[0,0]*p[0,0]*2-x[1,0]*p[1,0]
    b = 2*p[0,0]**2+p[1,0]**2
    r = a/b
    return x+r*p

def get_H_DFP(H,x,g):
    return H + (x*x.T)/(x.T*g) - (H*g*g.T*H)/(g.T*H*g)

def get_H_BFGS(H,x,g):
    xg = x.T*g
    b = 1+(g.T*H*g)/xg
    H_ans = H + b[0,0]*(x*x.T)/xg - (x*g.T*H+H*g*x.T)/xg
    return H_ans

def Quasi_Newton(get_H_method, method_name):
    e = 0.01
    n = 2
    H = np.eye(2)
    k = 0
    xk = np.mat([[2],[1]])
    gk = get_g(xk)
    p = (-1)*H*gk[0]
    xk_1 = get_x(xk,p)
    gk_1 = get_g(xk_1)

    while gk_1[1] >=e:
        if k==n-1:
            xk = xk_1
            gk = gk_1
            H = np.eye(2)
            k=0
        else:
            dg = gk_1[0]-gk[0]
            dx = xk_1-xk
            H = get_H_method(H,dx,dg)
            k = k+1
        p = (-1)*H*gk_1[0]
        xk_1 = get_x(xk_1, p)
        gk_1 = get_g(xk_1)
    print(method_name+"_answer:",[xk_1[0,0],xk_1[1,0]])

Quasi_Newton(get_H_DFP, "DFP")
Quasi_Newton(get_H_BFGS, "BFGS")

{\quad} 4.2结果