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

文章目录前言一、What is The Dogleg Method？信赖域算法原理Dogleg Method 方法信赖域算法流程二、How to use The Dogleg MethodQuestion代码实现总结前言最近在上王晓老师的最优化算法课程。课程偏硬核。记录作业中信赖域算法中狗腿（The Dogleg）方法的实现。一、What is The Dogleg Method？信赖域算法原理把

等待戈多。

8706人浏览 · 2021-10-17 19:27:57

等待戈多。 · 2021-10-17 19:27:57 发布

前言

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

一、What is The Dogleg Method？

信赖域算法原理

把 $m i n f (x, y)$ 转化为一维求步长 $s_k$ 问题。其中 $s_k=\min\limits_{p}m_{k}(p)$ 。通过不断求步长，来进一步迭代出新的x,y。其中函数 $m_k(p)$ 是将函数 $f (x, y)$ 在点 $f(x_k,y_k)$ 泰勒展开后，得到的一个近似函数，也就是其展开后的前三项，如下：
$\min\limits_{p\in R^n}m_k(p)=f_k+g_k^Tp+\frac{1}{2}p^TB_kp,\qquad s.t.||p||\leq \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 ^2f(x_k)$ 是二阶Hession。
且 $f(x_k+p)=f_k+g_k^Tp+\frac{1}{2}p^T\nabla ^2f(x_k+tp)p,\qquad t\in (0,1)$
可以看到 $f(x_k+p)$ 与 $m_k$ 之间有一个近似误差 $o(||p||^2)$ ，当 $p$ 很小时候，两者误差也就很小，反之亦然（近似函数准确性无法保证）。
则我们从之前求 $f(x_k+p)$ 的最小值转到求 $m_k$ 的最小值。及问题转换为：
$\min\limits_{p\in R^n}m_k(p)\qquad s.t.||p||\leq \Delta_k，\Delta_k>0$
而如何求上述约束方程（子问题）的解 $s_k$ ，就用到了The Dogleg Method。

Dogleg Method 方法

参考下图：

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

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

信赖域算法流程

附录（上述算法框架中出现的 $s_k,ρ_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 $m i n f (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 = ∇^2f(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')