信赖域算法-The Dogleg Method（含例题及Python实现）

---

前言

最近在上王晓老师的最优化算法课程。课程偏硬核。记录作业中信赖域算法中狗腿（The Dogleg）方法的实现。

一、What is The Dogleg Method？

信赖域算法原理

把$minf(x,y)$转化为一维求步长$s_k$问题。其中$s_k=\underset{p}{\min }{m}_k(p)$。通过不断求步长，来进一步迭代出新的x,y。其中函数$m_k(p)$是将函数$f(x,y)$在点$f(x_k,y_k)$泰勒展开后，得到的一个近似函数，也就是其展开后的前三项，如下：
$$\underset{p\in R^n}{\min }{m}_k(p)=f_k+g_k^{T}p+\frac{1}{2}{p}^T{B}_{k}p,s.t.||p||\le {\Delta }_k，{\Delta }_k>0$$

其中${\Delta }_k>0$是信赖域半径（trust-region radius）,$s_k$是上述方程的解，成为试探步（trial step）。
$g_k=\nabla f(x_k)$是一阶gradient。$B_k={\nabla }^{2}f(x_k)$是二阶Hession。
且$$f(x_k+p)=f_k+g_k^{T}p+\frac{1}{2}{p}^T{\nabla }^{2}f(x_k+tp)p,t\in (0,1)$$
可以看到$f(x_k+p)$与$m_k$之间有一个近似误差$o(||p|{|}^2)$，当$p$很小时候，两者误差也就很小，反之亦然（近似函数准确性无法保证）。
则我们从之前求$f(x_k+p)$的最小值转到求$m_k$的最小值。及问题转换为：
$$\underset{p\in R^n}{\min }{m}_k(p)s.t.||p||\le {\Delta }_k，{\Delta }_k>0$$
而如何求上述约束方程（子问题）的解$s_k$，就用到了The Dogleg Method。

Dogleg Method 方法

参考下图：

在上述子问题的情况下，有全局最优解：
$$P^B=-B^{-1}g$$
有$m_k$沿着负梯度方向的全局最优解：
$$P^U=-\frac{g^{T}g}{g^{T}Bg}g$$
而$P^B，P^U$可能在信赖域外，也可能在信赖域内，因此需要一个判断条件。注：$\tilde{p}(\tau )相当于s_k$
$$\tilde{p}(\tau )=\{\begin{array}{ll}\tau p^U,& 0\le \tau \le 1\\ p^U+(\tau -1)\left(p^B-p^U\right),& 1\le \tau \le 2\end{array}$$
当B在信赖域内部时（$\tau ＝2$），取方向为$p^B$方向，$x_{k+1}$为B点
当U在信赖域外部($0\le \tau \le 1$)，取$P^U$和边界的交点
当U在内部，B在外部($1\le \tau \le 2$)，取$UB$连线和信赖域边界的交点

这里$\tau$的计算过程有点恼火。代码中看

信赖域算法流程

附录（上述算法框架中出现的$s_k,{\rho }_k$的来历）：

二、How to use The Dogleg Method

Question

Give:
$$f(x,y)=100(y-x^2{)}^2+(1-x{)}^2$$
The initial point is $(-1.5,0.5)$,compute $minf(x,y)$
Note:
Write a program on trust region method with subproblems solved by the Dogleg method. Apply it to minimize this function. Choose $B_k={\nabla }^{2}f(x_k)$.
Experiment with the update rule of trust region. Give the first two iterates.

代码实现

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.float)
    B = B.astype(np.float)
    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)
    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')

Result presentation

可以看到，x1,x2都收敛到x1*,x2*。

总结

主要参考及感谢：
UCAS王晓老师PPT
信赖域算法+DogLeg[python实现]
信赖域算法+DogLeg[matlab实现]