设一元目标函数 f ( x ) f(x) f(x)在区间 [ a 0 , b 0 ] ⊆ R [a_0,b_0]\subseteq\text{R} [a0,b0]⊆R(其长度记为 λ \lambda λ)上为单峰函数,且在 ( a 0 , b 0 ) (a_0,b_0) (a0,b0)内连续可导,即其导函数 f ′ ( x ) f'(x) f′(x)在 ( a 0 , b 0 ) (a_0,b_0) (a0,b0)内连续。在此增强的条件下,可以加速迭代计算压缩区间的过程。仍然设置计算精度为 ε > 0 \varepsilon>0 ε>0。首次迭代,即 k = 1 k=1 k=1时,插入点取 a 1 ′ = a 0 + b 0 2 a'_1=\frac{a_0+b_0}{2} a1′=2a0+b0。若 f ′ ( a 1 ′ ) = 0 f'(a'_1)=0 f′(a1′)=0,即 a 1 ′ a'_1 a1′为 f ( x ) f(x) f(x)在 [ a 0 , b 0 ] [a_0,b_0] [a0,b0]的驻点,由单峰函数性质值, a 1 ′ a'_1 a1′为 f ( x ) f(x) f(x)在 ( a 0 , b 0 ) (a_0,b_0) (a0,b0)内维一的极小值点 x 0 x_0 x0,停止迭代。否则,若 f ′ ( a 1 ′ ) > 0 f'(a'_1)>0 f′(a1′)>0,意味着 f ( x ) f(x) f(x)沿 a 1 ′ a'_1 a1′的左边方向下降(见下图(a)),故 x 0 ∈ ( a 0 , a 1 ′ ) x_0\in(a_0,a'_1) x0∈(a0,a1′)。取压缩区间 [ a 1 , b 1 ] = [ a 0 , a 1 ′ ] [a_1,b_1]=[a_0,a'_1] [a1,b1]=[a0,a1′]。若 f ′ ( a 1 ′ ) < 0 f'(a'_1)<0 f′(a1′)<0,则 f ( x ) f(x) f(x)沿 a 1 ′ a'_1 a1′的右边方向下降(见下图(b)),必有 x 0 ∈ ( a 1 ′ , b 0 ) x_0\in(a'_1,b_0) x0∈(a1′,b0),取压缩区间为 [ a 1 , b 1 ] = [ a 1 ′ , b 0 ] [a_1,b_1]=[a'_1,b_0] [a1,b1]=[a1′,b0]。无论哪种情形, [ a 1 , b 1 ] [a_1,b_1] [a1,b1]的长度为 λ / 2 \lambda/2 λ/2,即压缩比为 0.5 0.5 0.5。用此策略,继续迭代,直至得到 a k ′ a'_k ak′,使得 f ′ ( a k ′ ) = 0 f'(a'_k)=0 f′(ak′)=0或 λ / 2 k < ε \lambda/2^k<\varepsilon λ/2k<ε,停止迭代, a k + b k 2 \frac{a_k+b_k}{2} 2ak+bk即为最小值点 x 0 x_0 x0的近似值。这一方法由于每次都将当前区间作对分,故称为“二分法”。回顾黄金分割法(详见博文《一元函数搜索算法——黄金分割法》),迭代计算压缩区间的压缩比为 1 − ρ = 1 − 0.382 = 0.618 > 0.5 1-\rho=1-0.382=0.618>0.5 1−ρ=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+=1if 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 ) f(x) f(x),bracket表示初始单峰区间 ( a 0 , b 0 ) (a_0,b_0) (a0,b0),gtol表示容错误差 ε \varepsilon ε,options用来使minimize_scalar将gtol等实际参数传递给bisection。第3行读取单峰区间(详见博文《连续函数的单峰区间计算》)左右端点a0和b0。第4行计算单峰区间的中点为a1。第5行调用der_scalar函数(详见博文《一元函数导数的数值计算》)计算 f ( x ) f(x) 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 f(x)=x^2+4\cos{x} f(x)=x2+4cosx在 x = 1 x=1 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(x) f(x)为f。第4行调用myBracket函数计算 f ( x ) f(x) f(x)在 x = 1 x=1 x=1近旁的单峰区间。第5行调用scipy.optimize的minimize_scalar函数,传递f、bracket和bisection函数,计算 f ( x ) f(x) f(x)在 x = 1 x=1 x=1近旁的局部最优解。运行程序,输出
fun: 2.3168084197882135
nit: 26
x: 1.895494265556336
bisection以容错误差 ε = 1.48 × 1 0 − 8 \varepsilon=1.48\times10^{-8} ε=1.48×10−8,迭代26次,算得最优解近似值为1.895494265556336,最优解处函数近似值为2.3168084197882135。
例2 物资需从位于陆地的城市 A A A运送到位于水中的海岛 B B B,假定各点间距离如题图中所示,且物资在陆地及水中的运输速度分别为1和 1 2 \frac{1}{2} 21。试确定海岸线上码头建造位置 x x x,使得物资运输时间最短。
解:根据题意,算得目标函数(即物资运送时间) f ( x ) = 1 + x 2 + 2 1 + ( 2 − x ) 2 f(x)=\sqrt{1+x^2}+2\sqrt{1+(2-x)^2} f(x)=1+x2+21+(2−x)2, 0 ≤ x ≤ 2 0\leq x\leq2 0≤x≤2。解析方法解决此问题,需算得其导数
f ′ ( x ) = x 1 + x 2 − 2 ( 2 − x ) 1 + ( 2 − x ) 2 . f'(x)=\frac{x}{\sqrt{1+x^2}}-\frac{2(2-x)}{\sqrt{1+(2-x)^2}}. f′(x)=1+x2x−1+(2−x)22(2−x).
令其为0,算得驻点 x 0 x_0 x0。然后根据二阶导数 f ′ ′ ( x 0 ) f''(x_0) f′′(x0)的符号,判断 x 0 x_0 x0是否为极小值点。为求驻点 x 0 x_0 x0,需解方程 x 1 + x 2 = 2 ( 2 − x ) 1 + ( 2 − x ) 2 \frac{x}{\sqrt{1+x^2}}=\frac{2(2-x)}{\sqrt{1+(2-x)^2}} 1+x2x=1+(2−x)22(2−x)。这将导致解高次方程
3 x 4 − 12 x 3 + 15 x 2 − 16 x + 16 = 0. 3x^4-12x^3+15x^2-16x+16=0. 3x4−12x3+15x2−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 ) f(x) f(x)。第4行调用myBracket函数计算0附近包围 f ( x ) f(x) f(x)最优解 x 0 x_0 x0的区间bracket。 第5~6行调用minimize_scalar函数,传递f、bracket、bisection以及提供给bisection的eps等参数,用二分法计算 f ( x ) f(x) f(x)的最优解近似值。运行程序,输出
fun: 4.037643276202614nit: 23x: 1.5382642555236816
意味着最优解 x 0 x_0 x0的近似值为1.5383, f ( x 0 ) f(x_0) f(x0)的近似值为4.0377。