设函数
f
(
x
)
f(x)
f(x),在
[
a
,
b
]
[a,b]
[a,b]上二阶连续可微且有唯一的最小值点
x
0
x_0
x0。由于
f
(
x
)
f(x)
f(x)是
[
a
,
b
]
[a,b]
[a,b]上的单峰函数,故
f
′
′
(
x
)
>
0
f''(x)>0
f′′(x)>0,
x
∈
(
a
,
b
)
x\in(a,b)
x∈(a,b)。对给定
x
k
∈
[
a
,
b
]
x_k\in[a,b]
xk∈[a,b],函数
f
(
x
)
f(x)
f(x)在
x
k
x_k
xk的近旁近似为
q
k
(
x
)
=
f
(
x
k
)
+
f
′
(
x
k
)
(
x
−
x
k
)
+
1
2
f
"
(
x
k
)
(
x
−
x
k
)
2
.
q_k(x)=f(x_k)+f'(x_k)(x-x_k)+\frac{1}{2}f"(x_k)(x-x_k)^2.
qk(x)=f(xk)+f′(xk)(x−xk)+21f"(xk)(x−xk)2.
由于
q
k
′
(
x
k
)
=
f
′
(
x
k
)
q'_k(x_k)=f'(x_k)
qk′(xk)=f′(xk),
q
k
′
′
(
x
k
)
=
f
′
′
(
x
k
)
>
0
q''_k(x_k)=f''(x_k)>0
qk′′(xk)=f′′(xk)>0,故
x
k
−
f
′
(
x
k
)
f
′
′
(
x
k
)
x_k-\frac{f'(x_k)}{f''(x_k)}
xk−f′′(xk)f′(xk)是
q
k
(
x
)
q_k(x)
qk(x)的驻点,也是唯一的极小值点。
我们用以下策略计算
f
(
x
)
f(x)
f(x)的最优解
x
0
x_0
x0的搜索序列:取
x
1
∈
[
a
,
b
]
x_1\in[a,b]
x1∈[a,b]充分接近
x
0
x_0
x0,从
k
=
1
k=1
k=1开始作迭代,若
f
′
(
x
k
)
=
0
f'(x_k)=0
f′(xk)=0,由
f
(
x
)
f(x)
f(x)在
[
a
,
b
]
[a,b]
[a,b]上的可微性及单峰性知,找到最优解
x
0
=
x
k
x_0=x_k
x0=xk,如上图(a)所示。否则,我们用
q
k
(
x
)
q_k(x)
qk(x)的唯一极小值点
x
k
−
f
′
(
x
k
)
f
′
′
(
x
k
)
x_k-\frac{f'(x_k)}{f''(x_k)}
xk−f′′(xk)f′(xk)作为第
k
+
1
k+1
k+1个迭代点
x
k
+
1
x_{k+1}
xk+1,即
x
k
+
1
=
x
k
−
f
′
(
x
k
)
f
′
′
(
x
k
)
.
x_{k+1}=x_k-\frac{f'(x_k)}{f''(x_k)}.
xk+1=xk−f′′(xk)f′(xk).
无论
x
k
x_k
xk是处于
f
(
x
)
f(x)
f(x)的上升区间如上图(b)或处于
f
(
x
)
f(x)
f(x)的下降区间如上图(c)所示,
x
k
+
1
x_{k+1}
xk+1都有望比
x
k
x_k
xk更接近
f
(
x
)
f(x)
f(x)的最小值点
x
0
x_0
x0。循环往复,直至遇到
f
′
(
x
k
)
=
0
f'(x_k)=0
f′(xk)=0。
按此策略的搜索算法称为牛顿方法,代码实现牛顿算法。
from scipy.optimize import OptimizeResult
def newton(fun, x1, gtol, **options):
xk=x1
f1,f2=der_scalar(fun,xk)
k=1
while abs(f1)>=gtol:
k+=1
xk-=f1/f2
f1,f2=der_scalar(fun,xk)
bestx=xk
besty=fun(bestx)
return OptimizeResult(fun=besty, x=bestx, nit=k)
程序的第2~12行定义实现牛顿算法的函数newton。参数fun表示函数
f
(
x
)
f(x)
f(x),x1表示初始点
x
1
x_1
x1,gtol表示容错误差
ε
\varepsilon
ε,options用来使minimize_scalar将x1和gtol等实际参数传递给newton。第3行将迭代点xk初始化为x1。第4行调用导数计算函数der_scalar详见博文《一元函数导数的数值计算》)计算
f
(
x
)
f(x)
f(x)的一、二阶导数置于f1和f2。第5行将迭代次数置为1。第6~9行的while循环执行迭代计算:第7行迭代次数k自增1,第8行计算迭代点
x
k
−
f
′
(
x
k
)
f
′
′
(
x
k
)
x_k-\frac{f'(x_k)}{f''(x_k)}
xk−f′′(xk)f′(xk)更新xk。第9行再次调用der_scalar计算更新后迭代点处的一、二阶导数。循环往复,直至
∣
f
′
(
x
k
)
∣
<
ε
|f'(x_k)|<\varepsilon
∣f′(xk)∣<ε。第10、11行分别计算最优解近似值bestx及最优解处函数近似besty。值12行返回用计算结果besty(最优解处函数值
f
(
x
0
)
f(x_0)
f(x0)的近似值)、bestx(最优解
x
0
x_0
x0的近似值)和k(迭代次数)构建OptimizeResult类型对象作为返回值。
例1 用上列程序定义的newton函数,计算函数
f
(
x
)
=
x
2
+
4
cos
x
f(x)=x^2+4\cos{x}
f(x)=x2+4cosx在
x
1
=
1.5
x_1=1.5
x1=1.5近旁的极小值点。
解:下列代码完成计算。
import numpy as np #导入numpy
from scipy.optimize import minimize_scalar #导入minimize_scalar
f=lambda x:x**2+4*np.cos(x) #设置目标函数
minimize_scalar(f,method=newton,options={'x1':1.5,'eps':1.48e-8}) #计算最优解
利用代码中的注释信息不难理解本程序。运行程序,输出
fun: 2.316808419788213
nit: 7
x: -1.8954942647118507
读者可将此运行结果与博文《一元函数搜索算法——二分法》中例1对同一函数使用二分法计算的结果相比。在相同的容错误差水平下,牛顿法的效率(仅用7次迭代)显然高于二分法(用27次迭代)。