拟牛顿法(BFGS)，结合精确一维搜索确定步长

拟牛顿法一般只确定方向，需要结合一维搜索方法来确定步长。

步长的精确确定的可以使用牛顿法等，不精确的常使用Goldstein或者Wolfe。

本代码是拟牛顿法(BFGS)结合精确一维搜索的实现。只需要修改函数和梯度函数即可使用。

如有错误，请大佬们指出。

import numpy as np

def BFGS(x0, fun, dfun, eps, kmax):
    '''
    拟牛顿法（BFGS）公式来更新方向，需要结合步长搜索算法来使用

    输入:
    x0    初始点 (列向量)
    fun   函数名
    dfun  梯度函数名
    eps   精度要求
    kmax  函数的最大迭代次数

    输出:
    f	  函数在极小值 xk 处的目标函数值
    x	  函数采用此方法求得的极小点(列向量)
    k	  求极小点算法迭代次数
    '''
    k = 0
    n = len(x0)
    H0 = np.eye(n)
    Hk = H0
    xk = x0
    gk = dfun(xk)
    while k <= kmax:
        if np.linalg.norm(gk) < eps:
            break
        dk = - Hk @ gk
        lamd = getStepLengthByNewton(xk, dk)
        xk_old = xk
        xk = xk + lamd * dk
        gk_old = gk
        gk = dfun(xk)
        sk = xk - xk_old
        qk = gk - gk_old
        if np.matmul(sk.T, qk) > 0:
            yk = Hk @ qk
            w = sk / (sk.T @ qk) - yk / (qk.T @ yk)
            Hk = Hk + (sk @ sk.T) / (sk.T @ qk) - (yk @ yk.T) / (qk.T @ yk) + qk.T @ yk * w @ w.T
        k += 1
    
    f = fun(xk)
    x = xk
    display(x, f, k)
    return x, f, k
   
def getStepLengthByNewton(x, d):
    '''
    采用牛顿法，精确线性搜索，确定移动步长

    输入：
    x:     拟牛顿法的迭代点（列向量）
    d:     拟牛顿法确定好的方向(列向量)

    输出：
    lamd:  精确的移动步长(标量)
    '''
    lamd = 1.0         #初始猜测值
    e0 = 1e-6          #退出搜索循环的条件
    delta_a = 1e-6     #对lamd作差分的小量
    while(1):
        new_lamd = x + lamd * d
        new_lame_1 = x + (lamd - delta_a) * d
        new_lamd_2 = x + (lamd + delta_a) * d
        diff_lamd = (fun(new_lamd_2) - fun(new_lame_1)) / (2.0 * delta_a)
        if np.abs(diff_lamd) < e0:
            break
        ddiff_lamd = (fun(new_lamd_2) + fun(new_lame_1) - 2.0 * fun(new_lamd)) / (delta_a * delta_a)
        lamd = lamd - diff_lamd / ddiff_lamd
    return lamd

def display(x, f, k): #显示答案
    x = np.around(x, 6)
    print("x是:\n{}\n函数值是{:.6f}\n迭代步数为{}".format(x, f, k))

def fun(x):
    return x[0,0]**2 + x[1,0]**2 - x[0,0] * x[1,0] - 10 * x[0,0] - 4 * x[1,0] + 60

def dfun(x):
    g = np.array([[2 * x[0,0] -  x[1,0] -10], [2 * x[1,0]] - x[0,0] - 4])
    return g

x0 = np.array([[0], [0]])
BFGS(x0, fun, dfun, 0.01, 100)