前言

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

一、What is The Dogleg Method？

信赖域算法原理

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

其中Δk>0\Delta_k>0Δk>0是信赖域半径（trust-region radius）,sks_ksk是上述方程的解，成为试探步（trial step）。
gk=∇f(xk)g_k=\nabla f(x_k)gk=∇f(xk)是一阶gradient。Bk=∇2f(xk)B_k=\nabla ^2f(x_k)Bk=∇2f(xk)是二阶Hession。
且f(xk+p)=fk+gkTp+12pT∇2f(xk+tp)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+21pT∇2f(xk+tp)p,t∈(0,1)
可以看到f(xk+p)f(x_k+p)f(xk+p)与mkm_kmk之间有一个近似误差o(∣∣p∣∣2)o(||p||^2)o(∣∣p∣∣2)，当ppp很小时候，两者误差也就很小，反之亦然（近似函数准确性无法保证）。
则我们从之前求f(xk+p)f(x_k+p)f(xk+p)的最小值转到求mkm_kmk的最小值。及问题转换为：
min⁡p∈Rnmk(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>0p∈Rnminmk(p)s.t.∣∣p∣∣≤Δk，Δk>0
而如何求上述约束方程（子问题）的解sks_ksk，就用到了The Dogleg Method。

Dogleg Method 方法

参考下图：

在上述子问题的情况下，有全局最优解：
PB=−B−1gP^B=-B^{-1}gPB=−B−1g
有mkm_kmk沿着负梯度方向的全局最优解：
PU=−gTggTBggP^U=-\frac{g^Tg}{g^TBg}gPU=−gTBggTgg
而PB，PUP^B，P^UPB，PU可能在信赖域外，也可能在信赖域内，因此需要一个判断条件。注：p~(τ)相当于sk\tilde{p}(\tau)相当于s_kp~(τ)相当于sk
p~(τ)={τpU,0≤τ≤1pU+(τ−1)(pB−pU),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)(pB−pU),0≤τ≤11≤τ≤2
当B在信赖域内部时（τ＝2τ＝2τ＝2），取方向为pBp^BpB方向，xk+1x_{k+1}xk+1为B点
当U在信赖域外部(0≤τ≤10\leq τ \leq10≤τ≤1)，取PUP^UPU和边界的交点
当U在内部，B在外部(1≤τ≤21\leq τ \leq21≤τ≤2)，取UBUBUB连线和信赖域边界的交点

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

信赖域算法流程

附录（上述算法框架中出现的sk,ρks_k,ρ_ksk,ρk的来历）：

二、How to use The Dogleg Method

Question

Give:
f(x,y)=100(y−x2)2+(1−x)2f(x,y)=100(y-x^2)^2+(1-x)^2f(x,y)=100(y−x2)2+(1−x)2
The initial point is (−1.5,0.5)(-1.5,0.5)(−1.5,0.5),compute minf(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 Bk=∇2f(xk)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 pltdef function(x1,x2):"""定义函数的表达式Args:x1 : 变量x1x2 : 变量x2Returns:函数表达式"""return 100*(x2-x1**2)**2+(1-x1)**2def gradient_function(x1,x2):"""定义函数的一阶梯度Args:x1 : 变量x1x2 : 变量x2Returns:函数的一阶梯度"""g=[[-400*(x1*x2-x1**3)+2*x1-2],[200*(x2-x1**2)]]g = np.array(g)return gdef Hessian_function(x1,x2):"""定义函数二阶Hessian矩阵Args:x1 : 变量x1x2 : 变量x2Returns:函数二阶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 : 变量x1x2 : 变量x2p : 下降的试探步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 mkdef Dogleg_Method(x1,x2,delta):"""Dogleg Method实现Args:x1 : 变量x1x2 : 变量x2delta : 信赖域半径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-PUPB_norm = np.linalg.norm(PB)PU_norm = np.linalg.norm(PU)PB_U_norm = np.linalg.norm(PB_U)#判断τif PB_norm <= delta:tao = 2elif PU_norm >= delta:tao = delta/PU_normelse: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*PUelif 1<tao<=2:s_k = PU+(tao-1)*(PB-PU)return s_kdef TrustRegion(x1,x2,delta_max):"""信赖域算法Args:x1 : 初始值x1x2 : 初始值x2delta_max : 最大信赖域半径Returns:x1 : 优化后x1x2 : 优化后x2"""delta = delta_maxk = 0#计算初始的函数梯度范数#终止判别条件中的epsilonepsilon = 1e-9maxk = 1000x1_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:breakif k > maxk:break#利用DogLeg_Method求解子问题迭代步长sksk = 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 / 4elif r > 0.75 and np.linalg.norm(sk) == delta:delta = np.min((2 * delta,delta_max))else:passif r <= 0:passelse:x1 = x1_newx2 = x2_newk = k + 1x1_log.append(x1)x2_log.append(x2)return x1_log,x2_log
if __name__ =='__main__':x1=0.5 x2=1.5delta_max = 20x1_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.clfplt.title('x2_convergence')plt.plot(x2_log)plt.savefig('x2.png')

Result presentation

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

总结

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

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

数值最优化—无约束问题信赖域算法
目录一.参考二.线性搜索与信赖域算法区别三.信赖域算法四.信赖域子问题的求解 1. 精确求解方法 2. 折线方法(Dogleg Method 狗腿法) 一.参考 <数值最优化算法与理论& ...
【优化算法】信赖域折线（狗腿）法（Trust Region Dogleg Method）
出发点刘红英老师的最优化课程需要用到信赖域方法,求解问题并画图,网上找到的中文内容的矩阵实验室(MATLAB)信赖域折线(狗腿)法求解代码,都很繁琐,而且心理感觉不可靠(手动狗头),而且我没心思一点 ...
matlab cg steihaug,信赖域（一）：Cauchy Point与Dogleg
信赖域(一):Cauchy Point与Dogleg 王金戈从本文开始介绍优化方法的另一大类--信赖域方法. 基本思想先回顾一下前面几篇文章介绍过的线搜索方法.在线搜索方法的每次迭代中,先确定一个 ...
信赖域狗腿（dogleg）方法
信赖域狗腿方法信赖域方法(Trust-region methods)又称为TR方法,它是一种最优化方法,能够保证最优化方法总体收敛.现今,信赖域算法广泛应用于应用数学.物理.化学.工程学.计算机科学 ...
信赖域(Trust Region)算法和L-M算法
近期在将样条数据进行公式化. 这时候发现多项高斯函数函数很好用,里面发现了两种优化高斯函数参数的算法一种是信赖域算法(Trust Region),另一种是L-M算法(Levenberg-Marquad ...
数学知识-- 信赖域(Trust Region)算法是怎么一回事
信赖域(Trust Region)算法是怎么一回事转载自: https://www.codelast.com/原创信赖域trust-region算法是怎么一回事/ 如果你关心最优化(Optimiza ...
信赖域(Trust Region)算法
文章来源:http://www.codelast.com/ 如果你关心最优化(Optimization),你一定听说过一类叫作"信赖域(Trust Region)"的算法.在本文中 ...
迭代求解最优化问题——信赖域方法
信赖域方法前面提到了Line Search算法分为两步,首先确定方向,然后确定步长,实际上是假设近似模型在某一方向上可以较好的代表真实模型.Trust region算法则在此基础上,假设在一个选定的 ...
数值优化 Ch.4 信赖域方法
目录信赖域方法信赖域方法概述基于柯西点的算法柯西点为什么要改良柯西法狗腿法二维子空间极小化全局收敛性柯西点收益平稳点的收敛性子问题的迭代解法困难的情况定理4.1的证明基于近 ...

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

文章目录

前言