最优化方法Python计算：一元函数搜索算法——二次插值法

已知连续函数$f(x)$在$x^{\ast }$近旁存在最优解$x_0$。对博文《最优化方法Python计算：连续函数的单峰区间计算》讨论的$f(x)$单峰区间的包围算法稍加修改，可算得$f(x)$包含$x_0$的单峰区间$[a,c]$及含于$(a,c)$内的满足$f(a)>f(b)<f(c)$的点$b$，如下图所示。

相应地修改实现该算法的myBracket函数如下：

def myBracket(f,xstar,s=1e-4,lamd=2):
   a=xstar						#a初始化为xstar
   ya=f(a)
   b=a+s						#c初始化为a+s
   yb=f(b)
   if yb>ya:					#s为上升方向
      a,b=b,a					#交换a,b
      ya,yb=yb,ya
      s=-s						#调整s为下降方向
   c=b+s						#c置为b+s
   yc=f(c)
   while yc<=yb:				#c同a、b处于同一下降区间
      a,ya=b,yb					#a置为b
      b,yb=c,yc					#b置为c
      s*=lamd					#扩大步长s
      c=b+s						#重置c为b+s
      yb=f(b)
   if a>c:						#若a大c小
      a,c=c,a					#交换a,c
   return a,b,c

令$y_a=f(a)$，$y_b=f(b)$及$y_c=f(c)$。构造通过点$(a,y_a)$，$(b,y_b)$及$(c,y_c)$的二次函数$q(x)=p_0+p_{1}x+p_2{x}^2$，即
$$\{\begin{array}{l}{y}_a=p_0+p_{1}a+p_2{a}^2\\ y_b=p_0+p_{1}b+p_2{b}^2\\ y_c=p_0+p_{1}c+p_2{c}^2\end{array}.$$
解出$p_0$，$p_1$和$p_2$:
$$\{\begin{array}{l}{p}_0=\frac{bc{y}_a}{(b-a)(c-a)}-\frac{ac{y}_b}{(b-a)(c-b)}+\frac{ab{y}_c}{(c-a)(c-b)}\\ p_1=-\frac{(b+c)y_a}{(b-a)(c-a)}+\frac{(a+c)y_c}{(b-a)(c-b)}-\frac{(a+b)y_c}{(c-a)(c-b)}\\ p_2=\frac{y_a}{(b-a)(c-a)}-\frac{y_b}{(b-a)(c-b)}+\frac{y_c}{(c-a)(c-b)}\end{array}.$$
二次式
$$q(x)=p_0+p_{1}x+p_2{x}^2$$
称为$f(x)$在$a,b,c$处的二次插值多项式。由于$q(a)=f(a)$，$q(b)=f(b)$，$q(c)=f(c)$，故$q(a)>q(b)<q(c)$。即$(a,b)$为$q(x)$的一个单峰区间。而二次函数仅有一个最小值，故$q(x)$的最优解$x_0^{{}^{\prime }}\in (a,c)$，且为其驻点。即
$$x_0^{{}^{\prime }}=-\frac{1}{2}\cdot \frac{p_1}{p_2}=\frac{1}{2}\left[\frac{\frac{(b+c)y_a}{(b-a)(c-a)}-\frac{(a+c)y_c}{(b-a)(c-b)}+\frac{(a+b)y_c}{(c-a)(c-b)}}{\frac{y_a}{(b-a)(c-a)}-\frac{y_b}{(b-a)(c-b)}+\frac{y_c}{(c-a)(c-b)}}\right].$$
直观地看，从$f(x)$的极小值点$x_0$的近旁$x^{\ast }$出发，运用包围算法计算$f(x)$的单峰区间$(a,c)$，及$b\in (a,c)$，满足$f(x)>f(b)<f(c)$,过三点$(a,f(a))$，$(b,f(b))$和$(c,f(c))$的插值二次函数$q(x)$的极小值点$x_0^{{}^{\prime }}\in (a,c)$与出发点$x^{\ast }$相比，离$x_0$更近，如下图所示。

于是，对给定的容错误差$\epsilon$，从邻近$f(x)$的极小值点$x_0$的点$x^{\ast }$出发，用包围算法算得含点$x_0$的$f(x)$的单峰区间$(a,c)$及$b\in (a,c)$，若其长度$|c-a|<\epsilon$，则即以二次插函数$q(x)$的极小值点$x_0^{{}^{\prime }}$作为$x_0$的近似值。否则，以$x_0^{{}^{\prime }}$作为包围算法新的出发点$x^{\ast }$计算$f(x)$的单峰区间$(a,c)$及含于其内的点$b$，重复上述计算二次插值函数$q(x)$的最小值点$x_0^{{}^{\prime }}$。循环往复，直至$|c-a|<\epsilon$。此时，$x_0^{{}^{\prime }}$即为$x_0$的近似值。这一方法称为二次插值法。实现为如下的Python函数

from scipy.optimize import OptimizeResult
def interpolation(fun,xstar,eps=0.0002,**options):
    k=1
    a,b,c=myBracket(fun,xstar)
    ya,yb,yc=fun(a),fun(b),fun(c)
    p1=-ya*(b+c)/((b-a)*(c-a))+yb*(a+c)/((b-a)*(c-b))-yc*(a+b)/((c-a)*(c-b))
    p2=ya/((b-a)*(c-a))-yb/((b-a)*(c-b))+yc/((c-a)*(c-b))
    x01=-0.5*p1/p2
    while c-a>eps:
        xstar=x01
        a,b,c=myBracket(fun,xstar)
        ya,yb,yc=fun(a),fun(b),fun(c)
        p1=-ya*(b+c)/((b-a)*(c-a))+yb*(a+c)/((b-a)*(c-b))-yc*(a+b)/((c-a)*(c-b))
        p2=ya/((b-a)*(c-a))-yb/((b-a)*(c-b))+yc/((c-a)*(c-b))
        x01=-0.5*p1/p2
        k+=1
    bestx=x01
    besty=fun(bestx)
    return OptimizeResult(fun=besty, x=bestx, nit=k)

程序的第2~19行定义函数interpolation，实现二次插值算法。参数fun、xstar和eps分别表示目标函数$f(x)$，起点$x^{\ast }$和容错误差$\epsilon$。options用来实现minimize_scalar向interpolation传递xstar和eps等实际参数的机制。
第3行将迭代次数k初始化为1。
第4~8行执行第一次迭代：第4行调用myBracket函数，从$x^{\ast }$起计算$f(x)$的单峰区间$(a,c)$及$b\in (a,c)$，满足$f(a)>f(b)<f(c)$。第5行计算$f(x)$在$a,b,c$处的函数值为ya，yb，yc。第6、7行分别计算二次插值函数$q(x)$的一次项系数和二次项系数p1，p2。第8行计算$q(x)$的驻点$x_0^{{}^{\prime }}=-\frac{1}{2}\frac{p_1}{p_2}$赋予x01。
第9~16行的while循环执行余下的迭代操作：第10行将x01赋予xstar，以更新起点$x^{\ast }$。第11~15行执行与第4~8行相同的操作。第16行将迭代次数k自增1。循环往复，直至区间$(a,c)$的长度$c-a$不超过$\epsilon$。
需要注意的是，之所以要对表示容错误差$\epsilon$的参数eps设置缺省值0.0002，是因为myBracket中我们将步长s的缺省值设置为0.0001，使得所计算的单峰区间$(a,b)$的长度最小为0.0002。这样，才不至于使得此处的{\bf{while}}语句陷入死循环。读者在使用interpolation时需注意eps参数值不要小于myBracket的s参数的初始值。
例1 用上述程序定义的interpolation函数，计算函数$f(x)=x^2+4\cos x$在$x_1=1.5$近旁的极小值点。
解：下列代码完成计算。

from scipy.optimize import minimize_scalar
f=lambda x:x**2+4*np.cos(x)
xstar=1.5
res=minimize_scalar(f, method=interpolation, options={'xstar':1.5,'eps':0.001})
print(res)

程序的第2行设置目标函数$f(x)=x^2+4\cos x$赋予f，第3行设置起始点$x^{\ast }$赋予xstar，第4行调用minimize_scalar函数（第1行导入），传递f给参数fun、interpolation给method并通过options向interpolation传递参数1.5给xstar，0.001给eps。运行程序，输出

 fun: 2.316808419788213
 nit: 3
   x: 1.895494265404134

与博文《一元函数搜索算法——牛顿法》中例1的计算结果相比较，对同一函数$f(x)$，相同起点$x^{\ast }=1.5$，牛顿方法在容错误差$\epsilon =1.4\cdot 10^{-8}$下迭代7次的结果（$x_0$的近似值为1.8954942711282465），用二次插值法在容错误差$\epsilon =0.001$下仅迭代3次就可算得与前者直至千万分位（$x_0$的近似值为1.895494265404134）均一致的结果。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！