我们知道,闭区间上的一元连续函数必在区间上取得最大值和最小值。实践中我们需要能数值地确定含有
f
(
x
)
f(x)
f(x)的唯一最优解
x
0
x_0
x0的区间
[
a
,
b
]
[a,b]
[a,b]。这里介绍寻求连续函数
f
(
x
)
f(x)
f(x)在一点
x
∗
x^*
x∗附近单峰区间的包围算法及其Python实现。
算法的思想是从
x
∗
x^*
x∗开始沿着
f
(
x
)
f(x)
f(x)下降的方向逐步探索,直至第一次遇到上升:
(1)设定初始步长
s
s
s和缩放系数
λ
\lambda
λ;
(2)设定
a
,
c
a,c
a,c为
x
∗
x^*
x∗及
x
∗
+
s
x^*+s
x∗+s,并确保
f
(
a
)
>
f
(
c
)
f(a)>f(c)
f(a)>f(c)。必要时需调整搜索方向;
(3)取
b
b
b为
c
+
s
c+s
c+s。比较
f
(
c
)
f(c)
f(c)与
f
(
b
)
f(b)
f(b),若
f
(
c
)
<
f
(
b
)
f(c)<f(b)
f(c)<f(b),意味着
f
(
a
)
>
f
(
c
)
<
f
(
b
)
f(a)>f(c)<f(b)
f(a)>f(c)<f(b),即可断定
x
0
∈
(
a
,
b
)
x_0\in(a,b)
x0∈(a,b),
(
a
,
b
)
(a,b)
(a,b)即为所求;
(4)否则,令
a
a
a为
c
c
c,
c
c
c为
b
b
b,并调整步长
s
s
s为
λ
×
s
\lambda\times s
λ×s,转(3)。
通常,将步长
s
s
s初始化为
∣
s
∣
=
0.01
|s|=0.01
∣s∣=0.01,
λ
\lambda
λ初始化为2。当
s
s
s初始化为正数时,首次探索自左向右。反之,
s
s
s取负数初始值则首次探索自右向左。上图给出了一个按此思想探寻包围函数
f
(
x
)
f(x)
f(x)极小值点的区间
(
a
,
b
)
(a,b)
(a,b)的实例。初始时,步长
s
<
0
s<0
s<0,以起始点
x
∗
x^*
x∗为
a
a
a,
c
=
a
+
s
<
a
c=a+s<a
c=a+s<a位于
a
a
a的左边。如图中(a)所示。由于
x
∗
x^*
x∗位于
f
(
x
)
f(x)
f(x)的下降区间,故
f
(
c
)
>
f
(
a
)
f(c)>f(a)
f(c)>f(a)。交换
a
a
a和
c
c
c,如图中(b)所示。令
s
=
−
s
s=-s
s=−s为下降方向,则此后
s
>
0
s>0
s>0。取
b
=
c
+
s
b=c+s
b=c+s,如图中(c )所示。由于
f
(
b
)
<
f
(
c
)
f(b)<f(c)
f(b)<f(c),故将
a
a
a置为
c
c
c,
c
c
c置为
b
b
b,如图中(d)所示。将步长扩大
s
=
λ
s
s=\lambda s
s=λs,并置
b
=
c
+
s
b=c+s
b=c+s,由于
f
(
b
)
f(b)
f(b)仍然小于
f
(
c
)
f(c)
f(c)(见图中(e)),故再次将
a
a
a和
c
c
c移至
c
c
c,
b
b
b,扩大步长
s
s
s为
λ
s
\lambda s
λs并置
b
=
c
+
s
b=c+s
b=c+s。此时,
f
(
a
)
>
f
(
c
)
<
f
(
b
)
f(a)>f(c)<f(b)
f(a)>f(c)<f(b),见图中(f)。可见
f
(
x
)
f(x)
f(x)的极小值点被区间
(
a
,
b
)
(a,b)
(a,b)所包围,即
f
(
x
)
f(x)
f(x)在
(
a
,
b
)
(a,b)
(a,b)为一单峰函数。故
(
a
,
b
)
(a,b)
(a,b)即为所求。
下列代码实现算法。
def myBracket(f,xstar,s=1e-2,lamd=2):
a=xstar #a初始化为xstar
ya=f(a)
c=a+s #c初始化为a+s
yc=f(c)
if yc>ya: #s为上升方向
a,c=c,a #交换a,c
ya,yc=yc,ya
s=-s #调整s为下降方向
b=c+s #b置为c+s
yb=f(b)
while yb<=yc: #b同a、c处于同一下降区间
a,ya=c,yc #a置为c
c,yc=b,yb #c置为b
s*=lamd #扩大步长s
b=c+s #重置b为c+s
yb=f(b)
if a>b: #若a大b小
a,b=b,a #交换a,b
return a,b
函数myBracket有4个参数:f表示函数
f
(
x
)
f(x)
f(x)。xstar表示初始点
x
∗
x^*
x∗。s表示步长
s
s
s,缺省值为0.01。lamd表示放大系数
λ
\lambda
λ,缺省值为2。函数体中第2~5行分别初始化点
a
a
a和
c
c
c以及对应的函数值
y
a
y_a
ya和
y
c
y_c
yc。第6~9行的if语句矫正a,c保证ya>yc,及步长方向与函数值的下降方向一致。第10~11行初始化点
b
b
b及对应的函数值
y
c
y_c
yc。第10~15行的while循环执行区间
(
a
,
b
)
(a,b)
(a,b)的迭代,直至
y
c
<
y
b
y_c<y_b
yc<yb。第16~17行的if语句确保
a
<
b
a<b
a<b。
例1 设函数
f
(
x
)
=
sin
x
f(x)=\sin{x}
f(x)=sinx,用myBracket分别对
x
∗
=
0
x^*=0
x∗=0、
x
∗
=
4
π
3
x^*=\frac{4\pi}{3}
x∗=34π,计算单峰区间。
解:下列代码完成本例计算。
import numpy as np #导入numpy
xstar=0 #设置xstar为0
print(myBracket(np.sin, xstar)) #计算附近的单峰区间
xstar=4*np.pi/3 #设置xstar为4pi/3
print(myBracket(np.sin, xstar)) #计算附近的单峰区间
利用行内注释信息,不难理解本程序。运行程序,输出
(-2.55, -0.6300000000000001)
(4.50879020478639, 5.46879020478639)
第1行输出包含距
x
∗
=
0
x^*=0
x∗=0最近的局部最小值点
x
0
=
−
π
2
x_0=-\frac{\pi}{2}
x0=−2π的区间
(
−
2.55
,
−
0.63
)
(-2.55,-0.63)
(−2.55,−0.63),如下图(a)所示。第2行输出包含距
x
∗
=
4
π
3
x^*=\frac{4\pi}{3}
x∗=34π最近的局部最小值点
x
0
=
3
π
2
x_0=\frac{3\pi}{2}
x0=23π的区间
(
4.51
,
5.47
)
(4.51,5.47)
(4.51,5.47),如下图(b)所示。