- 前言
- 一、What is The Dogleg Method？
- 信赖域算法原理
- Dogleg Method 方法
- 信赖域算法流程
- 二、How to use The Dogleg Method
- Question
- 代码实现
- Result presentation
- 总结
前言
最近在上王晓老师的最优化算法课程。课程偏硬核。记录作业中信赖域算法中狗腿（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=minlimits_{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
minlimits_{pin 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
p∈Rnminmk(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^Tnabla ^2f(x_k+tp)p,qquad tin (0,1)
f(xk+p)=fk+gkTp+21pT∇2f(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(∣∣p∣∣2)，当
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
minlimits_{pin R^n}m_k(p)qquad s.t.||p||leq Delta_k，Delta_k>0
p∈Rnminmk(p)s.t.∣∣p∣∣≤Δk，Δk>0
而如何求上述约束方程（子问题）的解
s
k
s_k
sk，就用到了The Dogleg Method。
参考下图：
在上述子问题的情况下，有全局最优解：
P
B
=
−
B
−
1
g
P^B=-B^{-1}g
PB=−B−1g
有
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
PB，PU可能在信赖域外，也可能在信赖域内，因此需要一个判断条件。注：
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 2end{cases}
p~(τ)={τpU,pU+(τ−1)(pB−pU),0≤τ≤11≤τ≤2
当B在信赖域内部时（
τ
＝
2
τ＝2
τ＝2），取方向为
p
B
p^B
pB方向，
x
k
+
1
x_{k+1}
xk+1为B点
当U在信赖域外部(
0
≤
τ
≤
1
0leq τ leq1
0≤τ≤1)，取
P
U
P^U
PU和边界的交点
当U在内部，B在外部(
1
≤
τ
≤
2
1leq τ leq2
1≤τ≤2)，取
U
B
UB
UB连线和信赖域边界的交点
这里 τ τ τ的计算过程有点恼火。代码中看
信赖域算法流程
附录（上述算法框架中出现的
s
k
,
ρ
k
s_k,ρ_k
sk,ρk的来历）：
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(y−x2)2+(1−x)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 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实现]
