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

设函数$f(x)$，在$[a,b]$上二阶连续可微且有唯一的最小值点$x_0$。由于$f(x)$是$[a,b]$上的单峰函数，故$f^{\prime \prime }(x)>0$，$x\in (a,b)$。对给定$x_k\in [a,b]$，函数$f(x)$在$x_k$的近旁近似为
$$q_k(x)=f(x_k)+f^{\prime }(x_k)(x-x_k)+\frac{1}{2}{f}^{\prime \prime }(x_k)(x-x_k{)}^2.$$
由于$q_k^{\prime }(x_k)=f^{\prime }(x_k)$，$q_k^{\prime \prime }(x_k)=f^{\prime \prime }(x_k)>0$，故$x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)}$是$q_k(x)$的驻点，也是唯一的极小值点。

我们用以下策略计算$f(x)$的最优解$x_0$的搜索序列：取$x_1\in [a,b]$充分接近$x_0$，从$k=1$开始作迭代，若$f^{\prime }(x_k)=0$，由$f(x)$在$[a,b]$上的可微性及单峰性知，找到最优解$x_0=x_k$，如上图(a)所示。否则，我们用$q_k(x)$的唯一极小值点$x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)}$作为第$k+1$个迭代点$x_{k+1}$，即
$$x_{k+1}=x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)}.$$
无论$x_k$是处于$f(x)$的上升区间如上图(b)或处于$f(x)$的下降区间如上图（c）所示，$x_{k+1}$都有望比$x_k$更接近$f(x)$的最小值点$x_0$。循环往复，直至遇到$f^{\prime }(x_k)=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)$，x1表示初始点$x_1$，gtol表示容错误差$\epsilon$，options用来使minimize_scalar将x1和gtol等实际参数传递给newton。第3行将迭代点xk初始化为x1。第4行调用导数计算函数der_scalar详见博文《一元函数导数的数值计算》）计算$f(x)$的一、二阶导数置于f1和f2。第5行将迭代次数置为1。第6~9行的while循环执行迭代计算：第7行迭代次数k自增1，第8行计算迭代点$x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)}$更新xk。第9行再次调用der_scalar计算更新后迭代点处的一、二阶导数。循环往复，直至$|f^{\prime }(x_k)|<\epsilon$。第10、11行分别计算最优解近似值bestx及最优解处函数近似besty。值12行返回用计算结果besty（最优解处函数值$f(x_0)$的近似值）、bestx（最优解$x_0$的近似值）和k（迭代次数）构建OptimizeResult类型对象作为返回值。
例1 用上列程序定义的newton函数，计算函数$f(x)=x^2+4\cos x$在$x_1=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次迭代）。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！