证明梯度的运算法则（Python数学编程）

7.6 用梯度上升法求全局最大值

有时候我们只是想求一个函数的全局最大值，而不是所有的局部和全局最大值和最小值。例如，我们希望找到一个抛射角度使球能在水平方向上飞出最远的距离。接下来我们学习一种新的更实用的方法来解决这一类问题。 这种方法仅使用一阶导数，所以它只适用于存在一阶导数的函数。 这种方法称为梯度上升法(gradient ascent method)，是一种迭代求全局最大值的方法。由于梯度上升法需要大量的计算，因此它非常适合程序化求解而不是手工求解。下面我们以求解抛射角度问题为例来说明该方法。在第2章，我们推导了表达式:

来计算以角度\theta和速度u抛射的抛物运动中的物体的飞行时间。抛射射程R是抛射物体经过的总的水平距离，由式子给出，其中是初始速度的水平分量，等于。将和的公式代入，得到表达式：

下图展示了在0到90度之间取值时，每个角度对应的射程（飞行距离）。

从图中可以看出，最大射程在投射角度为45度附近获得。现在学习如何使用梯度上升法对进行数值求解。

梯度上升法是一种迭代方法：从的一个初值开始，例如0.001或，并逐渐接近对应于最大射程距离的值。

逐步接近的方程如下：

其中，表示步长，是R关于θ的导数。设定后，按下述步骤进行:

(1)使用之前的公式计算

(2)若的绝对值大于设定值ε，则定义并返回到步骤1,否则执行步骤3。

(3)是R取最大值时对应的θ的近似值。

这里的epsilon (ε )值的大小决定算法何时停止迭代。我们将在后面“步长和epsilon的角色"一节继续讨论。 以下的grad_ascent()函数实现了梯度上升算法。参数x0是迭代开始时的变量初始值，f1x是我们所要求的最大值的函数的导数，x是与函数自变量对应的Symbol对象。

 '''
 Use gradient ascent to find the angle at which the projectile
 has maximum range for a fixed velocity, 25 m/s
 '''
 import math
 from sympy import Derivative, Symbol, sin
 
 def grad_ascent(x0, f1x, x):
     epsilon = 1e-6      # ①
     step_size = 1e-4    # ②
     x_old = x0      ·   # ③
     x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()   # ③
     while abs(x_old - x_new) > epsilon:         # ⑤
         x_old = x_new
         x_new = x_old + step_size * f1x.subs({x: x_old}).evalf()
     return x_new
 def find_max_theta(R, theta):   # ⑥
     # Calculate the first derivative
     R1theta = Derivative(R, theta).doit()
     theta0 = 1e-3
     theta_max = grad_ascent(theta0, R1theta, theta)
     return theta_max        # ⑦
 if __name__ == '__main__':
     g = 9.8
     # Assuse initial velocity
     u = 25
     # Expression for range
     theta = Symbol('theta')
     R = u**2*sin(2*theta)/g     # ⑧
     
     theta_max = find_max_theta(R, theta)    # ⑨
     print('Theta: {0}'.format(math.degrees(theta_max)))
     print('Maximum Range : {0}'.format(R.subs({theta:theta_max})))

我们在①处设置epsilon值为1e-6，在②处设置步长为1e-4，epsilon 值必须始终是一个接近于0的非常小的正数，所选步长应该会使变量在每一次迭代中小幅度增加。epsilon值和步长的选择将在后面“步长和epsilon的角色”一节中更进一步讨论。

我们在③处设置x_old 为x0，并在④处首次计算x_new。我们使用subs()函数，以实现用x_old的值替换变量，然后使用evalf()函数计算数值结果。如果绝对值abs(x_old-x_new) 大于epsilon，则⑤处的while循环将继续执行，并按照梯度上升算法中的第1步和第2步持续更新x_old和x_new的值。一旦退出了循环，即abs(x_old-x_ new)<= epsilon,程序返回x_new，该值即为与函数最大值相对应的 变量值。

我们在⑥处定义find_max_theta()函数， 在这个函数中，我们计算R的一阶导数；创建一个标签theta0，并将其设置为1e-3；调用grad_ascent()函数， 并以theta0和R1theta这两个值作为输入参数，以Symbol对象theta 为第三个输入参数。一旦得到对应于最大函数值(theta_max)的θ值，在⑦中返回。 最后我们在⑧处创建表示水平射程的表达式，设置初始速度u=25,以及与角度θ对应的Symbol对象theta。然后在⑨中以R和theta为输入参数调用find_max_theta0函数。 当运行这个程序时，将看到如下输出：

 Theta: 44.997815081691805
 Maximum Range : 63.7755100185965

7.6.1 梯度上升的通用程序

我们可以稍微修改上一节的程序，使其成为梯度上升法的通用程序：

 '''
 Use gradient ascent to find the maximum value of a
 single_variable function
 '''
 import math
 from sympy import Derivative, Symbol, sin
 from sympy import sympify
 from sympy.core import SympifyError
 
 def grad_ascent(x0, f1x, x):
     epsilon = 1e-6
     step_size = 1e-4
     x_old = x0
     x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()
     while abs(x_old - x_new) > epsilon:
         x_old = x_new
         x_new = x_old + step_size * f1x.subs({x: x_old}).evalf()
     return x_new
 if __name__ == '__main__':
     f = input('Enter a function in one variable: ')
     var = input('Enter the variable to fiffernetiate with respect to: ')
     var0 = float(input('Enter the initial value of the variable: '))
     try:
         f = sympify(f)
     except SympifyError:
         print('Invalid function enterd')
     else:
         var = Symbol(var)   # ①
         d = Derivative(f, var).doit()   # ②
         var_max = grad_ascent(var0,d, var)  # ③
         print('{0}: {1}'.format(var.name, var_max))
         print('Maximum value : {0}'.format(f.subs({var:var_max})))

在此程序中，grad_ascent()函 数保持不变，然而现在程序将提示用户输入一个函数、函数的自变量以及变量初始值(梯度上升算法从这里开始)。一旦确定 SymPy可以识别用户的输入，我们在①处创建一个与变量相对应的Symbol对象，在②处计算关于该变量的一阶导数，并将这三个量(var0， d, var)作为输入参数调用grad_ascent()函数，在③处返回最大值。

下面是一次运行结果：

 Enter a function in one variable: 25*25*sin(2*theta)/9.8
 Enter the variable to fiffernetiate with respect to: theta
 Enter the initial value of the variable: 0.001
 theta: 0.785360029379083
 Maximum value : 63.7755100185965

此处输入的函数与我们上一次实现梯度上升法时的相同，θ值以弧度为单位输出。

下面是程序的另一个结果，这次是计算cosy函数的最大值：

 Enter a function in one variable: cos(y)
 Enter the variable to fiffernetiate with respect to: y
 Enter the initial value of the variable: 0.01
 y: 0.00999900001666658
 Maximum value : 0.999950010415832

该程序也适用于cos(y)+k这类函数，其中k是常数：

 Enter a function in one variable: cos(y) + K
 Enter the variable to fiffernetiate with respect to: y
 Enter the initial value of the variable: 0.01
 y: 0.00999900001666658
 Maximum value : K + 0.999950010415832

然而，该程序不适用于cos(ky)这类函数，因为该函数的一阶导数-ksin(ky)仍包含常数k，而Sympy不知道其数值。因此，Sympy不能执行梯度上升算法中的一个关键步骤，abs(x_old-x_new)>epsilon的比较。

7.6.2 关于初始值的附加证明

开始迭代梯度上升法时的变量初始值在算法中起着非常重要的作用。考虑之前使用过的函数x^5-30x^3+50x，这里我们使用梯度上升通用程序来计算最大值：

 Enter a function in one variable: x**5 - 30*x**3 + 50*x
 Enter the variable to fiffernetiate with respect to: x
 Enter the initial value of the variable: -2
 x: -4.17445116397103
 Maximum value : 705.959460322318

当找到最近的峰值时，梯度上升算法就停止了，但最近的峰值并不总是全局最大值。在这个例子中，如果你从初始值-2开始出发，程序停止是所在的峰值恰好是对应于设置的定义域中的全局最大值（大约为706）.为进一步验证，我们尝试一个不同的初始值：

 Enter a function in one variable: x**5 - 30*x**3 + 50*x
 Enter the variable to fiffernetiate with respect to: x
 Enter the initial value of the variable: 0.5
 x: 0.757452532565767
 Maximum value : 25.0846622605419

在这种情形下，梯度上升算法停止时所得到的最近峰值并不是函数的全局最大值。

7.6.3 步长和epsilon的角色

在梯度上升算法中，变量的下一个值的计算公式为，其中λ表示步长。步长决定了下一步的距离，它应该被设置得很小，以避免越过峰值。换句话说，如果x的当前值接近于函数的最大值点，那么下一步不应该越过这个峰值，否则算法将失效。另外，非常小的步长值将花费更多的计算时间。之前我们使用了固定的步长，但这个值并不适用于所有函数。

epsilon的值(ε )决定何时停止算法的迭代，它应该是一个足够小的数，小到我们确信x的值不再变化。我们希望一阶导数f'(x)在最大值点处为0，并且在理想情况下(参见梯度上升算法中的第2步)。然而，由于数值计算的不精确性，我们并不能精确地得到差值0,因此应该选取一个接近于0的epsilon值，它从实际的角度告诉我们x的值不再改变。我们在之前所有的函数中将epsilon设定为，这个值虽然足够小，也适用于那些f'(x)= 0有解的函数，例如sin(x),但可能并不适用于其他函数。因此，最好在最后验证最大值，以保证其正确性，并在需要时相应地调整epsilon的值。

梯度上升算法中的第2步也意味着，要让算法终止，方程式f'(x)=0必须有一个解，而对于例如或log(x)这类函数则不是这样。如果你把这类函数输入到之前的程序中，程序将无法计算出解，程序将持续运行。针对这种情形，我们可以通过检查f'(x)=0是否有解来改进梯度上升算法使其更有效。以下是改进的程序:

 '''
 Use gradient ascent to find the maximum value of a
 single_variable function. This also checks for the existence
 of a solution for the equation f'(x) = 0
 '''
 import math
 from sympy import Derivative, Symbol, solve
 from sympy import sympify
 from sympy.core import SympifyError
 
 def grad_ascent(x0, f1x, x):
     # Check if f1x=0 has a solution
     if not solve(f1x):
         print('Cannot continue, solution for {0} = 0 does not exist'.format(f1x))
         return
     epsilon = 1e-6
     step_size = 1e-4
     x_old = x0
     x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()
     while abs(x_old - x_new) > epsilon:
         x_old = x_new
         x_new = x_old + step_size * f1x.subs({x: x_old}).evalf()
     return x_new
 if __name__ == '__main__':
     f = input('Enter a function in one variable: ')
     var = input('Enter the variable to fiffernetiate with respect to: ')
     var0 = float(input('Enter the initial value of the variable: '))
     try:
         f = sympify(f)
     except SympifyError:
         print('Invalid function enterd')
     else:
         var = Symbol(var)
         d = Derivative(f, var).doit()
         var_max = grad_ascent(var0,d, var)
         if var_max:
             print('{0}: {1}'.format(var.name, var_max))
             print('Maximum value : {0}'.format(f.subs({var:var_max})))
 

在这个该井好的grad_ascent()函数中，我们在①处调用SymPy的solve()函数来判断方程式f'(x) = 0是否有解，方程式在本例中为f1x。如果无解，则输出提示信息并返回。另一个改进在②处__main__模块。我们检查grad_ascent()函数是否成功的返回了结果，如果是，则继续输出函数的最大值以及相应的自变量值。

这些改进使程序可以处理例如e^x或log(x)这类函数：

 Enter a function in one variable: log(x)
 Enter the variable to fiffernetiate with respect to: x
 Enter the initial value of the variable: 0.1
 Cannot continue, solution for 1/x = 0 does not exist

对于，你将看到同样的结果。

梯度下降算法

梯度上升算法的逆运算是沿着梯度下降，这是一种求函数最小值的办法，它与梯度上升算法类似，但不是沿着函数“向上爬“，而是”向下爬“。本章的挑战就会讨论两种算法的区别，并实现梯度下降算法。