Python实现最速下降法(The steepest descent method)详细案例

最新推荐文章于 2023-10-30 11:21:57 发布

刘兴禄

最新推荐文章于 2023-10-30 11:21:57 发布

阅读量1.3w

点赞数 37

分类专栏：运筹优化理论与实践

本文链接：https://blog.csdn.net/HsinglukLiu/article/details/109524062

版权

运筹优化理论与实践专栏收录该内容

47 篇文章 164 订阅

订阅专栏

文章目录

最速下降法(The steepest descent method)

.[TOC].

最速下降法(The steepest descent method)

本文中的课件来自清华大学深圳国际研究生院，物流与交通学部张灿荣教授《高级运筹学》课程。

在无约束非线性函数的最优化中，最速下降法(The steepest descent method)是一个著名的基础算法。本文就来用一个实例来学习最速下降法。

最速下降法的原理

假设有一个多元非线性函数 $f(x_1, x_2, \cdots, x_n)$ ，其定义域为 $\in \mathbb{R}^n$ , 那么如何求该函数的最大值或者最小值呢？

当然，我们可以用求导的方法。但是有些时候直接求导并不方便，此时我们可以用最速下降法来求解。

最速下降法的基本思想就是：

从一个初始点开始，逐步沿着以当前点为基准，函数值变化最快的方向走，一直走到最优解为止。因此，在有了一个初始点以后，我们就需要决策以下两个事情：
（1）下一步要朝着什么方向走(方向);
（2）沿着该方向走多远(步长)。

具体如何选择方向和步长，见下面的PPT。
在这里插入图片描述
下图是一个很直观的例子，在当前点，沿着不同的方向走，函数值的变化速度是不同的，但是在等高线的切线的垂直方向，是变化最快的。

在这里插入图片描述
也就是说，我们得去找平行于该点梯度的方向，沿着这个方向(当为max问题)或者沿着这个方向的反方向(当为min问题)去更新当前位置。

就像上图一样，一步一步，最终走到最优解对应的点。

Python实现最速下降法实例

我们看下面的函数
$\begin{aligned} f(x_1,x_2)=&2x_1x_2 +2x_2-x_1^2-2x_2^2 , \\ &x_1, x_2 \in \mathbb{R^1} \end{aligned}$

给定初始点 $(0.5, 0.5)$ ，用最速下降法找到该函数的最大值。

我们用python中的sympy包来实现最速下降法求解上面的问题。

首先，我们来介绍我们将要用到的sympy包中的几个函数。结合jupyter notebook来给大家做个展示。ps.这个包还是挺好用的，结合jupyter notebook之后，可视化非常不错，所见即所得。

`sympy`包中用到的函数

构建符号变量和符号函数

from sympy import * 
x_1 = symbols('x_1')
x_2 = symbols('x_2') 
              
fun = 2 * x_1 * x_2 + 2 * x_2 - x_1**2 - 2 * x_2**2
fun

这个是用来构造两个符号变量 $x_1, x_2$ ，就像代数中用字母代替变量一样。然后可以定义出我们的函数
$\begin{aligned} f(x_1,x_2)=&2x_1x_2 +2x_2-x_1^2-2x_2^2 , \end{aligned}$

jupyter notebook中的显示效果是这样的

在这里插入图片描述
可以看到jupyter notebook中直接就显示出了数学公式格式的形式，这是因为jupyter notebook中内嵌了LaTeX相关支持包的缘故。总之这样可视化就非常不错。

对符号函数求导

下面我们来计算 $f(x_1,x_2)$ 的偏导数， $\frac{\partial f}{\partial x_1}$ ，代码很简单，用函数diff(函数, 变量)

grad_1 = diff(fun, x_1) 
grad_1

如下图

在这里插入图片描述

求函数值

有了符号函数，我们怎么知道自变量 $x_1, x_2$ 取具体值得时候，符号函数的取值呢？我们用函数函数.subs({变量1:变量1的取值, ....}).evalf()，具体代码为

fun_value = fun.subs({x_1:4, x_2: 2}).evalf()
fun_value

在这里插入图片描述
如果说只有一个变量已知，那就是下面的情况

fun_value = fun.subs({x_1:4, x_2: 2}).evalf()
fun_value

在这里插入图片描述
还是非常智能的。

求解方程的零点

为了寻找下降速度最快的方向，我们需要利用之前PPT中的方法去求解方程组。这里我们用到函数solve(函数，变量)

我们举一个例子，假设有函数
$\begin{aligned} g(x) = 4x + 5 \end{aligned}$
我们求解
$\begin{aligned} g(x) = 4x + 5 = 0 \end{aligned}$
在这里插入图片描述
OK，准备好了上述函数，我们就可以开心的把最速下降法写出来了。具体代码见接下来的部分。

Python实现最速下降法求解上述算例的完整代码

import math
from sympy import *

# define symbol variable
x_1 = symbols('x_1')
x_2 = symbols('x_2')

# define objective function 
fun = 2 * x_1 * x_2 + 2 * x_2 - x_1 ** 2 - 2 * x_2 ** 2
fun

# take derivative of x_1 and x_2 
grad_1 = diff(fun, x_1)
grad_2 = diff(fun, x_2)

# define parameters
MaxIter = 100  
epsilon = 0.0001

# define initial point 
x_1_value = 0.5
x_2_value = 0.5

iter_cnt = 0
current_step_size = 10000 

grad_1_value = (float)(grad_1.subs({x_1:x_1_value, x_2: x_2_value}).evalf())
grad_2_value = (float)(grad_2.subs({x_1:x_1_value, x_2: x_2_value}).evalf()) 

current_obj = fun.subs({x_1: x_1_value, x_2: x_2_value}).evalf()

print('itCnt: %2d  cur_point (%3.2f, %3.2f)   cur_Obj: %5.4f     grad_1: %5.4f     grad_2 : %5.4f     step_size : %5.4f' % (iter_cnt, x_1_value, x_2_value, current_obj, grad_1_value, grad_2_value, current_step_size)) 
    
# while (iter_cnt <= MaxIter and abs(grad_1_value) + abs(grad_2_value) >= epsilon):
while(abs(grad_1_value) + abs(grad_2_value) >= epsilon):  
    iter_cnt += 1
    # find the step size
    t = symbols('t')
    x_1_updated = x_1_value + grad_1_value * t
    x_2_updated = x_2_value + grad_2_value * t
    Fun_updated = fun.subs({x_1: x_1_updated, x_2: x_2_updated})
    grad_t = diff(Fun_updated, t)
    t_value = solve(grad_t, t)[0]  # solve grad_t == 0

    # update x_1_value and x_2_value
    grad_1_value = (float)(grad_1.subs({x_1: x_1_value, x_2: x_2_value}).evalf()) 
    grad_2_value = (float)(grad_2.subs({x_1: x_1_value, x_2: x_2_value}).evalf()) 

    x_1_value = (float)(x_1_value + t_value * grad_1_value)
    x_2_value = (float)(x_2_value + t_value * grad_2_value) 

    current_obj = fun.subs({x_1: x_1_value, x_2: x_2_value}).evalf()
    current_step_size = t_value

    print('itCnt: %2d  cur_point (%3.2f, %3.2f)   cur_Obj: %5.4f     grad_1: %5.4f     grad_2 : %5.4f     step_size : %5.4f' % (iter_cnt, x_1_value, x_2_value, current_obj, grad_1_value, grad_2_value, current_step_size))

运行结果如下

itCnt:  0  cur_point (0.50, 0.50)   cur_Obj: 0.7500     grad_1: 0.0000     grad_2 : 1.0000     step_size : 10000.0000
itCnt:  1  cur_point (0.50, 0.75)   cur_Obj: 0.8750     grad_1: 0.0000     grad_2 : 1.0000     step_size : 0.2500
itCnt:  2  cur_point (0.50, 0.75)   cur_Obj: 0.8750     grad_1: 0.5000     grad_2 : 0.0000     step_size : 0.0000
itCnt:  3  cur_point (0.75, 0.75)   cur_Obj: 0.9375     grad_1: 0.5000     grad_2 : 0.0000     step_size : 0.5000
itCnt:  4  cur_point (0.75, 0.75)   cur_Obj: 0.9375     grad_1: 0.0000     grad_2 : 0.5000     step_size : 0.0000
itCnt:  5  cur_point (0.75, 0.88)   cur_Obj: 0.9688     grad_1: 0.0000     grad_2 : 0.5000     step_size : 0.2500
itCnt:  6  cur_point (0.75, 0.88)   cur_Obj: 0.9688     grad_1: 0.2500     grad_2 : 0.0000     step_size : 0.0000
itCnt:  7  cur_point (0.88, 0.88)   cur_Obj: 0.9844     grad_1: 0.2500     grad_2 : 0.0000     step_size : 0.5000
itCnt:  8  cur_point (0.88, 0.88)   cur_Obj: 0.9844     grad_1: 0.0000     grad_2 : 0.2500     step_size : 0.0000
itCnt:  9  cur_point (0.88, 0.94)   cur_Obj: 0.9922     grad_1: 0.0000     grad_2 : 0.2500     step_size : 0.2500
itCnt: 10  cur_point (0.88, 0.94)   cur_Obj: 0.9922     grad_1: 0.1250     grad_2 : 0.0000     step_size : 0.0000
itCnt: 11  cur_point (0.94, 0.94)   cur_Obj: 0.9961     grad_1: 0.1250     grad_2 : 0.0000     step_size : 0.5000
itCnt: 12  cur_point (0.94, 0.94)   cur_Obj: 0.9961     grad_1: 0.0000     grad_2 : 0.1250     step_size : 0.0000
itCnt: 13  cur_point (0.94, 0.97)   cur_Obj: 0.9980     grad_1: 0.0000     grad_2 : 0.1250     step_size : 0.2500
itCnt: 14  cur_point (0.94, 0.97)   cur_Obj: 0.9980     grad_1: 0.0625     grad_2 : 0.0000     step_size : 0.0000
itCnt: 15  cur_point (0.97, 0.97)   cur_Obj: 0.9990     grad_1: 0.0625     grad_2 : 0.0000     step_size : 0.5000
itCnt: 16  cur_point (0.97, 0.97)   cur_Obj: 0.9990     grad_1: 0.0000     grad_2 : 0.0625     step_size : 0.0000
itCnt: 17  cur_point (0.97, 0.98)   cur_Obj: 0.9995     grad_1: 0.0000     grad_2 : 0.0625     step_size : 0.2500
itCnt: 18  cur_point (0.97, 0.98)   cur_Obj: 0.9995     grad_1: 0.0312     grad_2 : 0.0000     step_size : 0.0000
itCnt: 19  cur_point (0.98, 0.98)   cur_Obj: 0.9998     grad_1: 0.0312     grad_2 : 0.0000     step_size : 0.5000
itCnt: 20  cur_point (0.98, 0.98)   cur_Obj: 0.9998     grad_1: 0.0000     grad_2 : 0.0312     step_size : 0.0000
itCnt: 21  cur_point (0.98, 0.99)   cur_Obj: 0.9999     grad_1: 0.0000     grad_2 : 0.0312     step_size : 0.2500
itCnt: 22  cur_point (0.98, 0.99)   cur_Obj: 0.9999     grad_1: 0.0156     grad_2 : 0.0000     step_size : 0.0000
itCnt: 23  cur_point (0.99, 0.99)   cur_Obj: 0.9999     grad_1: 0.0156     grad_2 : 0.0000     step_size : 0.5000
itCnt: 24  cur_point (0.99, 0.99)   cur_Obj: 0.9999     grad_1: 0.0000     grad_2 : 0.0156     step_size : 0.0000
itCnt: 25  cur_point (0.99, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0156     step_size : 0.2500
itCnt: 26  cur_point (0.99, 1.00)   cur_Obj: 1.0000     grad_1: 0.0078     grad_2 : 0.0000     step_size : 0.0000
itCnt: 27  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0078     grad_2 : 0.0000     step_size : 0.5000
itCnt: 28  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0078     step_size : 0.0000
itCnt: 29  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0078     step_size : 0.2500
itCnt: 30  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0039     grad_2 : 0.0000     step_size : 0.0000
itCnt: 31  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0039     grad_2 : 0.0000     step_size : 0.5000
itCnt: 32  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0039     step_size : 0.0000
itCnt: 33  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0039     step_size : 0.2500
itCnt: 34  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0020     grad_2 : 0.0000     step_size : 0.0000
itCnt: 35  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0020     grad_2 : 0.0000     step_size : 0.5000
itCnt: 36  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0020     step_size : 0.0000
itCnt: 37  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0020     step_size : 0.2500
itCnt: 38  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0010     grad_2 : 0.0000     step_size : 0.0000
itCnt: 39  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0010     grad_2 : 0.0000     step_size : 0.5000
itCnt: 40  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0010     step_size : 0.0000
itCnt: 41  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0010     step_size : 0.2500
itCnt: 42  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0005     grad_2 : 0.0000     step_size : 0.0000
itCnt: 43  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0005     grad_2 : 0.0000     step_size : 0.5000
itCnt: 44  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0005     step_size : 0.0000
itCnt: 45  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0005     step_size : 0.2500
itCnt: 46  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0002     grad_2 : 0.0000     step_size : 0.0000
itCnt: 47  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0002     grad_2 : 0.0000     step_size : 0.5000
itCnt: 48  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0002     step_size : 0.0000
itCnt: 49  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0002     step_size : 0.2500
itCnt: 50  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0001     grad_2 : 0.0000     step_size : 0.0000
itCnt: 51  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0001     grad_2 : 0.0000     step_size : 0.5000
itCnt: 52  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0001     step_size : 0.0000
itCnt: 53  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0000     grad_2 : 0.0001     step_size : 0.2500
itCnt: 54  cur_point (1.00, 1.00)   cur_Obj: 1.0000     grad_1: 0.0001     grad_2 : 0.0000     step_size : 0.0000

这样就完成了。

当然了，在这个例子中，从第3步左右的迭代开始，后续的点就非常近了，因此，步长就需要动态的调整。具体文献之后补充。

作者：刘兴禄，清华大学，清华伯克利深圳学院 (博士在读)
邮箱：hsinglul@163.com

刘兴禄

关注

37
点赞
踩
192

收藏

觉得还不错? 一键收藏
10
评论
Python实现最速下降法(The steepest descent method)详细案例

文章目录最速下降法(The steepest descent method)最速下降法的过程Python实现最速下降法实例`sympy`包中用到的函数构建符号变量和符号函数对符号函数求导求函数值求解方程的零点Python实现最速下降法求解上述算例的完整代码最速下降法(The steepest descent method)本文中的课件来自清华大学深圳国际研究生院，物流与交通学部张灿荣教授《高级运筹学》课程。在无约束非线性函数的最优化中，最速下降法(The steepest descent meth
复制链接

扫一扫