最优化方法Python计算：连续函数的单峰区间计算

我们知道，闭区间上的一元连续函数必在区间上取得最大值和最小值。实践中我们需要能数值地确定含有$f(x)$的唯一最优解 x 0 x_0 x0​的区间$[a,b]$。这里介绍寻求连续函数$f(x)$在一点 x ∗ x^* x∗附近单峰区间的包围算法及其Python实现。
算法的思想是从 x ∗ x^* x∗开始沿着$f(x)$下降的方向逐步探索，直至第一次遇到上升：
(1)设定初始步长$s$和缩放系数$\lambda$；
(2)设定$a,c$为 x ∗ x^* x∗及 x ∗ + s x^*+s x∗+s，并确保$f(a)>f(c)$。必要时需调整搜索方向；
(3)取$b$为$c+s$。比较$f(c)$与$f(b)$，若$f(c)<f(b)$，意味着$f(a)>f(c)<f(b)$，即可断定 x 0 ∈ ( a , b ) x_0\\in(a,b) x0​∈(a,b)，$(a,b)$即为所求；
(4)否则，令$a$为$c$，$c$为$b$，并调整步长$s$为$\lambda \times s$，转(3)。

通常，将步长$s$初始化为$|s|=0.01$，$\lambda$初始化为2。当$s$初始化为正数时，首次探索自左向右。反之，$s$取负数初始值则首次探索自右向左。上图给出了一个按此思想探寻包围函数$f(x)$极小值点的区间$(a,b)$的实例。初始时，步长$s<0$，以起始点 x ∗ x^* x∗为$a$，$c=a+s<a$位于$a$的左边。如图中(a)所示。由于 x ∗ x^* x∗位于$f(x)$的下降区间，故$f(c)>f(a)$。交换$a$和$c$，如图中(b)所示。令$s=-s$为下降方向，则此后$s>0$。取$b=c+s$，如图中(c )所示。由于$f(b)<f(c)$，故将$a$置为$c$，$c$置为$b$，如图中(d)所示。将步长扩大$s=\lambda s$，并置$b=c+s$，由于$f(b)$仍然小于$f(c)$（见图中(e)），故再次将$a$和$c$移至$c$，$b$，扩大步长$s$为$\lambda s$并置$b=c+s$。此时，$f(a)>f(c)<f(b)$，见图中(f)。可见$f(x)$的极小值点被区间$(a,b)$所包围，即$f(x)$在$(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)$。xstar表示初始点 x ∗ x^* x∗。s表示步长$s$，缺省值为0.01。lamd表示放大系数$\lambda$，缺省值为2。函数体中第2～5行分别初始化点$a$和$c$以及对应的函数值 y a y_a ya​和 y c y_c yc​。第6～9行的if语句矫正a，c保证ya>yc，及步长方向与函数值的下降方向一致。第10～11行初始化点$b$及对应的函数值 y c y_c yc​。第10～15行的while循环执行区间$(a,b)$的迭代，直至 y c < y b y_c<y_b yc​<yb​。第16～17行的if语句确保$a<b$。
例1 设函数$f(x)=\sin x$，用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)$，如下图(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)$，如下图(b)所示。