我们知道,闭区间上的一元连续函数必在区间上取得最大值和最小值。实践中我们需要能数值地确定含有 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初始化为xstarya=f(a)c=a+s #c初始化为a+syc=f(c)if yc>ya: #s为上升方向a,c=c,a #交换a,cya,yc=yc,yas=-s #调整s为下降方向b=c+s #b置为c+syb=f(b)while yb<=yc: #b同a、c处于同一下降区间a,ya=c,yc #a置为cc,yc=b,yb #c置为bs*=lamd #扩大步长sb=c+s #重置b为c+syb=f(b)if a>b: #若a大b小a,b=b,a #交换a,breturn 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)所示。