Python科学计算(第2版)
上QQ阅读APP看书,第一时间看更新

3.2.3 计算函数局域最小值

optimize库还提供了许多求函数最小值的算法:Nelder-Mead、Powell、CG、BFGS、Newton-CG、L-BFGS-B等。下面我们用一个实例观察这些优化函数是如何找到函数的最小值的。在本例中,要计算最小值的函数f(x,y)为:

f(x,y)=(1-x)2 +100(y-x2 )2

这个函数叫作Rosenbrock函数,它经常用来测试最小化算法的收敛速度。它有一个十分平坦的山谷区域,收敛到此山谷区域比较容易,但是在山谷区域搜索到最小点则比较困难。根据函数的计算公式不难看出此函数的最小值是0,在(1,1)处。

为了提高运算速度和精度,有些算法带有一个fprime参数,它是计算目标函数f( )对各个自变量的偏导数的函数。f(x,y)对变量x和y的偏导函数为:

而Newton-CG算法则需要计算海森矩阵,它是一个由自变量为向量的实值函数的二阶偏导数构成的方块矩阵,对于函数f(x1 , x2 , …, xn ),其海森矩阵的定义如下:

对于本例来说,海森矩阵为一个二阶矩阵:

下面使用各种最小值优化算法计算f(x,y)的最小值,根据其输出可知有些算法需要较长的收敛时间,而有些算法则利用导数信息更快地找到最小点。

    def target_function(x, y):
        return (1-x)**
2 + 100*
(y-x**
2)**
2    
    
    class TargetFunction(object):
    
        def __init__(self):
            self.f_points = [ ]
            self.fprime_points = [ ]
            self.fhess_points = [ ]
    
        def f(self, p):
            x, y = p.tolist( )
            z = target_function(x, y)
            self.f_points.append((x, y))
            return z
    
        def fprime(self, p):
            x, y = p.tolist( )
            self.fprime_points.append((x, y))
            dx = -2 + 2*
x - 400*
x*
(y - x**
2)
            dy = 200*
y - 200*
x**
2
            return np.array([dx, dy])
    
        def fhess(self, p):
            x, y = p.tolist( )
            self.fhess_points.append((x, y))
            return np.array([[2*
(600*
x**
2 - 200*
y + 1), -400*
x],
    [-400*
x, 200]])
    
    def fmin_demo(method):
        target = TargetFunction( )
        init_point =(-1, -1)
        res = optimize.minimize(target.f, init_point, 
                          method=method,
                          jac=target.fprime,
    hess=target.fhess)
        return res, [np.array(points) for points in 
                    (target.f_points, target.fprime_points, target.fhess_points)]
    
    methods = ("Nelder-Mead", "Powell", "CG", "BFGS", "Newton-CG", "L-BFGS-B")
    for method in methods:
        res, (f_points, fprime_points, fhess_points) = fmin_demo(method)
        print "{:12s}: min={:12g}, f count={:3d}, fprime count={:3d}, "\
    "fhess count={:3d}".format(
    method, float(res["fun"]), len(f_points), 
    len(fprime_points), len(fhess_points))
    Nelder-Mead : min= 5.30934e-10, f count=125, fprime count=  0, fhess count=  0
    Powell      : min=           0, f count= 52, fprime count=  0, fhess count=  0
    CG          : min=  7.6345e-15, f count= 34, fprime count= 34, fhess count=  0
    BFGS        : min= 2.31605e-16, f count= 40, fprime count= 40, fhess count=  0
    Newton-CG   : min= 5.22666e-10, f count= 60, fprime count= 97, fhess count= 38
    L-BFGS-B    : min=  6.5215e-15, f count= 33, fprime count= 33, fhess count=  0

图3-4显示了各种优化算法的搜索路径,图中用圆点表示调用f( )时的坐标点,圆点的颜色表示调用顺序;叉点表示调用fprime( )时的坐标点。图中用图像表示二维函数的值,值越大则颜色越浅,值越小则颜色越深。为了更清晰地显示函数的山谷区域,图中显示的实际上是通过对数函数log10( )对f(x,y)进行处理之后的结果。

图3-4 各种优化算法的搜索路径