精华内容
下载资源
问答
  • 信赖域算法matlab实现

    2021-03-29 15:30:15
    使用matlab实现的信赖域算法,中文注释 使用matlab实现的信赖域算法,中文注释 使用matlab实现的信赖域算法,中文注释 使用matlab实现的信赖域算法,中文注释
  • 信赖域dogleg法

    2018-06-21 10:39:03
    无约束优化理论,信赖域法框架,使用dogleg法寻找搜索方向。给出dogleg法应用的一个实例
  • 域Newton 算法的新型ELM网络(TRON-ELM), 并采用信赖域Newton 算法求解ELM网络的输出权值. 该算法首先 构造一个ELM网络代价函数的Newton 方程, 并将其作为一个无约束优化问题, 采用共轭梯度法求解, 避免了求代价...
  • 用最优化方法中的狗腿法求解信赖域子问题的算法编程实现
  • Python实现最速下降法、共轭梯度法和信赖域狗腿法源代码。可以直接运行,同时将迭代分析绘图。配有详细注释
  • 解等式约束的信赖域子问题,经典方法,文件很小
  • 应用目前流行的信赖域算法,并用带有信赖域技巧的截断共扼梯度法来解信赖域子问题。该方法是一种迭代正则化方法,它具有自动调节正则化参数的功能,克服了常用的Tikhonov正则化方法中正则化参数求取的困难。本文给出了...
  • 本文提出一个新的非线性最小二乘的信赖域方法,在该方法中每个信赖域子问题只需要一次求解,而且每次迭代的一维搜索步长因子是给定的,避开一维搜索的环节,大大地提高了算法效率.文中证明了在一定的条件下算法的...
  • 改进的固定步长的自适应信赖域算法,杭丹,王晓燕,对于无约束优化问题,在HEILONG文章的启示下构造了一个新的R-函数,并提出了一种新的不重解子问题的非单调自适应信赖域算法。与其文
  • 信赖域算法中Dogleg算法的demo,被优化的函数很简单,不过算法框架还是正确的
  • 求解大规模无约束最小化问题,使用信赖域算法,其中信赖域子问题使用截断共轭梯度法
  • 最优化方法:第3章 线性搜索与信赖域方法.ppt
  • 提出一种既有界变量又有线性等式约束的非线性优化问题的信赖域内点算法,在合理的条件下所提供的算法不仅具有整体收敛性而且保持局部收敛速率。数值计算结果说明算法的有效性。
  • 无约束优化的Broyden族信赖域算法,耿玲玲,贺祖国,本文给出了一种信赖域方法与线搜索方法的结合,信赖域方法是近二十年发展起来的一类重要的数值计算方法,它与传统的线搜索方法并
  • 信赖域算法原理Dogleg Method 方法信赖域算法流程二、How to use The Dogleg MethodQuestion代码实现总结 前言 最近在上王晓老师的最优化算法课程。课程偏硬核。记录作业中信赖域算法中狗腿(The Dogleg)方法的...


    前言

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

    一、What is The Dogleg Method?

    信赖域算法原理

    m i n f ( x , y ) minf(x,y) minf(x,y)转化为一维求步长 s k s_k sk问题。其中 s k = min ⁡ p m k ( p ) s_k=\min\limits_{p}m_{k}(p) sk=pminmk(p)。通过不断求步长,来进一步迭代出新的x,y。其中函数 m k ( p ) m_k(p) mk(p)是将函数 f ( x , y ) f(x,y) f(x,y)在点 f ( x k , y k ) f(x_k,y_k) f(xk,yk)泰勒展开后,得到的一个近似函数,也就是其展开后的前三项,如下:
    min ⁡ p ∈ R n m k ( p ) = f k + g k T p + 1 2 p T B k p , s . t . ∣ ∣ p ∣ ∣ ≤ Δ k , Δ k > 0 \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 pRnminmk(p)=fk+gkTp+21pTBkp,s.t.pΔkΔk>0

    其中 Δ k > 0 \Delta_k>0 Δk>0是信赖域半径(trust-region radius), s k s_k sk是上述方程的解,成为试探步(trial step)。
    g k = ∇ f ( x k ) g_k=\nabla f(x_k) gk=f(xk)是一阶gradient。 B k = ∇ 2 f ( x k ) B_k=\nabla ^2f(x_k) Bk=2f(xk)是二阶Hession。
    f ( x k + p ) = f k + g k T p + 1 2 p T ∇ 2 f ( x k + t p ) p , t ∈ ( 0 , 1 ) 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(xk+p)=fk+gkTp+21pT2f(xk+tp)p,t(0,1)
    可以看到 f ( x k + p ) f(x_k+p) f(xk+p) m k m_k mk之间有一个近似误差 o ( ∣ ∣ p ∣ ∣ 2 ) o(||p||^2) o(p2),当 p p p很小时候,两者误差也就很小,反之亦然(近似函数准确性无法保证)。
    则我们从之前求 f ( x k + p ) f(x_k+p) f(xk+p)的最小值转到求 m k m_k mk的最小值。及问题转换为:
    min ⁡ p ∈ R n m k ( p ) s . t . ∣ ∣ p ∣ ∣ ≤ Δ k , Δ k > 0 \min\limits_{p\in R^n}m_k(p)\qquad s.t.||p||\leq \Delta_k,\Delta_k>0 pRnminmk(p)s.t.pΔkΔk>0
    而如何求上述约束方程(子问题)的解 s k s_k sk,就用到了The Dogleg Method

    Dogleg Method 方法

    参考下图:
    1

    在上述子问题的情况下,有全局最优解:
    P B = − B − 1 g P^B=-B^{-1}g PB=B1g
    m k m_k mk沿着负梯度方向的全局最优解:
    P U = − g T g g T B g g P^U=-\frac{g^Tg}{g^TBg}g PU=gTBggTgg
    P B , P U P^B,P^U PBPU可能在信赖域外,也可能在信赖域内,因此需要一个判断条件。注: p ~ ( τ ) 相 当 于 s k \tilde{p}(\tau)相当于s_k p~(τ)sk
    p ~ ( τ ) = { τ p U , 0 ≤ τ ≤ 1 p U + ( τ − 1 ) ( p B − p U ) , 1 ≤ τ ≤ 2 \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} p~(τ)={τpU,pU+(τ1)(pBpU),0τ11τ2
    当B在信赖域内部时( τ = 2 τ=2 τ2),取方向为 p B p^B pB方向, x k + 1 x_{k+1} xk+1为B点
    当U在信赖域外部( 0 ≤ τ ≤ 1 0\leq τ \leq1 0τ1),取 P U P^U PU和边界的交点
    当U在内部,B在外部( 1 ≤ τ ≤ 2 1\leq τ \leq2 1τ2),取 U B UB UB连线和信赖域边界的交点

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

    信赖域算法流程

    1
    附录(上述算法框架中出现的 s k , ρ k s_k,ρ_k sk,ρk的来历):
    1
    2

    二、How to use The Dogleg Method

    Question

    Give:
    f ( x , y ) = 100 ( y − x 2 ) 2 + ( 1 − x ) 2 f(x,y)=100(y-x^2)^2+(1-x)^2 f(x,y)=100(yx2)2+(1x)2
    The initial point is ( − 1.5 , 0.5 ) (-1.5,0.5) (1.5,0.5),compute m i n f ( x , y ) min f(x,y) 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 = ∇ 2 f ( x k ) B_k = ∇^2f(x_k) Bk=2f(xk).
    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实现]

    展开全文
  • 信赖域算法原理

    千次阅读 2019-09-29 13:31:16
    (5)若 (rho_k≥0.75) 且 (||d_k||=h_k) ,说明这一步已经迈到了信赖域半径的边缘,并且步子有点小,可以尝试扩大信赖域半径,令 (h_k+1=2h_k) ,可以走到下一点,即 (x_k+1=x_k+d_k) (6)若 (0.25) ,说明这一步迈...

    提到最优化方法,常见的有梯度下降法(衍生出来的有批梯度下降,随机梯度下降)、牛顿法(衍生出来了拟牛顿)等。我们知道,最优化在机器学习中,是为了优化损失函数,求得其最小值,即为(mathop {min }limits_theta f({x_theta })),其中 (theta) 为损失函数的参数,最优化的目的就是找到最佳的(theta)使得损失函数最小。梯度下降的方法是求出损失函数在某一点的梯度,然后沿着负梯度方向走一小步,然后继续求新点的梯度,继续的迭代,直到达到迭代限定的次数,或者梯度极小,则迭代结束,求得最小值。对于牛顿法的原理,这里简单推导下: 先用泰勒展开去逼近目标函数(f(x))即为 [varphi (x) = f({x_k}) + f’({x_k})(x - {x_k}) + frac{1}{2}f’’({x_k}){(x - {x_k})^2}] (phi (x))是二阶展开式,高阶项被略去,既然我们把其当做目标函数的逼近式,则我们对该函数求最值,那么采用的方法就是求出展开后的二次逼近式的导数,然后另其等于0,即如下: [varphi ‘(x) = 0,] [f’({x_k}) + f’’({x_k})(x - {x_k}) = 0] 然后得出 [x = {x_k} - frac{f’({x_k})} {f’’({x_k})}] 这里我们感性的理解下这种思路,我们知道,当函数在点(x_k)泰勒展开时,取其前几阶式(牛顿法取的是二阶)用于逼近原始的函数,那么只有在展开的那个点附近,才能用泰勒展开的式子作为近似原始函数,因为我们取的毕竟是有限的阶数(高阶的我们已省略),那么展开的函数离(x_k)点较远的点可能会有较大的误差,那么上述的式子中,我们对展开函数进行求极值,求出的只是我们的近似函数的极值,并非原始函数的极值,但跟原始的极值很接近,那么我们继续在新的点上展开,然后继续求极值,那么就会越来越接近原始函数的极值。想象一下,如果原始函数在展开的点较为平滑,那么其二次逼近式在那一点附近就会很接近原始的函数,那么每次求得的近似函数的极值会跟原始的函数的极值很接近,那么牛顿法收敛的就会很快。但是如果原始函数在站开点处是一个极其不平滑,坑坑洼洼的函数,那么我们在那一点展开的二次逼近式,就可能会跟原始函数相差很远,那么求得的极值可能会跟原始函数的极值不是很接近,甚至会变得更远,这也许就是牛顿法的一些局限性的原因吧。 在上面提到的方法中我们发现,梯度下降是先确定方向,然后沿着方向走,牛顿法也很类似,展开逼近函数后,沿着逼近函数的极值,其实就是沿着一个二阶的方向走。那么本文要说的是另外的一种思路,如果我们先不确定方向,而是先确定下一步我们要走的长度,从初始点开始,确定下一步要走的长度,然后以该长度为半径的空间区域内,寻找最小的一个点,走到该点,然后再从新的点确定长度,接着依次的按照这种方法进行迭代,直到找到函数的最小值点。接下来我们详述下该方法的原理。

    信赖域方法原理
    沿着上述的问题,我们提出几个问题,根据这几个问题,展开描述信赖域方法的实现原理:

    信赖域法在求得步长范围内空间的最小值点时,求的是什么函数的最小值点,是原始的损失函数的最小值点吗?
    2.确定完步长之后,在步长范围内空间是怎样求得最小值点的?

    3.步长是怎样确定的?

    下面一一展开。

    问题一
    信赖域法在求得步长范围内空间的最小值点时,求的是什么函数的最小值点,是原始的损失函数的最小值点吗?
    求的不是原始函数的极值点。信赖域法跟牛顿法有一点很相似,同样是在初始点进行泰勒展开,取的同样是二阶式,即如下

    [f(x) approx f({x_k}) + nabla f{({x_k})^{rm T}}(x - {x_k}) + frac{1}{2}{(x - {x_k})^{rm T}}{nabla ^2}f({x_k})(x - {x_k})]

    然后进行化简,令(d=x-x_k)得

    [varphi (d) = f({x_k})nabla f{({x_k})^{rm T}}d + frac{1}{2}{d^{rm T}}{nabla ^2}f({x_k})d]

    我们把原始的损失函数在(x_k)处展开取得二阶式用于逼近原始的损失函数,然后按照一定的规则,确定一个步长,把展开的逼近函数在该空间内求最小值,实际的数学公式如下,

    [ begin{cases} varphi (d) = f({x_k})nabla f{({x_k})^{rm T}}d + frac{1}{2}{d^{rm T}}{nabla ^2}f({x_k})d\ ||d|| le {h_k}\ end{cases} ]

    解释下,为什么是限制(||d||),我们知道(d=x-x_k),那么可以理解为我们要找到的最佳点为(x),距离我们当前展开点的长度。

    问题二
    确定完步长之后,在步长范围内空间是怎样求得最小值点的?
    在第一个问题的解答,我们知道信赖域法把损失函数在某一点展开后,取二阶式,然后再确定一个步长,然后求其最小值。很明显,这是一个带有约束条件的最优化问题,求带有约束条件尤其是不等式约束条件的最优化问题我们常用到的就是KKT条件,这里我们简单讲述下KKT条件,详细的请看我总结的另一篇:KKT条件推导(还未更新)。KKT条件专门用于求解带有不等式约束条件的最优化问题,它把待最优化问题中的目标函数、不等式以及等式约束条件,转化为可求解的方程组,也就是说如果原始问题存在极值解,那么极值解一定满足KKT的条件(即为必要条件),其中满足的条件可写为一个方程组,解这个方程组即可。

    问题三
    步长是怎样确定的?
    当我们求得在步长空间范围内的最小值点时,就要确定下一步的步长范围(注意,在算法的初始,我们会给出一个初始的步长和初始点)。确定下一步步长的原则是,看在求得的最小值点处,实际的损失函数的下降值与二次逼近式的预测下降值进行对比,如果实际的下降值满足预期的,则使得函数走到该点(更新(x_k+1)),并且继续使用上一步的步长,甚至扩大步长,如果实际的下降不满足预期,则缩小步长,如果非常不满足,甚至不会更新(x_k+1),下面使用公式描述下。 假设函数求得的最佳值为(d_k),即为我们要更新的向量,那么实际的下降值为 (f({x_k}) - f({x_k} + {d_k})) 我们预测的下降值为原始的值,减去我们逼近函数的最小值,即为 (f({x_k}) - varphi ({d_k})) 二者的比值为
    [{rho _k} = frac{f({x_k}) - f({x_k} + {d_k})}{f({x_k}) - varphi ({d_k})}] 我们根据这个比值可以来自定义步长更新规则,下面举例说明。 假设初始信赖域半径为 (h_1) ,初始点为 (x_1) ,求得的最佳的极值点为 (d_1) ,那么我们假设两个参数 (0 < mu < eta < 1) 更新迭代点如果比值 (rho_1 le mu) 时,则不更新 (x_2) ,即为 (x_2 = x_1) ,如果 (rho_1 > mu) ,则 (x_2 = x_1+d_1) 。 更新信赖域半径如果 (rho_1 le mu) 时,则缩小信赖域半径,(h_2=frac{1}{2}h_1) ,如果 (mu < rho_1 < eta) ,则令 (h_2=h_1) ,如果 (rho_1 ge eta) ,则令 (h_2=2h_1) 依次迭代下去。

    信赖域算法的步骤:
    (1)从初始点 (x_0) ,初始信赖域半径 (h_0) 开始迭代

    (2)到第 k 步时,求出二次逼近式

    (3)解信赖域模型,求出位移 (d_k) ,计算 (rho_k)

    (4)若 (rho_k≤0.25) ,说明步子迈得太大了,应缩小信赖域半径,令 (h_k+1=frac{1}{2}h_k) ,这时不应该走到下一点,而应“原地踏步”,即 (x_k+1=x_k)

    (5)若 (rho_k≥0.75) 且 (||d_k||=h_k) ,说明这一步已经迈到了信赖域半径的边缘,并且步子有点小,可以尝试扩大信赖域半径,令 (h_k+1=2h_k) ,可以走到下一点,即 (x_k+1=x_k+d_k)

    (6)若 (0.25<rho_k<0.75) ,说明这一步迈出去之后,处于“可信赖”和“不可信赖”之间,可以维持当前的信赖域半径,令 (h_k+1=h_k) ,且可以走到下一点,即 (x_k+1=x_k+d_k)

    感性的理解信赖域算法
    最优化算法中迭代的思想:
    由于目标函数过于复杂,无法一次性的求出函数的极值点,我们随机的找一个点,并且在该点处的附近用一个稍微简单的函数作为近似函数来代替目标函数,因为在一个点附近可以认为函数是平滑的,然后求出近似函数的极值,即为接近目标函数极值的点,然后更新该点,继续在该点处用近似函数代替目标函数,再次求出新的极值,一次次的逼近原始的目标函数的极值点。

    信赖域半径更新规则的原因
    当实际的下降量远小于预测的下降量时,说明了我们选的步长太大,为何这么说呢?因为我们知道逼近函数只有在展开点处附近接近目标函数,但是步长太大,我们的逼近函数在较远处与目标函数就会相差很大,这样逼近函数认为下降较大的点,实际目标函数并不会下降很多,所以我们要缩小步长,保证逼近函数与原始函数相差不大,这样保证下降的方向是对的。 同样的如果下降的很符合预期,说明什么,说明我们的逼近函数非常接近目标函数,两种可能,一种是目标函数在该点处的较大范围附近是平滑的,一种是,我们的步长选的太小了,在很小步长的范围内逼近函数非常接近目标函数,无论是哪种可能,我们都要扩大步长,来加快收敛。

    泰勒展开的重要性
    感觉泰勒展开就是这些最优化算法的基础,把复杂的目标函数在某一点处展开作为近似的目标函数是整个最优化算法精髓。

    展开全文
  • 研究了一类随机线性互补问题的解法,采川信赖域线搜索与拟牛顿方法相结合的方法对其进行求解,在适当的假设条件下进行收敛性分析,得到了算法的全局收敛性,表明了算法的可行性和有效性。
  • 我们需要比较不同信赖域算法的作用以及信赖域算法和线搜索型方法的比较。我们的数值结果需要包括函数的最优值,如果条件允许的话还需要提供最优点,函数调用次数和迭代次数以及程序的运行结果,如果结果不满足条件,...

    问题描述

    采用Hebden,More-Sorensen(MS)和二维子空间算法求解Wood function, Extended Powell singular function,Trigonometric function。我们需要比较不同信赖域算法的作用以及信赖域算法和线搜索型方法的比较。我们的数值结果需要包括函数的最优值,如果条件允许的话还需要提供最优点,函数调用次数和迭代次数以及程序的运行结果,如果结果不满足条件,那么需要分析出现这种情况的可能原因。终止条件是 ∥ f ( X k + 1 ) − f ( X k ) ∥ < 1 0 − 8 \parallel f(X_{k+1}) - f(X_{k}) \parallel < 10^{-8} f(Xk+1)f(Xk)<108
    Wood function
    { f 1 ( x ) = 10 ( x 2 − x 1 2 ) , f 2 ( x ) = 1 − x 1 , f 3 ( x ) = 9 0 1 2 ( x 4 − x 3 2 ) , f 4 ( x ) = 1 − x 3 , f 5 ( x ) = 1 0 1 2 ( x 2 + x 4 − 2 ) , f 6 ( x ) = 1 0 − 1 2 ( x 2 − x 4 ) . x i n i t = ( − 3 , − 1 , − 3 , − 1 ) , f = 0 , a t ( 1 , 1 , 1 , 1 ) \left \{ \begin{array}{r l} f_{1}(x) = 10(x_2 - x_{1}^{2}),\\ f_{2}(x) = 1 - x_{1},\\ f_{3}(x) = 90^{\frac{1}{2}}(x_4 - x_{3}^{2}),\\ f_{4}(x) = 1 - x_3 ,\\ f_{5}(x) = 10^{\frac{1}{2}}(x_2 + x_4 - 2),\\ f_{6}(x) = 10^{-\frac{1}{2}}(x_2 - x_4).\\ \\ x_{init} = (-3,-1,-3,-1),\\ f = 0,\quad at \quad(1,1,1,1) \end{array} \right. f1(x)=10(x2x12),f2(x)=1x1,f3(x)=9021(x4x32),f4(x)=1x3,f5(x)=1021(x2+x42),f6(x)=1021(x2x4).xinit=(3,1,3,1),f=0,at(1,1,1,1)
    Extended Powell singular function
    { n v a r i a b l e , n i s a m u l t i p l e o f 4 , f 4 l − 3 ( x ) = x 4 l − 3 + 10 x l − 2 , f 4 l − 2 ( x ) = 5 1 2 ( x 4 l − 1 − x 4 l ) , f 4 l − 1 ( x ) = ( x 4 l − 2 − 2 x 4 l − 1 ) 2 , f 4 l ( x ) = 1 0 1 2 ( x 4 l − 3 − x 4 l ) 2 . x i n i t = ( ξ J ) , ξ 4 J − 3 = 3 , ξ 4 J − 2 = − 1 , ξ 4 J − 1 = 0 , ξ 4 J = 1 , f = 0 a t x i n i t . \left \{ \begin{array}{r l} n \quad variable,n \quad is \quad a \quad multiple \quad of \quad 4,\\ f_{4l-3}(x) = x_{4l-3} + 10x_{l-2},\\ f_{4l-2}(x) = 5^{\frac{1}{2}}(x_{4l-1} - x_{4l}),\\ f_{4l-1}(x) = (x_{4l-2} - 2x_{4l-1})^{2},\\ f_{4l}(x) = 10^{\frac{1}{2}}(x_{4l-3} - x_{4l})^2 .\\ \\ x_{init} = (\xi_{J}),\\ \xi_{4J-3} = 3,\xi_{4J-2} = -1,\xi_{4J-1} = 0,\xi_{4J} = 1,\\ f = 0\quad at \quad x_{init}. \end{array} \right. nvariable,nisamultipleof4,f4l3(x)=x4l3+10xl2,f4l2(x)=521(x4l1x4l),f4l1(x)=(x4l22x4l1)2,f4l(x)=1021(x4l3x4l)2.xinit=(ξJ),ξ4J3=3,ξ4J2=1,ξ4J1=0,ξ4J=1,f=0atxinit.
    Trigonometric function
    { n v a r i a b l e , f i ( x ) = n + i ( 1 − c o s ( x i ) ) − s i n ( x i ) − ∑ j n c o s ( x j ) , x i n i t = ( 1 n , 1 n , … , 1 n ) , f = 0. \left \{ \begin{array}{r l} n \quad variable,\\ f_{i}(x) = n + i(1 - cos(x_i)) - sin(x_i) - \sum_{j}^{n} cos(x_j),\\ x_{init} = (\frac{1}{n},\frac{1}{n},\ldots,\frac{1}{n}),\\ f = 0. \end{array} \right. nvariable,fi(x)=n+i(1cos(xi))sin(xi)jncos(xj),xinit=(n1,n1,,n1),f=0.

    信赖域方法介绍

    { q k ( d ) = f k + g k T d + 1 2 d T G k d , γ k = f ( x k + d ) − f ( x k ) q ( d ) − q ( 0 ) . \left \{ \begin{array}{r l} q_{k}(d) = f_{k} + g_{k}^{T}d + \frac{1}{2}d^{T}G_{k}d,\\ \gamma_{k} = \frac{f(x_k + d) - f(x_k)}{q(d)-q(0)}. \end{array} \right. {qk(d)=fk+gkTd+21dTGkd,γk=q(d)q(0)f(xk+d)f(xk).

    1:给定初始点 X 0 , δ 0 > 0 X_0,\delta_0 > 0 X0,δ0>0,k = 0;
    2:是否满足条件梯度为0(当k=0时),是否有 ∥ f ( X k + 1 ) − f ( X k ) ∥ < 1 0 − 8 \parallel f(X_{k+1}) - f(X_{k}) \parallel < 10^{-8} f(Xk+1)f(Xk)<108。(当 k > 0 k > 0 k>0);
    3:求解子问题,给出下一个迭代点的方向 d k d_k dk
    4:计算 γ k \gamma_k γk,若 γ k ≤ 0.25 , δ k + 1 = δ k 4 \gamma_k \leq 0.25,\delta_{k+1} = \frac{\delta_k}{4} γk0.25,δk+1=4δk;若 γ k ≥ 0.25 \gamma_k \geq 0.25 γk0.25并且 ∥ d k ∥ = δ k \parallel d_k \parallel = \delta_k dk=δk,则 δ k + 1 = 2 δ k \delta_{k+1} = 2\delta_k δk+1=2δk;否则 δ k + 1 = δ k \delta_{k+1} = \delta_k δk+1=δk
    5:若 γ k ≤ 0 \gamma_k \leq 0 γk0,则 x k + 1 = x k x_{k+1}=x_k xk+1=xk,否则 x k + 1 = x k + d k x_{k+1}=x_k + d_k xk+1=xk+dk k = k + 1 k = k+1 k=k+1,转步2。

    以上步骤中的核心就是下降方向的寻找,基于此我们引入Hebden,More-Sorensen和二维子空间算法。
    ## Hebden
    考虑方程 ( G k + ν I ) d = g k (G_k + \nu I)d = g_{k} (Gk+νI)d=gk,我们知道 d = d ( ν ) d=d(\nu) d=d(ν),下面简要描述Hebden算法的过程:\par
    首先,如果 ∥ d ( 0 ) ∥ ≤ δ k \parallel d(0) \parallel \leq \delta_{k} d(0)δk,则 d k = d ( 0 ) d_k = d(0) dk=d(0),否则的话,我们需要寻找一个 ν \nu ν,使得 ∥ d ( ν ) ∥ = δ k \parallel d(\nu) \parallel = \delta_k d(ν)=δk,根据这个思想我们定义一个函数: φ ( ν ) = ∥ d ( ν ) ∥ − δ k \varphi (\nu) = \parallel d(\nu) \parallel - \delta_k φ(ν)=d(ν)δk,也就是说我们需要求出 φ ( ν ) \varphi(\nu) φ(ν)的零点。为了节约篇幅我们直接给出零点的迭代公式,第j步迭代的公式为:

    ν j + 1 = ν j − ( φ ( ν j ) + δ k ) ∗ φ ( ν j ) φ ′ ( ν j ) δ k \nu_{j+1}=\nu_j - \frac{(\varphi(\nu_j) + \delta_k)*\varphi(\nu_j)}{\varphi^{'}(\nu_j)\delta_k} νj+1=νjφ(νj)δk(φ(νj)+δk)φ(νj)
    求出 ν \nu ν以后,代入方程\把对应的方向求解出来即可。
    More-Sorensen
    这里修改函数 φ \varphi φ的定义为 φ ( ν ) = 1 δ k − 1 ∥ d ( ν ) ∥ \varphi(\nu)=\frac{1}{\delta_k} - \frac{1}{\parallel d(\nu) \parallel} φ(ν)=δk1d(ν)1,,第j步迭代的公式为:

    ν j + 1 = ν j − φ ( ν ) φ ′ ( ν ) \nu_{j+1}=\nu_j - \frac{\varphi(\nu)}{\varphi^{'}(\nu)} νj+1=νjφ(ν)φ(ν)
    二维子空间极小化方法
    G k G_k Gk正定的时候,求解

    min ⁡ d q ( d )  s.t.  ∥ d ∥ ≤ δ , p ∈ span ⁡ [ g , G − 1 g ] \min _{d} q(d) \quad \text { s.t. }\|d\| \leq \delta, p \in \operatorname{span}\left[g, G^{-1} g\right] mindq(d) s.t. dδ,pspan[g,G1g]
    G k G_k Gk不正定的时候,
    min ⁡ d q ( d )  s.t.  ∥ d ∥ ≤ δ , p ∈ span ⁡ [ g , ( G + ν I ) − 1 g ] \min _{d} q(d) \quad \text { s.t. }\|d\| \leq \delta, p \in \operatorname{span}\left[g, (G + \nu I)^{-1} g\right] mindq(d) s.t. dδ,pspan[g,(G+νI)1g]
    其中 ν ∈ ( − λ n , − 2 λ n ) \nu \in (-\lambda_n,-2\lambda_n) ν(λn,2λn) λ n \lambda_n λn G k G_{k} Gk最小的特征值。

    数值实验

    Extended Powell singular function
    在这里插入图片描述
    Trigonometric function
    在这里插入图片描述
    Wood function
    关于这个函数,本人发现如果选取初始点 X 0 = ( − 3 , − 1 , − 3 , − 1 ) T X_0 = (-3,-1,-3,-1)^T X0=(3,1,3,1)T,上面使用的信赖域方法全都不可用,针对这个问题,本人尝试换一个初始点,采用np.random.randn函数生成一个随机点来开始迭代,这个时候会得到正确的解,下面这个表格展示的就是针对随机初始点的结果。另外就是关于线搜索方法求解,即使使用随机初始点,线搜索仍然没有收敛到精确解,反而在几步之内就停止迭代,经过本人调试发现线搜索的迭代过程中没有找到正确的步长,导致后面迭代基本不动,这里说明线搜索过程还需要进一步完善。如果修改代码使得线搜索能够找到合适的步长就变成下一个需要解决的问题。
    四种方法得到最后的收敛解分别是: ( 0.999997150.999994281.000002851.00000572 ) T (0.99999715 0.99999428 1.00000285 1.00000572)^T (0.999997150.999994281.000002851.00000572)T,
    ( 0.999997150.999994281.000002851.00000572 ) T (0.99999715 0.99999428 1.00000285 1.00000572)^T (0.999997150.999994281.000002851.00000572)T,
    ( 0.999997150.999994281.000002851.00000572 ) T (0.99999715 0.99999428 1.00000285 1.00000572)^T (0.999997150.999994281.000002851.00000572)T,
    ( 0.466528390.039486521.394356131.64631839 ) T (0.46652839 0.03948652 1.39435613 1.64631839)^T (0.466528390.039486521.394356131.64631839)T
    在这里插入图片描述

    代码

    只放EPS函数的代码

    import numpy as np
    import time
    
    def FF(X):
        n = X.shape[1]
        x0 = X[:,0:n:4];x1 = X[:,1:n:4]
        x2 = X[:,2:n:4];x3 = X[:,3:n:4]
        return ((x0 + 10*x1)**2 + 5*(x2 - x3)**2 + (x1 - 2*x2)**4 + 10*(x0 - x3)**4).sum()
    def GG(X):
        n = X.shape[1]
        x0 = X[:,0:n:4];x1 = X[:,1:n:4]
        x2 = X[:,2:n:4];x3 = X[:,3:n:4]
        vec = np.zeros_like(X)
        vec[:,0:n:4] = 2*(x0 + 10*x1) + 40*(x0 - x3)**3
        vec[:,1:n:4] = 20*(x0 + 10*x1) + 4*(x1 - 2*x2)**3
        vec[:,2:n:4] = 10*(x2 - x3) - 8*(x1 - 2*x2)**3
        vec[:,3:n:4] = -10*(x2 - x3) - 40*(x0 - x3)**3
        return vec
    def HH(X):
        n = X.shape[1]
        x0 = X[:,0:n:4];x1 = X[:,1:n:4]
        x2 = X[:,2:n:4];x3 = X[:,3:n:4]
        diag = np.zeros(n)
        diag[0:n:4] = 2 + 120*(x0 - x3)**2;diag[1:n:4] = 200 + 12*(x1 - 2*x2)**2
        diag[2:n:4] = 10 + 48*(x1 - 2*x2)**2;diag[3:n:4] = 10 + 120*(x0 - x3)**2
        xie1 = np.zeros(n - 1);xie3 = np.zeros(n - 3)
        xie1[0:n - 1:4] = 20;xie1[1:n - 1:4] = -24*(x1 - 2*x2)**2;xie1[2:n - 1:4] = - 10
        xie3[0:n - 3:4] = -120*(x0  - x3)**2
        return np.diag(diag,0) + np.diag(xie1,1) + np.diag(xie1,-1) + np.diag(xie3,3) + np.diag(xie3,-3)
    def gold(FF,GG,d,X):
        rho = 0.3
        
        am = 0;aM = 1;al = 0.5*(am + aM)
        for i in range(100):
            rig1 = FF(X) + rho*al*GG(X)@d.T
            rig2 = FF(X) + (1 - rho)*al*GG(X)@d.T
            lef = FF(X + al*d)
            if lef <= rig1:
                if lef >= rig2:
                    break
                else:
                    am = al
                    al = am + 0.618*(aM - am)
            else:
                aM = al
                al = am + 0.618*(aM - am)
        al_k = al
        #print(i)
        return al_k,2*(i + 1)
    def QQ(X,d):#d:[n,1]
        return GG(X)@d + 0.5*d.T@HH(X)@d
    def CG(A,b):#只能求解对称正定矩阵
        x = np.zeros_like(b)
        eps = 1e-5
        k = 0
        r = b - A@x;rho = r.T@r
        p = np.zeros_like(r)
        while (np.sqrt(rho) > eps*np.linalg.norm(b)) and (k < 20*b.shape[0]):
            k = k + 1
            if k == 1:
                p = r
            else:
                beta = rho/rho_h;p = r + beta*p
            w = A@p;alpha = rho/(p.T@w);x = x + alpha*p
            r = r - alpha*w;rho_h = rho;rho = r.T@r
        #print(k)
        return x
    def nubest(delta,HH,GG,X,N):#针对Hebden方法求解nu
        n = X.shape[1]
        d = CG(HH(X),-GG(X).T)
        k = 0
        nu = 0
        while (np.linalg.norm(d) > delta) and (k < N):
            k = k + 1
            d = CG(HH(X) + nu*np.eye(n,n),-GG(X).T)#[n,1]
            d_ = CG(HH(X) + nu*np.eye(n,n),- d)#[n,1]
            nu = nu - (np.linalg.norm(d) - delta)*np.linalg.norm(d)**2/(delta*(d*d_).sum())
        return nu
    def nuMS(delta,HH,GG,X,N):#针对MS方法求解nu
        n = X.shape[1]
        d = CG(HH(X),-GG(X).T)
        k = 0
        nu = 0
        while (np.linalg.norm(d) > delta) and (k < N):
            k = k + 1
            d = CG(HH(X) + nu*np.eye(n,n),-GG(X).T)#[n,1]
            d_ = CG(HH(X) + nu*np.eye(n,n),- d)#[n,1]
            nu = nu + (delta - np.linalg.norm(d))*np.linalg.norm(d)**2/(delta*(d*d_).sum())
        return nu
    def solution(FF,GG,HH,X,N,mode):
        tic = time.time()
        X_new = X.copy()
        n = X.shape[1]
        delta = 10
        fun_iter = 0
        if mode == 'Hebden':
            for k in range(N):
                nu = nubest(delta,HH,GG,X,2*n)
                d = CG(HH(X) + nu*np.eye(n,n),-GG(X).T)
                gama = (FF(X + d.T) - FF(X)).sum()/(QQ(X,d).sum())
                fun_iter += 2
                #print(QQ(X,d))
                if gama <= 0.25:
                    delta = 0.25*delta
                else:
                    if gama >= 0.75 and abs(np.linalg.norm(d) - delta) < 1e-5:
                        delta = 2*delta
                    else:
                        delta = 1.0*delta
                if gama > 0:
                    X_new = X + d.T
                else:
                    X_new = X
                #print('gama = %.4f,delta=%.4e,d=%.4e,fenzi = %.4f'%(gama,delta,np.linalg.norm(d),FF(X) - FF(X + d.T)))
                if gama > 0 and abs(FF(X_new) - FF(X)) < 1e-8:
                    break
                else:
                    X = X_new
                if k%20 == 0:
                    print('the iteration:%d,gama = %.4f,fenzi = %.4f'%(k,gama,FF(X) - FF(X + d.T)))
            ela = time.time() - tic
            print('the end iteration:%d,the grad_2 = %.5e,the time:%.3f,number:%d'
                  %(k + 1,np.linalg.norm(GG(X_new)),ela,fun_iter))
        if mode == 'erwei':
            for k in range(N):
                
                if min(np.linalg.eig(HH(X))[0]) > 0:
                    nu = 0
                else:
                    nu = nubest(delta,HH,GG,X,2*n)
                if np.linalg.norm(CG(HH(X) + nu*np.eye(n,n),-GG(X).T)) <= delta:
                    d = CG(HH(X) + nu*np.eye(n,n),-GG(X).T)
                else:
                    d_old = -GG(X).T*GG(X)@GG(X).T/(GG(X)@(HH(X) + nu*np.eye(n,n))@GG(X).T)
                    x = d_old;y = CG(HH(X) + nu*np.eye(n,n),-GG(X).T) - d_old
                    if np.linalg.norm(d_old) >= delta:
                        d = d_old*delta/np.linalg.norm(d_old)
                    else:
                        beta = (- x.T@y + np.sqrt(delta*y.T@y))/(y.T@y)
                        
                        d = x + beta*y
                
                gama = (FF(X + d.T) - FF(X)).sum()/(QQ(X,d).sum())
                fun_iter += 2
                #print(QQ(X,d))
                if gama <= 0.25:
                    delta = 0.25*delta
                else:
                    if gama >= 0.75 and abs(np.linalg.norm(d) - delta) < 1e-5:
                        delta = 2*delta
                    else:
                        delta = 1.0*delta
                if gama > 0:
                    X_new = X + d.T
                else:
                    X_new = X
                #print('gama = %.4f,delta=%.4e,d=%.4e,fenzi = %.4f'%(gama,delta,np.linalg.norm(d),FF(X) - FF(X + d.T)))
                if gama > 0 and abs(FF(X_new) - FF(X)) < 1e-8:
                    break
                else:
                    X = X_new
                if k%20 == 0:
                    print('the iteration:%d,gama = %.4f,fenzi = %.4f'%(k,gama,FF(X) - FF(X + d.T)))
            ela = time.time() - tic
            print('the end iteration:%d,the grad_2 = %.5e,the time:%.3f,number:%d'
                  %(k + 1,np.linalg.norm(GG(X_new)),ela,fun_iter))
        if mode == 'MS':
            for k in range(N):
                nu = nuMS(delta,HH,GG,X,2*n)
                d = CG(HH(X) + nu*np.eye(n,n),-GG(X).T)
                gama = (FF(X + d.T) - FF(X)).sum()/(QQ(X,d).sum())
                fun_iter += 2
                #print(QQ(X,d))
                if gama <= 0.25:
                    delta = 0.25*delta
                else:
                    if gama >= 0.75 and abs(np.linalg.norm(d) - delta) < 1e-5:
                        delta = 2*delta
                    else:
                        delta = 1.0*delta
                if gama > 0:
                    X_new = X + d.T
                else:
                    X_new = X
                #print('gama = %.4f,delta=%.4e,d=%.4e,fenzi = %.4f'%(gama,delta,np.linalg.norm(d),FF(X) - FF(X + d.T)))
                if gama > 0 and abs(FF(X_new) - FF(X)) < 1e-8:
                    break
                else:
                    X = X_new
                if k%20 == 0:
                    print('the iteration:%d,gama = %.4f,fenzi = %.4f'%(k,gama,FF(X) - FF(X + d.T)))
            ela = time.time() - tic
            print('the end iteration:%d,the grad_2 = %.5e,the time:%.3f,number:%d'
                  %(k + 1,np.linalg.norm(GG(X_new)),ela,fun_iter))
        if mode == 'newton':
            for k in range(N):
                d = CG(HH(X),-GG(X).T).T
                al,ite = gold(FF,GG,d,X)
                X_new = X + al*d
                fun_iter += 2
                fun_iter += ite
                if abs(FF(X_new) - FF(X)) < 1e-8:
                    break
                else:
                    X = X_new
                if k%20 == 0:
                    print('the iteration:%d'%(k))
            ela = time.time() - tic
            print('the end iteration:%d,the grad_2 = %.5e,the time:%.3f,number:%d'
                  %(k + 1,np.linalg.norm(GG(X_new)),ela,fun_iter))
        return X_new
                
            
    n = 100
    X = np.zeros([1,n])
    for i in range(int(n/4)):
        X[:,4*i] = 3.0
        X[:,4*i + 1] = -1
        X[:,4*i + 3] = 1
    N = 300
    MODE = ['Hebden','erwei','MS','newton']
    for mode in MODE:
        X_new = solution(FF,GG,HH,X,N,mode)
        print('the mode=%s,the value = %.5e'%(mode,FF(X_new)))
        print('--------------------')
    
    展开全文
  • 信赖域算法

    千次阅读 2018-10-02 16:28:51
    如果你关心最优化(Optimization),你一定听说过一类叫作“信赖域(Trust Region)”的算法。在本文中,我将讲述一下信赖域算法与一维搜索的区别、联系,以及信赖域算法的数学思想,实现过程。 【1】信赖域算法与一...

    如果你关心最优化(Optimization),你一定听说过一类叫作“信赖域(Trust Region)”的算法。在本文中,我将讲述一下信赖域算法与一维搜索的区别、联系,以及信赖域算法的数学思想,实现过程。

    【1】信赖域算法与一维搜索算法的区别、联系
    最优化的目标是找到极小值点,在这个过程中,我们需要从一个初始点开始,先确定一个搜索方向 d,在这个方向上作一维搜索(line search),找到此方向上的可接受点(例如,按两个准则的判定)之后,通过一定的策略调整搜索方向,然后继续在新的方向上进行一维搜索,依此类推,直到我们认为目标函数已经收敛到了极小值点。
    这种通过不断调整搜索方向,再在搜索方向上进行一维搜索的技术被很多很多算法采用,也取得了很实际的工程意义,但是,我们非要这样做不可吗?有没有另外一种途径,可以不通过“调整搜索方向→进行一维搜索”的步骤,也能求得极小值点?当然有,这就是信赖域算法干的好事。

    为了说明这两种途径所实现的算法的区别和联系,请允许我做一个可能不太恰当,但是比较形象的比喻:

    上图表述的是:如果把求最优解的过程比喻为“造一个零件”的过程的话,那么,使用一维搜索的那些算法和信赖域算法就像是两种不同的工艺,它们分别使用不同的技术(一维搜索&信赖域方法)——即两种不同的材料作为达成最终目标的基础。
    作为一个了解最优化理论并不多的人,我从我看到过的书得到的感受就是:相比使用一维搜索的那一类算法,貌似信赖域算法们的应用还不够那么多。当然这仅仅是个人感觉,勿扔砖...

    【2】信赖域算法的基本思想
    信赖域和line search同为最优化算法的基础算法,但是,从“Trust Region”这个名字你就可以看出,它是没有line search过程的,它是直接在一个region中“search”。
    在一维搜索中,从 xk 点移动到下一个点的过程,可以描述为: xk+αkdk 
    此处 αkdk 就是在 dk 方向上的位移,可以记为 sk 
    而信赖域算法是根据一定的原则,直接确定位移 sk ,同时,与一维搜索不同的是,它并没有先确定搜索方向 dk 。如果根据“某种原则”确定的位移能使目标函数值充分下降,则扩大信赖域;若不能使目标函数值充分下降,则缩小信赖域。如此迭代下去,直到收敛。

    关于这种寻优的方法,我这里又有一个比喻,希望能帮助你理解:

    要从上海火车站去人民广场,有两种方法:
    ①可以先定一个方向,比如先向西走,走着走着发现方向有点不对(人民广场应该是时尚地标啊,怎么越走感觉越郊区了呢),就调整一下方向,变成向东南方向走,诸如此类。

    ②用信赖域算法,就比如,我先划一个圈,然后在这个圈里面找离人民广场可能最接近的点,如果我的圈划得太大了,一下子就划到了莘庄(不熟悉上海的同学可以查一下地图),我一步就走到了上海南站,那还得了,马上给我回来,把圈缩小到两个地铁站的距离之内,然后再在里面找离人民广场最近的点。

    【3】信赖域算法的数学模型
    前面说了,根据一定的原则,可以直接确定位移,那么,这个原则是什么呢?
    答:利用二次模型模拟目标函数 f(x) ,再用二次模型计算出位移 s 。根据位移 s 可以确定下一点 x+s ,从而可以计算出目标函数的下降量(下降是最优化的目标),再根据下降量来决定扩大信赖域或缩小信赖域。
    那么,我该如何判定要扩大还是缩小信赖域呢?为了说明这个问题,必须先描述信赖域算法的数学模型:


    第一个式子就是我们用于模拟目标函数的二次模型,其自变量为 s ,也就是我们要求的位移。 gk 为梯度, Gk 为Hesse矩阵,袁亚湘的书上说,如果Hesse矩阵不好计算,可以利用“有限差分”来近似 Gk (不好意思我不懂),或者用拟牛顿方法来构造Hesse矩阵的近似矩阵。
    第二个式子中的 hk 是第 k 次迭代的信赖域上界(或称为信赖域半径),因此第二个式子表示的就是位移要在信赖域上界范围内。此外,第二个式子中的范数是没有指定是什么范数的,例如,是2-范数还是∞-范数之类的(在实际中都有算法用这些范数)。

    现在又回到了上面的问题:我该如何判定要扩大还是缩小信赖域呢?通过衡量二次模型与目标函数的近似程度,可以作出判定:
    第 k 次迭代的实际下降量为: Δfk=fk−f(xk+sk) 
    第 k 次迭代的预测下降量为: Δmk=fk−m(sk) 
    定义比值: rk=ΔfkΔmk 
    这个比值可以用于衡量二次模型与目标函数的近似程度,显然 r 值越接近1越好。

     

    由此,我们就可以给出一个简单的信赖域算法了。

    【4】信赖域算法的步骤
    一个考虑周全的信赖域算法可能非常麻烦,为了说明其步骤,这里只说明基本的迭代步骤:

    • 从初始点 x0 ,初始信赖域半径 h0=∥g0∥ 开始迭代
    • 到第 k 步时,计算 gk 和 Gk
    • 解信赖域模型,求出位移 sk ,计算 rk
    • 若 rk≤0.25 ,说明步子迈得太大了,应缩小信赖域半径,令 hk+1=∥sk∥4
    • 若 rk≥0.75 且 ∥sk∥=hk ,说明这一步已经迈到了信赖域半径的边缘,并且步子有点小,可以尝试扩大信赖域半径,令 hk+1=2hk
    • 若 0.25<rk<0.75 ,说明这一步迈出去之后,处于“可信赖”和“不可信赖”之间,可以维持当前的信赖域半径,令 hk+1=hk
    • 若 rk≤0 ,说明函数值是向着上升而非下降的趋势变化了(与最优化的目标相反),这说明这一步迈得错得“离谱”了,这时不应该走到下一点,而应“原地踏步”,即 xk+1=xk ,并且和上面 rk≤0.25 的情况一样缩小信赖域。反之,在 rk>0 的情况下,都可以走到下一点,即 xk+1=xk+sk


    【5】最重要的一种信赖域算法:Levenberg-Marquardt算法
    当信赖域模型中的范数 ∥s∥≤hk 取2-范数时(即 ∥s∥2≤hk ),就得到了Levenberg-Marquardt算法(简称LM算法)的数学模型:

    具体请看这里

    【6】信赖域算法的收敛性
    信赖域算法具有整体收敛性。这个证明我没看(太长了),此处略。

    转载自:https://www.codelast.com/%E5%8E%9F%E5%88%9B%E4%BF%A1%E8%B5%96%E5%9F%9Ftrust-region%E7%AE%97%E6%B3%95%E6%98%AF%E6%80%8E%E4%B9%88%E4%B8%80%E5%9B%9E%E4%BA%8B/

    展开全文
  • 信赖域方法

    千次阅读 2019-06-11 21:31:42
    参考: 1: https://www.codelast.com/原创信赖域trust-region算法是怎么一回事/
  • 最优化理论——信赖域方法算法思想算法步骤子问题算法步骤代码示例原问题代码示例 算法思想 算法步骤 子问题 算法步骤 代码 function [d,val,lam,k]=trustq(gk,Bk,dta) % 功能: 求解信赖域子问题: min qk(d)...
  • 信赖域-Dogleg方法

    2021-04-24 15:45:12
    信赖域方法求极值点,初始值为[-1, 0],误差为1e-5,求Cauchy点时采用dogleg方法。 1. 代码实现: import numpy as np import matplotlib.pyplot as plt # Rosenbrock函数 def rosenbrock(x): return 100 * (x...
  • 提出并实现了用于广域阻尼控制器结合参数设计的新方法,结合非线性时域仿真技术与信赖域方法,解决了优化目标函数构造、控制参数初值获取和信赖域法子问题模型获取3个子问题,实现了多运行方式下、多个广域阻尼控制器的...
  • 提出利用信赖域狗腿法解决通用非球面光学元件抛光后置求解的强非线性问题。基于低序体表示和齐次变换方法建立适用于任意非球面面形和任意轨迹抛光的正向运动学模型;基于信赖域狗腿法建立强非线性正向运动学模型的...
  • 信赖域狗腿(dogleg)方法

    万次阅读 2018-06-03 18:14:53
    信赖域狗腿方法 信赖域方法(Trust-region methods)又称为TR方法,它是一种最优化方法,能够保证最优化方法总体收敛。现今,信赖域算法广泛应用于应用数学、物理、化学、工程学、计算机科学、生物学与医学等学科...
  • 目录 信赖域方法 信赖域方法概述 基于柯西点的算法 柯西点 为什么要改良柯西法 狗腿法 二维子空间极小化 全局收敛性 柯西点收益 平稳点的收敛性 子问...
  • 信赖域方法前面提到了Line Search算法分为两步,首先确定方向,然后确定步长,实际上是假设近似模型在某一方向上可以较好的代表真实模型。Trust region算法则在此基础上,假设在一个选定的可信赖区域中,可以选择一...
  • 关于信赖域方法的具体介绍,请参照博客添加链接描述
  • 提出了基于柔度的最小...采用信赖域方法求解该问题,使优化过程具有更强的鲁棒性和可靠性,同时解决了柔度灵敏度方法应用在复杂结构中的问题。最后,通过某导弹发射台的骨架损伤数值仿真,验证了该方法的可行性和有效性。

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 6,316
精华内容 2,526
关键字:

信赖域