数值优化 —— 信赖域算法（DogLeg算法）（python实现）

相关：

数值优化 —— 信赖域算法

https://zhuanlan.zhihu.com/p/364296114

根据https://zhuanlan.zhihu.com/p/364296114可以知道DogLeg算法的置信域算法的步骤如下：

不过，需要注意的是这个算法步骤可能存在描述的错误，原因https://zhuanlan.zhihu.com/p/364296114在随后给出的代码实现和这个算法步骤是有一定不同的地方的。

给出根据这边算法步骤（修正后的）所给出的DogLeg算法的代码实现：

import numpy as np
import matplotlib.pyplot as plt

def function(x1, x2):
    """定义函数的表达式

    Args:
        x1 : 变量x1
        x2 : 变量x2

    Returns:
        函数表达式
    """
    return 100*(x2-x1**2)**2+(1-x1)**2

def gradient_function(x1, x2):
    """定义函数的一阶梯度

    Args:
        x1 : 变量x1
        x2 : 变量x2

    Returns:
        函数的一阶梯度
    """
    g=[[-400*(x1*x2-x1**3)+2*x1-2], [200*(x2-x1**2)]]
    g = np.array(g)
    return g

def Hessian_function(x1, x2):
    """定义函数二阶Hessian矩阵

    Args:
        x1 : 变量x1
        x2 : 变量x2

    Returns:
        函数二阶Hessian矩阵
    """
    H = [[-400*(x2-3*x1**2)+2, -400*x1], [-400*x1, 200]]
    H = np.array(H)
    return H

#定义m_k函数
def mk_function(x1, x2, p):
    """近似函数m_k(p)

    Args:
        x1 : 变量x1
        x2 : 变量x2
        p : 下降的试探步

    Returns:
        mk : 近似函数m_k(p)
    """
    p = np.array(p)
    fk = function(x1, x2)
    gk = gradient_function(x1, x2)
    Bk = Hessian_function(x1, x2)
    mk = fk + np.dot(gk.T, p) + 0.5 * np.dot(np.dot(p.T, Bk), p) 
    return mk

def Dogleg_Method(x1,x2,delta):
    """Dogleg Method实现

    Args:
        x1 : 变量x1
        x2 : 变量x2
        delta : 信赖域半径
    Returns:
        s_k : 试探步
    """

    g = gradient_function(x1,x2)
    B = Hessian_function(x1,x2)
    g = g.astype(np.float32)
    B = B.astype(np.float32)

    inv_B = np.linalg.inv(B)
    PB = np.dot(-inv_B, g)
    PU = -(np.dot(g.T, g)/(np.dot(g.T, B).dot(g))) * g
    PB_U = PB-PU
    PB_norm = np.linalg.norm(PB)
    PU_norm = np.linalg.norm(PU)
    PB_U_norm = np.linalg.norm(PB_U)
    #判断τ
    if PB_norm <= delta:
        tao = 2
    elif PU_norm >= delta:
        tao = delta/PU_norm
    else:
        factor = np.dot(PU.T, PB_U) * np.dot(PU.T, PB_U)
        tao = -2 * np.dot(PU.T, PB_U) + 2 * np.math.sqrt(factor - PB_U_norm * PB_U_norm * (PU_norm * PU_norm - delta * delta))
        tao = tao / (2 * PB_U_norm * PB_U_norm) + 1

    #确定试探步
    if 0<=tao<=1:
        s_k = tao*PU
    elif 1<tao<=2:
        s_k = PU+(tao-1)*(PB-PU)
    print("tao: ", tao, "s_k: ", s_k)
    return s_k

def TrustRegion(x1, x2, delta_max):
    """信赖域算法

    Args:
        x1 : 初始值x1
        x2 : 初始值x2
        delta_max : 最大信赖域半径

    Returns:
        x1 : 优化后x1
        x2 : 优化后x2
        
    """
    delta = delta_max
    k = 0
    #计算初始的函数梯度范数
    #终止判别条件中的epsilon
    epsilon = 1e-9
    maxk = 1000
    x1_log=[]
    x2_log=[]
    x1_log.append(x1)
    x2_log.append(x2)
    
    #设置终止判断，判断函数fun的梯度的范数是不是比epsilon小
    while True:
        g_norm = np.linalg.norm(gradient_function(x1, x2))
        if g_norm < epsilon:
            break
        if k > maxk:
            break
        #利用DogLeg_Method求解子问题迭代步长sk
        sk = Dogleg_Method(x1, x2, delta)
        x1_new = x1 + sk[0][0]
        x2_new = x2 + sk[1][0]
        fun_k = function(x1, x2)
        fun_new = function(x1_new, x2_new)
        #计算下降比
        r = (fun_k - fun_new) / (mk_function(x1, x2, [[0],[0]]) - mk_function(x1, x2, sk))
        
        if r < 0.25:
            delta = delta / 4
        elif r > 0.75 and np.linalg.norm(sk) == delta:
            delta = np.min((2 * delta, delta_max))
        else:
            pass

        if r <= 0:
            pass
        else:
            x1 = x1_new
            x2 = x2_new
            k = k + 1
            x1_log.append(x1)
            x2_log.append(x2)

    return x1_log, x2_log


if __name__ =='__main__':
    x1=0.5 
    x2=1.5
    delta_max = 20
    x1_log, x2_log = TrustRegion(x1, x2, delta_max)
    print('x1迭代结果: ', x1_log, '\nx2迭代结果: ', x2_log)

    plt.figure()
    plt.title('x1_convergence')
    plt.plot(x1_log)
    plt.savefig('x1.png')

    plt.figure()
    plt.clf
    plt.title('x2_convergence')
    plt.plot(x2_log)
    plt.savefig('x2.png')