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

设一元目标函数$f(x)$在区间$[a_0,b_0]\subseteq \text{R}$（其长度记为$\lambda$）上为单峰函数，且在$(a_0,b_0)$内连续可导，即其导函数$f^{\prime }(x)$在$(a_0,b_0)$内连续。在此增强的条件下，可以加速迭代计算压缩区间的过程。仍然设置计算精度为$\epsilon >0$。首次迭代，即$k=1$时，插入点取$a_1^{\prime }=\frac{a_0+b_0}{2}$。若$f^{\prime }(a_1^{\prime })=0$，即$a_1^{\prime }$为$f(x)$在$[a_0,b_0]$的驻点，由单峰函数性质值，$a_1^{\prime }$为$f(x)$在$(a_0,b_0)$内维一的极小值点$x_0$，停止迭代。否则，若$f^{\prime }(a_1^{\prime })>0$，意味着$f(x)$沿$a_1^{\prime }$的左边方向下降（见下图(a)），故$x_0\in (a_0,a_1^{\prime })$。取压缩区间$[a_1,b_1]=[a_0,a_1^{\prime }]$。若$f^{\prime }(a_1^{\prime })<0$，则$f(x)$沿$a_1^{\prime }$的右边方向下降（见下图(b)），必有$x_0\in (a_1^{\prime },b_0)$，取压缩区间为$[a_1,b_1]=[a_1^{\prime },b_0]$。无论哪种情形，$[a_1,b_1]$的长度为$\lambda /2$，即压缩比为$0.5$。用此策略，继续迭代，直至得到$a_k^{\prime }$，使得$f^{\prime }(a_k^{\prime })=0$或$\lambda /2^k<\epsilon$，停止迭代，$\frac{a_k+b_k}{2}$即为最小值点$x_0$的近似值。这一方法由于每次都将当前区间作对分，故称为“二分法”。回顾黄金分割法，迭代计算压缩区间的压缩比为$1-\rho =1-0.382=0.618>0.5$，所以二分法比黄金分割法效率稍有提高。
将上述算法实现为如下的Python函数：

from scipy.optimize import OptimizeResult
bisection(fun,bracket,gtol,**options):
    a0,b0=bracket						#单峰区间端点
    a1=(a0+b0)/2						#中点
    f1,_=der_scalar(fun,a1)				#中点处导数
    k=1									#迭代次数
    while abs(f1)>gtol:					#重复迭代
        k+=1
        if f1<0:						#导数为负
            a0=a1						#修改左端点
        else:							#导数为正
            b0=a1						#修改右端点
        a1=(a0+b0)/2					#更新中点
        f1,_=der_scalar(fun,a1)			#更新中点处导数
    bestx=a1    
    besty=fun(bestx)
    return OptimizeResult(fun=besty, x=bestx, nit=k)

程序中第2-17行定义用二分策略计算目标函数局部最优解的bisection函数。参数fun表示目标函数$f(x)$，bracket表示初始单峰区间$(a_0,b_0)$，gtol表示容错误差$\epsilon$，options用来使minimize_scalar将gtol等实际参数传递给bisection。第3行读取单峰区间（详见博文《连续函数的单峰区间计算》）左右端点a0和b0。第4行计算单峰区间的中点为a1。第5行调用der_scalar函数（详见博文《一元函数导数的数值计算》）计算$f(x)$在区间中点处的导数为f1。第6行将迭代次数k初始化为1。第7-14行的while循环执行重复迭代，直至当前区间中点处的导数绝对值接近0。循环体中，第8行迭代次数k自增1，第9~12行的if-else分支根据当前中点处导数的符号，决定下一次迭代的单峰区间。第13行更新当前区间中点，第14行更新中点处导数。第15、16行用a1分别设置最优解bestx和最优解处函数值besty。第17行用besty、bestx、k构造OptimizeResult对象并返回。
例1 用bisection函数计算函数$f(x)=x^2+4\cos x$在$x=1$附近的局部最优解。
解：下列代码计算本例。

import numpy as np
from scipy.optimize import minimize_scalar
f=lambda x:x**2+4*np.cos(x)
bracket=myBracket(f,1)
res=minimize_scalar(f,bracket,method=bisection, options={'eps':1.48e-8})
print(res)

程序很简单，第3行定义目标函数$f(x)$为f。第4行调用myBracket函数计算$f(x)$在$x=1$近旁的单峰区间。第5行调用scipy.optimize的minimize_scalar函数，传递f、bracket和bisection函数，计算$f(x)$在$x=1$近旁的局部最优解。运行程序，输出

fun: 2.3168084197882135
nit: 26
x: 1.895494265556336

bisection以容错误差$\epsilon =1.48\times 10^{-8}$，迭代26次，算得最优解近似值为1.895494265556336，最优解处函数近似值为2.3168084197882135。
例2 物资需从位于陆地的城市$A$运送到位于水中的海岛$B$，假定各点间距离如题图中所示，且物资在陆地及水中的运输速度分别为1和$\frac{1}{2}$。试确定海岸线上码头建造位置$x$，使得物资运输时间最短。

解：根据题意，算得目标函数（即物资运送时间）$f(x)=\sqrt{1+x^2}+2\sqrt{1+(2-x{)}^2}$，$0\le x\le 2$。解析方法解决此问题，需算得其导数
$$f^{\prime }(x)=\frac{x}{\sqrt{1+x^2}}-\frac{2(2-x)}{\sqrt{1+(2-x{)}^2}}.$$
令其为0,算得驻点$x_0$。然后根据二阶导数$f^{\prime \prime }(x_0)$的符号，判断$x_0$是否为极小值点。为求驻点$x_0$，需解方程$\frac{x}{\sqrt{1+x^2}}=\frac{2(2-x)}{\sqrt{1+(2-x{)}^2}}$。这将导致解高次方程
$$3x^4-12x^3+15x^2-16x+16=0.$$
手算的工作量极大。下面我们用黄金分割法数值地计算这个问题，代码如下。

import numpy as np									#导入numpy
from scipy.optimize import minimize_scalar			#导入minimize_scalar
f=lambda x:np.sqrt(1+x**2)+2*np.sqrt(1+(2-x)**2)	#设置目标函数
bracket=myBracket(f, 0)								#计算包围最优解的区间
res=minimize_scalar(f, bracket,method=bisection,	#用二分法计算f(x)的最有接近似值
                    options={'eps':1.48e-8})
print(res)

程序第3行设置目标函数$f(x)$。第4行调用myBracket函数计算0附近包围$f(x)$最优解$x_0$的区间bracket。 第5~6行调用minimize_scalar函数，传递f、bracket、bisection以及提供给bisection的eps等参数，用二分法计算$f(x)$的最优解近似值。运行程序，输出

 fun: 4.037643276202614
 nit: 23
   x: 1.5382642555236816

意味着最优解$x_0$的近似值为1.5383，$f(x_0)$的近似值为4.0377。