【Python讲优化】S06E09 利用最速下降法求多元函数的极值

我们这一讲将要介绍的最速下降法,是梯度法的一种改进实现。

1.SymPy库的介绍

在具体进行算法介绍之前,我们先花一点时间来专门谈谈$python$的$SymPy$库,$SymPy$库是$Python$的数学符号计算库,用它可以进行数学表达式的符号推导和计算,可以很方便的进行公式推导、级数展开、积分、微分以及解方程等重要运算。

可能大家对于符号计算这个名字感觉有些陌生和奇怪,我们下面结合例子来慢慢熟悉他:

1.1.符号导入

我们举一个最典型的例子:欧拉公式。

$e^{i\pi}+1=0$,这里面有自然对数$e$,虚数$i$,圆周率$\pi$等数学符号,如果我们想用$python$来描述这个等式,运用之前的知识也可以办到,只不过相对而言比较麻烦,这里如果使用$SymPy$库就会非常简单,他可以一次性将这些符号都导入进来,并完成公式的计算:

代码片段:

from sympy import *

print("E**(I*pi)+1={}".format(E**(I*pi)+1))

运行结果:

E**(I*pi)+1=0

1.2.自定义符号

上面的欧拉公式中,我们所用的符号$e$、$i$、$\pi$都是已有的常量,如果我们的目标是实现一个带有自变量$x$的表达式,该如何处理?这里,我们需要自己定义新加入的变量$x$,我们看一个例子:验证$e^{ix}$这个表达式展开为实部和虚部的形式:$e^{ix}=isinx+cosx$。

代码片段:

from sympy import *

x = Symbol('x', real=True)
print(expand(E**(I*x), complex=True))

运行结果:

I*sin(x) + cos(x)

在这个例子中,我们利用$Symbol$函数定义了一个实变量符号$x$,并基于他定义了一个新的表达式$e^{ix}$,然后我们使用$expand$方法在复数的范围内将$e^{ix}$展开为了实部+虚部的形式:$isin(x) + cos(x)$。

1.3.任意阶数的泰勒展开

之前我们学习过函数的泰勒展开的有关知识,实际上利用$SymPy$库对指定函数进行任意阶数的泰勒展开是非常容易的,我们下面分别对$sin(x)$和$cos(x)$进行$10$阶泰勒展开:

代码片段:

from sympy import *

x = Symbol('x', real=True)
sin_s = series(sin(x), x, 0, 10)
cos_s = series(cos(x), x, 0, 10)
print('sin(x)={}'.format(sin_s))
print('cos(x)={}'.format(cos_s))

运行结果:

sin(x)=x - x**3/6 + x**5/120 - x**7/5040 + x**9/362880 + O(x**10)
cos(x)=1 - x**2/2 + x**4/24 - x**6/720 + x**8/40320 + O(x**10)

可以看出,得到的泰勒展开式的结果是非常准确的,还包含了无穷小量的表达。

1.4.自定义符号替换

在实际的程序当中,我们经常需要对表达式中的变量进行名称替换,最简单的,比如把变量名$x$替换成变量名$y$,至于说用处是什么,我们在后面的例子中会见到,我们先看看如何使用,试着把上面的$cos(x)$的级数展开式中的变量$x$替换成$y$:

代码片段:

from sympy import *

x = Symbol('x', real=True)
y = Symbol('y', real=True)
cos_s = series(cos(x), x, 0, 10)
print('cos(x)={}'.format(cos_s))
cos_s = cos_s.subs(x, y)
print('cos(y)={}'.format(cos_s))

运行结果:

cos(x)=1 - x**2/2 + x**4/24 - x**6/720 + x**8/40320 + O(x**10)
cos(y)=1 - y**2/2 + y**4/24 - y**6/720 + y**8/40320 + O(y**10)

这样就完成了变量名称的替换工作。

1.5.求导与微分

我们来看看如何利用$SymPy$库里的$diff$函数对函数进行求导运算。

这里,我们分别举三个例子,看看导数$\frac{d}{dx}sin(2x)$、二阶导$\frac{d^2}{dx^2}(x^2+2x+1)$、偏导数$\frac{\partial^2}{\partial x \partial y}(x^2y^2+2x^3+y^2)$的符号求解方法:

代码片段:

from sympy import *

x = Symbol('x', real=True)
y = Symbol('y', real=True)

print(diff(sin(2*x),x))
print(diff(x**2+2*x+1,x,2))
print(diff(x**2*y**2+2*x**3+y**2,x,1,y,1))

运行结果:

2*cos(2*x)
2
4*x*y

从结果中,我们很轻松的得到了三种导数结果的符号表达式:

$\frac{d}{dx}sin(2x)=2cos(2x)$

$\frac{d^2}{dx^2}x^2+2x+1=2$

$\frac{\partial^2}{\partial x \partial y}x^2y^2+2x^3+y^2=4xy$

1.6.解方程

$SymPy$库当中的$solve$方法可以用来解方程,这个就非常方便了,我们看两个简单的例子:

一次方程:$x+1=0$

二次方程:$x^2+3x+2=0$

代码片段:

from sympy import *

x = Symbol('x', real=True)

f1 = x + 1
f2 = x**2+3*x+2

print(solve(f1))
print(solve(f2))

运行结果:

[-1]
[-2, -1]

得到的结果是一个列表,里面包含了方程所有的解。

1.7.表达式求值

这个用法其实非常关键,前面我们用各种方式生成了表达式,如何对我们自定义的符号变量赋值,并计算出表达式的值,看下面的具体实现:

我们求当$x=2$的时候,函数$f(x)=x^2+3x+2$的取值:

很简单,表达式求值的本质也是符号变量的替换,就是把自定义的符号变量$x$替换成目标取值$2$。

代码片段:

from sympy import *
x = Symbol('x', real=True)
f = x**2+3*x+2

print(f.subs(x,2))

运行结果:

12.0000000000000

之所以花了这么多的篇幅来介绍一个新的$python$库用法,目的是在接下来的最速下降法的实现中使用他们,以简化我们的代码实现。

2.最速下降法的核心思想

最速下降法是梯度下降法的一种实现形式,每一步都是沿着负梯度的方向迭代前进,但是同上一讲里介绍的方法不同之处在于学习率的选择方式:

在基础的梯度下降法当中,学习率$\lambda_k$是固定的,即$p_{k+1}=p_k-\nabla f(p_k)\lambda_k$。

但是最速下降法则不然,搜索方向仍然是当前点$p_k$的负梯度方向$-\nabla f(p_k)$,但是迭代的目标点是需要满足在这个搜索方向上使得函数$f$取得该方向上的极小值,即:

$p_{k+1}=p_k-\alpha_k\nabla f(p_k)$,使得函数取值$f(p_{k+1})$在搜索的一维方向上取得极小值。

更进一步的,在$p_k$点进行迭代时,由于此时$p_k$和$\nabla f(p_k)$都是已知的,实际上的工作就是锁定$\alpha_k$的取值,使得函数$f(p_k-\alpha_k\nabla f(p_k))$取得极小值。此时在寻找函数极小值的过程中,仅仅只有一个未知数,那就是$\alpha_k$,因此我们面对的问题实际上就是一个一元函数极值点的搜索问题,方法有很多,这在我们前面的内容中已经详细讲解过了,相信对大家没有难度。

每次迭代的过程都是在当前迭代点的负梯度方向上寻找使得函数取得极小值点的$\alpha_k$值,然后依次迭代出下一个迭代点,最终直至满足预设阈值条件,迭代停止。

3.算法步骤

这里,我们来总结一下最速下降法的步骤:

首先,我们从一个初始的迭代点$p_0$开始;

每一轮迭代,对于当前点$p_k$,计算当前点的梯度$\nabla f(p_k)$,如果梯度的模长$|\nabla f(p_k)|$小于预设的阈值$\epsilon$,则停止迭代,当前点$p_k$即为函数的极小值点;

否则:寻找使得$f(p_k-\alpha_k\nabla f(p_k))$取得极小值的$\alpha_k$,然后迭代出新的取值点:$p_{k+1}=p_k-\alpha_k\nabla f(p_k)$;

循环往复的进行上述迭代过程。

4.代码演示

top Created with Sketch.